"""Check stored T615 intervals, analytic anchors and companion identities.
Run: sage -python readback_checks.py, with NUMBERDB_API_KEY for draft reads.
The campaign wrapper supports NUMBERDB_KEY_FROM_STDIN=1. Companion tables are
read through the public API with this account's draft authorization. The
repository check.stored helper requires Django and is not used on this host.
"""
import os,sys,json,pathlib
from decimal import Decimal,localcontext
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from generate import INDICES

def number(item):
    return str(item['number'] if isinstance(item,dict) else item)

def interval(text):
    d=Decimal(text)
    if '.' not in text and 'e' not in text.lower():return d,d,Decimal(0)
    ulp=Decimal(10)**d.as_tuple().exponent
    return d-ulp,d+ulp,ulp

def main():
    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('\"\'')
        os.environ['NUMBERDB_API_KEY']=token
    d=numberdb.table('T615');radii=numberdb.table('T613');masses=numberdb.table('T614')
    assert d['Title']=='Central-to-mean density ratio of Lane–Emden polytropes'
    assert d['Display properties']['number-header']==r'$C_n=\rho_c/\bar\rho$'
    assert d['Parameters']['n']['type']=='Q'
    assert list(d['Numbers'])==list(INDICES)
    p=pathlib.Path(__file__).parent
    computed=json.loads((p/'numerical-evidence.json').read_text())['values']
    independent=json.loads((p/'independent-evidence.json').read_text())['values']
    skill=json.loads((p/'skill-evidence.json').read_text())
    assert skill['match'] is True
    assert set(skill['versions'].values())=={'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199'}
    report={'table':'T615','skill_manifest':skill['manifest'],'rows':{}}
    B=RealBallField(700)
    with localcontext() as ctx:
        ctx.prec=250
        for n in INDICES:
            s=number(d['Numbers'][n]);lo,hi,ulp=interval(s);center=Decimal(s)
            if n=='0':
                assert s=='1'
            else:
                assert len(center.as_tuple().digits)>=100,(n,'too few digits')
                if n=='1':
                    target=B.pi()**2/3
                    assert target.is_finite() and B(target.rad())<B('1e-190')
                    assert B(str(lo))<target and target<B(str(hi))
                for work in ['150','210']:
                    assert lo<Decimal(computed[n][work]['ratio'])<hi,(n,work)
            error=abs(Decimal(independent[n]['ratio'])-center)/abs(center)
            assert error<Decimal('1e-40'),(n,error)
            rlo,rhi,_=interval(number(radii['Numbers'][n]))
            mlo,mhi,_=interval(number(masses['Numbers'][n]))
            assert rlo<Decimal(computed[n]['210']['radius'])<rhi,('radius',n)
            assert mlo<Decimal(computed[n]['210']['mass'])<mhi,('mass',n)
            # Positive interval arithmetic: the stored companion intervals must
            # give a concentration interval intersecting this row's interval.
            clo=rlo**3/(3*mhi);chi=rhi**3/(3*mlo)
            assert max(clo,lo)<=min(chi,hi),('companion identity',n)
            # A separate integral identity using the stored concentration to
            # infer mass from the independently integrated radius.
            xi=Decimal(independent[n]['radius']);J=Decimal(independent[n]['pressure_integral'])
            mu=xi**3/(3*center)
            nn=Decimal(n.split('/')[0])/(Decimal(n.split('/')[1]) if '/' in n else 1)
            virial=(nn+1)*mu*mu/((5-nn)*xi)
            verr=abs(J-virial)/abs(virial)
            assert verr<Decimal('1e-40'),('stored virial',n,verr)
            report['rows'][n]={'exact':n=='0','stored_interval_pass':True,
               'significant_digits':None if n=='0' else len(center.as_tuple().digits),
               'independent_relative_difference':str(error),'companion_identity_pass':True,
               'stored_virial_relative_error':str(verr)}
        anchor=number(d['Numbers']['1']);lo,hi,ulp=interval(anchor)
        assert not lo<Decimal(anchor)+100*ulp<hi
        report['perturbed_value_rejected']=True
    print('DATA '+json.dumps(report),flush=True)
if __name__=='__main__':main()