"""Independent numerical and graph checks for numberdb.org/T626.

Run: sage -python checks.py > numerical-evidence.json
Dependencies: numberdb, mpmath, numpy, and the attached generate.py.
Quadrature/finite tori are independent checks, not the source of certification.
For a database readback, run stored_checks.py after filling the draft.
"""
import json
import sys
from collections import Counter, deque
from decimal import Decimal
from functools import lru_cache
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
import mpmath as mp
import numpy as np
import generate as g


def key(p):
    return ','.join(str(p[k]) for k in g.Resistance.parameters)


def unbox(v):
    return v['number'] if isinstance(v,dict) else v


def endpoint_rationals(v):
    if hasattr(v,'lower'): return v.lower().exact_rational(),v.upper().exact_rational()
    return QQ(v),QQ(v)


def written_bounds(text):
    if '.' not in text and 'e' not in text.lower():
        q=QQ(text)
        return q,q
    d=Decimal(text)
    center=QQ(Fraction(d))
    unit=QQ(10)**d.as_tuple().exponent
    return center-unit,center+unit


def bfs(lattice,radius):
    start=('A',0,0)
    dist={start:0}
    queue=deque([start])
    while queue:
        s,m,n=queue.popleft()
        if dist[s,m,n] == radius: continue
        if lattice == 'triangular':
            neighbors=[('A',m+a,n+b) for a,b in [(1,0),(-1,0),(0,1),(0,-1),(1,1),(-1,-1)]]
        elif s == 'A': neighbors=[('B',m,n),('B',m-1,n),('B',m-1,n-1)]
        else: neighbors=[('A',m,n),('A',m+1,n),('A',m+1,n+1)]
        for v in neighbors:
            if v not in dist:
                dist[v]=dist[s,m,n]+1
                queue.append(v)
    return dist


@lru_cache(None)
def integral(m,n,dps):
    with mp.workdps(dps):
        diff=abs(m-n)
        total=m+n
        def f(y):
            if not y: return mp.mpf(diff)
            c=mp.cos(y)
            if y == mp.pi/2:
                return (1-mp.cos(total*y))/2 if diff == 0 else mp.mpf(1)/2
            x=2*mp.asinh(mp.sin(y/2)*mp.sqrt((c+2)/c))
            decay=mp.exp(-diff*x)
            numerator=-mp.expm1(-diff*x)+2*decay*mp.sin(total*y/2)**2
            return numerator/(mp.sin(y)*mp.sqrt(3+mp.sin(y)**2))
        return mp.quad(f,[0,mp.pi/8,mp.pi/4,3*mp.pi/8,mp.pi/2],method='gauss-legendre')/mp.pi


def quadrature(p,dps):
    m,n=int(p['m']),int(p['n'])
    with mp.workdps(dps):
        if p['lattice']=='triangular': return integral(m,n,dps)
        if p['sublattice']=='A': return 3*integral(m,n,dps)
        return integral(m,n,dps)+integral(m+1,n,dps)+integral(m+1,n+1,dps)


def torus(N):
    theta=2*np.pi*np.arange(N)/N
    x,y=theta[:,None],theta[None,:]
    # Build each network's Fourier Laplacian from its own declared edges.
    lam=6-2*np.cos(x)-2*np.cos(y)-2*np.cos(x+y)
    lam[0,0]=1
    inv=1/lam
    inv[0,0]=0
    green=np.fft.ifft2(inv).real
    t=2*(green[0,0]-green)
    coupling=1+np.exp(1j*x)+np.exp(1j*(x+y))
    det=9-abs(coupling)**2
    det[0,0]=1
    diag=3/det
    cross=coupling/det
    # At zero wave vector the nonconstant basis mode has eigenvalue six.
    diag[0,0]=1/12
    cross[0,0]=-1/12
    aa=np.fft.ifft2(diag).real
    ba=np.fft.ifft2(cross).real
    return t,2*(aa[0,0]-aa),2*(aa[0,0]-ba)


