#!/usr/bin/env python3
"""Experiments 1–4, 9 (empirical part), 11 (extraction part), 12, 16 — attribute extraction evaluation.

Reads data/processed/{images,text,props}_<tier>*.parquet and the frame; writes CSV tables to results/tables/ and a
machine-readable results/raw/exp_agreement.json. No API calls.
"""
import json
import sys
from pathlib import Path

import numpy as np
import pandas as pd

sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
from src.evaluation.metrics import agreement_table, cohen_kappa, cluster_bootstrap, krippendorff_alpha_nominal  # noqa: E402
from src.evaluation.reference import add_reference_columns  # noqa: E402
from src.extraction.aggregate import GRADE_ORD, NONINFO  # noqa: E402

ROOT = Path(__file__).resolve().parents[1]
P = ROOT / "data" / "processed"
T = ROOT / "results" / "tables"
R = {}
TIERS = [t for t in ["haiku", "sonnet", "opus"] if (P / f"props_{t}_K10.parquet").exists()]
KCOMMON = 10


def load_props(tier, k=KCOMMON):
    return pd.read_parquet(P / f"props_{tier}_K{k}.parquet").set_index("uid")


def load_text(tier):
    f = P / f"text_{tier}.parquet"
    return pd.read_parquet(f).set_index("uid") if f.exists() else None


smp = add_reference_columns(pd.read_parquet(P / "sample.parquet")).set_index("uid")


