"""Central-to-mean density ratio of Lane–Emden polytropes -- numberdb.org/T615.

Run with SageMath (install numberdb once with sage -pip install numberdb):
    sage -python generate.py             # recompute and verify every stored row
    sage -python generate.py --fill      # write entries to this existing draft
Download solver.py and the files listed in `files` alongside this script.
NUMBERDB_API_KEY permits draft reads and entry writes. The campaign wrapper
also supports NUMBERDB_KEY_FROM_STDIN=1 and NUMBERDB_PUBLISH=1.
Neither command offers, reviews, or publishes the draft.

The ratio is C_n = xi_1**3/(3*mu_n). The incompressible value C_0=1 is
exact. C_1=pi**2/3 uses Arb. Nonlinear models are freshly integrated with
150/210 working decimal digits, orders 160/210, step caps 1/12 and 1/18,
and central matching radii 1/2 and 2/3. Only 100 significant digits are
stored. Agreement does not certify ODE truncation or propagation errors.
The solver is the unchanged attached solver from T613 and T614; it computes
both radius and mass, rather than dividing their rounded stored entries.
The selection and its distinct roles are explained in argument-selection.md.
The mathematical domain is rational 0 <= n < 5; this implementation supports
only the selected integer/half-integer models. No precomputed numerical file
is used to generate values. independent_checks.py uses the radial coordinate
and a virial identity to check the surface-coordinate solver independently.
"""
import os
import sys
import json
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from solver import radius

INDICES = ('0','1/2','1','3/2','2','5/2','3','7/2','4','9/2')
CONFIG = {150:(160,12,'1/2'), 210:(210,18,'2/3')}

@lru_cache(maxsize=32)
def background(n, working):
    order, divisor, start = CONFIG[working]
    x,q,meta = radius(n,working,order,divisor,start,trace=True)
    ratio = x**3/(3*q)
    if os.environ.get('NUMBERDB_EVIDENCE') == '1':
        print('BACKGROUND '+json.dumps({'n':n,'working':working,'radius':str(x),
              'mass':str(q),'ratio':str(ratio),'meta':meta}),flush=True)
    return x,q,ratio

class LaneEmdenDensityRatio(numberdb.Generator):
    table = 'T615'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'heuristic (agreement-checked)'
    files = ('generate.py','solver.py','validate.py','check.py','dry_run.py',
             'argument-selection.md','numerical-evidence.json','independent_checks.py',
             'independent-evidence.json','readback_checks.py','skill-evidence.json',
             'source-manifest.json')

    def enumerate(self):
        for n in INDICES:
            yield {'n':n}

    def value(self, params, digits):
        n = str(QQ(params['n']))
        if n not in INDICES:
            raise ValueError('Use a selected integer or half-integer index in [0,5).')
        if digits > 100:
            raise ValueError('This precision policy was validated for at most 100 digits.')
        if n == '0':
            return {'number':QQ(1),'comment':r'The incompressible model has constant interior density, so $\rho_c=\bar\rho$.'}
        if n == '1':
            B=RealBallField(numberdb.bits(digits,losing=80))
            return {'number':B.pi()**2/3,'comment':r'The linear profile $\theta_1(\xi)=\sin(\xi)/\xi$ gives $C_1=\pi^2/3$.'}
        return numberdb.agreeing(lambda working:str(background(n,working)[2]),at=(150,210))

if __name__ == '__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
    generator=LaneEmdenDensityRatio()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        for params in generator.enumerate(): generator.value(params,generator.digits)
        print(generator.publish(message='Add ten checked Lane-Emden density ratios and reproducible sources.'))
    else:
        report=generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)