#!/usr/bin/env python3
"""Experiments 5–8, 10, 11, 13, 14, 15 — cost-approach indications under baselines B0–B4.

Reads data/processed/levels_<tier>.parquet (and the "sold" supplement), applies src/cost_model/cost.py, and writes
results/tables/exp*.csv + results/raw/exp_cost.json. Benchmarks: (i) assessor building value indexed with the StatCan
BCPI (assessor alignment), (ii) asking price (market alignment), (iii) recorded sale price after listing (supplement),
(iv) Vrai Prix production estimates. Nothing is fitted on prices.
"""
import json
import sys
from pathlib import Path

import numpy as np
import pandas as pd
import yaml

sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
from src.cost_model import cost  # noqa: E402
from src.evaluation.metrics import cluster_bootstrap, ratio_stats  # noqa: E402

ROOT = Path(__file__).resolve().parents[1]
P = ROOT / "data" / "processed"
T = ROOT / "results" / "tables"
CFG = cost.CFG
R = {}
TIERS = [t for t in ["haiku", "sonnet", "opus"] if (P / f"levels_{t}.parquet").exists()]
MAIN = "sonnet" if "sonnet" in TIERS else TIERS[0]
LEVELS = ["B0", "B1", "B2", "B3", "B4"]


def levels(tier):
    L = pd.read_parquet(P / f"levels_{tier}.parquet")
    return {lvl: g.set_index("uid") for lvl, g in L.groupby("level")}


def run_level(df, p=None):
    out = cost.cost_indication(df, p)
    out["price"] = df["price"]; out["weight"] = df["weight"]; out["muni"] = df["roll_municipalite"]
    out["age_band"] = df["age_band"]; out["price_tier"] = df["price_tier"]; out["cma"] = df["cma"]; out["source"] = df["source"]
    out["ratio_assessor"] = out["improvements_depreciated"] / out["assessor_building_indexed"]
    out["ratio_price"] = out["V_cost"] / out["price"]
    return out


def summarise(out, label):
    s = {"label": label, "n": int(len(out)),
         "V_median": float(out["V_cost"].median()), "L_share_median": float((out["L"] / out["V_cost"]).median()),
         "RCN_median": float(out["RCN"].median()), "D_median": float(out["D"].median()), "eff_age_median": float(out["effective_age"].median())}
    a = ratio_stats(out["improvements_depreciated"], out["assessor_building_indexed"])
    s.update({f"assessor_{k}": v for k, v in a.items()})
    m = ratio_stats(out["V_cost"], out["price"])
    s.update({f"price_{k}": v for k, v in m.items()})
    return s


