Wisconsin Diagnostic Breast Cancer
Benchmark & Calibration

EXPECTED_FILES = {
    "erl02_core": ("results/erl02/core_initial/benchmark_dev_core.json", "118307ed3880f3617947c898a796194466936b83c752ae874971be93b564e004"),
    "erl02_expanded": ("results/erl02/expanded_run01/benchmark_dev_expanded.json", "a4253e67f27c383211ba0df401202386ede305ab7f104bbe02cf4d4016f45212"),
    "erl03b_report": ("results/erl03b/expanded_oof_run01/benchmark_dev_expanded.json", "a4253e67f27c383211ba0df401202386ede305ab7f104bbe02cf4d4016f45212"),
    "erl03b_oof": ("results/erl03b/expanded_oof_run01/outer_oof_predictions_v1.csv", "fa2bb63447e0ee21c7796e3f19e4fbb3a160dec6c581d634b872671b3316f894"),
    "erl03b_manifest": ("results/erl03b/expanded_oof_run01/run_provenance.json", "88352b4f4e6ac412b2c2bd950631cc6b7fba8b2e92250befb5ed46dcc15a113a"),
    "erl03c_report": ("results/erl03c/calibration_run01/calibration_comparison_v0_3_1.json", "b62cda44133f998ebb705cd3ac7c87a17068f581bf29a815f4cb435d450708da"),
    "erl03c_oof": ("results/erl03c/calibration_run01/outer_oof_calibration_v1_1.csv", "ec25daa0253194d0bc2c0ff78e877a19e68b2c2e0399ce7982aed21e02cc985d"),
    "erl03c_manifest": ("results/erl03c/calibration_run01/run_provenance.json", "16c1ed341cec287a848f48641ff8629af553ec7c37126beeb96c188b25cdc1a9"),
}


def find_work_root() -> Path:
    override = os.environ.get("HEALTHCARE_MLOPS_WORK", "").strip()
    candidates = [Path(override)] if override else []
    for start in (Path.cwd(), Path.cwd().parent):
        candidates.extend([start, *start.parents])
    for candidate in candidates:
        if candidate.is_dir() and (candidate / "results" / "erl03c" / "calibration_run01").is_dir():
            return candidate.resolve()
    raise FileNotFoundError(
        "No se encontró healthcare-mlops-benchmark-work. "
        "Configure HEALTHCARE_MLOPS_WORK con la ruta de trabajo original; "
        "no copie ni vuelva a crear resultados congelados."
    )


