"""MP1 · paso 30: métricas (implementación 1) y aceptación pre-registrada (§7, §12) de un desenlace ou|des|mod.
Lee /w/interno/mp1_<h>/predicciones.csv (+ /w/interno/b2_<h>.csv para B2) y /w/salidas/20_<h>.json.
Escribe /w/salidas/30_<h>.json y /w/interno/mp1_<h>/eval_prueba.csv (las filas evaluadas, para la segunda implementación).
"""
import sys, json, os
import numpy as np, pandas as pd
import statsmodels.api as sm
from scipy.special import logit
from sklearn.metrics import average_precision_score

H = sys.argv[1]
SEM, B = 20261026, 2000
D = f"/w/interno/mp1_{H}"
r = pd.read_csv(f"{D}/predicciones.csv", dtype={"licitacion_id": "string", "lote_u": "string", "organo_clave": "string"},
                parse_dates=["fecha_pub"])
inf20 = json.load(open(f"/w/salidas/20_{H}.json"))
if os.path.exists(f"/w/interno/b2_{H}.csv"):
    b2 = pd.read_csv(f"/w/interno/b2_{H}.csv", dtype={"licitacion_id": "string", "lote_u": "string"})
    r = r.merge(b2, on=["licitacion_id", "lote_u"], how="left")
else:
    r["p_b2"] = np.nan

def grupos_score(s):
    u, g = np.unique(s, return_inverse=True)
    return g, len(u)

def auc_w(g, G, y, w):
    wp = np.bincount(g, weights=w * y, minlength=G); wn = np.bincount(g, weights=w * (1 - y), minlength=G)
    below = np.cumsum(wn) - wn
    den = wp.sum() * wn.sum()
    return float(np.sum(wp * (below + 0.5 * wn)) / den) if den > 0 else np.nan

def bootstrap(df, cols):
    """AUC puntual e IC 95 % por bootstrap de órganos (los lotes sin órgano, cada uno su grupo) de cada columna, y las
    réplicas para las diferencias emparejadas."""
    y = df["y"].to_numpy(float)
    cl = df["organo_clave"].astype(object).where(df["organo_clave"].notna(), "__lote__" + df.index.astype(str))
    cid, C = pd.factorize(cl)[0], cl.nunique()
    gs = {c: grupos_score(df[c].to_numpy(float)) for c in cols}
    w1 = np.ones(len(df))
    punto = {c: auc_w(*gs[c], y, w1) for c in cols}
    rng = np.random.default_rng(SEM)
    rep = {c: np.empty(B) for c in cols}
    for b in range(B):
        m = np.bincount(rng.integers(0, C, C), minlength=C).astype(float)
        w = m[cid]
        for c in cols:
            rep[c][b] = auc_w(*gs[c], y, w)
    ic = {c: [float(np.percentile(rep[c], 2.5)), float(np.percentile(rep[c], 97.5))] for c in cols}
    return punto, ic, rep

def deciles(p, y):
    q = pd.qcut(p, 10, duplicates="drop")
    t = pd.DataFrame({"p": p, "y": y, "q": q}).groupby("q", observed=True).agg(n=("y", "size"), predicha=("p", "mean"),
                                                                                 observada=("y", "mean"), E=("p", "sum"),
                                                                                 O=("y", "sum"),
                                                                                 V=("p", lambda s: float(np.sum(s * (1 - s)))))
    t["dif_puntos"] = 100 * (t["predicha"] - t["observada"])
    return t.reset_index(drop=True)

def pendiente(p, y):
    X = sm.add_constant(logit(np.clip(p, 1e-6, 1 - 1e-6)))
    f = sm.GLM(y, X, family=sm.families.Binomial()).fit()
    return float(f.params[1]), float(f.params[0])