# --------------------------------------------------------------------------- Exp 5/6: decomposition and level comparison
def exp5_6():
    rows, dec = [], []
    for tier in TIERS:
        lv = levels(tier)
        for lvl in LEVELS:
          if lvl not in lv:
              continue
          for anc in ["B", "A"]:
            out = run_level(lv[lvl], {"anchor": anc})
            if anc == "B":
                out.to_parquet(P / f"cost_{tier}_{lvl}.parquet")
            s = summarise(out, f"{tier}/{lvl}/{anc}"); s["tier"] = tier; s["level"] = lvl; s["anchor"] = anc
            # cluster bootstrap on the two headline metrics
            for metric, pred, ref in [("assessor_mdape", "improvements_depreciated", "assessor_building_indexed"), ("price_mdape", "V_cost", "price")]:
                pt, lo, hi = cluster_bootstrap(out, lambda d: float((d[pred] / d[ref] - 1).abs().median() * 100), "muni", B=499)
                s[f"{metric}_lo"], s[f"{metric}_hi"] = lo, hi
            rows.append(s)
            if anc != "B":
                continue
            for grp in ["price_tier", "age_band", "cma"]:
                g = out.groupby(grp, observed=True).agg(n=("V_cost", "size"), L=("L", "median"), RCN=("RCN", "median"), D=("D", "median"),
                                                        S=("S", "median"), improvements=("improvements_depreciated", "median"), V=("V_cost", "median"),
                                                        ratio_assessor_median=("ratio_assessor", "median"), ratio_price_median=("ratio_price", "median")).reset_index()
                g["tier"] = tier; g["level"] = lvl; g["group"] = grp; g = g.rename(columns={grp: "value"})
                dec.append(g)
    df = pd.DataFrame(rows); df.to_csv(T / "exp5_levels_summary.csv", index=False)
    pd.concat(dec).to_csv(T / "exp5_decomposition.csv", index=False)
    R["exp5"] = df.to_dict(orient="records")
    # Exp 6: effect of visual grade/condition — B2 vs B0 deltas by visual grade & condition tercile (main tier)
    lv = levels(MAIN)
    b0 = run_level(lv["B0"]); b2 = run_level(lv["B2"]); d2 = lv["B2"]
    eff = pd.DataFrame({"grade": d2["grade"], "grade_visible": d2.get("grade_visible", True), "condition": d2["condition"],
                        "delta_RCN_pct": (b2["RCN"] / b0["RCN"] - 1) * 100, "delta_D_pts": (b2["D"] - b0["D"]) * 100,
                        "delta_V_pct": (b2["V_cost"] / b0["V_cost"] - 1) * 100, "delta_impr_pct": (b2["improvements_depreciated"] / b0["improvements_depreciated"] - 1) * 100,
                        "ratio_assessor_B0": b0["ratio_assessor"], "ratio_assessor_B2": b2["ratio_assessor"]})
    eff["cond_tercile"] = pd.qcut(eff["condition"], 3, labels=["low", "mid", "high"], duplicates="drop")
    g1 = eff.groupby("grade", observed=True)[["delta_RCN_pct", "delta_impr_pct", "delta_V_pct", "ratio_assessor_B0", "ratio_assessor_B2"]].median().reset_index()
    g2 = eff.groupby("cond_tercile", observed=True)[["delta_D_pts", "delta_impr_pct", "delta_V_pct", "ratio_assessor_B0", "ratio_assessor_B2"]].median().reset_index()
    g1.to_csv(T / "exp6_by_grade.csv", index=False); g2.to_csv(T / "exp6_by_condition.csv", index=False)
    eff.to_parquet(P / f"exp6_effects_{MAIN}.parquet")
    R["exp6"] = {"by_grade": g1.to_dict(orient="records"), "by_condition": g2.to_dict(orient="records"),
                 "share_grade_visible": float(eff["grade_visible"].mean()) if "grade_visible" in eff else None}
    return df


# --------------------------------------------------------------------------- Exp 7/8: sensitivities
def exp7_8():
    lv = levels(MAIN)
    base = lv["B3"]
    rows = []
    nm = CFG["temporal_regional"]["non_metro_discount_scenario"]
    scenarios = {
        "baseline (anchor B, Québec 2026 guides)": {},
        "anchor A (Altus Montréal hard costs)": {"anchor": "A"},
        "non-metro regions −15 % (Dessins Drummond 2026)": {"regional_level_factor": base["cma"].map(lambda c: 1.0 if c in nm["metro_cmas"] else nm["factor_other"])},
        "storey multiplier flat (1.0 all)": {"storey_multiplier": {"1": 1.0, "1.5": 1.0, "2": 1.0, "2.5": 1.0, "3": 1.0}},
        "regional level 0.90": {"regional_level_factor": 0.90}, "regional level 1.10": {"regional_level_factor": 1.10},
        "relative growth since 2023 (CMA BCPI / Montréal BCPI)": {"regional_level_factor": "relative"},
        "grade costs: low band": {"grade_per_sqft": {k: v[0] for k, v in cost.anchor()["grade_sensitivity_per_sqft"].items()}},
        "grade costs: high band": {"grade_per_sqft": {k: v[1] for k, v in cost.anchor()["grade_sensitivity_per_sqft"].items()}},
        "soft costs 5 %": {"soft_cost_share": 0.05}, "soft costs 15 %": {"soft_cost_share": 0.15},
        "entrepreneurial incentive 0 %": {"entrepreneurial_incentive": 0.0}, "entrepreneurial incentive 15 %": {"entrepreneurial_incentive": 0.15},
        "add-ons −30 %": {"add_on_shift": -0.30}, "add-ons +30 %": {"add_on_shift": 0.30},
        "economic life 50 y": {"economic_life": 50}, "economic life 75 y": {"economic_life": 75},
        "condition→effective age: half strength": {"condition_strength": 0.5}, "condition→effective age: 1.5× strength": {"condition_strength": 1.5},
        "external obsolescence 10 %": {"external_obsolescence": 0.10},
        "land not indexed": {"land_index_share": 0.0}, "land half indexed": {"land_index_share": 0.5},
        "no depreciation (RCN + S + L)": {"economic_life": 1e9},
        "no depreciation, land not indexed": {"economic_life": 1e9, "land_index_share": 0.0},
    }
    for name, p in scenarios.items():
        p = dict(p)
        if isinstance(p.get("regional_level_factor"), str) and p["regional_level_factor"] == "relative":
            # ratio of the CMA's BCPI level (2023=100 base) to Montréal's at the valuation quarter
            b = cost.load_bcpi(); q = CFG["valuation_quarter"]
            lvl = b[b.refPer == q].set_index("series")["value"]
            series_of = CFG["temporal_regional"]["cma_to_series"]
            p["regional_level_factor"] = base["cma"].map(lambda c: lvl.get(series_of.get(c, "composite"), lvl["composite"]) / lvl["montreal"])
        out = run_level(base, p)
        s = summarise(out, name); s["scenario"] = name; s["block"] = "7" if any(k in name for k in ["regional", "grade costs", "soft", "entrepreneurial", "add-ons", "baseline", "anchor", "metro", "storey"]) else "8"
        rows.append(s)
    df = pd.DataFrame(rows); df.to_csv(T / "exp7_8_sensitivity.csv", index=False)
    R["exp7_8"] = df.to_dict(orient="records")
    return df


