"""Independent checks for the Thomson energies in numberdb.org/T478.

Run beside generate.py, check.py and sources.json:
    sage -python checks.py > numerical-checks.json
If stored-T478.json is present, the same checks also test the stored decimals.
No network or optimization is performed. Printed source energies are compared
with a direct sum over the independently supplied Cartesian coordinates.
"""
import hashlib
import json
from fractions import Fraction
from itertools import combinations, product
from pathlib import Path

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.version import version as sage_version
from numberdb._write import to_text
from generate import ThomsonEnergies, PROVEN, SOURCES, closed_form
from check import exactness, measure, names_its_rings


def rational(text):
    q = Fraction(text)
    return QQ(q.numerator) / QQ(q.denominator)


def interval(text, R):
    if '/' in text or ('.' not in text and 'e' not in text.lower()):
        return R(rational(text)), QQ(0)
    from decimal import Decimal
    d = Decimal(text)
    return R(rational(text)), QQ(10) ** int(d.as_tuple().exponent)


def contained(value, text, R):
    center, radius = interval(text, R)
    assert value.is_finite()
    return bool(abs(value - center) < R(radius)) if radius else bool(value == center)


def points_from_source(n, R):
    words = SOURCES['deeter'][str(n)]['coordinate_text'].split()
    assert len(words) == 3*n, (n, len(words))
    return [tuple(R(rational(w)) for w in words[k:k+3]) for k in range(0,3*n,3)]


def normalise(points, R):
    out=[]
    for point in points:
        norm2=sum((x*x for x in point),R(0))
        assert norm2 > 0
        out.append(tuple(x/norm2.sqrt() for x in point))
    return out


def distance_energy(points, R):
    total = R(0)
    for a,b in combinations(points,2):
        d2=sum(((x-y)**2 for x,y in zip(a,b)),R(0))
        assert d2 > 0
        total += 1/d2.sqrt()
    return total


def dot_energy(points, R):
    total=R(0)
    for i in range(len(points)):
        for j in range(i):
            d2=2-2*sum((points[i][k]*points[j][k] for k in range(3)),R(0))
            assert d2 > 0
            total += 1/d2.sqrt()
    return total


def exact_points(n,R):
    zero=R(0); one=R(1)
    if n==2:return [(zero,zero,one),(zero,zero,-one)]
    if n in (3,5):
        tri=[(one,zero,zero),(-one/2,R(3).sqrt()/2,zero),(-one/2,-R(3).sqrt()/2,zero)]
        return tri if n==3 else tri+[(zero,zero,one),(zero,zero,-one)]
    if n==4:
        return [tuple(R(x)/R(3).sqrt() for x in t) for t in [(1,1,1),(1,-1,-1),(-1,1,-1),(-1,-1,1)]]
    if n==6:
        return [tuple(R(sign) if j==axis else zero for j in range(3)) for axis in range(3) for sign in (-1,1)]
    if n==12:
        phi=(1+R(5).sqrt())/2
        return normalise([(zero,R(s),t*phi) for s,t in product((-1,1),repeat=2)]+[(R(s),t*phi,zero) for s,t in product((-1,1),repeat=2)]+[(t*phi,zero,R(s)) for s,t in product((-1,1),repeat=2)],R)
    raise ValueError(n)