# --------------------------------------------------------------------------- Experiment 1: agreement with reference fields
def exp1():
    rows = []
    for tier in TIERS:
        v = load_props(tier).reindex(smp.index)
        t = load_text(tier)
        pairs = [
            ("cladding_primary ∈ broker siding set", "visual", v["cladding_primary"],
             pd.Series([s if s else None for s in smp["ref_cladding_set"]], index=smp.index), "set"),
            ("roof_covering vs broker roofing", "visual", v["roof_covering"], smp["ref_roof"], "cat"),
            ("basement_finish vs broker basement", "visual", v["basement_finish"], smp["ref_basement"].replace({"no_basement": None}), "cat"),
            ("garage present (visual) where broker reports a garage", "visual",
             v["garage_visible"].map(lambda x: "present" if x in ("attached", "detached", "carport") else ("none" if x == "none_visible" else x)),
             smp["ref_garage"].map(lambda x: None if (x is None or (isinstance(x, float) and np.isnan(x))) else "present"), "cat"),
            ("garage type (attached/detached) vs broker", "visual", v["garage_visible"].where(v["garage_visible"].isin(["attached", "detached"]), "not_visible"),
             smp["ref_garage"].where(smp["ref_garage"].isin(["attached", "detached"]), None), "cat"),
            ("storeys_visible vs roll storeys (±0.5 tolerated)", "visual", v["storeys_visible"], smp["ref_storeys"], "storeys"),
            ("building_type_visual vs listing type", "visual", v["building_type_visual"], smp["ref_building_type"], "cat"),
            ("kitchen_counters vs broker kitchen field", "visual", v["kitchen_counters"], smp["ref_counters"], "cat"),
            ("kitchen_cabinets vs broker kitchen field", "visual", v["kitchen_cabinets"], smp["ref_cabinets"], "cat"),
        ]
        if t is not None:
            t = t.reindex(smp.index)
            pairs += [
                ("TEXT cladding ∈ broker siding set", "text", t["cladding"].map(lambda l: l[0] if isinstance(l, (list, np.ndarray)) and len(l) else None),
                 pd.Series([s if s else None for s in smp["ref_cladding_set"]], index=smp.index), "set"),
                ("TEXT roof vs broker roofing", "text", t["roof_covering"], smp["ref_roof"], "cat"),
                ("TEXT basement vs broker basement", "text", t["basement_finish"].replace({"no_basement": None}), smp["ref_basement"].replace({"no_basement": None}), "cat"),
                ("TEXT storeys vs roll storeys (±0.5)", "text", t["storeys"], smp["ref_storeys"], "storeys"),
                ("TEXT year built vs roll year (±1)", "text", t["year_built_claim"].replace({0: None}), smp["roll_year_built"], "year"),
                ("TEXT building footprint × storeys vs roll floor area (±15 %)", "text",
                 (t["building_footprint_claim_sqft"].replace({0: np.nan}) / 10.7639) * smp["roll_storeys"].fillna(1).clip(1, 3),
                 smp["roll_floor_area_m2"], "area"),
            ]
        for name, channel, pred, ref, kind in pairs:
            if kind == "set":
                inf = ~pred.isin(NONINFO) & pred.notna()
                refok = ref.notna()
                both = inf & refok
                acc = float(np.mean([p in s for p, s in zip(pred[both], ref[both])])) if both.sum() else np.nan
                row = {"n_ref": int(refok.sum()), "coverage": float(inf[refok].mean()) if refok.sum() else np.nan,
                       "n_compared": int(both.sum()), "accuracy": acc, "kappa": np.nan, "macro_f1": np.nan}
            elif kind == "storeys":
                num = pred.map({"1": 1, "1.5": 1.5, "2": 2, "2.5": 2.5, "3+": 3, "3": 3})
                refn = ref.map({"1": 1, "1.5": 1.5, "2": 2, "2.5": 2.5, "3+": 3})
                both = num.notna() & refn.notna()
                row = {"n_ref": int(refn.notna().sum()), "coverage": float(num.notna()[refn.notna()].mean()),
                       "n_compared": int(both.sum()),
                       "accuracy": float(((num - refn).abs() <= 0.5)[both].mean()) if both.sum() else np.nan,
                       "accuracy_strict": float((num == refn)[both].mean()) if both.sum() else np.nan,
                       "kappa": cohen_kappa(num[both].round().astype(int), refn[both].round().astype(int)) if both.sum() else np.nan,
                       "macro_f1": np.nan}
            elif kind == "year":
                p = pd.to_numeric(pred, errors="coerce"); both = p.notna() & ref.notna()
                row = {"n_ref": int(ref.notna().sum()), "coverage": float(p.notna().mean()), "n_compared": int(both.sum()),
                       "accuracy": float(((p - ref).abs() <= 1)[both].mean()) if both.sum() else np.nan,
                       "mae_years": float((p - ref).abs()[both].median()) if both.sum() else np.nan, "kappa": np.nan, "macro_f1": np.nan}
            elif kind == "area":
                both = pred.notna() & ref.notna()
                ratio = (pred / ref)[both]
                row = {"n_ref": int(ref.notna().sum()), "coverage": float(pred.notna().mean()), "n_compared": int(both.sum()),
                       "accuracy": float(((ratio - 1).abs() <= 0.15).mean()) if both.sum() else np.nan,
                       "median_ratio": float(ratio.median()) if both.sum() else np.nan,
                       "mdape": float((ratio - 1).abs().median() * 100) if both.sum() else np.nan, "kappa": np.nan, "macro_f1": np.nan}
            else:
                row = agreement_table(pred, ref)
            rows.append({"tier": tier, "channel": channel, "attribute": name, **row})
    df = pd.DataFrame(rows)
    df.to_csv(T / "exp1_agreement.csv", index=False)
    R["exp1"] = df.to_dict(orient="records")
    # confusion matrices for the key visual attributes (Sonnet, common K)
    tier = "sonnet" if "sonnet" in TIERS else TIERS[0]
    v = load_props(tier).reindex(smp.index)
    conf = {}
    for name, pred, ref in [("roof", v["roof_covering"], smp["ref_roof"]),
                            ("basement", v["basement_finish"], smp["ref_basement"].replace({"no_basement": None})),
                            ("storeys", v["storeys_visible"], smp["ref_storeys"]),
                            ("garage", v["garage_visible"], smp["ref_garage"])]:
        m = pred.notna() & ref.notna() & ~pred.isin(["not_applicable"])
        ct = pd.crosstab(ref[m], pred[m])
        ct.to_csv(T / f"exp1_confusion_{name}_{tier}.csv")
        conf[name] = ct.to_dict()
    R["exp1_confusion_tier"] = tier
    return df


