"""Independent checks for external simplex angles, numberdb.org/T597.

    sage -python checks.py                 # compute and check locally
    sage -python checks.py --stored        # also check every stored interval
For a private draft set NUMBERDB_API_KEY. Campaign wrappers can instead set
NUMBERDB_KEY_FROM_STDIN=1 and pipe the key on stdin. This file writes nothing.

The independent numeric path uses mpmath tanh-sinh integration after the
change of variables x = sqrt(2*(k+1))*t. It shares the published integral,
not the generator's code or quadrature library. It is an agreement check,
not the basis for certification. The geometric wedge and spherical-triangle
checks, normal Gram matrices and vertex normal-fan partition provide checks
outside the integral implementation. Controls run before the full checks.
"""
import sys
import os
import json
import time
from math import comb
from decimal import Decimal
from fractions import Fraction

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
import mpmath as mp
from generate import ExternalSimplexAngles, quadrature


def independent(d, k, digits):
    with mp.workdps(digits):
        m = d-k
        q = mp.mpf(k+1)
        def integrand(t):
            # The pair +/-t avoids subtracting nearly equal erfs in the tails.
            u=mp.erfc(-t)/2
            v=mp.erfc(t)/2
            return mp.exp(-q*t*t)*(u**m+v**m)
        value=mp.sqrt(q/mp.pi)*mp.quad(integrand,[0,1,2,4,8,mp.inf],method='tanh-sinh')
        assert mp.isfinite(value)
        return +value


def mp_of(x):
    if x in QQ:
        z=QQ(x)
        return mp.mpf(str(z.numerator()))/mp.mpf(str(z.denominator()))
    return mp.mpf(str(x))


def rational_decimal(s):
    z=Fraction(Decimal(s))
    return QQ(z.numerator)/QQ(z.denominator)


def stored_records(document):
    # Reuse the client normalizer: API tables may be nested by parameter.
    from numberdb._generate import _stored_entries
    assert list(document['Parameters']) == ['d', 'k']
    return [{'params': dict(zip(('d', 'k'), identity.split(','))),
             'number': value}
            for identity, value in _stored_entries(document).items()]


def written_interval(text):
    text=str(text)
    if '/' in text or ('.' not in text and 'e' not in text.lower()):
        x=QQ(text); return x,x
    v=Decimal(text); ulp=Decimal(1).scaleb(v.as_tuple().exponent)
    center=rational_decimal(text)
    radius=rational_decimal(str(ulp))
    return center-radius,center+radius


def geometry(values):
    R=RealBallField(500)
    counts={'wedge':0,'spherical_triangle':0,'normal_gram':0,'vertices':0,'facets':0,'reflection':0}
    for (d,k), v in values.items():
        m=d-k
        # Normals centroid-e_i have Gram I_m-J_m/(d+1).
        # The inverse is I_m+J_m/(k+1), exactly the claimed covariance.
        for i in range(m):
            for j in range(m):
                h=sum(((QQ(1) if i==l else QQ(0))-QQ(1)/(d+1))*
                      ((QQ(1) if l==j else QQ(0))+QQ(1)/(k+1)) for l in range(m))
                assert h==(1 if i==j else 0)
        counts['normal_gram']+=1
        v=R(v)
        if k==0:
            assert QQ(1)/(d+1) in v; counts['vertices']+=1
        if m==1:
            assert QQ(1)/2 in v; counts['facets']+=1
        if m==2:
            wedge=(-R(1)/d).arccos()/(2*R.pi())
            assert v.overlaps(wedge); counts['wedge']+=1
        if m==3:
            triangle=(3*(-R(1)/(d-1)).arccos()-R.pi())/(4*R.pi())
            assert v.overlaps(triangle); counts['spherical_triangle']+=1
        if m>=3 and m%2:
            def g(j):
                if j==0: return R(1)
                if j==1: return R(QQ(1)/2)
                return R(values[(k+j,k)])
            reflected=sum(((-1)**j)*comb(m,j)*g(j) for j in range(m))-2*v
            assert reflected.contains_zero(); counts['reflection']+=1
    return counts