def run():
    R=RealBallField(numberdb.bits(60,losing=64))
    S=RealBallField(numberdb.bits(90,losing=64))
    H=RealBallField(numberdb.bits(200,losing=64))
    # Known positive controls first: these must not all silently return zero.
    assert distance_energy([(R(0),R(0),R(1)),(R(0),R(0),R(-1))],R)==R(QQ(1)/2)
    control=distance_energy(exact_points(3,R),R)
    assert abs(control-R(3).sqrt()) < R(10)**(-55)
    assert control > 1
    assert not contained(R(2),'1.000000',R)
    values={int(p['n']):ThomsonEnergies().value(p,100) for p in ThomsonEnergies().enumerate()}
    assert not exactness(values),exactness(values)
    assert not names_its_rings(Path(__file__).with_name('generate.py'))
    written={n:to_text(v['number'],v.get('digits',100)) for n,v in values.items()}
    report={'sage_version':sage_version,'controls':'antipodal energy 1/2, triangle energy sqrt(3), and false comparison rejected','entries':len(values),'measurement':measure({n:{**values[n],'number':written[n]} for n in values}),'closed_forms':{},'coordinate_checks':{}}
    max_norm=R(0); max_source_ulp=R(0); max_ccd_ulp=R(0); max_stored_ulp=R(0)
    for n in PROVEN:
        point_energy=distance_energy(exact_points(n,H),H)
        assert abs(2*distance_energy([tuple(2*x for x in p) for p in exact_points(n,H)],H)-point_energy)<H(10)**(-190)
        direct=closed_form(n,100)
        assert contained(point_energy,written[n],H),(n,written[n])
        if n!=2:
            assert direct.is_finite() and R(direct.rad()) < R(10)**(-110)
        if n>=4:
            assert contained(point_energy,SOURCES['deeter'][str(n)]['energy'],H),n
        report['closed_forms'][str(n)]={'stored':written[n],'accuracy_bits':None if n==2 else int(direct.accuracy()),'independent_coordinate_sum':str(point_energy)}
    for n in range(4,101):
        raw=points_from_source(n,R)
        dev=max(abs(sum((x*x for x in p),R(0))-1) for p in raw)
        assert dev<R(10)**(-10),(n,str(dev))
        max_norm=max(max_norm,dev)
        pts=normalise(raw,R)
        energy=distance_energy(pts,R)
        other=dot_energy(normalise(points_from_source(n,S),S),S)
        assert energy.is_finite() and other.is_finite()
        assert R(energy.rad())<R(10)**(-50) and S(other.rad())<S(10)**(-80)
        assert abs(S(energy)-other)<S(10)**(-50),(n,'different formulas disagree')
        quoted=SOURCES['deeter'][str(n)]['energy']
        center,ulp=interval(quoted,R)
        discrepancy=abs(energy-center)/ulp
        assert discrepancy < 2, (n, str(discrepancy))
        max_source_ulp=max(max_source_ulp,discrepancy)
        # Source coordinate precision limits the check, not the exact rows.
        if n not in PROVEN:
            assert contained(energy,written[n],R),n
            sc,su=interval(written[n],R)
            max_stored_ulp=max(max_stored_ulp,abs(energy-sc)/su)
        if n>=10:
            ctext=SOURCES['cambridge'][str(n)]
            assert contained(energy,ctext,R),(n,'CCD',str(energy),ctext)
            ccenter,culp=interval(ctext,R)
            max_ccd_ulp=max(max_ccd_ulp,abs(energy-ccenter)/culp)
        assert abs(distance_energy([tuple(2*x for x in p) for p in pts],R)*2-energy)<R(10)**(-50)
        report['coordinate_checks'][str(n)]={'energy':str(energy),'dot_energy_90_digits':str(other),'source':quoted,'source_error_in_last_places':str(discrepancy),'raw_norm2_error':str(dev)}
    # Jensen's lower bound holds for every configuration, with equality for simplices.
    for n,text in written.items():
        center,ulp=interval(text,H)
        bound=H(QQ(n*(n-1))/2)*H(QQ(n-1)/(2*n)).sqrt()
        if n in (2,3,4):assert abs(center-bound)<=H(ulp)+H(10)**(-190)
        else: assert center-H(ulp)>bound,n
    report['max_stored_error_in_last_places']=str(max_stored_ulp)
    report['max_source_error_in_last_places']=str(max_source_ulp)
    report['max_cambridge_error_in_last_places']=str(max_ccd_ulp)
    report['max_raw_norm2_error']=str(max_norm)
    stored_path=Path(__file__).with_name('stored-T478.json')
    if stored_path.exists():
        doc=json.loads(stored_path.read_text())
        block=doc['Numbers']
        if isinstance(block, dict):
            stored={int(k):(v['number'] if isinstance(v,dict) else v) for k,v in block.items()}
        else:
            stored={int(e['params']['n']):e['number'] for e in block}
        assert set(stored)==set(range(2,101))
        assert list(stored)==list(range(2,101))
        for n,text in stored.items():
            if n in PROVEN: value=distance_energy(exact_points(n,H),H); field=H
            else:value=distance_energy(normalise(points_from_source(n,R),R),R);field=R
            assert contained(value,text,field),(n,text)
            assert text==written[n],(n,'different stored spelling')
        report['stored_rows_checked']=len(stored)
    report['ok']=True
    return report


if __name__=='__main__':
    print(json.dumps(run(),indent=2))