# --------------------------------------------------------------------------- Experiment 2: channel information content
def exp2():
    rows = []
    for tier in TIERS:
        v = load_props(tier).reindex(smp.index)
        t = load_text(tier)
        if t is None:
            continue
        t = t.reindex(smp.index)
        for attr, vcol, tcol in [("building_type", "building_type_visual", "building_type"), ("storeys", "storeys_visible", "storeys"),
                                 ("cladding", "cladding_primary", None), ("roof_covering", "roof_covering", "roof_covering"),
                                 ("garage", "garage_visible", "garage"), ("basement_finish", "basement_finish", "basement_finish"),
                                 ("kitchen_grade", "kitchen_grade", "kitchen_grade"), ("bathroom_grade", "bathroom_grade", "bathroom_grade"),
                                 ("kitchen_counters", "kitchen_counters", "kitchen_counters"), ("flooring_main", "flooring_main", "flooring_main")]:
            vi = ~v[vcol].isin(NONINFO) & v[vcol].notna()
            if tcol is None:
                ti = t["cladding"].map(lambda l: isinstance(l, (list, np.ndarray)) and len(l) > 0 and l[0] not in NONINFO)
                tv = t["cladding"].map(lambda l: l[0] if isinstance(l, (list, np.ndarray)) and len(l) else None)
            else:
                ti = ~t[tcol].isin(NONINFO) & t[tcol].notna(); tv = t[tcol]
            vv = v[vcol].replace({"3+": "3", "none_visible": "none"})
            both = vi & ti
            agree = float((vv[both] == tv[both]).mean()) if both.sum() else np.nan
            rows.append({"tier": tier, "attribute": attr, "visual_coverage": float(vi.mean()), "text_coverage": float(ti.mean()),
                         "either_coverage": float((vi | ti).mean()), "both": int(both.sum()), "agreement_when_both": agree,
                         "visual_only_share": float((vi & ~ti).mean()), "text_only_share": float((ti & ~vi).mean())})
    df = pd.DataFrame(rows)
    df.to_csv(T / "exp2_channels.csv", index=False)
    R["exp2"] = df.to_dict(orient="records")
    return df


# --------------------------------------------------------------------------- Experiment 3: image count sensitivity
def exp3():
    rows = []
    for tier in TIERS:
        ks = [1, 2, 4, 8, 10, 12] if tier == "haiku" else [1, 2, 4, 8, 10]
        ref = load_props(tier, ks[-1])
        for k in ks:
            p = load_props(tier, k).reindex(ref.index)
            row = {"tier": tier, "K": k, "n": int(p["n_images_ok"].notna().sum()),
                   "share_with_kitchen": float((p["n_kitchen"] > 0).mean()), "share_with_bathroom": float((p["n_bathroom"] > 0).mean()),
                   "share_with_basement": float((p["n_basement"] > 0).mean()), "share_with_exterior": float((p["n_exterior"] > 0).mean()),
                   "grade_available": float(p["grade_composite"].notna().mean()), "condition_available": float(p["condition_composite"].notna().mean()),
                   "grade_stable_vs_full": float((p["grade_composite"] == ref["grade_composite"])[ref["grade_composite"].notna()].mean()),
                   "condition_mad_vs_full": float((p["condition_composite"] - ref["condition_composite"]).abs().median()),
                   "cladding_stable_vs_full": float((p["cladding_primary"] == ref["cladding_primary"])[~ref["cladding_primary"].isin(NONINFO)].mean())}
            rows.append(row)
    df = pd.DataFrame(rows)
    df.to_csv(T / "exp3_image_count.csv", index=False)
    R["exp3"] = df.to_dict(orient="records")
    # photo quality moderator: agreement with broker roofing/siding by photo-quality tercile (Sonnet)
    tier = "sonnet" if "sonnet" in TIERS else TIERS[0]
    v = load_props(tier).reindex(smp.index)
    q = pd.qcut(v["photo_quality_mean"], 3, labels=["low", "mid", "high"])
    rows = []
    for lab in ["low", "mid", "high"]:
        m = q == lab
        rows.append({"tier": tier, "photo_quality_tercile": lab, "n": int(m.sum()),
                     "roof_agreement": agreement_table(v.loc[m, "roof_covering"], smp.loc[m, "ref_roof"])["accuracy"],
                     "cladding_in_set": float(np.mean([p in s for p, s in zip(v.loc[m, "cladding_primary"], smp.loc[m, "ref_cladding_set"]) if s and p not in NONINFO])) if m.sum() else np.nan,
                     "confidence_mean": float(v.loc[m, "confidence_mean"].mean())})
    pd.DataFrame(rows).to_csv(T / "exp3_photo_quality.csv", index=False)
    R["exp3_quality"] = rows
    return df