def metricas(df, base_rate):
    y = df["y"].to_numpy(float); p = df["p"].to_numpy(float)
    pc = np.clip(p, 1e-6, 1 - 1e-6)
    brier = float(np.mean((p - y) ** 2)); brier0 = float(np.mean((base_rate - y) ** 2))
    out = {"n": int(len(df)), "positivos": int(y.sum()), "tasa": float(y.mean()), "brier": brier,
           "brier_skill": float(1 - brier / brier0), "log_loss": float(-np.mean(y * np.log(pc) + (1 - y) * np.log(1 - pc))),
           "ap": float(average_precision_score(y, p)), "cal_global_puntos": float(100 * (p.mean() - y.mean())),
           "o_e": float(y.sum() / p.sum())}
    out["pendiente"], out["intercepto_cal"] = pendiente(p, y)
    dec = deciles(p, y); out["deciles"] = dec.to_dict(orient="records")
    return out, dec

rc = r[r["ventana"] == "recal"].copy(); te = r[r["ventana"] == "test"].copy().reset_index(drop=True)
base_rc = float(rc["y"].mean())
inf = {"desenlace": H, "pre_registro_md5": "3da2be07b085b96c3f80dd0a826748bd", "eleccion": inf20["eleccion"],
       "calibrador": inf20["calibrador"], "avisados": {"n": 0, "nota": "ningún anuncio avisado todavía: con y sin avisados es lo mismo"}}
m_te, dec_te = metricas(te, base_rc)
m_rc, dec_rc = metricas(rc, base_rc)
inf["prueba"] = m_te; inf["recal"] = m_rc
cols = ["p", "p_raw", "b1"] + (["b3"] if H in ("mod", "ser") else [])
pt, ic, rep = bootstrap(te, cols)
inf["prueba"]["auc"] = pt; inf["prueba"]["auc_ic"] = ic
dif = {}
for c in [x for x in cols if x not in ("p", "p_raw")]:
    dd = rep["p"] - rep[c]
    dif[f"p-{c}"] = {"delta": pt["p"] - pt[c], "ic": [float(np.percentile(dd, 2.5)), float(np.percentile(dd, 97.5))]}
# B2 (E1/E1c) en los lotes de la prueba que tienen su predicción
if H in ("ou", "des"):
    tb = te[te["p_b2"].notna()].reset_index(drop=True)
    inf["prueba"]["b2_cobertura"] = {"con_b2": int(len(tb)), "sin_b2": int(len(te) - len(tb))}
    if len(tb):
        pt2, ic2, rep2 = bootstrap(tb, ["p", "p_b2"])
        dd = rep2["p"] - rep2["p_b2"]
        dif["p-b2"] = {"auc_p": pt2["p"], "auc_b2": pt2["p_b2"], "delta": pt2["p"] - pt2["p_b2"],
                       "ic": [float(np.percentile(dd, 2.5)), float(np.percentile(dd, 97.5))],
                       "brier_p": float(np.mean((tb["p"] - tb["y"]) ** 2)), "brier_b2": float(np.mean((tb["p_b2"] - tb["y"]) ** 2))}
inf["prueba"]["diferencias"] = dif
# AUC de la ventana de recalibración (Anexo B, «año anterior»)
ptr, icr, _ = bootstrap(rc.reset_index(drop=True), ["p", "b1"])
inf["recal"]["auc"] = ptr; inf["recal"]["auc_ic"] = icr
inf["intervalos"] = inf20.get("intervalos")
inf["deriva"] = {"alarmas_psi_ge_025": inf20["deriva"]["alarmas_psi_ge_025"], "psi_p_bruta": inf20["deriva"]["psi"].get("p_bruta"),
                 "calibracion_por_mes": inf20["deriva"]["calibracion_por_mes"]}