def main():
    start=time.monotonic()
    mp.mp.dps=165
    # A normalization control returning a known nonzero answer, then a
    # nontrivial vertex value, with no generator shortcuts.
    for d,k,exact in [(2,1,QQ(1)/2),(2,0,QQ(1)/3),(7,0,QQ(1)/8)]:
        assert abs(independent(d,k,125)-mp_of(exact))<mp.mpf('1e-120')
        assert exact in quadrature(d,k)
    print('CONTROLS PASS',flush=True)
    generator=ExternalSimplexAngles()
    params=list(generator.enumerate())
    assert [(p['d'],p['k']) for p in params]==[(d,k) for d in range(2,13) for k in range(d)]
    values={}; independent_values={}; records=[]
    min_bits=10**9; max_relative=mp.mpf(0); max_difference=mp.mpf(0)
    for p in params:
        d,k=p['d'],p['k']; key=(d,k)
        value=generator.value(p,100); values[key]=value
        a=independent(d,k,125); b=independent(d,k,155)
        assert abs(a-b)<mp.mpf('1e-119'), (key,'precision agreement')
        # Independent midpoint differences need not lie within the much
        # narrower certification radius: allow 110 digits for this check.
        center=value.center() if hasattr(value,'center') else value
        discrepancy=abs(mp_of(center)-b)
        assert discrepancy<mp.mpf('1e-110'), (key,'independent comparison')
        max_difference=max(max_difference,discrepancy)
        if hasattr(value,'rad'):
            relative=mp_of(value.rad())/abs(mp_of(value.center()))
            assert relative<mp.mpf('1e-110')
            max_relative=max(max_relative,relative)
            min_bits=min(min_bits,int(value.accuracy()))
        independent_values[key]=b
        records.append({'d':d,'k':k,'certified':str(value),'independent_155':mp.nstr(b,155)})
        if k==d-1: print('CHECKED dimension',d,'elapsed',round(time.monotonic()-start,2),flush=True)
    counts=geometry(values)
    stored_count=0
    if '--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
            token=sys.stdin.read().strip()
            if '=' in token and token.split('=',1)[0].replace('_','').isalnum():
                token=token.split('=',1)[1].strip().strip('\"\'')
            numberdb.configure(api_key=token)
        report=generator.verify(sample=None)
        assert report.ok,str(report)
        print('GENERATOR VERIFY',str(report),flush=True)
        document=numberdb.table('T597')
        data=stored_records(document)
        assert isinstance(data,list) and len(data)==77
        stored_values={}
        for entry in data:
            p=entry['params']; key=(int(p['d']),int(p['k']))
            low,high=written_interval(entry['number'])
            expected=independent_values[key]
            if low==high:
                assert abs(mp_of(low)-expected)<mp.mpf('1e-145')
            else:
                assert mp_of(low)<expected<mp_of(high), (key,'stored interval')
            R=RealBallField(500)
            # Exact endpoint arithmetic; convert each endpoint separately,
            # then enclose their interval in a ball.
            stored_values[key]=R((low+high)/2,R((high-low)/2).upper())
            stored_count+=1
        assert list(stored_values)==list(values)
        geometry(stored_values)
        print('STORED INTERVALS AND GEOMETRY PASS',stored_count,flush=True)
    evidence={'table':'T597','entries':77,'exact_entries':22,'ball_entries':55,
              'min_relative_accuracy_bits':min_bits,
              'max_relative_radius':mp.nstr(max_relative,12),
              'max_independent_difference':mp.nstr(max_difference,12),
              'checks':counts,'stored_checked':stored_count,
              'environment':generator.environment(),'mpmath_version':mp.__version__,
              'seconds':round(time.monotonic()-start,2),'records':records}
    print('EVIDENCE_JSON '+json.dumps(evidence),flush=True)


if __name__=='__main__':
    main()