# --------------------------------------------------------------------------- Experiment 4: cross-model agreement
def exp4():
    if len(TIERS) < 2:
        return None
    props = {t: load_props(t) for t in TIERS}
    common = set.intersection(*[set(p.index) for p in props.values()])
    rows = []
    attrs = ["building_type_visual", "storeys_visible", "cladding_primary", "roof_covering", "garage_visible", "basement_finish",
             "kitchen_grade", "bathroom_grade", "kitchen_counters", "kitchen_cabinets", "flooring_main", "renovation_evidence", "grade_composite"]
    for a in attrs:
        lab = np.array([[props[t].loc[u, a] if (props[t].loc[u, a] not in NONINFO) else None for u in sorted(common)] for t in TIERS], dtype=object)
        row = {"attribute": a, "n_listings": len(common), "alpha_nominal": krippendorff_alpha_nominal(lab)}
        for i in range(len(TIERS)):
            for j in range(i + 1, len(TIERS)):
                a1 = pd.Series(lab[i]); a2 = pd.Series(lab[j])
                m = a1.notna() & a2.notna()
                row[f"agree_{TIERS[i]}_{TIERS[j]}"] = float((a1[m] == a2[m]).mean()) if m.sum() else np.nan
                row[f"kappa_{TIERS[i]}_{TIERS[j]}"] = cohen_kappa(a1[m], a2[m]) if m.sum() else np.nan
                row[f"n_{TIERS[i]}_{TIERS[j]}"] = int(m.sum())
        rows.append(row)
    # continuous: condition composite
    for a in ["condition_composite", "exterior_condition", "interior_condition", "photo_quality_mean"]:
        row = {"attribute": a, "n_listings": len(common), "alpha_nominal": np.nan}
        for i in range(len(TIERS)):
            for j in range(i + 1, len(TIERS)):
                x = props[TIERS[i]].loc[sorted(common), a]; y = props[TIERS[j]].loc[sorted(common), a]
                m = x.notna() & y.notna()
                row[f"corr_{TIERS[i]}_{TIERS[j]}"] = float(np.corrcoef(x[m], y[m])[0, 1]) if m.sum() > 3 else np.nan
                row[f"mad_{TIERS[i]}_{TIERS[j]}"] = float((x - y).abs()[m].median()) if m.sum() else np.nan
                row[f"n_{TIERS[i]}_{TIERS[j]}"] = int(m.sum())
        rows.append(row)
    df = pd.DataFrame(rows)
    df.to_csv(T / "exp4_cross_model.csv", index=False)
    R["exp4"] = df.to_dict(orient="records")
    # per-image scene agreement
    imgs = {t: pd.read_parquet(P / f"images_{t}.parquet") for t in TIERS}
    key = ["uid", "image_index"]
    m = imgs[TIERS[0]][key + ["scene", "confidence"]].merge(imgs[TIERS[1]][key + ["scene", "confidence"]], on=key, suffixes=("_a", "_b"))
    R["exp4_scene_agreement"] = {"tiers": TIERS[:2], "n_images": int(len(m)), "agreement": float((m.scene_a == m.scene_b).mean()),
                                 "kappa": cohen_kappa(m.scene_a, m.scene_b)}
    pd.crosstab(m.scene_a, m.scene_b).to_csv(T / f"exp4_scene_confusion_{TIERS[0]}_{TIERS[1]}.csv")
    return df