# --------------------------------------------------------------------------- Exp 10/11: holdout CMA, roll vintage, subgroups
def exp10_calibration():
    """Calibrate the single depreciation parameter (economic life) so that the median ratio of depreciated improvements to
    the assessor's indexed building value equals 1 on the calibration CMAs; evaluate on the geographic holdout (Québec CMA).
    Prices are never used. Reported for B0 and B3 (main tier)."""
    hold = yaml.safe_load(open(ROOT / "configs" / "data.yaml"))["sample"]["holdout_cma_for_geo_oos"]
    lv = levels(MAIN)
    rows = []
    for anc, lvl in [(a, l) for a in ["B", "A"] for l in ["B0", "B3", "B4"]]:
        d = lv[lvl]
        cal, test = d[d["cma"] != hold], d[d["cma"] == hold]
        # (a) life-only calibration (reported even if it hits the grid bound), (b) RCN level-only calibration at life 60,
        # (c) joint grid. The implied hard cost per sq ft for a standard-grade house is reported for (b)/(c).
        def med_ratio(part, p):
            return float(run_level(part, {**p, "anchor": anc})["ratio_assessor"].median())
        life_only = min(range(40, 201, 5), key=lambda L: abs(med_ratio(cal, {"economic_life": L}) - 1))
        level_only = min([round(0.8 + 0.05 * i, 2) for i in range(29)], key=lambda k: abs(med_ratio(cal, {"regional_level_factor": k}) - 1))
        joint = min([(k, L) for k in [round(0.8 + 0.1 * i, 2) for i in range(15)] for L in (50, 60, 75)],
                    key=lambda kl: abs(med_ratio(cal, {"regional_level_factor": kl[0], "economic_life": kl[1]}) - 1))
        std = cost.anchor({"anchor": anc})["grade_per_sqft"]["standard"]
        variants = [("public parameters", {}), (f"life-only calibration ({life_only} y)", {"economic_life": life_only}),
                    (f"RCN level-only calibration (×{level_only} → standard ≈ ${std * level_only:.0f}/sq ft hard cost)", {"regional_level_factor": level_only}),
                    (f"joint calibration (×{joint[0]}, life {joint[1]} y)", {"regional_level_factor": joint[0], "economic_life": joint[1]})]
        for name, part in [("calibration", cal), ("holdout", test)]:
            for tag, p in variants:
                out = run_level(part, {**p, "anchor": anc})
                a = ratio_stats(out["improvements_depreciated"], out["assessor_building_indexed"]); m = ratio_stats(out["V_cost"], out["price"])
                rows.append({"anchor": anc, "level": lvl, "set": name, "parameters": tag, "economic_life": p.get("economic_life", CFG["depreciation"]["economic_life_years"]),
                             "level_factor": p.get("regional_level_factor", 1.0), "implied_standard_per_sqft": std * p.get("regional_level_factor", 1.0),
                             "n": len(out), "D_median": float(out["D"].median()), "assessor_median_ratio": a["median_ratio"], "assessor_mdape": a["mdape"],
                             "assessor_cod": a["cod"], "price_median_ratio": m["median_ratio"], "price_mdape": m["mdape"], "price_within_20pct": m["within_20pct"]})
    df = pd.DataFrame(rows); df.to_csv(T / "exp10_calibration.csv", index=False)
    R["exp10_calibration"] = df.to_dict(orient="records")
    return df


