"""
Career Planner — Graph Suite
Generates ~20 charts from the Monte Carlo simulation output.
"""
import pandas as pd
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.ticker as mticker
import seaborn as sns
from pathlib import Path

# ── Setup ─────────────────────────────────────────────────────────────────────
import sys
DATA   = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(__file__).parent / "career_permutations_v2.csv"
OUTDIR = Path(__file__).parent / "graphs"
OUTDIR.mkdir(exist_ok=True)

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.05)
plt.rcParams.update({
    "figure.dpi": 150,
    "savefig.bbox": "tight",
    "savefig.dpi": 150,
})

M = lambda x: x / 1_000_000  # scale to millions

# Inflation deflators: convert nominal future-dollar figures to 2026 real dollars.
# Sim runs age 18-50 (year 1-33). Age 50 = year 33; start year = 2026, so values are in 2058-ish $$.
# Using CPI = 3%/yr (COST_INFLATION in the simulation).
DEFL = {
    35: (1.03 ** 18),   # age 35 = year 18
    40: (1.03 ** 23),   # age 40 = year 23
    45: (1.03 ** 28),   # age 45 = year 28
    50: (1.03 ** 33),   # age 50 = year 33
}

df = pd.read_csv(DATA)

# Clean mixed-type columns
df["median_actual_buy_age"]  = pd.to_numeric(df["median_actual_buy_age"],  errors="coerce")
df["median_debt_free_age"]   = pd.to_numeric(df["median_debt_free_age"],   errors="coerce")
df["requested_buy_age"]      = pd.to_numeric(df["requested_buy_age"],      errors="coerce")

# Real (2026) dollar equivalents — deflate each nominal milestone by the CPI compounding factor.
# These remove the distortion from 30 years of 3% inflation: $5M nominal ≈ $2M in 2026 dollars.
for age, defl in DEFL.items():
    col = f"net_worth_p50_age{age}"
    if col in df.columns:
        df[f"real_nw_p50_age{age}"] = df[col] / defl

# Accessible net worth at 50: total minus super (super locked until age 60, preservation age).
# Liquid savings + property equity is what someone can actually use at age 50.
df["accessible_nw_p50_age50"] = (
    df["net_worth_p50_age50"] - df["median_super_age50"]
)

CAREER_ORDER = [
    "FIFO E&I",
    "FIFO — Supervisor track",
    "Mining Engineering (FIFO)",
    "Petroleum Engineering (FIFO)",
    "FIFO Electrical Engineer",
    "FIFO I&C Engineer",
    "Residential Mining Electrician",
    "Local Trade — stay Geraldton",
    "Trade → Bridge → Engineering",
    "RF/Satellite Engineering",
    "Engineering — Canberra defence",
    "Cloud/Solutions Architect",
    "Security Architect",
    "Aerospace Engineering",
    "Mechatronics/Robotics Engineering",
    "Cybersecurity",
    "Software Engineering",
    "Data Science/AI Engineering",
    "RAAF Technical Officer",
]

PALETTE = dict(zip(CAREER_ORDER, sns.color_palette("tab20", len(CAREER_ORDER))))

def fmt_m(x, _=None):
    return f"${x:.1f}M"

def save(name):
    p = OUTDIR / f"{name}.png"
    plt.savefig(p)
    plt.close()
    print(f"  saved {name}.png")

# ── 1. Career P10/P50/P90 range — best single strategy per career ─────────────
print("[1] Career wealth range (best scenario per path)...")

# best scenario = highest P50 across all parameter combos per career
best = (df.sort_values("net_worth_p50_age50", ascending=False)
         .drop_duplicates("career_path")
         .set_index("career_path")
         .reindex(CAREER_ORDER)
         .dropna(how="all"))

fig, ax = plt.subplots(figsize=(13, 7))
y = np.arange(len(best))
ax.barh(y, M(best["net_worth_p90_age50"] - best["net_worth_p10_age50"]),
        left=M(best["net_worth_p10_age50"]),
        height=0.6, color=[PALETTE[c] for c in best.index], alpha=0.35, label="P10–P90 range")
ax.scatter(M(best["net_worth_p50_age50"]), y, color=[PALETTE[c] for c in best.index],
           zorder=5, s=60, label="P50 (median)")
ax.scatter(M(best["net_worth_p10_age50"]), y, color=[PALETTE[c] for c in best.index],
           marker="|", s=120, zorder=5)
ax.scatter(M(best["net_worth_p90_age50"]), y, color=[PALETTE[c] for c in best.index],
           marker="|", s=120, zorder=5)
