"""Independent checks of T470. Run: sage -python checks.py.

With NUMBERDB_API_KEY (or NUMBERDB_KEY_FROM_STDIN=1), also check all stored
entries and the n=3 relation with T468. The exact polygon integration below
uses arbitrary convex polygon vertices and shares no formula or helper with
the generator. Disk quadrature integrates signed hull-edge areas, rather than
counting vertices and subtracting missing area. It is diagnostic; Arb in the
generator supplies the certified bounds. Simulations use a fixed random seed.
"""
import json
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from check import exactness, measure, agrees_with, names_its_rings
from generate import ExpectedHullVolume, expected_ratio, REGIONS, MAX_N
import mpmath as mp


def cross(a, b):
    return a[0]*b[1] - a[1]*b[0]


def polygon_ratio(vertices, n):
    """Integrate the signed area contributed by every possible hull edge.

    Boundary chord A(s)->B(t) intersects two different polygon sides. If
    D=B-A, the endpoint change of variables has factor
    -cross(D,A')*cross(D,B')*(v-u) for 0<u<v<1. The signed cross product of
    the two sampled points is (v-u)*cross(A,B); their remaining n-2 points
    are on the left with probability (left_area/area)^(n-2).
    Integral of (v-u)^2 on the unit triangle is 1/12. The shoelace formula
    then supplies n(n-1)/2. Every other integration is an exact polynomial
    integral over 0<s,t<1.
    """
    P = PolynomialRing(QQ, names=('s', 't'))
    s,t = P.gens()
    vertices = [tuple(map(QQ,v)) for v in vertices]
    m = len(vertices)
    area = sum(cross(vertices[i],vertices[(i+1)%m]) for i in range(m))/2
    assert area > 0
    total = QQ(0)
    for i in range(m):
        for j in range(m):
            if i == j:
                continue
            a0, a1 = vertices[i], vertices[(i+1)%m]
            b0, b1 = vertices[j], vertices[(j+1)%m]
            da = (a1[0]-a0[0],a1[1]-a0[1])
            db = (b1[0]-b0[0],b1[1]-b0[1])
            a = (a0[0]+s*da[0],a0[1]+s*da[1])
            b = (b0[0]+t*db[0],b0[1]+t*db[1])
            d = (b[0]-a[0],b[1]-a[1])
            chain = [a,b]
            k=(j+1)%m
            while k != (i+1)%m:
                chain.append(vertices[k]); k=(k+1)%m
            left = sum(cross(chain[k],chain[(k+1)%len(chain)]) for k in range(len(chain)))/2
            jacobian = -cross(d,da)*cross(d,db)
            poly = cross(a,b)*jacobian*(left/area)**(n-2)
            total += sum(c/((e[0]+1)*(e[1]+1)) for e,c in poly.dict().items())
    return QQ(n*(n-1))/24 * total/area**3


def disk_area(n, precision):
    """Direct signed shoelace-area integral in line-normal coordinates."""
    with mp.workdps(precision):
        def f(theta):
            sn,cs=mp.sin(theta),mp.cos(theta)
            fleft=(theta-sn*cs)/mp.pi
            return -cs*sn**5*fleft**(n-2)
        return + (mp.mpf(4)*n*(n-1)/(3*mp.pi**2)
                  *mp.quad(f,[0,mp.pi/4,mp.pi/2,3*mp.pi/4,mp.pi]))


def decimal_ball_text(x):
    with mp.workdps(180):
        return mp.mpf(str(x.center()))


def stored_entries(tid):
    doc=numberdb.table(tid)
    out={}
    # Numbers is a nested parameter dictionary in website documents.
    for region, rows in doc['Numbers'].items():
        for n, entry in rows.items():
            out[(region,int(n))]=entry['number'] if isinstance(entry,dict) else entry
    return out


def text_bounds(text):
    from decimal import Decimal
    from fractions import Fraction
    s=str(text)
    if '/' in s:
        q=Fraction(s);return QQ(q.numerator)/q.denominator, QQ(q.numerator)/q.denominator
    if '.' not in s and 'e' not in s.lower():return QQ(s),QQ(s)
    d=Decimal(s); q=Fraction(d); ulp=Fraction(10)**d.as_tuple().exponent
    return QQ((q-ulp).numerator)/(q-ulp).denominator, QQ((q+ulp).numerator)/(q+ulp).denominator


