"""Monte-Carlo identifiability stress test for the coupling estimator (honest reporting, no tuning to truth)."""
import sys, json, numpy as np
sys.path.insert(0, __file__.rsplit('/', 1)[0])
import pranapd as P
truth = {k: .8 * v for k, v in P.PAIRS.items()}
out = {}
for n in (60, 200):
    ests = {k: [] for k in P.PAIRS}
    for seed in range(1, 9):
        c, e = P.synth_cohort(n=n, seed=seed)
        est = P.estimate_coupling(c, e)
        for k in P.PAIRS: ests[k].append(est[k])
    out[n] = {f"{P.DOMAINS[i]}-{P.DOMAINS[j]}": dict(true=round(truth[(i, j)], 3), mean=round(float(np.mean(v)), 3),
              sd=round(float(np.std(v)), 3), bias=round(float(np.mean(v)) - truth[(i, j)], 3),
              sign_ok=int(np.sum(np.array(v) > 0))) for (i, j), v in ests.items()}
print(json.dumps(out, indent=1))
json.dump(out, open('/home/user/webapp/python/identifiability_results.json', 'w'), indent=1)