# --------------------------------------------------------------------------- Experiment 9 (empirical): what sellers do not show
def exp9():
    rows = []
    for tier in TIERS:
        v = load_props(tier).reindex(smp.index)
        for scene, col in [("basement", "n_basement"), ("kitchen", "n_kitchen"), ("bathroom", "n_bathroom"), ("exterior", "n_exterior")]:
            shown = v[col] > 0
            row = {"tier": tier, "scene": scene, "share_shown": float(shown.mean())}
            # association of showing with broker basement claim / age / price / condition
            if scene == "basement":
                rb = smp["ref_basement"]
                for lab in ["finished", "partially_finished", "unfinished"]:
                    m = rb == lab
                    row[f"shown_given_{lab}"] = float(shown[m].mean()) if m.sum() else np.nan
                    row[f"n_{lab}"] = int(m.sum())
            row["shown_median_price"] = float(smp.loc[shown, "price"].median()); row["notshown_median_price"] = float(smp.loc[~shown, "price"].median())
            row["shown_mean_age"] = float(smp.loc[shown, "age_2026"].mean()); row["notshown_mean_age"] = float(smp.loc[~shown, "age_2026"].mean())
            row["shown_cond"] = float(v.loc[shown, "condition_composite"].mean()); row["notshown_cond"] = float(v.loc[~shown, "condition_composite"].mean())
            rows.append(row)
    df = pd.DataFrame(rows)
    df.to_csv(T / "exp9_scene_presence.csv", index=False)
    R["exp9"] = df.to_dict(orient="records")
    # position of first basement / kitchen image in seller order (Haiku, 12 images)
    tier = "haiku" if "haiku" in TIERS else TIERS[0]
    img = pd.read_parquet(P / f"images_{tier}.parquet")
    pos = img[img.ok].groupby(["uid", "scene"])["image_index"].min().reset_index()
    R["exp9_first_position"] = pos.groupby("scene")["image_index"].describe()[["count", "mean", "50%"]].to_dict()
    pos.groupby("scene")["image_index"].describe().to_csv(T / "exp9_first_position.csv")
    return df


# --------------------------------------------------------------------------- Experiment 9 (simulation): strategic omission of poor rooms
def exp9_sim(n_draws=40, seed=20261020):
    """Counterfactual MNAR: sellers omit kitchen/bathroom/basement photographs with a probability that falls with the room's
    condition (p = 0.8 at condition 1 → 0.05 at condition 5). We delete such images from the observed set, re-aggregate, and
    measure how the property-level condition, grade and cost indication shift. This bounds the bias from selective
    presentation that is *not* already in the data (which already contains the sellers' actual omissions)."""
    from src.extraction.aggregate import aggregate_images
    import sys
    sys.path.insert(0, str(ROOT))
    from src.cost_model import cost
    rng = np.random.default_rng(seed)
    tier = "sonnet" if "sonnet" in TIERS else TIERS[0]
    img = pd.read_parquet(P / f"images_{tier}.parquet"); img = img[img.ok & img.image_index.lt(KCOMMON)]
    full = aggregate_images(img).set_index("uid")
    lv = pd.read_parquet(P / f"levels_{tier}.parquet"); b2 = lv[lv.level == "B2"].set_index("uid").reindex(full.index)
    base_cost = cost.cost_indication(b2)["V_cost"]
    rows = []
    interior = img.scene.isin(["kitchen", "bathroom", "basement"])
    cond = img["interior_condition"].astype(float).where(img["interior_condition"].astype(float) > 0, 3.0)
    p_drop = np.clip(0.8 - 0.1875 * (cond - 1), 0.05, 0.8)  # 0.8 at 1 → 0.05 at 5
    for d in range(n_draws):
        keep = ~(interior & (rng.uniform(size=len(img)) < p_drop))
        agg = aggregate_images(img[keep]).set_index("uid").reindex(full.index)
        b = b2.copy()
        b["grade"] = agg["grade_composite"].fillna("standard"); b["condition"] = agg["condition_composite"].fillna(3.0)
        b["basement"] = agg["basement_finish"].where(~agg["basement_finish"].isin(NONINFO), "unfinished")
        v = cost.cost_indication(b)["V_cost"]
        rows.append({"draw": d, "share_images_dropped": float((~keep).mean()),
                     "cond_shift": float((agg["condition_composite"] - full["condition_composite"]).mean()),
                     "grade_available": float(agg["grade_composite"].notna().mean()), "grade_available_full": float(full["grade_composite"].notna().mean()),
                     "share_grade_up": float((agg["grade_composite"].map(GRADE_ORD) > full["grade_composite"].map(GRADE_ORD)).mean()),
                     "share_grade_down": float((agg["grade_composite"].map(GRADE_ORD) < full["grade_composite"].map(GRADE_ORD)).mean()),
                     "median_V_ratio": float((v / base_cost).median()), "mean_V_ratio": float((v / base_cost).mean())})
    df = pd.DataFrame(rows); df.to_csv(T / "exp9_simulation_draws.csv", index=False)
    summ = df.drop(columns="draw").mean().to_dict(); summ["n_draws"] = n_draws; summ["tier"] = tier
    pd.DataFrame([summ]).to_csv(T / "exp9_simulation.csv", index=False)
    R["exp9_sim"] = summ
    return df