# --- aceptación (§12) -------------------------------------------------------------------------------------------------
A = {}
cov = (inf["intervalos"] or {}).get("cobertura_test")
if H in ("ou", "des", "ser"):
    d9 = [abs(x["dif_puntos"]) for x in m_te["deciles"][:-1]]; d10 = abs(m_te["deciles"][-1]["dif_puntos"])
    A["1_auc"] = pt["p"] >= 0.68 and ptr["p"] >= 0.66
    A["2_b1"] = dif["p-b1"]["delta"] >= 0.015 and dif["p-b1"]["ic"][0] > 0
    if H == "ser":
        A["3_b3"] = dif["p-b3"]["delta"] >= 0
    else:
        A["3_b2"] = ("p-b2" in dif) and dif["p-b2"]["delta"] >= 0 and dif["p-b2"]["brier_p"] <= dif["p-b2"]["brier_b2"]
    A["4_calibracion"] = (max(d9) <= 5 and d10 <= 8 and abs(m_te["cal_global_puntos"]) <= 2 and 0.85 <= m_te["pendiente"] <= 1.15)
    A["5_cobertura"] = cov is not None and 0.85 <= cov <= 0.95
    A["detalle_calibracion"] = {"max_dif_deciles_1_9": max(d9), "dif_decil_10": d10, "global": m_te["cal_global_puntos"],
                                "pendiente": m_te["pendiente"]}
    duros, blandos = ["1_auc", "2_b1", "4_calibracion"], (["3_b3"] if H == "ser" else ["3_b2", "5_cobertura"])
    potencia = (m_te["positivos"] >= 300) if H == "ser" else (m_te["n"] >= 5000)
else:
    tol = [abs(x["O"] - x["E"]) <= 2.58 * np.sqrt(x["V"]) + 0.10 * x["E"] for x in m_te["deciles"]]
    A["1_auc"] = pt["p"] >= 0.75
    A["2_ap"] = m_te["ap"] >= 3 * m_te["tasa"]
    A["3_b1_b3"] = (dif["p-b1"]["delta"] >= 0.03 and dif["p-b1"]["ic"][0] > 0 and dif["p-b3"]["delta"] >= 0.03 and dif["p-b3"]["ic"][0] > 0)
    A["4_calibracion"] = (0.8 <= m_te["o_e"] <= 1.25 and 0.8 <= m_te["pendiente"] <= 1.2 and all(tol))
    A["5_cobertura"] = cov is not None and 0.85 <= cov <= 0.95
    A["detalle_calibracion"] = {"o_e": m_te["o_e"], "pendiente": m_te["pendiente"], "deciles_dentro": int(sum(tol)), "deciles": len(tol)}
    duros, blandos = ["1_auc", "3_b1_b3", "4_calibracion"], ["2_ap", "5_cobertura"]
    potencia = m_te["positivos"] >= 100
A["pendientes_fuera_de_este_script"] = ["6_ex_ante (test dinámico)", "7_segunda_implementacion", "8_cotejo"]
if not potencia:
    ver = "NO CONFIRMADA (sin potencia)"
elif not all(A[k] for k in duros):
    ver = "NO PASA"
elif not all(A[k] for k in blandos):
    ver = "NO CONFIRMADA (pasan los criterios duros; falla " + ", ".join(k for k in blandos if not A[k]) + ")"
else:
    ver = "PASA a falta de 6-8"
A["veredicto_provisional"] = ver
inf["aceptacion"] = A
te[["licitacion_id", "lote_u", "organo_clave"] + (["cpv_cel"] if "cpv_cel" in te else []) + ["y", "p", "p_raw", "b1"] + (["b3"] if H in ("mod", "ser") else []) + ["p_b2"]].to_csv(
    f"{D}/eval_prueba.csv", index=False)
json.dump(inf, open(f"/w/salidas/30_{H}.json", "w"), indent=1, ensure_ascii=False, default=float)
print(json.dumps({"auc": pt, "ic": ic, "dif": dif, "recal_auc": ptr, "cal": A["detalle_calibracion"], "cobertura": cov,
                  "veredicto": ver}, ensure_ascii=False, default=float, indent=1))