"""Growth constants of polyominoes and lattice animals -- numberdb.org/T484.

Run with SageMath (or a compatible passagemath installation):
    sage -pip install numberdb
    sage -python generate.py              # verify all stored values
    sage -python checks.py                # offline numerical and convention checks

Reading the unpublished draft requires NUMBERDB_API_KEY. This script has no
write or publication mode. The accompanying fill script was used once to put
the checked document in the draft.

This is a reproducible transcription and uncertainty conversion, not a new
Monte Carlo simulation. Jensen, J. Stat. Phys. 102 (2001), 865-881, supplies
4.062570(8). Hsu, Nadler and Grassberger, J. Phys. A 38 (2005), 775-806,
Table 1 supplies log(gamma)=2.1218588(25), partially constrained fit; section
5.2 supplies 5.20789(4) and 10.61539(6) for bond animals counted by edges.
The source estimates are heuristic. Arb only certifies the conversion of
that assumed input interval, never the unknown limiting constant itself.

For the one exponentiation we use 192 bits, then repeat at 384 bits in the
checks. The decimal midpoint is rounded to six places and the radius upwards
to six places, enclosing the complete exponentiated source interval. The
radius, not the midpoint's printed length, states the supported precision.
The three direct transcriptions keep the uncertainty written in the paper.
"""
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_BITS = 192
PARAMETERS = ('class', 'lattice', 'd', 'size', 'labelling')


def decimal_rational(text):
    """Exact conversion, also on modular Sage lacking QQ(decimal_string)."""
    text = str(text)
    if '.' not in text:
        return QQ(int(text))
    whole, fraction = text.split('.')
    sign = -1 if text.startswith('-') else 1
    return QQ(int(whole)) + sign * QQ(int(fraction)) / QQ(10**len(fraction))


def fixed(value, places):
    """Write an integer multiple of 10^-places without float conversion."""
    scale = 10**places
    n = int(value * scale)
    assert QQ(n) / QQ(scale) == value
    sign = '-' if n < 0 else ''
    n = abs(n)
    return '%s%d.%0*d' % (sign, n // scale, places, n % scale)


def cubic_site_text(bits=WORKING_BITS):
    field = RealBallField(bits)
    a = decimal_rational('2.1218588')
    error = decimal_rational('0.0000025')
    lower = field(a-error).exp()
    upper = field(a+error).exp()
    center = field(a).exp()
    scale = QQ(10**6)
    rounded = (center * scale + field(QQ(1)/2)).floor()
    midpoint = QQ(rounded) / scale
    # Use outward Arb endpoints, then outward decimal rounding.
    gap = max(field(midpoint)-lower, upper-field(midpoint))
    radius = QQ((gap.upper()*scale).ceil()) / scale
    assert lower > field(midpoint-radius)
    assert upper < field(midpoint+radius)
    return fixed(midpoint,6) + ' +/- ' + fixed(radius,6)


class AnimalGrowth(numberdb.Generator):
    table = 'T484'
    parameters = PARAMETERS
    type = 'R'
    digits = 5
    rigour = 'heuristic'
    files = ('generate.py', 'checks.py', 'source-data.json', 'evidence.md')

    def enumerate(self):
        for animal_class, size in [('site','cells'), ('bond','edges')]:
            for lattice, dimension in [('square','2'), ('sc','3')]:
                yield dict(zip(PARAMETERS,
                    (animal_class,lattice,dimension,size,'unlabelled')))

    def value(self, params, digits):
        assert params in list(self.enumerate()), 'outside this generator range'
        key = params['class'], params['lattice']
        if key == ('site','square'):
            return {'number':'4.062570 +/- 0.000008', 'comment':
                r"Klarner's constant. The series estimate is $4.062570(8)$ CITE{Jensen}. The rigorous bounds $4.00253<\gamma\leq4.5238$ are due to CITE{BRS} and CITE{Bui}, respectively."}
        if key == ('site','sc'):
            return {'number':cubic_site_text(), 'comment':
                r'The growth constant of fixed polycubes. Table 1 of CITE{HNG} gives $\log\gamma=2.1218588(25)$ for the partially constrained fit. The number shown encloses its exponentiated uncertainty interval, with outward decimal rounding. A rigorous upper bound is $\gamma\leq9.3835$ CITE{BS}.'}
        if key == ('bond','square'):
            return {'number':'5.20789 +/- 0.00004', 'comment':
                r'The Monte Carlo estimate is $5.20789(4)$ in Section 5.2 of CITE{HNG}, with size equal to the number of chosen edges.'}
        return {'number':'10.61539 +/- 0.00006', 'comment':
                r'The Monte Carlo estimate is $10.61539(6)$ in Section 5.2 of CITE{HNG}, with size equal to the number of chosen edges.'}


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    result = AnimalGrowth().verify(sample=None)
    print(result)
    sys.exit(0 if result.ok else 1)