# --------------------------------------------------------------------------- Error taxonomy from evidence strings (roof and cladding disagreements)
def exp1_taxonomy():
    tier = "sonnet" if "sonnet" in TIERS else TIERS[0]
    img = pd.read_parquet(P / f"images_{tier}.parquet"); img = img[img.ok & img.image_index.lt(KCOMMON) & img.scene.str.startswith("exterior", na=False)]
    img = img.merge(smp[["ref_roof", "ref_cladding_set"]], left_on="uid", right_index=True)
    roof = img[img.ref_roof.notna() & ~img.roof_covering.isin(NONINFO)]
    dis = roof[roof.roof_covering != roof.ref_roof]
    cats = {"roof not fully visible / angle": r"not (fully )?visible|partial|angle|cannot see|obscur|hidden|only|edge",
            "snow or weather": r"snow|neige|winter|wet",
            "flat / low-slope confusion": r"flat|low[- ]slope|membrane|gravel",
            "metal vs shingle look-alike": r"metal|standing seam|tôle|tole|steel",
            "distance / resolution": r"distan|far|aerial|drone|small|resolution|blur",
            "other": r"."}
    rows = []
    for name, pat in cats.items():
        m = dis.evidence.fillna("").str.contains(pat, case=False, regex=True)
        rows.append({"attribute": "roof covering", "category": name, "n": int(m.sum()), "share_of_disagreements": float(m.mean()) if len(dis) else np.nan})
        dis = dis[~m] if name != "other" else dis
    ndis = int((roof.roof_covering != roof.ref_roof).sum())
    for r in rows:
        r["share_of_disagreements"] = r["n"] / ndis if ndis else np.nan
    clad = img[img.ref_cladding_set.map(lambda s: bool(s)) & img.cladding.map(lambda l: isinstance(l, (list, np.ndarray)) and len(l) > 0 and l[0] not in NONINFO)]
    cl_dis = clad[[not any(x in s for x in l) for s, l in zip(clad.ref_cladding_set, clad.cladding)]]
    pairs = pd.Series([f"{l[0]} (image) vs {'/'.join(sorted(s))} (broker)" for s, l in zip(cl_dis.ref_cladding_set, cl_dis.cladding)]).value_counts().head(8)
    for k, v in pairs.items():
        rows.append({"attribute": "cladding", "category": k, "n": int(v), "share_of_disagreements": float(v / len(cl_dis)) if len(cl_dis) else np.nan})
    df = pd.DataFrame(rows); df.to_csv(T / "exp1_error_taxonomy.csv", index=False)
    R["exp1_taxonomy"] = {"roof_disagreements": ndis, "roof_compared": int(len(roof)), "cladding_disagreements": int(len(cl_dis)), "cladding_compared": int(len(clad))}
    return df