def main():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    gen=ExpectedHullVolume()
    values={(p['region'],int(p['n'])):gen.value(p,100) for p in gen.enumerate()}
    bare={k:v['number'] if isinstance(v,dict) else v for k,v in values.items()}
    assert not exactness(values)
    assert not names_its_rings(os.path.join(os.path.dirname(__file__),'generate.py'))
    triangle=[(0,0),(1,0),(0,1)]
    square=[(0,0),(1,0),(1,1),(0,1)]
    # Controls with answers known before the integration is trusted.
    assert polygon_ratio(triangle,3)==QQ(1)/12
    assert polygon_ratio(square,3)==QQ(11)/144
    assert polygon_ratio(triangle,4)==QQ(1)/6
    assert polygon_ratio(square,4)==QQ(11)/72
    computed={k:v for k,v in bare.items() if k[0]!='unit-disk'}
    independent={k:polygon_ratio(triangle if k[0]=='equilateral-triangle' else square,k[1]) for k in computed}
    assert not agrees_with(computed,lambda k:independent[k])
    print('Exact polygon edge integrals: 36/36, including normalization controls.',file=sys.stderr)
    # An invertible, non-similarity affine map and a translation.
    def affine(v):return (2*v[0]+v[1]+3,v[0]+3*v[1]-2)
    for k in computed:
        vs=triangle if k[0]=='equilateral-triangle' else square
        assert polygon_ratio([affine(v) for v in vs],k[1])==independent[k]
    print('Exact affine-transformed polygon integrals: 36/36.',file=sys.stderr)
    R=RealBallField(600)
    pi=R.pi()
    special={3:R(35)/(48*pi**2),4:R(35)/(24*pi**2),5:R(175)/(72*pi**2)-R(23023)/(6912*pi**4)}
    controls={}
    with mp.workdps(180):
        # A quadrature control for the same direct edge-area integrator.
        ctrl=disk_area(3,170)
        assert abs(ctrl-mp.mpf(35)/(48*mp.pi**2))<mp.mpf('1e-160')
        max_difference=mp.mpf(0)
        for n in range(3,MAX_N+1):
            a=disk_area(n,135); b=disk_area(n,175)
            assert abs(a-b)<mp.mpf('1e-130')
            x=bare[('unit-disk',n)]
            discrepancy=abs(decimal_ball_text(x)-b)
            assert discrepancy<mp.mpf('1e-118')
            max_difference=max(max_difference,discrepancy)
            controls[str(n)]={'independent_area_integral':mp.nstr(b,130),'difference':mp.nstr(discrepancy,6)}
            if n in special:
                assert special[n].is_finite() and R(special[n].rad())<R(10)**-170
                assert x.overlaps(special[n])
    print('Independent disk signed-area quadratures: 18/18 at 135 and 175 digits.',file=sys.stderr)
    for region in REGIONS:
        for n in range(3,MAX_N):
            assert bare[(region,n)] < bare[(region,n+1)]
    for region in REGIONS:
        assert (R(bare[(region,4)])-2*R(bare[(region,3)])).contains_zero()
    relative_radii=[float(v.rad()/abs(v.center())) for k,v in bare.items() if k[0]=='unit-disk']
    report={'entries':len(values),'exact_rationals':len(computed),'certified_disk_enclosures':18,
            'measure':measure(values),'max_relative_radius':max(relative_radii),
            'polygon_geometric_integrals':36,'affine_checks':36,
            'disk_quadrature_precisions':[135,175],'disk_checks':controls,
            'known_disk_closed_forms_checked':[3,4,5],
            'monotonicity_checks':51,'v4_equals_2v3_checks':3}
    if os.environ.get('NUMBERDB_API_KEY'):
        stored=stored_entries('T470')
        assert set(stored)==set(bare)
        with mp.workdps(180):
            for k,v in stored.items():
                lo,hi=text_bounds(v)
                if k in independent: assert lo==hi==independent[k]
                else:
                    b=mp.mpf(controls[str(k[1])]['independent_area_integral'])
                    assert mp.mpf(str(lo.numerator()))/mp.mpf(str(lo.denominator())) < b < mp.mpf(str(hi.numerator()))/mp.mpf(str(hi.denominator()))
                    x=bare[k];assert R(lo)<=R(R(x).lower()) and R(R(x).upper())<=R(hi)
        # Check the neighboring table's actual stored convex-position probabilities.
        neighbor=stored_entries('T468')
        for region in REGIONS:
            plo,phi=text_bounds(neighbor[(region,4)])
            vlo,vhi=text_bounds(stored[(region,3)])
            assert max(vlo,(1-phi)/4)<=min(vhi,(1-plo)/4)
        report['stored_independent_checks']=54
        report['stored_sylvester_relation_checks']=3
        for region in REGIONS:
            for n in range(3,MAX_N):
                assert text_bounds(stored[(region,n)])[1] < text_bounds(stored[(region,n+1)])[0]
        report['stored_monotonicity_checks']=51
        print('Stored entries and T468 relation checked.',file=sys.stderr)
    print(json.dumps(report,indent=2))


if __name__=='__main__':
    main()