ax.set_yticks(y)
ax.set_yticklabels(best.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Net Worth at 50")
ax.set_title("Net Worth at 50 — Best Scenario per Career Path\n(P10 / Median / P90 across 200 Monte Carlo runs)", fontweight="bold")
ax.legend(loc="lower right")
save("01_career_wealth_range")

# ── 2. P50 wealth by career, strategy, and location (heatmap) ────────────────
print("[2] P50 heatmap: career × location (best savings/deposit per cell)...")

pivot_data = (
    df.groupby(["career_path", "location"])["net_worth_p50_age50"]
    .max()
    .unstack("location")
    .reindex(CAREER_ORDER)
)

fig, ax = plt.subplots(figsize=(11, 8))
sns.heatmap(
    M(pivot_data), ax=ax,
    fmt=".2f", annot=True, cmap="YlOrRd",
    cbar_kws={"label": "Median Net Worth at 50 ($M)"},
    linewidths=0.4,
)
ax.set_title("Median Net Worth at 50 by Career × Property Location\n(best savings/deposit combo per cell)", fontweight="bold")
ax.set_xlabel("Property Location / Strategy")
ax.set_ylabel("")
plt.xticks(rotation=25, ha="right")
save("02_heatmap_career_location")

# ── 3. Strategy comparison per career ─────────────────────────────────────────
print("[3] Strategy comparison bar chart...")

strat_best = (
    df.groupby(["career_path", "strategy"])["net_worth_p50_age50"]
    .max()
    .unstack("strategy")
    .reindex(CAREER_ORDER)
)
strat_colors = {"BUY_HOLD": "#2196F3", "BUY_LAND_BUILD": "#4CAF50", "RENT_FOREVER": "#FF5722"}

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(strat_best))
w = 0.26
for i, (strat, col) in enumerate(strat_colors.items()):
    vals = strat_best[strat] if strat in strat_best.columns else pd.Series([np.nan]*len(strat_best))
    ax.bar(x + (i-1)*w, M(vals), width=w, label=strat.replace("_", " ").title(),
           color=col, alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(strat_best.index, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Property Strategy Impact on Median Net Worth at 50\n(best deposit/savings/buy-age per strategy)", fontweight="bold")
ax.legend()
save("03_strategy_comparison")

# ── 4. Savings rate sensitivity ───────────────────────────────────────────────
print("[4] Savings rate sensitivity...")

sr_data = (
    df.groupby(["career_path", "savings_rate"])["net_worth_p50_age50"]
    .max()
    .unstack("savings_rate")
    .reindex(CAREER_ORDER)
)

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(sr_data))
sr_colors = {0.55: "#90A4AE", 0.65: "#B0BEC5", 0.75: "#42A5F5", 0.85: "#1565C0", 0.90: "#0D47A1"}
sr_present = {sr: col for sr, col in sr_colors.items() if sr in sr_data.columns}
n_sr = len(sr_present)
w = 0.80 / n_sr
for i, (sr, col) in enumerate(sr_present.items()):
    offset = (i - (n_sr - 1) / 2) * w
    ax.bar(x + offset, M(sr_data[sr]), width=w,
           label=f"{int(sr*100)}% savings rate", color=col, alpha=0.9)
ax.set_xticks(x)
ax.set_xticklabels(sr_data.index, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Best-case Median Net Worth at 50")
ax.set_title("Impact of Savings Rate on Net Worth at 50 by Career\n(best property/deposit combo)", fontweight="bold")
ax.legend()
save("04_savings_rate_sensitivity")

# ── 5. Deposit size impact ────────────────────────────────────────────────────
print("[5] 20% vs 50% deposit comparison...")

dep_data = df[df["strategy"] != "RENT_FOREVER"].groupby(
    ["career_path", "deposit_pct"])["net_worth_p50_age50"].max().unstack()

fig, ax = plt.subplots(figsize=(13, 6))
careers = [c for c in CAREER_ORDER if c in dep_data.index]
x = np.arange(len(careers))
dep_cols = {0.10: ("10% deposit (+LMI)", "#EF9A9A"), 0.20: ("20% deposit", "#FF7043"), 0.50: ("50% deposit", "#26A69A")}
dep_present = {d: (lbl, col) for d, (lbl, col) in dep_cols.items() if d in dep_data.columns}
n_dep = len(dep_present)
w = 0.80 / n_dep
for i, (d, (lbl, col)) in enumerate(dep_present.items()):
    offset = (i - (n_dep - 1) / 2) * w
    ax.bar(x + offset, M(dep_data.reindex(careers)[d]), width=w, label=lbl, color=col, alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(careers, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Deposit Size Comparison — Impact on Median Net Worth at 50\n(10% + LMI / 20% / 50%; best savings rate/buy-age per column)", fontweight="bold")
ax.legend()
save("05_deposit_comparison")

# ── 6. Net worth components breakdown (stacked bar, best scenario per career) ─
print("[6] Wealth component breakdown...")

fig, ax = plt.subplots(figsize=(13, 7))
y = np.arange(len(best))
liq   = M(best["median_liquid_age50"].clip(lower=0))
sup_  = M(best["median_super_age50"])
prop_ = M(best["median_property_value_age50"])
mort_ = M(best["median_mortgage_remaining_age50"])

ax.barh(y, liq,          height=0.6, label="Liquid savings",      color="#4CAF50")
ax.barh(y, sup_,         height=0.6, left=liq, label="Super",     color="#2196F3")
ax.barh(y, prop_ - mort_,height=0.6, left=liq+sup_, label="Property equity", color="#FF9800")
ax.barh(y, mort_,        height=0.6, left=liq+sup_+prop_-mort_,
        label="Mortgage remaining", color="#F44336", alpha=0.6)
ax.set_yticks(y)
ax.set_yticklabels(best.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Median value at 50")
ax.set_title("Wealth Component Breakdown at 50 — Best Scenario per Career", fontweight="bold")
ax.legend(loc="lower right")
save("06_wealth_components")

# ── 7. Variance (P90–P10) vs Median wealth scatter — risk/reward ──────────────
print("[7] Risk vs reward scatter...")

best_all = (df.sort_values("net_worth_p50_age50", ascending=False)
              .drop_duplicates("career_path"))

fig, ax = plt.subplots(figsize=(11, 8))
for _, row in best_all.iterrows():
    cp = row["career_path"]
    ax.scatter(M(row["net_worth_p50_age50"]), M(row["net_worth_range_age50"]),
               color=PALETTE.get(cp, "grey"), s=120, zorder=5)
    ax.annotate(cp.replace(" Engineering","Eng").replace("Trade → Bridge → ","→"),
                (M(row["net_worth_p50_age50"]), M(row["net_worth_range_age50"])),
                fontsize=8, textcoords="offset points", xytext=(6, 3))
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Median Net Worth at 50 (P50)  →  Higher is better")
ax.set_ylabel("P90–P10 Outcome Range  →  Lower is more predictable")
ax.set_title("Risk vs Reward — Best Scenario per Career\n(lower-right quadrant = high wealth, low variance)", fontweight="bold")
ax.axvline(M(best_all["net_worth_p50_age50"].median()), color="grey", linestyle="--", alpha=0.4)
ax.axhline(M(best_all["net_worth_range_age50"].median()), color="grey", linestyle="--", alpha=0.4)
ax.text(0.02, 0.98, "← Low wealth\nHigh risk", transform=ax.transAxes, va="top",
        fontsize=8, color="grey")
ax.text(0.98, 0.02, "High wealth →\nLow risk", transform=ax.transAxes, ha="right",
        fontsize=8, color="grey")
save("07_risk_reward_scatter")

# ── 8. Probability of being debt-free by 50 ───────────────────────────────────
print("[8] % debt-free by 50...")

df_buy = df[df["strategy"] == "BUY_HOLD"].copy()
debt_free = (
    df_buy.groupby("career_path")["pct_debt_free_by_50"]
    .mean()
    .reindex(CAREER_ORDER)
    .dropna()
)

fig, ax = plt.subplots(figsize=(11, 6))
colors = ["#E53935" if v < 30 else "#FB8C00" if v < 70 else "#43A047" for v in debt_free]
bars = ax.barh(range(len(debt_free)), debt_free.values, color=colors, alpha=0.85)
ax.set_yticks(range(len(debt_free)))
ax.set_yticklabels(debt_free.index, fontsize=10)
ax.set_xlabel("% of Monte Carlo runs debt-free by age 50")
ax.axvline(50, color="grey", linestyle="--", alpha=0.5, label="50% threshold")
ax.set_title("Probability of Being Mortgage-Free by Age 50\n(average across all BUY_HOLD scenarios per career)", fontweight="bold")
for bar, val in zip(bars, debt_free.values):
    ax.text(val + 0.5, bar.get_y() + bar.get_height()/2, f"{val:.0f}%",
            va="center", fontsize=9)
patches = [mpatches.Patch(color=c, label=l) for c, l in
           [("#43A047","≥70%"),("#FB8C00","30–70%"),("#E53935","<30%")]]
ax.legend(handles=patches, loc="lower right")
save("08_debt_free_probability")

# ── 9. Sell-age impact (BUY_HOLD — when is the best time to move?) ───────────
print("[9] Sell-age timing impact...")

hold = df[df["strategy"] == "BUY_HOLD"].copy()
hold["sell_label"] = hold["sell_age"].apply(lambda x: "Never sell" if pd.isna(x) else f"Sell at {int(x)}")
sell_pivot = (
    hold.groupby(["career_path", "sell_label"])["net_worth_p50_age50"]
    .max()
    .unstack("sell_label")
    .reindex(CAREER_ORDER)
)
sell_cols_all = ["Never sell", "Sell at 30", "Sell at 32", "Sell at 35", "Sell at 40", "Sell at 45", "Sell at 50"]
sell_colors_all = ["#0D1B2A", "#1A237E", "#283593", "#1976D2", "#42A5F5", "#90CAF9", "#BBDEFB"]
sell_cols = [c for c in sell_cols_all if c in sell_pivot.columns]
sell_colors = [sell_colors_all[sell_cols_all.index(c)] for c in sell_cols]

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(sell_pivot))
w = 0.85 / len(sell_cols)
for i, (col, color) in enumerate(zip(sell_cols, sell_colors)):
    offset = (i - len(sell_cols)/2 + 0.5) * w
    ax.bar(x + offset, M(sell_pivot[col]), width=w, label=col, color=color, alpha=0.88)
ax.set_xticks(x)
ax.set_xticklabels(sell_pivot.index, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Effect of Sell Timing on Net Worth at 50 (BUY_HOLD strategies)\n(sell PPOR and rebuild in Geraldton at that age)", fontweight="bold")
ax.legend()
save("09_sell_age_timing")

# ── 10. Buy age impact (23 vs 28) ────────────────────────────────────────────
print("[10] Buy age comparison...")

buy_data = df[df["strategy"] == "BUY_HOLD"].groupby(
    ["career_path", "requested_buy_age"])["net_worth_p50_age50"].max().unstack()

fig, ax = plt.subplots(figsize=(13, 6))
buy_age_cols = {23.0: ("Buy at 23", "#7B1FA2"), 25.0: ("Buy at 25", "#9C27B0"), 28.0: ("Buy at 28", "#CE93D8")}
available = {age: label for age, (label, _) in buy_age_cols.items() if age in buy_data.columns}
careers = [c for c in CAREER_ORDER if c in buy_data.index]
x = np.arange(len(careers))
n = len(available)
w = 0.8 / n
for i, (age, (label, color)) in enumerate(buy_age_cols.items()):
    if age not in buy_data.columns:
        continue
    offset = (i - (n - 1) / 2) * w
    ax.bar(x + offset, M(buy_data.reindex(careers)[age]), width=w,
           label=label, color=color, alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(careers, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Buy Age Comparison — Impact on Net Worth at 50\n(BUY_HOLD, best savings rate/deposit per bar)", fontweight="bold")
ax.legend()
save("10_buy_age_comparison")

# ── 11. Super balance by career ───────────────────────────────────────────────
print("[11] Super at 50 by career...")

super_data = df.groupby("career_path")["median_super_age50"].max().reindex(CAREER_ORDER).dropna()

fig, ax = plt.subplots(figsize=(11, 6))
colors_s = [PALETTE[c] for c in super_data.index]
bars = ax.barh(range(len(super_data)), M(super_data), color=colors_s, alpha=0.85)
ax.set_yticks(range(len(super_data)))
ax.set_yticklabels(super_data.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Median Super Balance at 50")
ax.set_title("Superannuation Balance at Age 50 by Career\n(best scenario — note RAAF benefits from 16.4% employer rate)", fontweight="bold")
for bar, val in zip(bars, M(super_data)):
    ax.text(val + 0.005, bar.get_y() + bar.get_height()/2, f"${val:.2f}M",
            va="center", fontsize=8.5)
save("11_super_at_50")

# ── 12. Property value vs mortgage remaining (equity picture) ─────────────────
print("[12] Property equity picture...")

buy_scenarios = df[(df["strategy"] != "RENT_FOREVER") &
                   (df["median_property_value_age50"] > 0)].copy()
equity_best = (buy_scenarios.sort_values("net_worth_p50_age50", ascending=False)
               .drop_duplicates("career_path")
               .set_index("career_path")
               .reindex(CAREER_ORDER)
               .dropna(how="all"))

fig, ax = plt.subplots(figsize=(13, 7))
y = np.arange(len(equity_best))
prop_m = M(equity_best["median_property_value_age50"])
mort_m = M(equity_best["median_mortgage_remaining_age50"])
equity_m = prop_m - mort_m

ax.barh(y, equity_m,  height=0.55, label="Property equity", color="#FF8F00", alpha=0.85)
ax.barh(y, mort_m,    height=0.55, left=equity_m, label="Remaining mortgage", color="#EF5350", alpha=0.7)
ax.set_yticks(y)
ax.set_yticklabels(equity_best.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Median value at 50")
ax.set_title("Property Equity vs Remaining Mortgage at 50\n(best scenario per career, property-owning paths only)", fontweight="bold")
ax.legend()
save("12_property_equity")

# ── 13. P50 vs P10 "floor" — downside protection ─────────────────────────────
print("[13] Downside floor (P10)...")

floor = (df.sort_values("net_worth_p10_age50", ascending=False)
          .drop_duplicates("career_path")
          .set_index("career_path")
          .reindex(CAREER_ORDER)
          .dropna(how="all"))

fig, ax = plt.subplots(figsize=(11, 6))
y = np.arange(len(floor))
ax.barh(y, M(floor["net_worth_p10_age50"]),  height=0.45,
        label="P10 (bad-luck floor)", color="#EF5350", alpha=0.85)
ax.scatter(M(floor["net_worth_p50_age50"]), y, color=[PALETTE[c] for c in floor.index],
           s=70, zorder=5, label="P50 (median)", marker="D")
ax.set_yticks(y)
ax.set_yticklabels(floor.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Net Worth at 50")
ax.set_title("Worst-Case Floor (P10) vs Median Net Worth at 50\n(best scenario per career)", fontweight="bold")
ax.legend()
save("13_downside_floor")

# ── 14. RENT_FOREVER vs best BUY scenario per career ─────────────────────────
print("[14] Rent vs buy comparison...")

rent_nw  = df[df["strategy"] == "RENT_FOREVER"].groupby("career_path")["net_worth_p50_age50"].max()
buy_nw   = df[df["strategy"] != "RENT_FOREVER"].groupby("career_path")["net_worth_p50_age50"].max()
compare  = pd.DataFrame({"Rent forever": rent_nw, "Best buy scenario": buy_nw}).reindex(CAREER_ORDER).dropna()

fig, ax = plt.subplots(figsize=(13, 6))
x = np.arange(len(compare))
w = 0.35
ax.bar(x - w/2, M(compare["Rent forever"]),     width=w, label="Rent forever", color="#78909C", alpha=0.85)
ax.bar(x + w/2, M(compare["Best buy scenario"]), width=w, label="Best buy scenario", color="#F57C00", alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(compare.index, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Renting Forever vs Best Property Strategy — Net Worth at 50", fontweight="bold")
ax.legend()
save("14_rent_vs_buy")

# ── 15. Perth vs Geraldton vs Canberra median NW by career ───────────────────
print("[15] Location comparison per career...")

loc_best = (
    df[df["strategy"] == "BUY_HOLD"]
    .groupby(["career_path", "location"])["net_worth_p50_age50"]
    .max()
    .unstack("location")
    .reindex(CAREER_ORDER)
)
loc_cols  = ["GERALDTON_HOUSE", "PERTH", "CANBERRA"]
loc_names = {"GERALDTON_HOUSE": "Geraldton", "PERTH": "Perth", "CANBERRA": "Canberra"}
loc_colors = {"GERALDTON_HOUSE": "#8D6E63", "PERTH": "#1E88E5", "CANBERRA": "#43A047"}

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(loc_best))
w = 0.26
for i, col in enumerate(loc_cols):
    if col in loc_best.columns:
        ax.bar(x + (i-1)*w, M(loc_best[col]), width=w,
               label=loc_names[col], color=loc_colors[col], alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(loc_best.index, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Where You Buy Matters — Net Worth at 50 by Location\n(BUY_HOLD, best savings/deposit/buy-age per bar)", fontweight="bold")
ax.legend()
save("15_location_comparison")

# ── 16. Geraldton land & build vs Geraldton house buy ─────────────────────────
print("[16] Land+build vs buy existing (Geraldton)...")

gld_hold  = df[(df["location"] == "GERALDTON_HOUSE") & (df["strategy"] == "BUY_HOLD")]
gld_build = df[(df["location"] == "GERALDTON_LAND")  & (df["strategy"] == "BUY_LAND_BUILD")]

gld_h = gld_hold.groupby("career_path")["net_worth_p50_age50"].max().reindex(CAREER_ORDER).dropna()
gld_b = gld_build.groupby("career_path")["net_worth_p50_age50"].max().reindex(CAREER_ORDER).dropna()
careers_gld = [c for c in CAREER_ORDER if c in gld_h.index or c in gld_b.index]

fig, ax = plt.subplots(figsize=(13, 6))
x = np.arange(len(careers_gld))
w = 0.35
ax.bar(x - w/2, M(gld_h.reindex(careers_gld)), width=w,
       label="Buy existing house", color="#8D6E63", alpha=0.85)
ax.bar(x + w/2, M(gld_b.reindex(careers_gld)), width=w,
       label="Buy land + build", color="#D7CCC8", alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(careers_gld, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Geraldton: Buy Existing House vs Buy Land and Build\n(best savings/deposit/buy-age per bar)", fontweight="bold")
ax.legend()
save("16_geraldton_build_vs_buy")

# ── 17. Violin — full distribution of P50 outcomes by strategy ───────────────
print("[17] Distribution of P50 by strategy...")

fig, ax = plt.subplots(figsize=(10, 6))
strat_map = {"BUY_HOLD": "Buy & Hold", "BUY_LAND_BUILD": "Land + Build", "RENT_FOREVER": "Rent Forever"}
plot_df = df.copy()
plot_df["strategy_label"] = plot_df["strategy"].map(strat_map)
plot_df["nw_p50_m"] = M(plot_df["net_worth_p50_age50"])
order = ["Buy & Hold", "Land + Build", "Rent Forever"]
pal = {"Buy & Hold": "#1E88E5", "Land + Build": "#43A047", "Rent Forever": "#FB8C00"}
sns.violinplot(data=plot_df, x="strategy_label", y="nw_p50_m",
               order=order, palette=pal, ax=ax, cut=0, inner="quartile")
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Strategy")
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Distribution of Median Net Worth at 50 by Property Strategy\n(all career paths, all parameter combos)", fontweight="bold")
save("17_violin_strategy_nw")

# ── 18. Scatter: savings rate vs NW, coloured by career ──────────────────────
print("[18] Savings rate vs NW scatter...")

fig, ax = plt.subplots(figsize=(11, 7))
for career in CAREER_ORDER:
    sub = df[df["career_path"] == career]
    ax.scatter(sub["savings_rate"], M(sub["net_worth_p50_age50"]),
               color=PALETTE[career], alpha=0.25, s=12)
# Overlay means
means = df.groupby(["career_path", "savings_rate"])["net_worth_p50_age50"].mean().reset_index()
for career in CAREER_ORDER:
    sub = means[means["career_path"] == career].sort_values("savings_rate")
    ax.plot(sub["savings_rate"], M(sub["net_worth_p50_age50"]),
            color=PALETTE[career], linewidth=2, label=career)
sr_ticks = sorted(df["savings_rate"].unique())
ax.set_xticks(sr_ticks)
ax.set_xticklabels([f"{int(s*100)}%" for s in sr_ticks])
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Savings Rate")
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Savings Rate vs Net Worth at 50 — All Careers\n(dots = individual scenarios, lines = average per savings rate)", fontweight="bold")
ax.legend(fontsize=7, ncol=2, loc="upper left")
save("18_savings_rate_vs_nw_scatter")

# ── 19. Median actual buy age by career + location ───────────────────────────
print("[19] When do they actually buy?...")

buy_age_data = (
    df[df["strategy"] != "RENT_FOREVER"]
    .dropna(subset=["median_actual_buy_age"])
    .groupby(["career_path", "location"])["median_actual_buy_age"]
    .min()
    .unstack("location")
    .reindex(CAREER_ORDER)
)
loc_cols3 = ["GERALDTON_HOUSE", "PERTH", "CANBERRA"]

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(buy_age_data))
w = 0.26
for i, col in enumerate(loc_cols3):
    if col in buy_age_data.columns:
        ax.bar(x + (i-1)*w, buy_age_data[col], width=w,
               label=loc_names[col], color=loc_colors[col], alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(buy_age_data.index, rotation=35, ha="right", fontsize=9)
ax.axhline(23, color="green",  linestyle="--", alpha=0.4, label="Age 23 target")
ax.axhline(25, color="purple", linestyle="--", alpha=0.4, label="Age 25 target")
ax.axhline(28, color="orange", linestyle="--", alpha=0.4, label="Age 28 target")
ax.set_ylabel("Median Actual Buy Age")
ax.set_title("When Can Each Career Actually Afford to Buy?\n(earliest feasible median buy age across deposit/savings combos)", fontweight="bold")
ax.legend()
save("19_actual_buy_age")

# ── 20. "Top 10 all-round" comparison — composite score ──────────────────────
print("[20] Top 10 composite leaderboard...")

# Composite score: normalise P50, P10, -range, debt-free pct, then average
scored = df.copy()
scored["pct_df_fill"] = scored["pct_debt_free_by_50"].fillna(0)
for col, weight in [("net_worth_p50_age50", 0.40),
                    ("net_worth_p10_age50",  0.25),
                    ("pct_df_fill",           0.20)]:
    mn, mx = scored[col].min(), scored[col].max()
    scored[f"norm_{col}"] = (scored[col] - mn) / (mx - mn) * weight
# Penalise high variance (lower is better)
mn, mx = scored["net_worth_range_age50"].min(), scored["net_worth_range_age50"].max()
scored["norm_range"] = (1 - (scored["net_worth_range_age50"] - mn) / (mx - mn)) * 0.15
scored["composite"] = (scored["norm_net_worth_p50_age50"] +
                        scored["norm_net_worth_p10_age50"] +
                        scored["norm_pct_df_fill"] +
                        scored["norm_range"])

top10 = scored.nlargest(10, "composite")[
    ["career_path","location","strategy","savings_rate","deposit_pct",
     "net_worth_p50_age50","net_worth_p10_age50","pct_debt_free_by_50","composite"]
].reset_index(drop=True)

fig, ax = plt.subplots(figsize=(13, 5))
ax.axis("off")
table_data = []
for _, row in top10.iterrows():
    table_data.append([
        row["career_path"],
        row["location"],
        row["strategy"].replace("_", " ").title(),
        f"{int(row['savings_rate']*100)}%",
        f"{int(row['deposit_pct']*100)}%",
        f"${row['net_worth_p50_age50']/1e6:.2f}M",
        f"${row['net_worth_p10_age50']/1e6:.2f}M",
        f"{row['pct_debt_free_by_50']:.0f}%",
        f"{row['composite']:.3f}",
    ])
table = ax.table(
    cellText=table_data,
    colLabels=["Career","Location","Strategy","Save%","Dep%","P50 NW","P10 NW","Debt-free%","Score"],
    loc="center", cellLoc="center",
)
table.auto_set_font_size(False)
table.set_fontsize(8.5)
table.scale(1, 1.6)
for j in range(9):
    table[(0, j)].set_facecolor("#1565C0")
    table[(0, j)].set_text_props(color="white", fontweight="bold")
for i in range(1, 11):
    col = "#E3F2FD" if i % 2 == 0 else "white"
    for j in range(9):
        table[(i, j)].set_facecolor(col)
ax.set_title("Top 10 All-Round Scenarios — Composite Score\n(40% P50 wealth + 25% P10 floor + 20% debt-free chance + 15% low-variance)",
             fontweight="bold", pad=20, y=0.98)
save("20_top10_composite")

# ── 21. Variable importance — tornado chart ──────────────────────────────────
print("[21] Variable importance (tornado)...")

# For each lever compute: (mean NW at best value) - (mean NW at worst value)
# This shows which knob moves the needle most.
prop_df = df[df["strategy"] != "RENT_FOREVER"].copy()

lever_swings = {}

# Career choice
career_means = df.groupby("career_path")["net_worth_p50_age50"].mean()
lever_swings["Career path"] = (career_means.max() - career_means.min(),
                               career_means.idxmax(), career_means.idxmin())

# Savings rate
sr_means = df.groupby("savings_rate")["net_worth_p50_age50"].mean()
best_sr = sr_means.idxmax(); worst_sr = sr_means.idxmin()
lever_swings["Savings rate"] = (sr_means.max() - sr_means.min(),
                                f"{int(best_sr*100)}%", f"{int(worst_sr*100)}%")

# Location (property buyers only)
loc_means = prop_df.groupby("location")["net_worth_p50_age50"].mean()
lever_swings["Location"] = (loc_means.max() - loc_means.min(),
                            loc_means.idxmax(), loc_means.idxmin())

# Deposit size
dep_means = prop_df.groupby("deposit_pct")["net_worth_p50_age50"].mean()
best_dep = dep_means.idxmax(); worst_dep = dep_means.idxmin()
lever_swings["Deposit size"] = (dep_means.max() - dep_means.min(),
                                f"{int(best_dep*100)}%", f"{int(worst_dep*100)}%")

# Buy age
buyage_means = prop_df.groupby("requested_buy_age")["net_worth_p50_age50"].mean()
best_ba = buyage_means.idxmax(); worst_ba = buyage_means.idxmin()
lever_swings["Buy age"] = (buyage_means.max() - buyage_means.min(),
                           f"Age {int(best_ba)}", f"Age {int(worst_ba)}")

levers_sorted = sorted(lever_swings.items(), key=lambda x: x[1][0])
labels  = [k for k, _ in levers_sorted]
swings  = [M(v[0]) for _, v in levers_sorted]
best_labels  = [v[1] for _, v in levers_sorted]
worst_labels = [v[2] for _, v in levers_sorted]

fig, ax = plt.subplots(figsize=(11, 5))
colors_t = ["#1E88E5" if s == max(swings) else "#64B5F6" for s in swings]
bars = ax.barh(range(len(labels)), swings, color=colors_t, alpha=0.85)
ax.set_yticks(range(len(labels)))
ax.set_yticklabels(labels, fontsize=11)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Swing in Mean Net Worth at 50 (best vs worst value of that lever)")
ax.set_title("Which Lever Matters Most?\nVariability in Net Worth Attributable to Each Decision", fontweight="bold")
for i, (bar, swing, best, worst) in enumerate(zip(bars, swings, best_labels, worst_labels)):
    ax.text(bar.get_width() + 0.01, bar.get_y() + bar.get_height()/2,
            f"Best: {best} / Worst: {worst}", va="center", fontsize=9, color="#333")
ax.margins(x=0.32)
save("21_variable_importance_tornado")

# ── 22. Wealth trajectory (age 35 → 50) — top careers ────────────────────────
print("[22] Wealth trajectory over time...")

milestone_ages = [35, 40, 45, 50]
milestone_cols = ["net_worth_p50_age35", "net_worth_p50_age40",
                  "net_worth_p50_age45", "net_worth_p50_age50"]

# Best scenario per career (highest P50 at 50)
traj = (df.sort_values("net_worth_p50_age50", ascending=False)
         .drop_duplicates("career_path")
         .set_index("career_path")
         .reindex(CAREER_ORDER)
         .dropna(how="all"))

fig, ax = plt.subplots(figsize=(12, 7))
for career in traj.index:
    row = traj.loc[career]
    vals = [M(row[c]) for c in milestone_cols]
    color = PALETTE.get(career, "grey")
    ax.plot(milestone_ages, vals, color=color, linewidth=2, marker="o", markersize=5, label=career)
    ax.annotate(career.replace(" Engineering", " Eng").replace("(FIFO)", ""),
                (50, vals[-1]), fontsize=7.5,
                textcoords="offset points", xytext=(5, 0), va="center", color=color)

ax.set_xlim(34, 58)
ax.set_xticks(milestone_ages)
ax.set_xticklabels([f"Age {a}" for a in milestone_ages])
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth (P50)")
ax.set_title("Wealth Trajectory — Age 35 to 50\n(best scenario per career: top savings rate, optimal deposit and buy age)", fontweight="bold")
ax.axhline(0, color="grey", linestyle="--", alpha=0.3)
save("22_wealth_trajectory")

# ── 23. Wealth at milestone ages (grouped bar) ────────────────────────────────
print("[23] Wealth at milestone ages...")

# Show median outcome (75% savings, 20% deposit, buy_age 25 or nearest) per career
# Use best scenario per career to keep it optimistic but show the growth arc
milestone_labels = ["Age 35", "Age 40", "Age 45", "Age 50"]
milestone_colors = ["#B3E5FC", "#4FC3F7", "#0288D1", "#01579B"]

top_n = 10
top_careers = (traj["net_worth_p50_age50"]
               .dropna()
               .sort_values(ascending=False)
               .head(top_n)
               .index.tolist())

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(top_careers))
n_m = len(milestone_labels)
w = 0.78 / n_m
for i, (col, label, color) in enumerate(zip(milestone_cols, milestone_labels, milestone_colors)):
    offset = (i - (n_m - 1) / 2) * w
    vals = [M(traj.loc[c, col]) if c in traj.index else 0 for c in top_careers]
    ax.bar(x + offset, vals, width=w, label=label, color=color, alpha=0.9)
ax.set_xticks(x)
ax.set_xticklabels(top_careers, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth (P50)")
ax.set_title("Wealth Milestones — Net Worth at Ages 35, 40, 45 and 50\n(top 10 careers by final NW; best scenario per career)", fontweight="bold")
ax.legend()
save("23_wealth_milestones")

# ── 24. Opportunity cost vs best career ───────────────────────────────────────
print("[24] Opportunity cost chart...")

best_nw = traj["net_worth_p50_age50"].dropna()
top_career = best_nw.idxmax()
top_val    = best_nw.max()
gaps = (top_val - best_nw).sort_values(ascending=True)  # sorted: smallest gap at top

fig, ax = plt.subplots(figsize=(12, 7))
y = np.arange(len(gaps))
colors_oc = [PALETTE.get(c, "grey") for c in gaps.index]
bars = ax.barh(y, M(gaps), color=colors_oc, alpha=0.85)
ax.set_yticks(y)
ax.set_yticklabels(gaps.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel(f"Median NW at 50 foregone vs {top_career}")
ax.set_title(f"Opportunity Cost — What Each Career Path Costs vs the Best\n(gap in median net worth at 50 relative to {top_career})", fontweight="bold")
for bar, val in zip(bars, M(gaps)):
    if val > 0.01:
        ax.text(val + 0.01, bar.get_y() + bar.get_height()/2,
                f"−${val:.2f}M", va="center", fontsize=8.5, color="#555")
    else:
        ax.text(0.01, bar.get_y() + bar.get_height()/2,
                "← Best career", va="center", fontsize=8.5, color="#1565C0", fontweight="bold")
save("24_opportunity_cost")

# ── 25. Realistic median vs best case ────────────────────────────────────────
print("[25] Realistic median vs best case...")

# "Realistic median": 75% savings rate, 20% deposit, buy_age 25 (or 23 if 25 not available)
# vs "Best case": highest P50 per career across all combos
realistic_mask = (
    (df["savings_rate"] == 0.75) &
    (df["deposit_pct"] == 0.20) &
    (df["requested_buy_age"].isin([25.0, 23.0]))
)
realistic = (df[realistic_mask]
             .groupby("career_path")["net_worth_p50_age50"]
             .max()
             .reindex(CAREER_ORDER)
             .dropna())

best_case = (df.groupby("career_path")["net_worth_p50_age50"]
             .max()
             .reindex(CAREER_ORDER)
             .dropna())

compare2 = pd.DataFrame({"Realistic (75% save, 20% dep)": realistic,
                          "Best case": best_case}).dropna()

fig, ax = plt.subplots(figsize=(13, 7))
x = np.arange(len(compare2))
w = 0.35
ax.bar(x - w/2, M(compare2["Realistic (75% save, 20% dep)"]), width=w,
       label="Realistic (75% savings, 20% deposit)", color="#66BB6A", alpha=0.85)
ax.bar(x + w/2, M(compare2["Best case"]), width=w,
       label="Best case (optimal savings/deposit)", color="#FFA726", alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(compare2.index, rotation=35, ha="right", fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_ylabel("Median Net Worth at 50")
ax.set_title("Realistic vs Best-Case Scenario — Median Net Worth at 50\n(realistic = 75% savings rate, 20% deposit, buy at 25)", fontweight="bold")
ax.legend()
save("25_realistic_vs_best_case")

# ── 26. Real vs nominal NW — stripping 30 years of inflation ─────────────────
print("[26] Real vs nominal net worth comparison...")

real_best = (df.sort_values("net_worth_p50_age50", ascending=False)
               .drop_duplicates("career_path")
               .set_index("career_path")
               .reindex(CAREER_ORDER)
               .dropna(how="all"))

fig, ax = plt.subplots(figsize=(13, 7))
y = np.arange(len(real_best))
nom = M(real_best["net_worth_p50_age50"])
real = M(real_best["real_nw_p50_age50"])

ax.barh(y - 0.18, nom,  height=0.32, color="#90CAF9", alpha=0.9, label=f"Nominal (future $)")
ax.barh(y + 0.18, real, height=0.32, color="#1565C0", alpha=0.9, label=f"Real (2026 $, ÷{DEFL[50]:.2f})")
ax.set_yticks(y)
ax.set_yticklabels(real_best.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Median Net Worth at Age 50")
ax.set_title(
    f"Nominal vs Real (2026) Net Worth at 50 — Best Scenario per Career\n"
    f"30 years of 3% CPI inflates nominal values by ×{DEFL[50]:.2f}  ·  divide by that to get today's purchasing power",
    fontweight="bold"
)
ax.legend()
save("26_real_vs_nominal_nw")

# ── 27. Accessible NW vs total NW (super is locked until 60) ─────────────────
print("[27] Accessible vs total net worth at 50...")

acc_best = real_best.copy()
acc_best["accessible"] = M(acc_best["accessible_nw_p50_age50"].clip(lower=0))
acc_best["locked_super"] = M(acc_best["median_super_age50"])

fig, ax = plt.subplots(figsize=(13, 7))
y = np.arange(len(acc_best))
ax.barh(y, acc_best["accessible"],   height=0.55, color="#43A047", alpha=0.9, label="Accessible (liquid + property equity)")
ax.barh(y, acc_best["locked_super"], height=0.55, left=acc_best["accessible"],
        color="#B0BEC5", alpha=0.7, label="Super (locked until age 60)")
ax.set_yticks(y)
ax.set_yticklabels(acc_best.index, fontsize=10)
ax.xaxis.set_major_formatter(mticker.FuncFormatter(fmt_m))
ax.set_xlabel("Median Net Worth at Age 50")
ax.set_title(
    "Accessible vs Locked Wealth at Age 50\n"
    "Super cannot be touched until preservation age (60) — this is the real liquidity picture",
    fontweight="bold"
)
ax.legend()
save("27_accessible_vs_total_nw")

# ── Done ──────────────────────────────────────────────────────────────────────
files = sorted(OUTDIR.glob("*.png"))
print(f"\nDone — {len(files)} charts saved to {OUTDIR}/")
for f in files:
    size_kb = f.stat().st_size // 1024
    print(f"  {f.name}  ({size_kb} KB)")
