"""Independent checks for Growth constants of the classes of planar graphs, T481.

sage -python checks.py                  # local computations only
sage -python checks.py --stored         # also check the stored draft
The latter needs NUMBERDB_API_KEY while the table is private. The campaign
may instead set NUMBERDB_KEY_FROM_STDIN=1 and pipe its key to stdin.

Uses the planar-map characteristic equations and the network equation
(Giménez--Noy, equations 2.2, 2.3, 2.5) independently of the generator's Y.
The all/connected radius is recomputed from the appendix's B2 coefficient,
rather than the generator's explicit r(t). mpmath calculations at 160 and
200 decimal places are numerical cross-checks; the generator's Arb brackets
and ball evaluation are what certify its digits.

The maximal-planar value is checked using Tutte's factorial count for rooted
simple triangulations; the 3-connected value from the quartic map equation's
critical point. These test the graph/map and vertex/edge conventions.
"""
import json
import os
import sys
from math import factorial
import mpmath as mp
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from check import exactness, measure, names_its_rings, agrees_with
from generate import PlanarGrowth, constants, root_bracket, radii, y_of_t


def independent(dps):
    with mp.workdps(dps):
        # Control for root finding: sqrt(2), evaluated by an independent path.
        control = mp.findroot(lambda x: x*x-2, (mp.mpf(1),mp.mpf(2)))
        assert abs(control-mp.sqrt(2)) < mp.mpf(10)**(-dps+5)

        def network(u):
            v = (u+1)/(3*u-1)
            z = v/(1+u)**2
            x = u/(z*(1+v)**2)
            m = x*x*z*z*(1/(1+x*z)+1/(1+z)-1-(1+u)**2*(1+v)**2/(1+u+v)**3)
            return m/(2*x*x*z)-mp.log((1+z)/2)+x*z*z/(1+x*z)

        u = mp.findroot(network,(mp.mpf('0.53'),mp.mpf('0.54')))
        assert abs(network(u)) < mp.mpf(10)**(-dps+5)
        v = (u+1)/(3*u-1)
        d = v/(1+u)**2
        r = u/(d*(1+v)**2)
        assert abs(4*r*d*d*(1+u)*(1+v)-1) < mp.mpf(10)**(-dps+5)
        t = 1/(3*u)
        b2 = (1-t)**3*(3*t-1)*(1+3*t)*(1+t)**3*mp.log(1+t)/(256*t**6)
        b2 -= (1-t)**3*(1+3*t)*mp.log(1+3*t)/(32*t**3)
        b2 += (1+3*t)**2*(1-t)**6*mp.log(1+2*t)/(512*t**6)
        b2 += (1-t)**4*(185*t**4+698*t**3-217*t**2-160*t+6)/(1024*t**4*(1+3*t)*(3+t))
        gamma = mp.exp(-b2/r)/r
        # Map critical point with edge variable z=1.
        u3 = mp.findroot(lambda w: (3*w-1)*(1+w)-1,(mp.mpf('0.54'),mp.mpf('0.56')))
        r3 = u3/(1+(1+u3)**2)**2
        answers = {'all':gamma,'connected':gamma,'2-connected':1/r,
                   '3-connected':1/r3,'maximal':mp.mpf(256)/27}
        return {k:mp.nstr(v,dps-5) for k,v in answers.items()}


def as_decimal_interval(text):
    """NumberDB decimal convention: exact rational or +/- one final place."""
    if '/' in text or ('.' not in text and 'e' not in text.lower()):
        q = QQ(text)
        return q,q
    mantissa,_,exp = text.lower().partition('e')
    exponent = int(exp or '0')
    sign = -1 if mantissa.startswith('-') else 1
    mantissa = mantissa.lstrip('+-')
    before,_,after = mantissa.partition('.')
    scale = QQ(10)**(exponent-len(after))
    middle = sign*QQ(before+after)*scale
    return middle-scale,middle+scale


def check_rows(rows, high):
    field = RealBallField(700)
    assert set(rows)==set(high)
    for key,text in rows.items():
        lo,hi = as_decimal_interval(text)
        hlo,hhi = as_decimal_interval(high[key])
        v = field((hlo+hhi)/2)
        if key == 'maximal':
            assert lo == hi == QQ(256)/27
        else:
            assert field(lo) < v < field(hi),key
            assert hi-lo <= QQ(2)/10**98,key
    assert rows['all']==rows['connected']
    intervals={k:as_decimal_interval(v) for k,v in rows.items()}
    for smaller,larger in [('maximal','3-connected'),('3-connected','2-connected'),('2-connected','connected')]:
        assert intervals[smaller][1] < intervals[larger][0]
    # The exact algebraic polynomial has opposite signs at the stored bounds.
    lo,hi=intervals['3-connected']
    f=lambda x:27*x*x-544*x-512
    assert f(lo)<0<f(hi)
    print('Every written interval contains the independent value; equality, ordering, and algebraic bracket pass.')


def main():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    gen=PlanarGrowth()
    vals=constants(100)
    assert not exactness(vals)
    assert not names_its_rings(__file__.replace('checks.py','generate.py'))
    field,lo,hi=root_bracket(100)
    assert y_of_t(field(lo))<1<y_of_t(field(hi))
    assert hi-lo < QQ(1)/10**110
    for key,val in vals.items():
        if key!='maximal':
            assert val.is_finite() and QQ(val.rad()) < QQ(1)/10**108
    print('Arb root bracket, finite balls, and radii below 1e-108: PASS')
    low,high=independent(160),independent(200)
    with mp.workdps(210):
        for k in high:
            assert abs(mp.mpf(low[k])-mp.mpf(high[k]))<mp.mpf('1e-150')
    rows={k:to_text(v,100) for k,v in vals.items()}
    print(json.dumps(rows,indent=2))
    check_rows(rows,high)
    # Exact Tutte counts for n+2 vertices; check several known small cases.
    tri=lambda n: QQ(2*factorial(4*n-3))/(factorial(n)*factorial(3*n-1))
    assert [tri(n) for n in range(1,7)]==[1,1,3,13,68,399]
    ratio=lambda n: QQ((4*n-2)*(4*n-1)*(4*n)*(4*n+1))/((n+1)*(3*n)*(3*n+1)*(3*n+2))
    assert not agrees_with({n:tri(n+1)/tri(n) for n in range(1,51)},ratio)
    # Leading coefficients of the rational ratio prove the limit, exactly.
    assert QQ(4**4)/(3**3)==vals['maximal']
    print('Tutte counts, 50 exact ratios, and ratio limit 256/27: PASS')
    print('Independent network/map computations agree at 160 and 200 digits: PASS')
    print('Measurement:',measure({p['class']:gen.value(p,100) for p in gen.enumerate()}))
    if '--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        report=gen.verify(sample=None)
        print(report);assert report.ok
        doc=numberdb.table('T481')
        stored_rows={k.split(',')[0]:v for k,v in _stored_entries(doc).items()}
        check_rows(stored_rows,high)
        print('Stored values pass all independent interval checks.')


if __name__=='__main__':
    main()