"""Independent checks for numberdb.org/T474.

Run with: sage -python generate.py --check
This checks the transcription and published independent estimates, not a
new bootstrap or Monte Carlo calculation. All interval endpoints are exact
rational numbers, so no float conversion enters an uncertainty comparison.
"""
import itertools
import re
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from check import agrees_with, exactness, measure, prose

# Transcribed from the source's parenthesised uncertainties, separately from
# the explicit centre/radius strings in generate.py.
PUBLISHED = {
    ('sigma', 'sigma', 'epsilon'): '1.05185373(11)',
    ('epsilon', 'epsilon', 'epsilon'): '1.53244304(58)',
}
MONTE_CARLO = {
    ('sigma', 'sigma', 'epsilon'): '1.051(1)',
    ('epsilon', 'epsilon', 'epsilon'): '1.533(5)',
}
PREVIOUS_BOOTSTRAP = {
    ('sigma', 'sigma', 'epsilon'): '1.0518537(41)',
    ('epsilon', 'epsilon', 'epsilon'): '1.532435(19)',
}


def parenthesised(text):
    match = re.fullmatch(r'(-?\d+\.\d+)\((\d+)\)', text)
    assert match, text
    centre = Fraction(match[1])
    radius = Fraction(int(match[2]), 10 ** len(match[1].split('.')[1]))
    return centre - radius, centre + radius


def interval(text):
    centre, radius = (Fraction(x.strip()) for x in text.split('+/-'))
    assert radius > 0
    return centre - radius, centre + radius


def overlaps(a, b):
    return max(a[0], b[0]) <= min(a[1], b[1])


def contains(outer, inner):
    return outer[0] <= inner[0] <= inner[1] <= outer[1]


def run_checks(entries):
    # Controls must return known nonzero answers, and reject a non-overlap.
    assert parenthesised('1.234(56)') == (Fraction(589, 500), Fraction(129, 100))
    assert interval('2 +/- 0.5') == (Fraction(3, 2), Fraction(5, 2))
    assert overlaps(interval('2 +/- 0.5'), interval('3 +/- 0.5'))
    assert not overlaps(interval('2 +/- 0.1'), interval('3 +/- 0.1'))
    assert not contains(interval('2 +/- 0.1'), interval('2 +/- 0.2'))
    # A shifted centre and an incorrectly parsed radius must both fail.
    assert interval('1.05185374 +/- 0.00000011') != parenthesised('1.05185373(11)')
    assert interval('1.05185373 +/- 0.000000011') != parenthesised('1.05185373(11)')
    assert not exactness(entries), exactness(entries)
    assert not prose(entries), prose(entries)
    assert set(entries) == set(PUBLISHED)
    intervals = {k: interval(v['number']) for k, v in entries.items()}
    assert not agrees_with(intervals, lambda k: parenthesised(PUBLISHED[k]))

    parity = {'sigma': -1, 'epsilon': 1}
    allowed = set()
    forbidden = set()
    for triple in itertools.combinations_with_replacement(('sigma', 'epsilon'), 3):
        if parity[triple[0]] * parity[triple[1]] * parity[triple[2]] == 1:
            allowed.add(triple)
        else:
            forbidden.add(triple)
    assert allowed == set(entries)
    assert forbidden == {('sigma', 'sigma', 'sigma'), ('sigma', 'epsilon', 'epsilon')}

    # Check the coordinate exponents as identities, not rounded evaluations:
    # total dilation weight and each insertion's weight under inversion.
    ring = PolynomialRing(QQ, names=('d_i', 'd_j', 'd_k'))
    dimensions = ring.gens()
    for _key in entries:
        a, b, c = dimensions
        exponents = (a + b - c, b + c - a, c + a - b)
        assert sum(exponents) == a + b + c
        assert exponents[0] + exponents[2] == 2 * a
        assert exponents[0] + exponents[1] == 2 * b
        assert exponents[1] + exponents[2] == 2 * c
        assert sum((exponents[0] + 1, *exponents[1:])) != a + b + c
        for perm in itertools.permutations(range(3)):
            remapped = {(min(perm[i], perm[j]), max(perm[i], perm[j])):
                        dimensions[perm[i]] + dimensions[perm[j]] - dimensions[perm[k]]
                        for i, j, k in [(0, 1, 2), (1, 2, 0), (2, 0, 1)]}
            assert remapped == {(0, 1): exponents[0], (1, 2): exponents[1],
                                (0, 2): exponents[2]}

    comparisons = {}
    for key, box in intervals.items():
        assert 0 < box[0] < box[1]
        assert box[1] - box[0] < Fraction(1, 100000)
        for perm in set(itertools.permutations(key)):
            assert tuple(sorted(perm, key=lambda x: ('sigma', 'epsilon').index(x))) == key
        # Reversing sigma's sign preserves either stored coefficient.
        assert (-1) ** key.count('sigma') == 1
        # Reversing epsilon changes either sign, so its phase must be fixed.
        assert (-1) ** key.count('epsilon') == -1
        mc = parenthesised(MONTE_CARLO[key])
        old = parenthesised(PREVIOUS_BOOTSTRAP[key])
        assert contains(mc, box), (key, 'Monte Carlo discrepancy')
        assert contains(old, box), (key, '2016 bootstrap discrepancy')
        centre = sum(box, Fraction(0)) / 2
        mc_centre = sum(mc, Fraction(0)) / 2
        mc_radius = (mc[1] - mc[0]) / 2
        comparisons[','.join(key)] = {
            'source_endpoints': [str(x) for x in box],
            'inside_Monte_Carlo_error_interval': True,
            'inside_2016_bootstrap_error_interval': True,
            'centre_difference_in_MC_error_units': str(abs(centre - mc_centre) / mc_radius),
        }
    return {'entries': len(entries), 'controls': 'passed',
            'source_uncertainty_transcription': 'exact agreement',
            'formula_checks': 'dilation, inversion, all insertion permutations, spin-flip selection and phase changes',
            'measurement': measure(entries), 'comparisons': comparisons,
            'limitation': 'No conformal-bootstrap optimization or Monte Carlo simulation is rerun.'}


def from_document(document):
    numbers = document['Numbers']
    return {(i, j, k): leaf if isinstance(leaf, dict) else {'number': leaf}
            for i, second in numbers.items()
            for j, third in second.items()
            for k, leaf in third.items()}