def exp10_11():
    rows = []
    for tier in TIERS:
        lv = levels(tier)
        for lvl in ["B0", "B3"]:
            if lvl not in lv:
                continue
            out = run_level(lv[lvl])
            df0 = pd.read_parquet(P / "sample.parquet").set_index("uid").reindex(out.index)
            out["roll_vintage"] = df0["roll_market_date"].astype(str).str[:4]
            out["source"] = df0["source"]
            out["n_images"] = df0["n_images_used"]
            for grp in ["cma", "roll_vintage", "age_band", "price_tier", "source"]:
                for val, g in out.groupby(grp, observed=True):
                    if len(g) < 15:
                        continue
                    a = ratio_stats(g["improvements_depreciated"], g["assessor_building_indexed"]); m = ratio_stats(g["V_cost"], g["price"])
                    rows.append({"tier": tier, "level": lvl, "group": grp, "value": str(val), "n": len(g),
                                 "assessor_median_ratio": a["median_ratio"], "assessor_mdape": a["mdape"], "assessor_cod": a["cod"],
                                 "price_median_ratio": m["median_ratio"], "price_mdape": m["mdape"]})
    df = pd.DataFrame(rows); df.to_csv(T / "exp10_11_subgroups.csv", index=False)
    R["exp10_11"] = df.to_dict(orient="records")
    return df


# --------------------------------------------------------------------------- Exp 13: counterfactual scenarios
def exp13():
    lv = levels(MAIN); b3 = lv["B3"]
    base = run_level(b3)
    rows = []
    order = ["basic", "standard", "upgraded", "luxury"]
    for name, mod in [
        ("grade +1 level", lambda d: d.assign(grade=d["grade"].map(lambda g: order[min(3, order.index(g) + 1)] if g in order else g))),
        ("grade −1 level", lambda d: d.assign(grade=d["grade"].map(lambda g: order[max(0, order.index(g) - 1)] if g in order else g))),
        ("condition +1", lambda d: d.assign(condition=(d["condition"] + 1).clip(upper=5))),
        ("condition −1", lambda d: d.assign(condition=(d["condition"] - 1).clip(lower=1))),
        ("extensive renovation evidence", lambda d: d.assign(renovation="extensive")),
        ("basement finished", lambda d: d.assign(basement="finished")),
        ("attached garage added", lambda d: d.assign(garage="attached")),
    ]:
        out = run_level(mod(b3))
        dv = (out["V_cost"] / base["V_cost"]) - 1
        aff = dv.abs() > 1e-9
        rows.append({"scenario": name, "share_affected": float(aff.mean()),
                     "median_delta_V_pct": float(dv.median() * 100), "median_delta_V_pct_affected": float(dv[aff].median() * 100) if aff.any() else 0.0,
                     "median_delta_improvements_pct": float(((out["improvements_depreciated"] / base["improvements_depreciated"]) - 1).median() * 100),
                     "median_delta_improvements_pct_affected": float(((out["improvements_depreciated"] / base["improvements_depreciated"]) - 1)[aff].median() * 100) if aff.any() else 0.0,
                     "median_delta_V_cad": float((out["V_cost"] - base["V_cost"]).median()),
                     "median_delta_V_cad_affected": float((out["V_cost"] - base["V_cost"])[aff].median()) if aff.any() else 0.0})
    df = pd.DataFrame(rows); df.to_csv(T / "exp13_scenarios.csv", index=False)
    R["exp13"] = df.to_dict(orient="records")
    return df