# --------------------------------------------------------------------------- Experiment 12: confidence calibration
def exp12():
    rows = []
    for tier in TIERS:
        img = pd.read_parquet(P / f"images_{tier}.parquet")
        img = img[img.ok & img.image_index.lt(KCOMMON)]
        ext = img[img.scene.isin(["exterior_front", "exterior_rear_or_side", "exterior_aerial_or_street"])].copy()
        ext = ext.merge(smp[["ref_roof", "ref_cladding_set"]], left_on="uid", right_index=True)
        ext["roof_ok"] = np.where(ext.ref_roof.notna() & ~ext.roof_covering.isin(NONINFO), ext.roof_covering == ext.ref_roof, np.nan)
        ext["clad_ok"] = [np.nan if (not s or c in NONINFO or not isinstance(l, (list, np.ndarray)) or len(l) == 0) else float(any(x in s for x in l))
                          for s, c, l in zip(ext.ref_cladding_set, ext.cladding.map(lambda l: l[0] if isinstance(l, (list, np.ndarray)) and len(l) else None), ext.cladding)]
        ext["conf_bin"] = pd.cut(ext.confidence.astype(float), [0, .6, .7, .8, .9, 1.0], include_lowest=True)
        g = ext.groupby("conf_bin", observed=True).agg(n=("uid", "size"), conf_mean=("confidence", "mean"),
                                                        roof_acc=("roof_ok", "mean"), roof_n=("roof_ok", "count"),
                                                        clad_acc=("clad_ok", "mean"), clad_n=("clad_ok", "count")).reset_index()
        g["tier"] = tier; g["conf_bin"] = g["conf_bin"].astype(str)
        rows.append(g)
        # abstention: share of not_visible by scene for scene-relevant fields
    df = pd.concat(rows)
    df.to_csv(T / "exp12_calibration.csv", index=False)
    R["exp12"] = df.to_dict(orient="records")
    # selective prediction: accuracy of roof when keeping only images above a confidence threshold
    sel = []
    for tier in TIERS:
        img = pd.read_parquet(P / f"images_{tier}.parquet")
        img = img[img.ok & img.image_index.lt(KCOMMON) & img.scene.str.startswith("exterior", na=False)].merge(smp[["ref_roof"]], left_on="uid", right_index=True)
        img = img[img.ref_roof.notna() & ~img.roof_covering.isin(NONINFO)]
        for thr in [0.0, 0.6, 0.7, 0.8, 0.9]:
            m = img.confidence.astype(float) >= thr
            sel.append({"tier": tier, "threshold": thr, "kept_share": float(m.mean()), "n": int(m.sum()),
                        "roof_accuracy": float((img.roof_covering == img.ref_roof)[m].mean()) if m.sum() else np.nan})
    pd.DataFrame(sel).to_csv(T / "exp12_selective.csv", index=False)
    R["exp12_selective"] = sel
    return df


# --------------------------------------------------------------------------- Experiment 16: review-cost proxy
def exp16():
    rows = []
    for tier in TIERS:
        v = load_props(tier).reindex(smp.index)
        for thr in [0.7, 0.75, 0.8, 0.85]:
            flag = (v["confidence_mean"] < thr) | v["grade_composite"].isna() | (v["n_defects"] > 0)
            rows.append({"tier": tier, "rule": f"conf<{thr} or no grade or defects", "share_flagged": float(flag.mean()),
                         "share_no_grade": float(v["grade_composite"].isna().mean()), "share_defects": float((v["n_defects"] > 0).mean())})
    df = pd.DataFrame(rows)
    df.to_csv(T / "exp16_review_load.csv", index=False)
    R["exp16"] = df.to_dict(orient="records")
    return df


def main():
    R["tiers"] = TIERS; R["K_common"] = KCOMMON
    R["n_sample"] = int(len(smp))
    R["reference_coverage"] = {c: float(smp[c].map(lambda x: bool(x) if isinstance(x, set) else x is not None and x == x).mean())
                               for c in ["ref_cladding_set", "ref_roof", "ref_basement", "ref_garage", "ref_counters", "ref_cabinets", "ref_storeys", "ref_building_type"]}
    exp1(); exp1_taxonomy(); exp2(); exp3(); exp4(); exp9(); exp9_sim(); exp12(); exp16()
    json.dump(R, open(ROOT / "results" / "raw" / "exp_agreement.json", "w"), indent=1, default=str)
    print(json.dumps({k: v for k, v in R.items() if k in ("tiers", "reference_coverage", "exp4_scene_agreement", "exp9_first_position")}, indent=1, default=str))
    print(pd.read_csv(T / "exp1_agreement.csv").round(3).to_string())


if __name__ == "__main__":
    main()