"""Re-run the certificates of Table 2 of the paper from the saved profiles.

usage:  python3 reproduce.py [m ...]          (default: 2 3 4 5)

For each m this prints one JSON line with
  Q_lower : rigorous lower bound for the quotient Q(F_psi) of F_psi = psi(t_1+...+t_k) prod g(t_i)  (hcert.py)
  rc_loss : rigorous upper bound for the loss from restricting F_psi to R_k(u0, kappa)            (rc_loss_rig*.py)
  target  : m / (1/4 + omega),  omega = 0.012466, delta = 1e-4   (240 omega + 80 delta = 2.99984 < 3)
  net     : Q_lower - rc_loss - target     (net > 0 certifies M_k^R > m / (1/4 + omega))
Requirements: python3, numpy, scipy, python-flint (Arb; we used python-flint 0.9.0 with FLINT 3.6.0,
numpy 2.5.3, scipy 1.18.1, Python 3.13).  Running time: about 5, 10, 11 and 10 minutes for m = 2, 3, 4, 5 on one
core.  Our output is in results.json and ../logs/reproduce.log."""
import os, sys, json, time
import numpy as np
HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)
from trial_io import load_trial
from hcert import hcert
from fcore import CDF
from rc_loss_rig import rc_loss_rig
from rc_loss_rig2 import rc_loss_rig2
OMEGA, DELTA = 0.012466, 1e-4
TRIALS = {2: "A_m2_k2355", 3: "A_m3_k111300", 4: "A_m4_k5129500", 5: "A_m5_k234360000"}

def run(m):
    t0 = time.time()
    d = load_trial(os.path.join(HERE, "trials", TRIALS[m] + ".json"))
    k, gamma, beta, dl = int(d["k"]), float(d["theta"]), float(d["beta"]), float(d["delta"])
    cdf = CDF(k, gamma, beta, ns=(k - 1, k))
    c = hcert(k, gamma, beta, d["E"], d["h"], cdf=cdf, verbose=False, DELTA=dl)
    f = rc_loss_rig2 if m == 2 else rc_loss_rig      # first / second inequality of Lemma 4.3
    rc = f(k, gamma, beta, DELTA, OMEGA, c["den_lower"], cv=1.5)
    target = m / (0.25 + OMEGA)
    out = dict(m=m, k=k, gamma=gamma, beta=beta, Delta=dl, Xi1=c["Xi1"], Q_lower=c["Q_lower"],
               E_psi2_lower=c["den_lower"], E_psi2_upper=c["den_upper"], max_cdf_width=c["max_cdf_width"],
               rc_loss=rc["loss_upper"], rc_method=f.__name__, target=target,
               net=c["Q_lower"] - rc["loss_upper"] - target, secs=round(time.time() - t0, 1))
    print(json.dumps(out), flush=True)
    return out

if __name__ == "__main__":
    for m in ([int(a) for a in sys.argv[1:]] or [2, 3, 4, 5]):
        run(m)