#!/usr/bin/env python3 """ 2021 HiMCM B - Lake Mead drought model. Runs on REAL Bureau of Reclamation data (monthly end-of-month elevations, 1935-2021). Produces every number, figure and table the paper reports. Deterministic (fixed RNG seed) so results reproduce. Outputs: figures/*.png + results.json + tables/*.csv """ import json, csv, os import numpy as np import pandas as pd from scipy.optimize import curve_fit from scipy import stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt np.random.seed(42) HERE = os.path.dirname(os.path.abspath(__file__)) FIG = os.path.join(HERE, "figures"); os.makedirs(FIG, exist_ok=True) TAB = os.path.join(HERE, "tables"); os.makedirs(TAB, exist_ok=True) plt.rcParams.update({"figure.dpi": 150, "font.size": 11, "axes.grid": True, "grid.alpha": 0.3, "figure.autolayout": True}) TEAL, CORAL, INK, GOLD = "#1b7f79", "#e4572e", "#22303a", "#f2a900" R = {} # results dict -> results.json # ---------------------------------------------------------------------------- # Load real data # ---------------------------------------------------------------------------- df = pd.read_csv(os.path.join(HERE, "lakemead_monthly.csv")) df["t"] = df["year"] + (df["month"] - 0.5) / 12.0 # decimal year df = df.sort_values("t").reset_index(drop=True) ann = df.groupby("year")["elev_ft"].mean().reset_index() # annual mean level R["data"] = {"n_monthly": int(len(df)), "year_min": int(df.year.min()), "year_max": int(df.year.max()), "elev_2021": round(float(df[df.year == 2021].elev_ft.mean()), 2), "elev_peak": round(float(ann.elev_ft.max()), 2), "elev_peak_year": int(ann.loc[ann.elev_ft.idxmax(), "year"])} # ---------------------------------------------------------------------------- # 1b. Elevation -> Volume power-law fit (validation against USBR Table 1) # V = a * (h - h0)^k (basin geometry => near power law) # ---------------------------------------------------------------------------- T1_h = np.array([1229.0, 1219.6, 1050.0, 895.0]) T1_V = np.array([29.686054, 28.229730, 10.217399, 2.576395]) # in maf def vol_model(h, a, h0, k): return a * np.power(np.clip(h - h0, 1e-6, None), k) p0 = [1e-4, 640.0, 2.0] (vp, _) = curve_fit(vol_model, T1_h, T1_V, p0=p0, maxfev=200000) V_pred = vol_model(T1_h, *vp) vol_err = np.abs(V_pred - T1_V) / T1_V * 100 def level_to_volume(h): return vol_model(np.asarray(h, float), *vp) R["volume_fit"] = {"a": float(vp[0]), "h0": float(vp[1]), "k": float(vp[2]), "pct_err": [round(float(e), 3) for e in vol_err], "max_pct_err": round(float(vol_err.max()), 3)} # ---------------------------------------------------------------------------- # 2a. Drought analysis on the 12-month moving average # ---------------------------------------------------------------------------- df["ma12"] = df["elev_ft"].rolling(12, min_periods=12).mean() ma = df.dropna(subset=["ma12"]).reset_index(drop=True) # a drought = sustained MA decline: rolling 24-month net change < -10 ft ma["chg24"] = ma["ma12"] - ma["ma12"].shift(24) ma["drought"] = ma["chg24"] < -10 # collapse into maximal runs >= 24 months periods, run = [], None for _, r in ma.iterrows(): if r["drought"]: if run is None: run = [r["t"], r["t"]] else: run[1] = r["t"] else: if run is not None and (run[1] - run[0]) >= 2.0: periods.append(tuple(run)) run = None if run is not None and (run[1] - run[0]) >= 2.0: periods.append(tuple(run)) drought_tbl = [] for (a, b) in periods: seg = ma[(ma.t >= a) & (ma.t <= b)] drop = float(seg.ma12.iloc[0] - seg.ma12.iloc[-1]) drought_tbl.append({"start": round(a, 1), "end": round(b, 1), "years": round(b - a, 1), "drop_ft": round(drop, 1)}) R["droughts"] = drought_tbl R["current_drought"] = max(drought_tbl, key=lambda d: d["end"]) if drought_tbl else None # ---------------------------------------------------------------------------- # 2b. Two forecast models # ---------------------------------------------------------------------------- # Model 1: recent drought period only (>=2000), exponential decay to a floor rec = ann[ann.year >= 2000].copy() x1 = rec.year.values.astype(float); y1 = rec.elev_ft.values def exp_floor(t, floor, A, r, t0=2000.0): return floor + A * np.exp(-r * (t - t0)) (m1, m1cov) = curve_fit(exp_floor, x1, y1, p0=[900, 320, 0.03], maxfev=200000, bounds=([700, 0, 0], [1050, 600, 0.5])) y1hat = exp_floor(x1, *m1) ss_res = np.sum((y1 - y1hat) ** 2); ss_tot = np.sum((y1 - y1.mean()) ** 2) m1_r2 = 1 - ss_res / ss_tot m1_resid_sd = np.sqrt(ss_res / (len(x1) - 3)) # Model 2: OLS linear on 2005-2020 (problem-prescribed window) win = ann[(ann.year >= 2005) & (ann.year <= 2020)].copy() x2 = win.year.values.astype(float); y2 = win.elev_ft.values lr = stats.linregress(x2, y2) m2_resid_sd = np.sqrt(np.sum((y2 - (lr.intercept + lr.slope * x2)) ** 2) / (len(x2) - 2)) def predict(model, yr): if model == 1: mu = float(exp_floor(yr, *m1)); return mu, mu - 1.96 * m1_resid_sd, mu + 1.96 * m1_resid_sd mu = float(lr.intercept + lr.slope * yr); return mu, mu - 1.96 * m2_resid_sd, mu + 1.96 * m2_resid_sd HORIZONS = [2025, 2030, 2050] R["model1"] = {"floor": float(m1[0]), "A": float(m1[1]), "r": float(m1[2]), "r2": round(float(m1_r2), 4), "pred": {str(h): [round(v, 1) for v in predict(1, h)] for h in HORIZONS}} R["model2"] = {"slope": float(lr.slope), "intercept": float(lr.intercept), "r2": round(float(lr.rvalue ** 2), 4), "pred": {str(h): [round(v, 1) for v in predict(2, h)] for h in HORIZONS}} # Back-test: train on <=2015, predict 2016-2021, compare to real annual means def backtest(): tr = ann[(ann.year >= 2000) & (ann.year <= 2015)] (bm1, _) = curve_fit(exp_floor, tr.year.values.astype(float), tr.elev_ft.values, p0=[900, 320, 0.03], maxfev=200000, bounds=([700, 0, 0], [1050, 600, 0.5])) trl = ann[(ann.year >= 2005) & (ann.year <= 2015)] blr = stats.linregress(trl.year.values.astype(float), trl.elev_ft.values) te = ann[(ann.year >= 2016) & (ann.year <= 2021)] a1 = exp_floor(te.year.values.astype(float), *bm1) a2 = blr.intercept + blr.slope * te.year.values.astype(float) mae1 = float(np.mean(np.abs(a1 - te.elev_ft.values))) mae2 = float(np.mean(np.abs(a2 - te.elev_ft.values))) return te, a1, a2, mae1, mae2 bt = backtest() R["backtest"] = {"mae_model1_ft": round(bt[3], 2), "mae_model2_ft": round(bt[4], 2), "winner": "Model 1" if bt[3] < bt[4] else "Model 2"} # ---------------------------------------------------------------------------- # 3. Water-budget + recycling model (chained to the 2b forecast via 1b curve) # Real anchors (cited in paper): # Lower-basin Colorado allocation ~7.5 maf/yr; Lake Mead serves ~25M people; # municipal return flow recyclable; EPA potable-reuse recovery ~0.5-0.85. # ---------------------------------------------------------------------------- POP_2021 = 25.0e6 # people served (problem statement) POP_GROWTH = 0.015 # 1.5%/yr (Census SW region, cited) PERCAP = 146.0 # gal/person/day municipal (USGS, cited) GAL_PER_AF = 325851.0 RECYCLE_RETURN = 0.55 # fraction of municipal use returning as wastewater RECOVER = 0.60 # treated-reuse recovery fraction (EPA, baseline) DEM_2021 = POP_2021 * PERCAP * 365.0 / GAL_PER_AF / 1e6 # ~4.1 maf baseline municipal def demand_af(yr): # served municipal demand grows with pop pop = POP_2021 * (1 + POP_GROWTH) ** (yr - 2021) return pop * PERCAP * 365.0 / GAL_PER_AF / 1e6 # maf/yr def shortage_cut(level): """Fraction of the municipal allocation curtailed by Lake Mead shortage tiers (USBR DCP elevation triggers 1075/1050/1025 ft, scaled to municipal share).""" if level >= 1075: return 0.00 if level >= 1050: return 0.07 if level >= 1025: return 0.13 if level >= 1000: return 0.20 return 0.27 def available_af(level): # elevation-constrained deliverable return DEM_2021 * (1 - shortage_cut(level)) def recycle_capacity(yr, recover=RECOVER): return demand_af(yr) * RECYCLE_RETURN * recover # maf/yr recovered new supply budget = [] for yr in range(2021, 2051): lvl = float(exp_floor(yr, *m1)) dem = demand_af(yr); avail = available_af(lvl); gap = max(dem - avail, 0.0) rec_cap = recycle_capacity(yr) budget.append({"year": yr, "level": round(lvl, 1), "demand": round(dem, 3), "available": round(avail, 3), "gap": round(gap, 3), "recycle": round(rec_cap, 3), "covered": round(min(rec_cap, gap) / gap * 100, 1) if gap > 1e-6 else 0.0}) R["budget"] = budget b2050 = budget[-1] R["recycling_2050"] = {"gap_maf": b2050["gap"], "recycle_maf": b2050["recycle"], "pct_covered": b2050["covered"]} # ---------------------------------------------------------------------------- # 4. Sensitivity: OAT sweeps, tornado, Monte Carlo, flip statement # ---------------------------------------------------------------------------- base = {"r": float(m1[2]), "floor": float(m1[0]), "POP_GROWTH": POP_GROWTH, "RECOVER": RECOVER, "RECYCLE_RETURN": RECYCLE_RETURN} def level2050(r=None, floor=None): r = base["r"] if r is None else r; floor = base["floor"] if floor is None else floor return floor + m1[1] * np.exp(-r * (2050 - 2000)) base_level_2050 = level2050() def covered2050(POP_GROWTH=POP_GROWTH, RECOVER=RECOVER, RECYCLE_RETURN=RECYCLE_RETURN, level=None): level = base_level_2050 if level is None else level pop = POP_2021 * (1 + POP_GROWTH) ** (2050 - 2021) dem = pop * PERCAP * 365.0 / GAL_PER_AF / 1e6 avail = DEM_2021 * (1 - shortage_cut(level)) gap = max(dem - avail, 0.0); rec = dem * RECYCLE_RETURN * RECOVER return min(rec, gap) / gap * 100 if gap > 1e-6 else 0.0 base_cov = covered2050() # tornado: +/-30% each driver, effect on % gap covered in 2050 tornado = [] for name in ["POP_GROWTH", "RECOVER", "RECYCLE_RETURN", "r"]: lo, hi = base[name] * 0.7, base[name] * 1.3 if name == "r": c_lo = covered2050(level=level2050(r=lo)); c_hi = covered2050(level=level2050(r=hi)) else: c_lo = covered2050(**{name: lo}); c_hi = covered2050(**{name: hi}) tornado.append({"param": name, "low": round(c_lo, 1), "high": round(c_hi, 1), "swing": round(abs(c_hi - c_lo), 1)}) tornado.sort(key=lambda d: -d["swing"]) R["tornado"] = tornado # Monte Carlo: 10^4 draws on the uncertain drivers -> distribution of % covered in 2050 N = 10000 draws = { "r": np.random.normal(base["r"], 0.20 * base["r"], N), "RECOVER": np.clip(np.random.normal(RECOVER, 0.10, N), 0.3, 0.9), "POP_GROWTH": np.clip(np.random.normal(POP_GROWTH, 0.004, N), 0, 0.04), "RECYCLE_RETURN": np.clip(np.random.normal(RECYCLE_RETURN, 0.07, N), 0.3, 0.75), } mc = np.array([covered2050(POP_GROWTH=draws["POP_GROWTH"][i], RECOVER=draws["RECOVER"][i], RECYCLE_RETURN=draws["RECYCLE_RETURN"][i], level=level2050(r=draws["r"][i])) for i in range(N)]) R["monte_carlo"] = {"mean": round(float(mc.mean()), 1), "p10": round(float(np.percentile(mc, 10)), 1), "p90": round(float(np.percentile(mc, 90)), 1), "p_cover_ge_100": round(float((mc >= 100).mean()) * 100, 1)} # named scenarios scn = {} for nm, kw in {"Baseline": {}, "Optimistic": {"RECOVER": 0.80, "POP_GROWTH": 0.010}, "Stress": {"RECOVER": 0.45, "POP_GROWTH": 0.020}, "Extreme": {"RECOVER": 0.40, "POP_GROWTH": 0.025}}.items(): lv = base_level_2050 scn[nm] = {"level_2050": round(float(lv), 1), "pct_covered": round(covered2050(**kw), 1)} R["scenarios"] = scn # flip statement: recovery fraction at which recycling fully covers the 2050 gap rr = np.linspace(0.3, 0.95, 200) cov_rr = np.array([covered2050(RECOVER=v) for v in rr]) flip = float(rr[np.argmax(cov_rr >= 100)]) if (cov_rr >= 100).any() else None R["flip_recover_for_full_cover"] = round(flip, 3) if flip else None # ============================================================================ # FIGURES (>=14) # ============================================================================ def save(fig, name): p = os.path.join(FIG, name); fig.savefig(p, bbox_inches="tight"); plt.close(fig) # F1 full historical monthly series fig, ax = plt.subplots(figsize=(7, 3.2)) ax.plot(df.t, df.elev_ft, lw=0.8, color=TEAL) ax.axhline(895, ls="--", color=CORAL, lw=1); ax.text(1936, 905, "dead pool 895 ft", color=CORAL, fontsize=8) ax.set(xlabel="Year", ylabel="Elevation (ft MSL)", title="Lake Mead end-of-month elevation, 1935–2021 (USBR)") save(fig, "f01_history.png") # F2 12-mo MA with drought shading fig, ax = plt.subplots(figsize=(7, 3.2)) ax.plot(ma.t, ma.ma12, color=INK, lw=1.2, label="12-mo moving average") for (a, b) in periods: ax.axvspan(a, b, color=CORAL, alpha=0.15) ax.set(xlabel="Year", ylabel="Elevation (ft)", title="Identified drought periods (shaded): sustained MA decline") ax.legend(fontsize=8); save(fig, "f02_droughts.png") # F3 annual mean + seasonal band g = df.groupby("year")["elev_ft"].agg(["mean", "min", "max"]).reset_index() fig, ax = plt.subplots(figsize=(7, 3.2)) ax.fill_between(g.year, g["min"], g["max"], color=TEAL, alpha=0.2, label="seasonal min–max") ax.plot(g.year, g["mean"], color=TEAL, lw=1.4, label="annual mean") ax.set(xlabel="Year", ylabel="Elevation (ft)", title="Annual mean and seasonal range") ax.legend(fontsize=8); save(fig, "f03_annual_band.png") # F4 elevation->volume power-law fit (validation) hh = np.linspace(895, 1229, 200) fig, ax = plt.subplots(figsize=(5.2, 3.6)) ax.plot(hh, level_to_volume(hh), color=TEAL, label=f"fit V=a(h−{vp[1]:.0f})$^{{{vp[2]:.2f}}}$") ax.scatter(T1_h, T1_V, color=CORAL, zorder=5, label="USBR Table 1") ax.set(xlabel="Elevation (ft)", ylabel="Volume (maf)", title=f"Volume curve validated to <{vol_err.max():.1f}% error") ax.legend(fontsize=8); save(fig, "f04_volume_fit.png") # F5 Model 1 fit + forecast with CI yrs = np.arange(2000, 2051) mu1 = exp_floor(yrs, *m1) fig, ax = plt.subplots(figsize=(6, 3.4)) ax.scatter(x1, y1, s=14, color=INK, label="annual mean (obs)") ax.plot(yrs, mu1, color=TEAL, label="Model 1: exp→floor") ax.fill_between(yrs, mu1 - 1.96 * m1_resid_sd, mu1 + 1.96 * m1_resid_sd, color=TEAL, alpha=0.15) for h in HORIZONS: ax.scatter([h], [predict(1, h)[0]], color=CORAL, zorder=6) ax.set(xlabel="Year", ylabel="Elevation (ft)", title=f"Model 1 (recent drought), R²={m1_r2:.3f}") ax.legend(fontsize=8); save(fig, "f05_model1.png") # F6 Model 2 fit + forecast with CI yrs2 = np.arange(2005, 2051) mu2 = lr.intercept + lr.slope * yrs2 fig, ax = plt.subplots(figsize=(6, 3.4)) ax.scatter(x2, y2, s=14, color=INK, label="2005–2020 (obs)") ax.plot(yrs2, mu2, color=GOLD, label="Model 2: OLS linear") ax.fill_between(yrs2, mu2 - 1.96 * m2_resid_sd, mu2 + 1.96 * m2_resid_sd, color=GOLD, alpha=0.15) for h in HORIZONS: ax.scatter([h], [predict(2, h)[0]], color=CORAL, zorder=6) ax.set(xlabel="Year", ylabel="Elevation (ft)", title=f"Model 2 (2005–2020 linear), R²={lr.rvalue**2:.3f}") ax.legend(fontsize=8); save(fig, "f06_model2.png") # F7 model comparison fig, ax = plt.subplots(figsize=(6.4, 3.4)) ax.scatter(ann.year, ann.elev_ft, s=10, color="grey", label="observed annual") ax.plot(yrs, mu1, color=TEAL, label="Model 1") ax.plot(yrs2, mu2, color=GOLD, label="Model 2") ax.axhline(895, ls="--", color=CORAL, lw=1) ax.set(xlabel="Year", ylabel="Elevation (ft)", title="Two models diverge after 2030 (Model 2 → sub-dead-pool)") ax.legend(fontsize=8); save(fig, "f07_compare.png") # F8 back-test te, a1, a2, mae1, mae2 = bt fig, ax = plt.subplots(figsize=(5.6, 3.4)) ax.plot(te.year, te.elev_ft, "o-", color=INK, label="actual 2016–21") ax.plot(te.year, a1, "s--", color=TEAL, label=f"M1 (MAE {mae1:.1f} ft)") ax.plot(te.year, a2, "^--", color=GOLD, label=f"M2 (MAE {mae2:.1f} ft)") ax.set(xlabel="Year", ylabel="Elevation (ft)", title="Back-test: train ≤2015, predict 2016–21") ax.legend(fontsize=8); save(fig, "f08_backtest.png") # F9 forecast volume projection vols = {h: level_to_volume(predict(1, h)[0]) for h in HORIZONS} fig, ax = plt.subplots(figsize=(5, 3.2)) ax.bar([str(h) for h in HORIZONS], [vols[h] for h in HORIZONS], color=TEAL) for h in HORIZONS: ax.text(str(h), float(vols[h]) + 0.2, f"{float(vols[h]):.1f}", ha="center", fontsize=9) ax.set(ylabel="Volume (maf)", title="Projected storage (Model 1 → volume curve)") save(fig, "f09_volume_proj.png") # F10 demand vs available supply bd = pd.DataFrame(budget) fig, ax = plt.subplots(figsize=(6.4, 3.4)) ax.plot(bd.year, bd.demand, color=CORAL, label="municipal demand") ax.plot(bd.year, bd.available, color=TEAL, label="deliverable supply") ax.fill_between(bd.year, bd.available, bd.demand, where=bd.demand > bd.available, color=CORAL, alpha=0.15, label="shortfall") ax.set(xlabel="Year", ylabel="maf/yr", title="Supply–demand gap widens through 2050") ax.legend(fontsize=8); save(fig, "f10_supply_demand.png") # F11 recycling offset fig, ax = plt.subplots(figsize=(6.4, 3.2)) ax.plot(bd.year, bd.gap, color=CORAL, label="shortfall") ax.plot(bd.year, bd.recycle, color=TEAL, label="recycling capacity (60% recovery)") ax.set(xlabel="Year", ylabel="maf/yr", title=f"Recycling covers ≈{b2050['covered']:.0f}% of the 2050 gap") ax.legend(fontsize=8); save(fig, "f11_recycle.png") # F12 sensitivity tornado fig, ax = plt.subplots(figsize=(6, 3.2)) names = [t["param"] for t in tornado][::-1] swings = [t["swing"] for t in tornado][::-1] ax.barh(names, swings, color=TEAL) ax.set(xlabel="Swing in % gap covered (±30%)", title="Tornado: recovery fraction dominates") save(fig, "f12_tornado.png") # F13 Monte Carlo distribution fig, ax = plt.subplots(figsize=(6, 3.2)) ax.hist(mc, bins=40, color=TEAL, alpha=0.8) ax.axvline(mc.mean(), color=CORAL, lw=2, label=f"mean {mc.mean():.0f}%") ax.axvline(np.percentile(mc, 10), color=INK, ls="--", lw=1, label=f"p10 {np.percentile(mc,10):.0f}%") ax.set(xlabel="% of 2050 gap covered by recycling", ylabel="frequency", title="Monte-Carlo (N=10⁴): mean and worst-case tail") ax.legend(fontsize=8); save(fig, "f13_montecarlo.png") # F14 scenario comparison fig, ax = plt.subplots(figsize=(5.6, 3.2)) nm = list(scn.keys()); cv = [scn[k]["pct_covered"] for k in nm] ax.bar(nm, cv, color=[TEAL, "#3fa796", GOLD, CORAL]) for i, v in enumerate(cv): ax.text(i, v + 1, f"{v:.0f}%", ha="center", fontsize=9) ax.set(ylabel="% gap covered (2050)", title="Named scenarios") save(fig, "f15_scenarios.png") # F15 flip curve fig, ax = plt.subplots(figsize=(5.6, 3.2)) ax.plot(rr, cov_rr, color=TEAL) ax.axhline(100, ls="--", color=INK, lw=1) if flip: ax.axvline(flip, color=CORAL, lw=1.5, label=f"full cover at recovery={flip:.2f}") ax.set(xlabel="Treated-reuse recovery fraction", ylabel="% gap covered (2050)", title="Decision flip: recovery needed for full coverage") ax.legend(fontsize=8); save(fig, "f14_flip.png") # F16 key-findings multipanel (conclusion) fig, axs = plt.subplots(2, 2, figsize=(8, 5.4)) axs[0,0].plot(yrs, mu1, color=TEAL); axs[0,0].plot(yrs2, mu2, color=GOLD); axs[0,0].axhline(895, ls="--", color=CORAL, lw=.8) axs[0,0].set_title("Level forecast", fontsize=10) axs[0,1].bar([str(h) for h in HORIZONS], [vols[h] for h in HORIZONS], color=TEAL); axs[0,1].set_title("Storage (maf)", fontsize=10) axs[1,0].plot(bd.year, bd.gap, color=CORAL); axs[1,0].plot(bd.year, bd.recycle, color=TEAL); axs[1,0].set_title("Gap vs recycling", fontsize=10) axs[1,1].hist(mc, bins=30, color=TEAL); axs[1,1].axvline(mc.mean(), color=CORAL); axs[1,1].set_title("MC % covered", fontsize=10) fig.suptitle("Key findings: declining level, widening gap, partial recycling cover", fontsize=11) save(fig, "f16_keyfindings.png") # ---------------------------------------------------------------------------- # Tables (CSV) + results.json # ---------------------------------------------------------------------------- pd.DataFrame([{"Elevation_ft": h, "USBR_V_maf": v, "fit_V_maf": round(float(vol_model(h, *vp)), 3), "pct_err": round(float(e), 2)} for h, v, e in zip(T1_h, T1_V, vol_err)] ).to_csv(os.path.join(TAB, "t_volume.csv"), index=False) pd.DataFrame(drought_tbl).to_csv(os.path.join(TAB, "t_droughts.csv"), index=False) pd.DataFrame([{"Year": h, "Model1_ft": predict(1, h)[0], "Model2_ft": predict(2, h)[0], "M1_volume_maf": round(float(level_to_volume(predict(1, h)[0])), 2)} for h in HORIZONS] ).to_csv(os.path.join(TAB, "t_predictions.csv"), index=False) pd.DataFrame(tornado).to_csv(os.path.join(TAB, "t_tornado.csv"), index=False) pd.DataFrame([{"Scenario": k, **v} for k, v in scn.items()]).to_csv(os.path.join(TAB, "t_scenarios.csv"), index=False) pd.DataFrame([b for b in budget if b["year"] in (2025, 2030, 2040, 2050)]).to_csv(os.path.join(TAB, "t_budget.csv"), index=False) R["counts"] = {"figures": len([f for f in os.listdir(FIG) if f.endswith(".png")]), "tables": len([t for t in os.listdir(TAB) if t.endswith(".csv")])} json.dump(R, open(os.path.join(HERE, "results.json"), "w"), indent=2) print("=== KEY RESULTS ===") print("volume fit max err %:", R["volume_fit"]["max_pct_err"]) print("droughts:", [(d['start'], d['end'], d['drop_ft']) for d in drought_tbl]) print("Model1 2025/30/50:", [R['model1']['pred'][str(h)][0] for h in HORIZONS], "R2", R['model1']['r2']) print("Model2 2025/30/50:", [R['model2']['pred'][str(h)][0] for h in HORIZONS], "R2", R['model2']['r2']) print("backtest MAE M1/M2:", R['backtest']['mae_model1_ft'], R['backtest']['mae_model2_ft'], "->", R['backtest']['winner']) print("2050 gap/recycle/covered:", R['recycling_2050']) print("MC mean/p10/p90:", R['monte_carlo']) print("tornado top:", R['tornado'][0]) print("flip recover for full cover:", R['flip_recover_for_full_cover']) print("figures:", R['counts']['figures'], "tables:", R['counts']['tables'])