# --------------------------------------------------------------------------- Exp 14: Monte Carlo propagation
def exp14(n_draws=400, seed=20261020):
    rng = np.random.default_rng(seed)
    lv = levels(MAIN); b4 = lv["B4"].copy()
    props = pd.read_parquet(P / f"props_{MAIN}_K10.parquet").set_index("uid").reindex(b4.index)
    order = ["basic", "standard", "upgraded", "luxury"]
    gs = cost.anchor()["grade_sensitivity_per_sqft"]
    base = run_level(b4)
    draws = np.zeros((n_draws, len(b4)))
    for i in range(n_draws):
        d = b4.copy()
        # attribute uncertainty: grade jitter scaled by (1 − confidence) and by missing grade; condition noise sd 0.5 (1.0 if missing)
        conf = props["confidence_mean"].fillna(0.75).values
        miss_g = d["uncertain_grade"].fillna(False).astype(float).values if "uncertain_grade" in d else np.zeros(len(d))
        shift = rng.normal(0, 0.35 + 0.6 * (1 - conf) + 0.6 * miss_g)
        gi = d["grade"].map(lambda g: order.index(g) if g in order else 1).values + shift
        d["grade"] = [order[int(min(3, max(0, round(x))))] for x in gi]
        miss_c = d["uncertain_condition"].fillna(False).astype(float).values if "uncertain_condition" in d else np.zeros(len(d))
        d["condition"] = (d["condition"] + rng.normal(0, 0.5 + 0.5 * miss_c)).clip(1, 5)
        # parameter uncertainty: unit cost within band (uniform), life ~ U(50,75), soft ~ U(.05,.15), EI ~ U(0,.15), land share ~ U(0,1), regional ~ U(.9,1.1)
        u = rng.uniform()
        p = {"grade_per_sqft": {k: v[0] + u * (v[1] - v[0]) for k, v in gs.items()},
             "economic_life": rng.uniform(50, 75), "soft_cost_share": rng.uniform(0.05, 0.15),
             "entrepreneurial_incentive": rng.uniform(0.0, 0.15), "land_index_share": rng.uniform(0, 1),
             "regional_level_factor": rng.uniform(0.9, 1.1), "add_on_shift": rng.uniform(-0.3, 0.3)}
        draws[i] = cost.cost_indication(d, p)["V_cost"].values
    lo, med, hi = np.percentile(draws, [10, 50, 90], axis=0)
    mc = pd.DataFrame({"V_point": base["V_cost"].values, "V_mc_median": med, "V_p10": lo, "V_p90": hi, "price": base["price"].values,
                       "assessor_total_indexed": (base["L"] + base["assessor_building_indexed"]).values}, index=b4.index)
    mc["rel_width"] = (mc["V_p90"] - mc["V_p10"]) / mc["V_mc_median"]
    mc["price_in_interval"] = (mc["price"] >= mc["V_p10"]) & (mc["price"] <= mc["V_p90"])
    mc["assessor_in_interval"] = (mc["assessor_total_indexed"] >= mc["V_p10"]) & (mc["assessor_total_indexed"] <= mc["V_p90"])
    mc.to_parquet(P / f"exp14_mc_{MAIN}.parquet")
    # variance decomposition: one-at-a-time, holding others at baseline
    oat = {}
    for name, p in [("unit cost band", {"grade_per_sqft": {k: v[1] for k, v in gs.items()}}), ("economic life 50", {"economic_life": 50}),
                    ("soft+EI high", {"soft_cost_share": 0.15, "entrepreneurial_incentive": 0.15}), ("land not indexed", {"land_index_share": 0.0}),
                    ("regional 1.10", {"regional_level_factor": 1.10})]:
        oat[name] = float(((cost.cost_indication(b4, p)["V_cost"] / base["V_cost"]) - 1).abs().median() * 100)
    d_g = b4.copy(); d_g["grade"] = d_g["grade"].map(lambda g: order[min(3, order.index(g) + 1)] if g in order else g)
    oat["grade +1"] = float(((cost.cost_indication(d_g)["V_cost"] / base["V_cost"]) - 1).abs().median() * 100)
    d_c = b4.copy(); d_c["condition"] = (d_c["condition"] + 1).clip(upper=5)
    oat["condition +1"] = float(((cost.cost_indication(d_c)["V_cost"] / base["V_cost"]) - 1).abs().median() * 100)
    R["exp14"] = {"n_draws": n_draws, "rel_width_median": float(mc["rel_width"].median()), "rel_width_p25_p75": [float(mc["rel_width"].quantile(.25)), float(mc["rel_width"].quantile(.75))],
                  "price_coverage_80": float(mc["price_in_interval"].mean()), "assessor_total_coverage_80": float(mc["assessor_in_interval"].mean()),
                  "one_at_a_time_median_abs_pct": oat}
    pd.DataFrame([R["exp14"]]).to_csv(T / "exp14_mc_summary.csv", index=False)
    pd.Series(oat, name="median_abs_delta_pct").to_csv(T / "exp14_oat.csv")
    return mc