def main():
    gen=g.Resistance()
    params=list(gen.enumerate())
    assert len(params)==len({key(p) for p in params})==84
    counts=Counter((p['lattice'],p['sublattice']) for p in params)
    assert counts=={('triangular','A'):48,('honeycomb','A'):15,('honeycomb','B'):21}
    graph_summary={}
    for lattice in ['triangular','honeycomb']:
        dist=bfs(lattice,12)
        reps=set()
        for (s,m,n),d in dist.items():
            if not d: continue
            expected=g.distance(m,n) if lattice=='triangular' else g.honey_distance(s,m,n)
            assert d==expected
            representative=g.canonical(m,n) if s=='A' else max(g.b_orbit(m,n))
            reps.add((s,*representative))
        selected={(p['sublattice'],int(p['m']),int(p['n'])) for p in params if p['lattice']==lattice}
        assert reps==selected
        shells=dict(sorted(Counter(dist.values()).items()))
        assert all(shells[d]==(6 if lattice=='triangular' else 3)*d for d in range(1,13))
        graph_summary[lattice]={'vertices_including_origin':len(dist),'shell_counts':shells,'orbits':len(reps)}
    # Published small values, evaluated before trusting quadrature.
    mp.mp.dps=160
    controls={(0,0):mp.mpf(0),(1,0):mp.mpf(1)/3,(2,0):mp.mpf(8)/3-4*mp.sqrt(3)/mp.pi,
              (2,1):-mp.mpf(2)/3+2*mp.sqrt(3)/mp.pi,(3,1):-5+10*mp.sqrt(3)/mp.pi}
    for mn,expected in controls.items():
        assert abs(integral(*mn,125)-expected)<mp.mpf('1e-118'),(mn,'quadrature control')
    print('Quadrature controls passed.',file=sys.stderr,flush=True)
    records=[]
    min_bits=10000
    max_error=mp.mpf(0)
    max_precision_difference=mp.mpf(0)
    max_cancellation=mp.mpf(0)
    for i,p in enumerate(params):
        poly=g.coefficients(p)
        assert all(c.parent() is QQ for c in poly.list())
        v=unbox(gen.value(p,100))
        lo,hi=endpoint_rationals(v)
        text=to_text(v,100)
        wlo,whi=written_bounds(text)
        assert wlo<=lo<=hi<=whi,(key(p),'serialization')
        if hi!=lo:
            assert v.is_finite() and lo>0
            relative=(hi-lo)/lo
            accuracy=relative.denominator().nbits()-relative.numerator().nbits()-1
            assert relative < QQ(2)**(-accuracy)
            assert accuracy>=350,(key(p),accuracy)
            min_bits=min(min_bits,int(accuracy))
        else: accuracy=None
        # Independent integral has no access to the recurrence coefficients.
        q1,q2=quadrature(p,125),quadrature(p,150)
        a,b=mp.mpf(str(poly[0])),mp.mpf(str(poly[1]))
        expected=a+b*mp.sqrt(3)/mp.pi
        error=abs(q2-expected)
        difference=abs(q1-q2)
        assert error<mp.mpf('1e-140'),(key(p),str(error))
        assert difference<mp.mpf('1e-118'),(key(p),str(difference))
        max_error=max(max_error,error)
        max_precision_difference=max(max_precision_difference,difference)
        cancellation=mp.log((abs(a)+abs(b*mp.sqrt(3)/mp.pi))/abs(expected),2)
        max_cancellation=max(max_cancellation,cancellation)
        records.append({'params':p,'number':text,'exact_coefficients':[str(poly[0]),str(poly[1])],
                        'relative_accuracy_bits_lower_bound':accuracy,'integral_150_digits':mp.nstr(q2,145),
                        'quadrature_absolute_error':mp.nstr(error,8)})
        if (i+1)%12==0: print('Checked %d/84 entries.'%(i+1),file=sys.stderr,flush=True)
    # Harmonicity at every selected triangular displacement, including the
    # additional outer neighbors needed to check the boundary shell.
    for p in params:
        m,n=int(p['m']),int(p['n'])
        if p['lattice']=='triangular':
            assert 6*g.triangular(m,n)==sum(g.triangular(m+a,n+b) for a,b in g.STEPS)
        else:
            s=p['sublattice']
            if s=='A': neighbors=[('B',m,n),('B',m-1,n),('B',m-1,n-1)]
            else: neighbors=[('A',m,n),('A',m+1,n),('A',m+1,n+1)]
            assert 3*g.coefficients(p)==sum(g.coefficients(dict(lattice='honeycomb',sublattice=s1,m=str(a),n=str(b))) for s1,a,b in neighbors)
    torus_summary=[]
    previous=None
    for N in [64,128,256]:
        t,ha,hb=torus(N)
        assert abs(t[1,0]-(N*N-1)/(3*N*N))<1e-12
        assert abs(hb[0,0]-(2*N*N-1)/(3*N*N))<1e-12
        errors=[]
        for p in params:
            array=t if p['lattice']=='triangular' else ha if p['sublattice']=='A' else hb
            observed=array[int(p['m'])%N,int(p['n'])%N]
            target=float(unbox(gen.value(p,100)))
            errors.append(abs(observed-target))
        if previous is not None: assert all(b<a for a,b in zip(previous,errors))
        previous=errors
        torus_summary.append({'side':N,'maximum_absolute_error':max(errors),'adjacent_controls_passed':True})
    evidence={'table':'T626','skill_manifest':'agents/runs/20261030T125612Z-build-skills/manifest.json',
       'skill_sha256': '28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
       'entries':84,'group_counts':{','.join(k):v for k,v in counts.items()},'graph_checks':graph_summary,
       'working_precision_bits':numberdb.bits(100,losing=g.WORKING_GUARD),'minimum_relative_accuracy_bits':min_bits,
       'maximum_cancellation_bits':mp.nstr(max_cancellation,12),'quadrature_working_digits':[125,150],
       'maximum_quadrature_error':mp.nstr(max_error,12),'maximum_precision_difference':mp.nstr(max_precision_difference,12),
       'exact_harmonicity_checks':84,'serialized_enclosures_checked':84,'torus_checks':torus_summary,
       'versions':{'numberdb':getattr(numberdb,'__version__','unknown'),'mpmath':mp.__version__,'numpy':np.__version__},
       'records':records}
    print(json.dumps(evidence,indent=2))


if __name__=='__main__': main()