"""Independent finite-source science replay; Python only, no interactive Engine."""
from pathlib import Path
import argparse
import json
import subprocess
import tempfile
import numpy as np
from numpy.polynomial.legendre import leggauss
from binary_reference import finite_source_image_area


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--materials', type=Path, default=Path(__file__).resolve().parent)
    parser.add_argument('--output', type=Path, default=Path('science-summary.json'))
    arguments = parser.parse_args()
    folder = arguments.materials.resolve()
    states = json.loads((folder / 'synthetic-domain-reference.json').read_text())['cases']
    results = []
    for state in states:
        value = finite_source_image_area(complex(*state['source']), state['rho'], state['s'], state['q'], limb_gamma=state['gamma'],
                                         use_curvature_bounds=True, relative_width=.01, maximum_cells=5000000)
        difference = abs(value['estimate'] / state['independent_contour'] - 1)
        assert value['status'] == 'bounded' and difference < .001
        assert value['lower'] <= state['independent_contour'] <= value['upper']
        results.append(difference)
    for gamma in [0., .538, 1.]:
        nodes, weights = leggauss(64); mu = (nodes + 1) / 2
        annulus_weights = 1.5 * gamma * weights / 2 * (1 - mu ** 2)
        assert np.min(annulus_weights) >= 0 and 1 - gamma >= 0
        assert abs(1 - gamma + np.sum(annulus_weights) - 1) < 2e-14
    # Independently compiled uniform contours, combined with positive annulus weights.
    course = json.loads((folder / 'synthetic-source-size-curves.json').read_text())
    with tempfile.TemporaryDirectory(prefix='ac23-uniform-reference-') as temporary:
        executable = Path(temporary) / 'reference'
        subprocess.run(['c++', '-O3', '-std=c++11', '-I' + str(folder / 'reference'), str(folder / 'reference_runner.cpp'),
                        str(folder / 'reference/VBBinaryLensingLibrary.cpp'), '-o', str(executable)], check=True)
        examples = []
        for row in course['rows']:
            parameters = row['parameters']; radius = row['sourceRadius']
            point = min(row['points'], key=lambda item: abs(item['time'] - .02))
            source = np.exp(1j * parameters['trajectoryAngleRadians']) * (point['time'] + 1j * parameters['impactParameter'])
            estimates = []
            for order in [64, 128]:
                nodes, weights = leggauss(order); mu = (nodes + 1) / 2
                radii = np.r_[radius, radius * np.sqrt(1 - mu ** 2)]
                lines = ''.join(f"{parameters['separation']} {parameters['massRatio']} {source.real} {source.imag} {rho} 0\n" for rho in radii)
                values = np.fromstring(subprocess.check_output([str(executable), '1e-7'], input=lines, text=True, timeout=90), sep=' ')
                gamma = parameters['limbGamma']
                estimates.append(float((1 - gamma) * values[0] + 1.5 * gamma * np.sum(weights / 2 * (1 - mu ** 2) * values[1:])))
            assert abs(estimates[-1] / estimates[-2] - 1) < .0001
            assert abs(estimates[-1] / point['magnification'] - 1) < .001
            examples.append({'radius': radius, 'native': point['magnification'], 'independent_positive_annulus': estimates[-1], 'last_relative_refinement': abs(estimates[-1] / estimates[-2] - 1)})
    report = {'status': 'PASS independent science replay', 'synthetic_states': len(states), 'maximum_relative_flux_difference': max(results),
              'source_examples': examples, 'qualification': 'Float64 geometric envelopes, not formal outward-rounded certificates; finite probes do not establish a continuum-wide bound.'}
    arguments.output.write_text(json.dumps(report, indent=2) + '\n')
    print(json.dumps(report, indent=2))


if __name__ == '__main__': main()