# --------------------------------------------------------------------------- Exp 15: market alignment incl. recorded sales
def exp15():
    rows = []
    lv = levels(MAIN)
    smp = pd.read_parquet(P / "sample.parquet").set_index("uid")
    for lvl in LEVELS:
        out = run_level(lv[lvl])
        s = {"level": lvl, "n": len(out)}
        s.update({f"price_{k}": v for k, v in ratio_stats(out["V_cost"], out["price"]).items()})
        # Vrai Prix production estimates as external benchmarks (unknown formula for cost_est)
        for col in ["vp_est", "vp_cost_est", "vp_model_est"]:
            s.update({f"{col}_{k}": v for k, v in ratio_stats(out["V_cost"], smp.loc[out.index, col]).items() if k in ("n", "median_ratio", "mdape")})
        rows.append(s)
    # production benchmarks against price (for context)
    for col in ["vp_est", "vp_cost_est", "roll_model_est"]:
        s = {"level": f"benchmark:{col}", "n": int(smp[col].notna().sum())}
        s.update({f"price_{k}": v for k, v in ratio_stats(smp[col], smp["price"]).items()})
        rows.append(s)
    # recorded sales after listing (supplement + sample units with a later sale)
    sold_rows = []
    sales = pd.read_parquet(P / "frame_sales_after_listing.parquet").sort_values("date").drop_duplicates("uid", keep="last").set_index("uid")
    for tag, lvfile, sfile in [("sold_supplement", P / f"levels_{MAIN}_sold.parquet", P / "sample_sold.parquet"), ("sample", P / f"levels_{MAIN}.parquet", P / "sample.parquet")]:
        if not lvfile.exists():
            continue
        L = pd.read_parquet(lvfile)
        for lvl in ["B0", "B3"]:
            d = L[L.level == lvl].set_index("uid")
            d = d[d.index.isin(sales.index)]
            if len(d) == 0:
                continue
            out = run_level(d)
            sp = sales.loc[out.index, "amount"]
            r = ratio_stats(out["V_cost"], sp); r.update({"set": tag, "level": lvl, "ask_vs_sale_median": float((out["price"] / sp).median())})
            sold_rows.append(r)
    df = pd.DataFrame(rows); df.to_csv(T / "exp15_market_alignment.csv", index=False)
    ds = pd.DataFrame(sold_rows); ds.to_csv(T / "exp15_recorded_sales.csv", index=False)
    R["exp15"] = {"market": df.to_dict(orient="records"), "sales": ds.to_dict(orient="records")}
    return df


def main():
    R["tiers"] = TIERS; R["main_tier"] = MAIN; R["schedule_version"] = CFG["version"]
    exp5_6(); exp7_8(); exp10_calibration(); exp10_11(); exp13(); exp14(); exp15()
    json.dump(R, open(ROOT / "results" / "raw" / "exp_cost.json", "w"), indent=1, default=str)
    print(pd.read_csv(T / "exp5_levels_summary.csv")[["tier", "level", "n", "V_median", "L_share_median", "D_median", "assessor_median_ratio", "assessor_mdape", "assessor_cod", "price_median_ratio", "price_mdape", "price_within_20pct"]].round(3).to_string())
    print(pd.read_csv(T / "exp7_8_sensitivity.csv")[["scenario", "V_median", "assessor_median_ratio", "assessor_mdape", "price_median_ratio", "price_mdape"]].round(3).to_string())
    print(json.dumps(R["exp14"], indent=1)); print(pd.read_csv(T / "exp15_recorded_sales.csv").round(3).to_string())
    print(pd.read_csv(T / "exp10_calibration.csv").round(3).to_string())
    print(pd.read_csv(T / "exp15_market_alignment.csv")[[c for c in pd.read_csv(T / "exp15_market_alignment.csv").columns if c in ("level", "n", "price_median_ratio", "price_mdape", "vp_est_median_ratio", "vp_cost_est_median_ratio", "vp_cost_est_mdape", "vp_model_est_median_ratio")]].round(3).to_string())


if __name__ == "__main__":
    main()