def sha256_file(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as stream:
        for block in iter(lambda: stream.read(1024 * 1024), b""):
            digest.update(block)
    return digest.hexdigest()


WORK_ROOT = find_work_root()
VERIFIED = {}
for name, (relative, expected_hash) in EXPECTED_FILES.items():
    source = (WORK_ROOT / relative).resolve(strict=True)
    results_root = (WORK_ROOT / "results").resolve(strict=True)
    if results_root not in source.parents:
        raise RuntimeError(f"Artefacto fuera de results: {name}")
    observed_hash = sha256_file(source)
    if observed_hash != expected_hash:
        raise RuntimeError(f"Integridad SHA256 incumplida: {name}")
    VERIFIED[name] = source

assert len(VERIFIED) == 8
print("ERL GATE 01 · SHA256 8/8 PASS · Modo de lectura exclusivamente")

def load_json(key: str) -> dict:
    with VERIFIED[key].open("r", encoding="utf-8") as stream:
        return json.load(stream)


def load_oof(key: str) -> list[dict]:
    with VERIFIED[key].open("r", encoding="utf-8", newline="") as stream:
        raw_rows = list(csv.DictReader(stream))
    fields = (
        "fold_id", "development_row_id", "y_true_malignant", "predicted_malignant",
        "selected_variant_index",
    )
    float_fields = ("p_malignant", "threshold_train_only")
    rows = []
    for item in raw_rows:
        row = dict(item)
        for field in fields:
            row[field] = int(row[field])
        for field in float_fields:
            row[field] = float(row[field])
        rows.append(row)
    return rows


report_02 = load_json("erl02_expanded")
report_03c = load_json("erl03c_report")
manifest_03b = load_json("erl03b_manifest")
manifest_03c = load_json("erl03c_manifest")
oof_03b = load_oof("erl03b_oof")
oof_03c = load_oof("erl03c_oof")

assert manifest_03b["oof"]["sha256"] == EXPECTED_FILES["erl03b_oof"][1]
assert manifest_03b["report_sha256"] == EXPECTED_FILES["erl03b_report"][1]
assert manifest_03c["output"]["outer_oof_calibration_v1_1.csv"]["sha256"] == EXPECTED_FILES["erl03c_oof"][1]
assert manifest_03c["output"]["calibration_comparison_v0_3_1.json"]["sha256"] == EXPECTED_FILES["erl03c_report"][1]
assert manifest_03b["historical_holdout_evaluated_in_this_run"] is False
assert manifest_03c["historical_holdout_evaluated_in_this_run"] is False
assert report_03c["provenance"]["historical_holdout_evaluated_in_this_run"] is False
assert report_02["selected_model_external_validation"] is False
assert report_03c["selected_model_external_validation"] is False
assert report_03c["config"]["ensemble"] is False
assert report_03c["config"]["calibration_method"] == "sigmoid"
assert report_03c["config"]["seed"] == 42
assert manifest_03c["run_id"] == report_03c["run_id"]
assert {x["run_id"] for x in oof_03c} == {manifest_03c["run_id"]}
assert {x["run_id"] for x in oof_03b} == {manifest_03b["run_id"]}
assert sha256_file(VERIFIED["erl02_expanded"]) == sha256_file(VERIFIED["erl03b_report"])


def kpi(title: str, value: str) -> str:
    return '<div><strong>' + html.escape(value) + '</strong><span>' + html.escape(title) + '</span></div>'


cards = (kpi("Datos de desarrollo", "455") + kpi("Folds externos", "5")
         + kpi("Modelos benchmark", "7") + kpi("OOF calibración", "1.820"))
display(HTML('<div class="erl-kpi">' + cards + '</div>'))
print("ERL GATE 02 · Manifiestos y contrato experimental PASS")

LABELS = {
    "dummy_prior": "Referencia por prevalencia",
    "logistic_regression": "Regresión logística",
    "logistic_elastic_net": "Logística ElasticNet",
    "svm_rbf": "SVM · RBF",
    "random_forest": "Random Forest",
    "xgboost": "XGBoost",
    "catboost": "CatBoost",
}
FAMILIES_03C = ("logistic_regression", "catboost")
METHODS = ("raw", "sigmoid")


def group_rows(rows: list[dict]) -> dict[tuple[str, str], list[dict]]:
    groups = {}
    for row in rows:
        groups.setdefault((row["model_family"], row["calibration_method"]), []).append(row)
    return groups


def check_partition(rows: list[dict], n_per_group: int) -> dict:
    groups = group_rows(rows)
    for key, obs in groups.items():
        ids = [r["development_row_id"] for r in obs]
        if len(obs) != n_per_group or set(ids) != set(range(n_per_group)):
            raise ValueError(f"Cobertura OOF incorrecta: {key}")
        if len(ids) != len(set(ids)):
            raise ValueError(f"Identificador OOF repetido: {key}")
        y = [r["y_true_malignant"] for r in obs]
        if Counter(y) != {0: 285, 1: 170}:
            raise ValueError(f"Clases inconsistentes: {key}")
        folds = Counter(r["fold_id"] for r in obs)
        if folds != {i: 91 for i in range(5)}:
            raise ValueError(f"Folds inconsistentes: {key}")
        for row in obs:
            p, cutoff = row["p_malignant"], row["threshold_train_only"]
            if not (math.isfinite(p) and math.isfinite(cutoff)
                    and 0 <= p <= 1 and 0 <= cutoff <= 1):
                raise ValueError(f"Score o umbral inválido: {key}")
            if row["predicted_malignant"] != int(p >= cutoff):
                raise ValueError(f"Decisión inconsistente: {key}")
    return groups


def metrics(obs: list[dict]) -> dict[str, float | int]:
    y = np.array([r["y_true_malignant"] for r in obs], dtype=int)
    p = np.array([r["p_malignant"] for r in obs], dtype=float)
    pred = np.array([r["predicted_malignant"] for r in obs], dtype=int)
    tp = int(((pred == 1) & (y == 1)).sum())
    fn = int(((pred == 0) & (y == 1)).sum())
    tn = int(((pred == 0) & (y == 0)).sum())
    fp = int(((pred == 1) & (y == 0)).sum())
    pclip = np.clip(p, 1e-15, 1 - 1e-15)
    return {
        "sensitivity_malignant": tp / (tp + fn),
        "specificity_benign": tn / (tn + fp),
        "false_negatives_malignant": fn,
        "false_positives_benign": fp,
        "average_precision": float(average_precision_score(y, p)),
        "roc_auc": float(roc_auc_score(y, p)),
        "brier": float(np.mean((p-y) ** 2)),
        "log_loss": float(np.mean(-(y * np.log(pclip) + (1-y) * np.log(1-pclip)))),
    }


def verify_metrics(actual: dict, recorded: dict, label: str) -> None:
    for name, observed in actual.items():
        if not math.isclose(observed, recorded[name], rel_tol=1e-10, abs_tol=1e-12):
            raise ValueError(f"Métrica no reproducible {label} / {name}")


pairs_03b = check_partition(oof_03b, 455)
pairs_03c = check_partition(oof_03c, 455)
assert set(pairs_03b) == {(k, "raw") for k in LABELS}
assert set(pairs_03c) == {(k, m) for k in FAMILIES_03C for m in METHODS}
assert len(oof_03b) == 3185 and len(oof_03c) == 1820

# Consistencia de etiquetas, folds e identificadores para todas las familias/métodos.
reference = {r["development_row_id"]: (r["y_true_malignant"], r["fold_id"])
             for r in pairs_03b[("dummy_prior", "raw")]}
for observation in oof_03b + oof_03c:
    k = observation["development_row_id"]
    assert reference[k] == (observation["y_true_malignant"], observation["fold_id"])

for (family, method), obs in pairs_03c.items():
    verify_metrics(metrics(obs), report_03c["comparison"][f"{family}/{method}"]["metrics"],
                   f"{family}/{method}")
    for fold in range(5):
        part = [r for r in obs if r["fold_id"] == fold]
        assert Counter(r["y_true_malignant"] for r in part) == {0: 57, 1: 34}
        recorded = report_03c["fold_details"][fold]["methods"][f"{family}/{method}"]["metrics"]
        verify_metrics(metrics(part), recorded, f"{family}/{method}/fold{fold}")

for (family, method), obs in pairs_03b.items():
    recorded = report_02["candidate_results"][family]["nested_oof_metrics"]
    verify_metrics(metrics(obs), recorded, f"ERL03B/{family}")

for family in FAMILIES_03C:
    prev = {r["development_row_id"]: r for r in pairs_03b[(family, "raw")]}
    for row in pairs_03c[(family, "raw")]:
        other = prev[row["development_row_id"]]
        assert row["p_malignant"] == other["p_malignant"]

print("ERL GATE 03 · Cobertura, clases, 5 folds, métricas y paridad RAW PASS")

BENCHMARK_ORDER = (
    "logistic_regression", "logistic_elastic_net", "svm_rbf", "catboost",
    "xgboost", "random_forest", "dummy_prior",
)
bench = {k: metrics(pairs_03b[(k, "raw")]) for k in BENCHMARK_ORDER}
fig, ax = plt.subplots(figsize=(11.7, 5.7))
ypos = np.arange(len(BENCHMARK_ORDER))
values = [bench[k]["brier"] for k in BENCHMARK_ORDER]
bar_colors = [
    COLORS["lr"] if k == "logistic_regression" else
    COLORS["cb"] if k == "catboost" else
    "#BBC6CB" if k == "dummy_prior" else "#80949D"
    for k in BENCHMARK_ORDER
]
ax.barh(ypos, values, color=bar_colors, height=.62)
ax.set_yticks(ypos, labels=[LABELS[k] for k in BENCHMARK_ORDER])
ax.invert_yaxis()
ax.set_xlim(0, max(values)*1.22)
ax.set_xlabel("Brier score · menor es mejor", labelpad=11)
ax.set_title("Siete familias, una referencia común", loc="left", pad=18)
ax.xaxis.grid(True, color=COLORS["grid"], zorder=0)
ax.set_axisbelow(True)
for y, value in zip(ypos, values):
    ax.text(value + max(values)*.016, y, f"{value:.4f}", va="center", fontsize=10,
            color=COLORS["ink"], fontweight="semibold")
fig.subplots_adjust(left=.245, right=.96, top=.88, bottom=.17)
plt.show()
print("Lectura: Brier compara calidad probabilística; el baseline evita interpretar scores sin referencia.")

def table(headers: list[str], rows: list[list[str]], caption: str) -> None:
    head = "".join(f"<th scope='col'>{html.escape(str(s))}</th>" for s in headers)
    body = "".join("<tr>" + "".join(f"<td>{html.escape(str(c))}</td>" for c in row)
                   + "</tr>" for row in rows)
    display(HTML("<div class='erl-table-wrap'><table class='erl-table'>"
                 f"<caption>{html.escape(caption)}</caption><thead><tr>{head}</tr></thead>"
                 f"<tbody>{body}</tbody></table></div>"))


table(
    ["Familia", "AP ↑", "AUROC ↑", "Sens. ↑", "Esp. ↑", "FN", "FP"],
    [[LABELS[k], f"{bench[k]['average_precision']:.4f}",
      f"{bench[k]['roc_auc']:.4f}",
      f"{100*bench[k]['sensitivity_malignant']:.2f} %",
      f"{100*bench[k]['specificity_benign']:.2f} %",
      str(bench[k]["false_negatives_malignant"]),
      str(bench[k]["false_positives_benign"])] for k in BENCHMARK_ORDER],
    "Tabla 1. Métricas exploratorias de siete familias sobre desarrollo OOF",
)

comparison = {f"{family}/{method}": metrics(pairs_03c[(family, method)])
              for family in FAMILIES_03C for method in METHODS}
rows = []
for family in FAMILIES_03C:
    for method in METHODS:
        m = comparison[f"{family}/{method}"]
        rows.append([
            LABELS[family], method.upper(), f"{m['brier']:.6f}",
            f"{m['log_loss']:.6f}", f"{100*m['sensitivity_malignant']:.2f} %",
            f"{100*m['specificity_benign']:.2f} %",
            str(m["false_negatives_malignant"]), str(m["false_positives_benign"]),
        ])
table(["Familia", "Método", "Brier ↓", "Log-loss ↓", "Sensibilidad", "Especificidad", "FN", "FP"], rows, "Tabla 2. Comparación pareada RAW / SIGMOID de dos familias")
for family in FAMILIES_03C:
    raw, sig = comparison[f"{family}/raw"], comparison[f"{family}/sigmoid"]
    for metric, report_key in (("brier", "delta_brier"), ("log_loss", "delta_log_loss")):
        delta = sig[metric] - raw[metric]
        expected = report_03c["comparison"][f"{family}/sigmoid_minus_raw"][report_key]
        assert math.isclose(delta, expected, rel_tol=1e-10, abs_tol=1e-12)
print("Los deltas RAW/SIGMOID coinciden con el protocolo congelado.")

def paired_effect_chart(metric: str, multiplier: float, units: str, title: str) -> None:
    arr = [multiplier * (comparison[f"{k}/sigmoid"][metric] -
                         comparison[f"{k}/raw"][metric]) for k in FAMILIES_03C]
    extent = max(abs(x) for x in arr) * 1.32
    fig, ax = plt.subplots(figsize=(11.3, 4.0))
    for index, (value, family) in enumerate(zip(arr, FAMILIES_03C)):
        color = COLORS["lr"] if family == "logistic_regression" else COLORS["cb"]
        ax.barh(index, value, color=color, height=.43)
        direction = "MEJORA" if value < 0 else "EMPEORA" if value > 0 else "SIN CAMBIO"
        offset = extent * .024 * (1 if value >= 0 else -1)
        ax.text(value + offset, index, f"{value:+.3f} · {direction}",
                va="center", ha="left" if value >= 0 else "right",
                color=COLORS["ink"], fontsize=10, fontweight="semibold")
    ax.set_yticks(range(len(FAMILIES_03C)), [LABELS[k] for k in FAMILIES_03C])
    ax.invert_yaxis()
    ax.axvline(0, color=COLORS["ink"], lw=1.2)
    ax.set_xlim(-extent*1.22, extent*1.22)
    ax.set_xlabel(f"Cambio {units} · SIGMOID − RAW (negativo = mejora)", labelpad=9)
    ax.set_title(title, loc="left", pad=15)
    ax.xaxis.grid(True, color=COLORS["grid"])
    ax.set_axisbelow(True)
    fig.subplots_adjust(left=.24, right=.97, top=.80, bottom=.27)
    plt.show()


paired_effect_chart("brier", 1000, "Brier × 1.000", "Calibración · efecto sobre Brier")
paired_effect_chart("log_loss", 1000, "log-loss × 1.000", "Calibración · efecto sobre log-loss")
print("Interpretación: LR mejora ligeramente ambas pérdidas; CatBoost empeora ambas.")

def reliability_bins(obs: list[dict]) -> tuple[np.ndarray, np.ndarray]:
    sorted_obs = sorted(obs, key=lambda r: r["p_malignant"])  # estable
    parts = np.array_split(np.arange(len(sorted_obs)), 5)
    pred = np.array([np.mean([sorted_obs[j]["p_malignant"] for j in part])
                     for part in parts])
    actual = np.array([np.mean([sorted_obs[j]["y_true_malignant"] for j in part])
                       for part in parts])
    return pred, actual


fig, ax = plt.subplots(figsize=(11.5, 6.1))
ax.plot([0, 1], [0, 1], color=COLORS["reference"], linestyle=(0, (5, 4)),
        lw=1.2, label="Calibración ideal")
styles = {
    ("logistic_regression", "raw"): (COLORS["lr"], "o", "-"),
    ("logistic_regression", "sigmoid"): (COLORS["lr"], "s", "--"),
    ("catboost", "raw"): (COLORS["cb"], "^", "-"),
    ("catboost", "sigmoid"): (COLORS["cb"], "D", "--"),
}
for key, (color, marker, linestyle) in styles.items():
    pred, observed = reliability_bins(pairs_03c[key])
    recorded = report_03c["comparison"][f"{key[0]}/{key[1]}"]["calibration"]["reliability_equal_count_bins"]
    assert np.allclose(pred, [q["mean_predicted"] for q in recorded], atol=1e-12)
    assert np.allclose(observed, [q["mean_observed"] for q in recorded], atol=1e-12)
    ax.plot(pred, observed, color=color, marker=marker, linestyle=linestyle,
            markersize=7, lw=2.2, label=f"{LABELS[key[0]]} · {key[1].upper()}")
ax.set(xlim=(-.02, 1.02), ylim=(-.02, 1.02),
       xlabel="Probabilidad pronosticada media de malignidad",
       ylabel="Proporción observada de malignidad")
ax.set_aspect("equal", adjustable="box")
ax.set_title("Fiabilidad de la probabilidad pronosticada", loc="left", pad=18)
ax.legend(loc="center left", bbox_to_anchor=(1.02, .5), fontsize=10, frameon=False)
ax.grid(color=COLORS["grid"], alpha=.95)
fig.subplots_adjust(left=.10, right=.66, top=.88, bottom=.15)
plt.show()
print("Grupos de 91 observaciones; proximidad a la diagonal no equivale a validación externa.")

fold_rows = []
for family in FAMILIES_03C:
    for method in METHODS:
        obs = pairs_03c[(family, method)]
        for fold in range(5):
            fold_m = metrics([r for r in obs if r["fold_id"] == fold])
            fold_rows.append([LABELS[family], method.upper(), str(fold),
                              f"{100*fold_m['sensitivity_malignant']:.2f} %",
                              f"{100*fold_m['specificity_benign']:.2f} %",
                              str(fold_m["false_negatives_malignant"]),
                              f"{fold_m['brier']:.4f}"])
table(["Familia", "Método", "Fold", "Sensibilidad", "Especificidad", "FN", "Brier"], fold_rows, "Tabla 3. Heterogeneidad descriptiva en las cinco particiones externas")

fig, ax = plt.subplots(figsize=(11.5, 5.6))
for key in styles:
    scores = [100 * metrics([r for r in pairs_03c[key] if r["fold_id"] == fold])
              ["sensitivity_malignant"] for fold in range(5)]
    color, marker, ls = styles[key]
    ax.plot(range(5), scores, linestyle=ls, marker=marker, color=color, lw=2.15,
            markersize=8, label=f"{LABELS[key[0]]} · {key[1].upper()}")
ax.axhline(95, color="#886A55", ls=(0, (4, 3)), lw=1.4,
           label="Objetivo interno ilustrativo: 95 %")
ax.set_xticks(range(5), labels=[f"Fold {i}" for i in range(5)])
ax.set_ylim(84, 102)
ax.set_ylabel("Sensibilidad de malignidad (%)")
ax.set_title("Variabilidad entre folds: el promedio no garantiza 95 %", loc="left", pad=20)
ax.legend(loc="upper left", bbox_to_anchor=(1.015, 1.0), frameon=False, fontsize=9.5)
ax.grid(axis="y", color=COLORS["grid"])
fig.subplots_adjust(left=.09, right=.69, top=.87, bottom=.16)
plt.show()

source_rows = []
for name, (relative, expected_hash) in EXPECTED_FILES.items():
    source_rows.append([name, relative.split("/")[-1], expected_hash[:16] + "…", "PASS"])
table(["Evidencia", "Archivo", "SHA256 (prefijo)", "Integridad"], source_rows, "Tabla 4. Fuentes verificadas antes del análisis")

# Control mínimo de la propia evidencia de no uso del holdout histórico.
assert manifest_03c["run_status"] == "COMPLETE"
assert manifest_03c["development_n"] == 455
assert report_03c["development_n"] == 455
assert report_02["development_n"] == 455
assert not manifest_03c["historical_holdout_evaluated_in_this_run"]
assert all(event["forbidden_intersection"] == 0 for event in report_03c["fit_audit"])
assert len(report_03c["fit_audit"]) == 320

print("ERL GATE 04 · Trazabilidad: 8 hashes, 320 eventos registrados, holdout excluido")
print("ERL GATE 05 · Sin fit, sin descarga, sin escritura de artefactos y sin exportar modelos")
print("DICTAMEN: evidencia descriptiva de investigación, no validación clínica")