model-3/scripts/fr_robust.py

"""Risk analysis for raising France's threshold.

The p**g shift model saturated at g=6 without fully explaining France's deficit, so g is
NOT identified. The decision must therefore be robust across the whole plausible range of
g - including g=1 (no shift at all), which is the downside case. Reports, for each g:
the optimal France threshold, and what threshold 0.90 costs or gains against 0.70.

Also propagates the uncertainty in the val US/India estimate into France's implied score,
since France is inferred by subtraction and the lever arm is 0.85/0.15 = 5.7x.
"""
import sys, glob
import numpy as np, pandas as pd, joblib
sys.path.insert(0, "src")
import config as C, decide as DC

b = joblib.load("work/results_A/model_bundle.joblib")
MM, THR = b["miss_mass"], b["param"]
s1 = pd.read_parquet(C.PATHS["test"]["source1"], columns=["entity_id", "country_norm"])
ctry = s1["country_norm"].astype(str).values
N = len(s1)
nq = {c: int((ctry == c).sum()) for c in np.unique(ctry)}
wfr = nq["france"] / N

d = pd.concat([pd.read_parquet(f) for f in
               sorted(glob.glob("work/results_A/test_chunks/chunk_*.parquet"))],
              ignore_index=True)
d = DC.resolve_conflicts(d, "p")
d["ctry"] = ctry[d["q"].values]

fr = d[d.ctry == "france"][["q", "p"]].sort_values("q", kind="stable")
fq, fp = fr["q"].values, fr["p"].values.astype(np.float64)
uq, starts = np.unique(fq, return_index=True)
idx = np.r_[starts, len(fq)]
n_fr, n_nocand = nq["france"], nq["france"] - len(uq)

def f05(g, t):
    pg = fp ** g
    sel = fp >= t
    tot = float(n_nocand)
    for i in range(len(uq)):
        a, z = idx[i], idx[i + 1]
        m = sel[a:z]
        k = int(m.sum())
        tp = float(pg[a:z][m].sum())
        T = tp + float(pg[a:z][~m].sum()) + MM
        tot += (np.exp(-T) if T > 0 else 1.0) if k == 0 else 1.25 * tp / (0.25 * T + k)
    return tot / n_fr

print("=" * 78)
print("ROBUSTNESS: France threshold choice across the whole plausible shift range")
print("=" * 78)
grid = [0.70, 0.75, 0.80, 0.85, 0.88, 0.90, 0.92, 0.95]
print(f"{'g':>5}{'F@0.70':>10}" + "".join(f"{'d@'+format(t,'.2f'):>10}" for t in grid[1:])
      + f"{'  best thr':>11}")
for g in [1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0]:
    base = f05(g, 0.70)
    vals = [(t, f05(g, t)) for t in grid]
    bt, bv = max(vals, key=lambda x: x[1])
    print(f"{g:>5.1f}{base:>10.5f}" + "".join(f"{v-base:>+10.5f}" for _, v in vals[1:])
          + f"{bt:>11.2f}")

print("\n" + "=" * 78)
print("DOWNSIDE AT g=1 (calibration transfers perfectly - the pessimistic case for")
print("this change) vs UPSIDE at g>=3, expressed on the LEADERBOARD")
print("=" * 78)
for g in [1.0, 2.0, 3.0, 4.0, 6.0]:
    base = f05(g, 0.70)
    for t in (0.85, 0.90):
        print(f"  g={g:<4} thr {t:.2f}: France {f05(g,t)-base:+.6f}  "
              f"-> leaderboard {wfr*(f05(g,t)-base):+.6f}")

print("\n" + "=" * 78)
print("HOW MUCH OF FRANCE'S DEFICIT CAN MISCALIBRATION EXPLAIN AT ALL?")
print("=" * 78)
print(f"  at threshold 0.70:  g=1 -> {f05(1.0,0.70):.6f}   g=6 -> {f05(6.0,0.70):.6f}")
print("  implied from leaderboard: 0.946062")
print("  => even an extreme shift leaves a residual, so part of France's deficit is")
print("     NOT the operating point: it is true pairs never retrieved or never scored.")

print("\n" + "=" * 78)
print("SENSITIVITY of France's implied score to the val US/India estimate")
print("=" * 78)
wus, win = nq["us"] / N, nq["india"] / N
for dv in (-0.003, -0.002, -0.001, 0.0, 0.001, 0.002, 0.003):
    usin = wus * (0.9832 + dv) + win * (0.9806 + dv)
    print(f"  val US/India {0.98177+dv:.5f}  ->  France {(0.976423-usin)/wfr:.5f}")