#!/usr/bin/env python3
"""Independent deterministic Laplace-integral check; no Monte Carlo estimates."""
import argparse
import json
from pathlib import Path
import subprocess
import mpmath as mp

parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--model', type=Path, default=Path(__file__).resolve().parent.parent / 'assets' / 'sensor-models.js', help='Path to the downloaded or repository sensor-models.js')
parser.add_argument('--output', type=Path, help='Optional full JSON report; the summary is always printed')
args = parser.parse_args()
model = args.model.resolve()
if not model.is_file():
    parser.error('The production model was not found; pass --model /path/to/sensor-models.js')
mp.mp.dps = 65
cases = sorted(set(
    [(n + gap, n, v) for n in [1, 2, 5, 20] for gap in [2, 3, 5] for v in [.5, 4, 16]]
    + [(p, n, v) for p in [100, 1000, 10000] for n, v in [(1, 4), (20, .5), (20, 4)]]
))
js = "const M=require(process.argv[1]);const c=" + json.dumps(cases) + ";console.log(JSON.stringify(c.map(([p,n,v])=>M.interpolatorAnalyticBenchmark(p,n,v))));"
computed = json.loads(subprocess.check_output(['node', '-e', js, str(model)], text=True))
records = []
for (p, n, variance), result in zip(cases, computed):
    v = mp.mpf(str(variance)); mean = mp.mpf(n) * p / (2 * v)
    # Rescale the Laplace integration variable: narrow endpoint peaks at large
    # p or noncentrality cannot be missed by ordinary quadrature on [0,1].
    knots = sorted(set([mp.mpf(0), mean] + [mp.mpf(t) for t in [1, 4, 16, 64, 256, 1024] if t < mean]))
    h = mp.quad(lambda t: (1-t/mean)**(mp.mpf(p)/2-2)*mp.exp(-t), knots) / (2*mean)
    k = p-n+1
    expected = mp.mpf(p-2)/(2*(k-2))*((k-3)+(2*(p-1)-(k-3)*(p-2+2*mean))*h)
    difference = abs(mp.mpf(result['test']) - expected)
    assert difference <= result['errorBound'], (p, n, variance, difference, result)
    assert result['errorBound'] < 1e-8
    records.append(dict(p=p, n=n, variance=variance, reference=mp.nstr(expected, 30),
                        javascript=result['test'], absolute_difference=float(difference),
                        reported_error_budget=result['errorBound']))

out = dict(method='65-digit Laplace integration, independently evaluated from the bounded Poisson implementation',
           cases=len(records), maximum_difference=max(r['absolute_difference'] for r in records),
           maximum_error_budget=max(r['reported_error_budget'] for r in records), results=records)
if args.output:
    args.output.write_text(json.dumps(out, indent=2)+'\n')
print(json.dumps({k:v for k,v in out.items() if k!='results'}, indent=2))