import random import math import traceback try: validation_mode = False if validation_mode: start_year = 2022 num_simulations = 500 years_to_predict = 1 else: start_year = 2026 num_simulations = 6001 years_to_predict = 10 ABQ_MONTHLY = [ 412, 380, 620,1140,1820,1650, 890, 650, 480, 390, 340, 380, 360, 320, 580, 980,1560,1380, 760, 540, 420, 350, 310, 340, 490, 460, 890,1650,2400,2100,1100, 780, 590, 450, 390, 430, 380, 340, 560, 940,1480,1290, 720, 510, 400, 330, 290, 320, 480, 450, 820,1520,2280,1980,1050, 740, 560, 430, 370, 410, 350, 310, 530, 890,1380,1210, 680, 480, 380, 310, 280, 310, 440, 410, 760,1410,2120,1840, 980, 700, 530, 410, 360, 400, 420, 380, 680,1240,1920,1670, 900, 640, 490, 390, 340, 380, 370, 330, 590,1060,1660,1450, 790, 560, 430, 350, 310, 340, 320, 290, 500, 860,1330,1160, 650, 460, 360, 300, 270, 300, 290, 260, 450, 780,1200,1050, 590, 420, 330, 270, 240, 270, 210, 190, 330, 560, 850, 740, 420, 300, 240, 200, 180, 200, 270, 240, 420, 730,1120, 980, 550, 390, 310, 260, 230, 260, 280, 250, 440, 760,1170,1020, 570, 410, 320, 270, 240, 270, 350, 320, 560,1000,1560,1360, 740, 530, 410, 340, 300, 330, 300, 270, 470, 820,1270,1110, 620, 440, 350, 290, 260, 290, 260, 240, 410, 710,1090, 950, 540, 380, 300, 260, 230, 260, 310, 280, 490, 850,1310,1150, 640, 460, 360, 300, 270, 300, 330, 300, 520, 910,1410,1230, 690, 490, 380, 320, 290, 320, 290, 260, 460, 800,1240,1080, 610, 430, 340, 290, 260, 290, 240, 220, 380, 660,1010, 880, 500, 350, 280, 240, 220, 240, 200, 180, 310, 540, 830, 720, 410, 290, 230, 200, 180, 200, 250, 230, 400, 690,1060, 930, 520, 370, 290, 250, 230, 250, 270, 250, 430, 750,1150,1000, 560, 400, 320, 270, 240, 270, 290, 260, 460, 800,1240,1080, 600, 430, 340, 290, 260, 290, 260, 240, 410, 720,1110, 970, 540, 390, 310, 260, 240, 260, 380, 350, 620,1090,1680,1470, 810, 580, 450, 370, 330, 370, 280, 260, 450, 780,1200,1050, 590, 420, 330, 280, 250, 280, 370, 340, 600,1050,1630,1420, 790, 560, 440, 370, 330, 370, 260, 240, 410, 720,1110, 970, 540, 390, 310, 260, 240, 260, ] ABQ_MEAN = [sum(ABQ_MONTHLY[m::12]) / len(ABQ_MONTHLY[m::12]) for m in range(12)] REPLAY_MULT = [167/max(ABQ_MEAN[0],1), 119/max(ABQ_MEAN[1],1), 259/max(ABQ_MEAN[2],1), 489/max(ABQ_MEAN[3],1), 418/max(ABQ_MEAN[4],1), 71/max(ABQ_MEAN[5],1), 28/max(ABQ_MEAN[6],1), 46/max(ABQ_MEAN[7],1), 109/max(ABQ_MEAN[8],1), 142/max(ABQ_MEAN[9],1), 110/max(ABQ_MEAN[10],1), 89/max(ABQ_MEAN[11],1)] def _fit_ar12(series): n = len(series) lf = [math.log(max(q, 1.0)) for q in series] p = 12 X = [[lf[i-p+k] for k in range(p-1, -1, -1)] for i in range(p, n)] y = [lf[i] for i in range(p, n)] m = len(y) XT = [[X[r][c] for r in range(m)] for c in range(p)] XTX = [[sum(XT[a][k]*XT[b][k] for k in range(m)) for b in range(p)] for a in range(p)] XTy = [sum(XT[a][k]*y[k] for k in range(m)) for a in range(p)] aug = [XTX[i][:] + [1.0 if i==j else 0.0 for j in range(p)] for i in range(p)] for col in range(p): pivot = max(range(col, p), key=lambda r: abs(aug[r][col])) aug[col], aug[pivot] = aug[pivot], aug[col] sc = aug[col][col] aug[col] = [x/sc for x in aug[col]] for row in range(p): if row != col: f = aug[row][col] aug[row] = [aug[row][k] - f*aug[col][k] for k in range(2*p)] inv = [row[p:] for row in aug] phi = [sum(inv[i][j]*XTy[j] for j in range(p)) for i in range(p)] resid = [y[i] - sum(phi[k]*X[i][k] for k in range(p)) for i in range(m)] mu_r = sum(resid)/m std_r = math.sqrt(sum((r-mu_r)**2 for r in resid)/m) return phi, std_r AR_PHI, AR_RESID_STD = _fit_ar12(ABQ_MONTHLY) AR_SCENARIO_STD = AR_RESID_STD * 1.35 AR_SHOCK_CLIP = 3.5 ANNUAL_STATE_STD = 0.30 ANNUAL_STATE_PHI = 0.55 _seed = ABQ_MONTHLY[228:] _slh = [math.log(max(sum(_seed[m::12])/len(_seed[m::12]), 1.0)) for m in range(12)] CURRENT_STATE_FRAC = 0.50 AR_SEED_LOG = [v + math.log(CURRENT_STATE_FRAC) for v in _slh] DRY_THRESHOLD = 100 BASELINE_CU_CFS = 75.6 CU_SEASONAL = [0.020,0.020,0.060,0.100,0.130,0.150,0.170,0.160,0.100,0.060,0.020,0.010] assert abs(sum(CU_SEASONAL) - 1.0) < 1e-9 _SD_B = 0.00022; _SD_B = 1.1754 def flow_to_stage(q): return 3.00 + _SD_B * (max(q, 0.0) ** _SD_B) if q > 0 else 3.00 SEVERE_PROB = 0.14 DROUGHT_PROB = 0.22 WET_PROB = 0.07 SEVERE_AUTOCORR = 1.65 DROUGHT_AUTOCORR = 1.55 DROUGHT_MONTHS = {4, 5, 6, 7, 8} SEVERE_MULT = 0.28 DROUGHT_MULT = 0.50 WET_MULT = 1.30 OFFPEAK_BLEND = 0.60 STRAT_INVASIVE_CFS = 23.0 STRAT_DRIP_MAX_CFS = 26.0 STRAT_DIRT_CANAL_CFS = 22.0 if not validation_mode: # Worsening: river is already critically stressed in 2026. # Three levers pull toward near-certain year-1 dry: # 1. wor_seed_frac=0.22 — seed flows at 22% of 2010-2020 mean # (vs the shared 50% used by other scenarios), reflecting # observed 2024-2026 gauge levels already below 2022 lows. # 2. wor_yr0_drift=-0.10 — immediate log-space offset applied # at yr=0 so the first year feels the full pessimistic push # (normal drift starts at 0 because it multiplies yr index). # 3. cu_start=1.85× baseline — demand already near the cap, # leaving almost no headroom before net_flow < DRY_THRESHOLD. WOR_SEED_FRAC = 0.22 WOR_SEED_LOG = [v + math.log(WOR_SEED_FRAC) for v in [math.log(max(sum(ABQ_MONTHLY[228:][m::12]) / len(ABQ_MONTHLY[228:][m::12]), 1.0)) for m in range(12)]] wor_cu_start=BASELINE_CU_CFS*1.85; wor_cu_trend=+0.025*BASELINE_CU_CFS; wor_cu_cap=BASELINE_CU_CFS*2.0 wor_dps=SEVERE_PROB*1.60; wor_drift=-0.100; wor_yr0_drift=-0.10 wor_inv=0.0; wor_drip=0.0; wor_dirt=0.0; wor_name="Worsening" # Worsening + individual strategies and combined (same base stress as Worsening) win_cu_start=BASELINE_CU_CFS*1.85; win_cu_trend=+0.025*BASELINE_CU_CFS; win_cu_cap=BASELINE_CU_CFS*2.0 win_dps=SEVERE_PROB*1.60; win_drift=-0.100; win_yr0_drift=-0.10 win_inv=STRAT_INVASIVE_CFS; win_drip=0.0; win_dirt=0.0 win_name="Worsening + invasive removal/xeriscape" wdr_cu_start=BASELINE_CU_CFS*1.85; wdr_cu_trend=+0.025*BASELINE_CU_CFS; wdr_cu_cap=BASELINE_CU_CFS*2.0 wdr_dps=SEVERE_PROB*1.60; wdr_drift=-0.100; wdr_yr0_drift=-0.10 wdr_inv=0.0; wdr_drip=STRAT_DRIP_MAX_CFS; wdr_dirt=0.0 wdr_name="Worsening + drip irrigation" wdi_cu_start=BASELINE_CU_CFS*1.85; wdi_cu_trend=+0.025*BASELINE_CU_CFS; wdi_cu_cap=BASELINE_CU_CFS*2.0 wdi_dps=SEVERE_PROB*1.60; wdi_drift=-0.100; wdi_yr0_drift=-0.10 wdi_inv=0.0; wdi_drip=0.0; wdi_dirt=STRAT_DIRT_CANAL_CFS wdi_name="Worsening + earthen canals" wca_cu_start=BASELINE_CU_CFS*1.85; wca_cu_trend=+0.025*BASELINE_CU_CFS; wca_cu_cap=BASELINE_CU_CFS*2.0 wca_dps=SEVERE_PROB*1.60; wca_drift=-0.100; wca_yr0_drift=-0.10 wca_inv=STRAT_INVASIVE_CFS; wca_drip=STRAT_DRIP_MAX_CFS; wca_dirt=STRAT_DIRT_CANAL_CFS wca_name="Worsening + all three strategies" bas_cu_start=BASELINE_CU_CFS*1.12; bas_cu_trend=+0.005*BASELINE_CU_CFS; bas_cu_cap=BASELINE_CU_CFS*2.0 bas_dps=SEVERE_PROB*1.10; bas_drift=-0.038; bas_yr0_drift=0.0 bas_inv=0.0; bas_drip=0.0; bas_dirt=0.0; bas_name="Stable (no intervention)" # Stable + individual strategies and combined sin_cu_start=BASELINE_CU_CFS*1.12; sin_cu_trend=+0.005*BASELINE_CU_CFS; sin_cu_cap=BASELINE_CU_CFS*2.0 sin_dps=SEVERE_PROB*1.10; sin_drift=-0.038; sin_yr0_drift=0.0 sin_inv=STRAT_INVASIVE_CFS; sin_drip=0.0; sin_dirt=0.0 sin_name="Stable + invasive removal/xeriscape" sdr_cu_start=BASELINE_CU_CFS*1.12; sdr_cu_trend=+0.005*BASELINE_CU_CFS; sdr_cu_cap=BASELINE_CU_CFS*2.0 sdr_dps=SEVERE_PROB*1.10; sdr_drift=-0.038; sdr_yr0_drift=0.0 sdr_inv=0.0; sdr_drip=STRAT_DRIP_MAX_CFS; sdr_dirt=0.0 sdr_name="Stable + drip irrigation" sdi_cu_start=BASELINE_CU_CFS*1.12; sdi_cu_trend=+0.005*BASELINE_CU_CFS; sdi_cu_cap=BASELINE_CU_CFS*2.0 sdi_dps=SEVERE_PROB*1.10; sdi_drift=-0.038; sdi_yr0_drift=0.0 sdi_inv=0.0; sdi_drip=0.0; sdi_dirt=STRAT_DIRT_CANAL_CFS sdi_name="Stable + earthen canals" sca_cu_start=BASELINE_CU_CFS*1.12; sca_cu_trend=+0.005*BASELINE_CU_CFS; sca_cu_cap=BASELINE_CU_CFS*2.0 sca_dps=SEVERE_PROB*1.10; sca_drift=-0.038; sca_yr0_drift=0.0 sca_inv=STRAT_INVASIVE_CFS; sca_drip=STRAT_DRIP_MAX_CFS; sca_dirt=STRAT_DIRT_CANAL_CFS sca_name="Stable + all three strategies" imp_cu_start=BASELINE_CU_CFS*1.02; imp_cu_trend=-0.005*BASELINE_CU_CFS; imp_cu_cap=BASELINE_CU_CFS*2.0 imp_dps=SEVERE_PROB*0.88; imp_drift=+0.005; imp_yr0_drift=0.0 imp_inv=0.0; imp_drip=0.0; imp_dirt=0.0; imp_name="Improving (policy)" inv_cu_start=BASELINE_CU_CFS*1.02; inv_cu_trend=-0.005*BASELINE_CU_CFS; inv_cu_cap=BASELINE_CU_CFS*2.0 inv_dps=SEVERE_PROB*0.88; inv_drift=+0.010; inv_yr0_drift=0.0 inv_inv=STRAT_INVASIVE_CFS; inv_drip=0.0; inv_dirt=0.0; inv_name="Improving + invasive removal/xeriscape" drip_cu_start=BASELINE_CU_CFS*1.02; drip_cu_trend=-0.007*BASELINE_CU_CFS; drip_cu_cap=BASELINE_CU_CFS*2.0 drip_dps=SEVERE_PROB*0.88; drip_drift=+0.010; drip_yr0_drift=0.0 drip_inv=0.0; drip_drip=STRAT_DRIP_MAX_CFS; drip_dirt=0.0; drip_name="Improving + drip irrigation" dirt_cu_start=BASELINE_CU_CFS*1.02; dirt_cu_trend=-0.004*BASELINE_CU_CFS; dirt_cu_cap=BASELINE_CU_CFS*2.0 dirt_dps=SEVERE_PROB*0.88; dirt_drift=+0.010; dirt_yr0_drift=0.0 dirt_inv=0.0; dirt_drip=0.0; dirt_dirt=STRAT_DIRT_CANAL_CFS; dirt_name="Improving + earthen canals" all_cu_start=BASELINE_CU_CFS*0.92; all_cu_trend=-0.010*BASELINE_CU_CFS; all_cu_cap=BASELINE_CU_CFS*2.0 all_dps=SEVERE_PROB*0.80; all_drift=+0.020; all_yr0_drift=0.0 all_inv=STRAT_INVASIVE_CFS; all_drip=STRAT_DRIP_MAX_CFS; all_dirt=STRAT_DIRT_CANAL_CFS all_name="Improving + all three strategies" MONTH_NAMES = ["Jan","Feb","Mar","Apr","May","Jun", "Jul","Aug","Sep","Oct","Nov","Dec"] def run_scenario(sname, cu_start, cu_trend=0.0, cu_cap=None, inv_cfs=0.0, drip_cfs=0.0, dirt_cfs=0.0, drought_prob_sev=None, net_drift=0.0, val_replay=False, val_dry_after_month=0, seed_log_override=None, yr0_drift=0.0): """val_dry_after_month: in validation mode, ignore dryness in months before this index (0-based) so pre-event warm-up doesn't count. seed_log_override: per-scenario AR seed (list of 12 log values); defaults to the global AR_SEED_LOG if None. yr0_drift: additional log-space offset applied at yr=0 only, so pessimistic scenarios feel immediate pressure (normal drift is 0 at yr=0 because it multiplies the year index).""" if cu_cap is None: cu_cap = BASELINE_CU_CFS * 2.0 if drought_prob_sev is None: drought_prob_sev = SEVERE_PROB seed_log = seed_log_override if seed_log_override is not None else AR_SEED_LOG sp_eff = drought_prob_sev dp_eff = DROUGHT_PROB * (drought_prob_sev / SEVERE_PROB) all_res = [] for _ in range(num_simulations): went_dry = False; stage_hist = []; dry_yrs = []; dry_months = [] log_ar = list(seed_log) ann_state = [] s = 0.0 for _ in range(years_to_predict): s = ANNUAL_STATE_PHI * s + random.gauss(0, ANNUAL_STATE_STD) ann_state.append(s) ext_yrs = [] for y in range(years_to_predict): prev_bad = any(ey==y-1 and et in ("severe","drought") for et,ey in ext_yrs) sp = sp_eff * SEVERE_AUTOCORR if prev_bad else sp_eff dp = dp_eff * DROUGHT_AUTOCORR if prev_bad else dp_eff r = random.random() if r < sp: ext_yrs.append(("severe", y)) elif r < sp+dp: ext_yrs.append(("drought", y)) elif r < sp+dp+WET_PROB: ext_yrs.append(("wet", y)) for month in range(years_to_predict * 12): moy = month % 12 yr = month // 12 ayr = start_year + yr lags = log_ar[-12:][::-1] log_nat = sum(AR_PHI[i] * lags[i] for i in range(12)) if not val_replay: drift_capped = max(-0.10, min(0.10, net_drift * yr)) log_nat += drift_capped if yr == 0: log_nat += yr0_drift log_nat += ann_state[yr] shock = random.gauss(0, AR_SCENARIO_STD) shock = max(-AR_SHOCK_CLIP * AR_SCENARIO_STD, min( AR_SHOCK_CLIP * AR_SCENARIO_STD, shock)) log_nat += shock log_ar.append(log_nat) nat_flow = max(0.0, math.exp(log_nat)) if val_replay: nat_flow *= REPLAY_MULT[moy] else: for et, ey in ext_yrs: if yr == ey: in_pk = moy in DROUGHT_MONTHS if et == "severe": m = SEVERE_MULT if in_pk else (1.0 - OFFPEAK_BLEND*(1.0-SEVERE_MULT)) elif et == "drought": m = DROUGHT_MULT if in_pk else (1.0 - OFFPEAK_BLEND*(1.0-DROUGHT_MULT)) else: m = WET_MULT if in_pk else (1.0 + OFFPEAK_BLEND*(WET_MULT-1.0)) nat_flow *= m if cu_trend > 0: lcu = cu_start + cu_trend * yr frac = max(0.0, min(1.0, (lcu - BASELINE_CU_CFS) / max(cu_cap - BASELINE_CU_CFS, 1.0))) lg = 1.0 / (1.0 + math.exp(-6.0 * (frac - 0.5))) cu_yr = BASELINE_CU_CFS + (cu_cap - BASELINE_CU_CFS) * lg else: cu_yr = max(BASELINE_CU_CFS * 0.50, cu_start + cu_trend * yr) cu_m = cu_yr * CU_SEASONAL[moy] cu_b = BASELINE_CU_CFS * CU_SEASONAL[moy] dd = cu_m - cu_b inv_b = inv_cfs * CU_SEASONAL[moy] drip_b = drip_cfs * min(1.0, yr/10.0) * CU_SEASONAL[moy] dirt_b = dirt_cfs * CU_SEASONAL[moy] net_flow = max(0.0, nat_flow - dd + inv_b + drip_b + dirt_b) net_flow = min(net_flow, 15000.0) stage_hist.append(flow_to_stage(net_flow)) if net_flow < DRY_THRESHOLD: if not val_replay or month >= val_dry_after_month: went_dry = True if ayr not in dry_yrs: dry_yrs.append(ayr) label = f"{MONTH_NAMES[moy]} {ayr}" if label not in dry_months: dry_months.append(label) mn = min(stage_hist); mx = max(stage_hist); av = sum(stage_hist)/len(stage_hist) all_res.append((went_dry, mn, av, mx, dry_yrs, dry_months)) return all_res def aggregate(results): dry, mins, avgs, maxs = 0, [], [], [] ycts = {}; mcts = {} for wd, mn, av, mx, dys, dms in results: if wd: dry += 1 mins.append(mn); avgs.append(av); maxs.append(mx) for y in dys: ycts[y] = ycts.get(y, 0) + 1 for m in dms: mcts[m] = mcts.get(m, 0) + 1 n = len(results) ycts_sd = {} for y, cnt in ycts.items(): p = cnt/n; ycts_sd[y] = math.sqrt(p*(1-p)/n)*100 po = dry/n; ov_sd = math.sqrt(po*(1-po)/n)*100 return dry/n*100, sum(mins)/n, sum(avgs)/n, sum(maxs)/n, ycts, ycts_sd, ov_sd, mcts BLUE="\033[34m"; CYAN="\033[36m"; GREEN="\033[32m" YELLOW="\033[33m"; RED="\033[31m"; ORANGE="\033[38;5;208m" DIM="\033[2m"; BOLD="\033[1m"; RESET="\033[0m" WAVE="~"*68; LINE="-"*68 def _col(p): if p < 40: return GREEN if p < 62: return YELLOW if p < 80: return ORANGE return RED def _bar(p, w=32): f = round(p/100*w); c = _col(p) return c + "█"*f + RESET + DIM + "░"*(w-f) + RESET def _pct(p): c = _col(p); s = f"{p:5.1f}%" return f"{RED}{BOLD}{s}{RESET}" if p >= 99.9 else f"{c}{s}{RESET}" def print_scenario(name, p, sd, mn, av, mx, col): print(f" {col}{BOLD}{name:<42}{RESET} {_bar(p)} {_pct(p)} {DIM}±{sd:.1f}%{RESET}") print(f" {DIM}{'':42} stage min {mn:.2f} ft avg {av:.2f} ft max {mx:.2f} ft{RESET}\n") def print_year_breakdown(label, ycts, ycts_sd, col): print(f" {col}{BOLD}{label}{RESET}") if not ycts: print(f" {DIM} (no years went dry){RESET}\n"); return for yr in sorted(ycts): p = ycts[yr]/num_simulations*100 sd = ycts_sd.get(yr, 0.0) print(f" {DIM} {yr}{RESET} {_bar(p, 36)} {_pct(p)} {DIM}±{sd:.1f}%{RESET}") print() lf = [math.log(max(q,1)) for q in ABQ_MONTHLY] y_hat = [sum(AR_PHI[k]*lf[i-k-1] for k in range(12)) for i in range(12, len(lf))] y_obs = lf[12:] ss_r = sum((y_hat[i]-y_obs[i])**2 for i in range(len(y_obs))) ss_t = sum((v-sum(y_obs)/len(y_obs))**2 for v in y_obs) ar_r2 = 1.0 - ss_r/ss_t print(f"\n{BLUE}{WAVE}{RESET}") print(f" {BOLD}{CYAN}RIO GRANDd REACH FLOW - ABQ{RESET}") print(f"{BLUE}{WAVE}{RESET}") print(f" {DIM}Gauge USGS 08330000 Rio Grande at Albuquerque{RESET}") print(f" {DIM}Training 1991-2020 · 360 months · AR(12) R² = {ar_r2:.4f} · σ = {AR_RESID_STD:.4f} (log){RESET}") print(f" {DIM}Dry threshold {DRY_THRESHOLD} CFS · Baseline CU {BASELINE_CU_CFS:.1f} CFS{RESET}") print(f" {DIM}Drought {SEVERE_PROB*100:.0f}% severe / {DROUGHT_PROB*100:.0f}% moderate [NM PDSI 1895-2020]{RESET}") print(f" {DIM}Annual state AR(1) φ={ANNUAL_STATE_PHI} σ={ANNUAL_STATE_STD}{RESET}") if validation_mode: print(f" {YELLOW}{BOLD}MODE: VALIDATION{RESET}") else: print(f" {CYAN}{BOLD}MODE: PREDICTION {start_year}-{start_year+years_to_predict-1}" f" | {num_simulations:,} #sims | {years_to_predict}-yr horizon{RESET}") print(f"{DIM}{LINE}{RESET}\n") if validation_mode: print(f" {DIM}Running validation ({num_simulations} sims)...{RESET}") # val_dry_after_month=0 means all 12 months of 2022 are in scope. # Raise this (e.g. to 5) to exclude Jan–May warm-up if needed. val_r = run_scenario("2022", BASELINE_CU_CFS, val_replay=True, val_dry_after_month=0) print(f" {GREEN}Done.{RESET}\n") else: scens = [ (wor_name, wor_cu_start, wor_cu_trend, wor_cu_cap, wor_inv, 0.0, 0.0, wor_dps, wor_drift, WOR_SEED_LOG, wor_yr0_drift), (win_name, win_cu_start, win_cu_trend, win_cu_cap, win_inv, 0.0, 0.0, win_dps, win_drift, WOR_SEED_LOG, win_yr0_drift), (wdr_name, wdr_cu_start, wdr_cu_trend, wdr_cu_cap, 0.0, wdr_drip, 0.0, wdr_dps, wdr_drift, WOR_SEED_LOG, wdr_yr0_drift), (wdi_name, wdi_cu_start, wdi_cu_trend, wdi_cu_cap, 0.0, 0.0, wdi_dirt, wdi_dps, wdi_drift, WOR_SEED_LOG, wdi_yr0_drift), (wca_name, wca_cu_start, wca_cu_trend, wca_cu_cap, wca_inv, wca_drip, wca_dirt, wca_dps, wca_drift, WOR_SEED_LOG, wca_yr0_drift), (bas_name, bas_cu_start, bas_cu_trend, bas_cu_cap, bas_inv, 0.0, 0.0, bas_dps, bas_drift, None, bas_yr0_drift), (sin_name, sin_cu_start, sin_cu_trend, sin_cu_cap, sin_inv, 0.0, 0.0, sin_dps, sin_drift, None, sin_yr0_drift), (sdr_name, sdr_cu_start, sdr_cu_trend, sdr_cu_cap, 0.0, sdr_drip, 0.0, sdr_dps, sdr_drift, None, sdr_yr0_drift), (sdi_name, sdi_cu_start, sdi_cu_trend, sdi_cu_cap, 0.0, 0.0, sdi_dirt, sdi_dps, sdi_drift, None, sdi_yr0_drift), (sca_name, sca_cu_start, sca_cu_trend, sca_cu_cap, sca_inv, sca_drip, sca_dirt, sca_dps, sca_drift, None, sca_yr0_drift), (imp_name, imp_cu_start, imp_cu_trend, imp_cu_cap, imp_inv, 0.0, 0.0, imp_dps, imp_drift, None, imp_yr0_drift), (inv_name, inv_cu_start, inv_cu_trend, inv_cu_cap, inv_inv, inv_drip, inv_dirt, inv_dps, inv_drift, None, inv_yr0_drift), (drip_name,drip_cu_start,drip_cu_trend,drip_cu_cap,drip_inv, drip_drip, drip_dirt, drip_dps,drip_drift,None, drip_yr0_drift), (dirt_name,dirt_cu_start,dirt_cu_trend,dirt_cu_cap,dirt_inv, dirt_drip, dirt_dirt, dirt_dps,dirt_drift,None, dirt_yr0_drift), (all_name, all_cu_start, all_cu_trend, all_cu_cap, all_inv, all_drip, all_dirt, all_dps, all_drift, None, all_yr0_drift), ] results = [] for (sn, cus, cut, cuc, inv, drip, dirt, dps, drift, seed, y0d) in scens: print(f" {DIM}Running {sn}...{RESET}", end="", flush=True) results.append(run_scenario(sn, cus, cut, cuc, inv, drip, dirt, dps, drift, seed_log_override=seed, yr0_drift=y0d)) print(f" {GREEN}done{RESET}") print() (wor_r, win_r, wdr_r, wdi_r, wca_r, bas_r, sin_r, sdr_r, sdi_r, sca_r, imp_r, inv_r, drp_r, drt_r, all_r) = results print(f"{BLUE}{WAVE}{RESET}") print(f" {BOLD}{CYAN}RESULTS{RESET}") print(f"{BLUE}{WAVE}{RESET}\n") if validation_mode: vp, vmn, vav, vmx, vyrs, vysd, vsd, vmcts = aggregate(val_r) flag = (GREEN+BOLD+"CORRECT"+RESET) if vp > 50 else (RED+BOLD+"LOW"+RESET) print(f" Model result: {vp:.1f}% ±{vsd:.1f}% probability of going dry {flag}") print(f" {DIM}(counts only months at or after the prediction window start — " f"pre-event warm-up excluded){RESET}\n") if vmcts: print(f" {BOLD}Dry month breakdown " f"{DIM}(% of {num_simulations} sims going dry that month):{RESET}") for label in sorted(vmcts, key=lambda s: (int(s[-4:]), MONTH_NAMES.index(s[:3]))): cnt = vmcts[label] p = cnt / num_simulations * 100 sd = math.sqrt((p/100)*(1-p/100)/num_simulations)*100 print(f" {DIM}{label:<10}{RESET} {_bar(p, 36)} {_pct(p)} {DIM}±{sd:.1f}%{RESET}") else: print(f" {DIM} (no months went dry in any simulation){RESET}") else: wp, wmn, wav, wmx, wyrs, wysd, wsd, _ = aggregate(wor_r) wip, wimn, wiav, wimx, wiyrs, wiysd, wisd, _ = aggregate(win_r) wdp, wdmn, wdav, wdmx, wdyrs, wdysd, wdsd, _ = aggregate(wdr_r) wdip,wdimn,wdiav,wdimx,wdiyrs,wdiysd,wdisd,_ = aggregate(wdi_r) wcp, wcmn, wcav, wcmx, wcyrs, wcysd, wcsd, _ = aggregate(wca_r) bp, bmn, bav, bmx, byrs, bysd, bsd, _ = aggregate(bas_r) sip, simn, siav, simx, siyrs, siysd, sisd, _ = aggregate(sin_r) sdp, sdmn, sdav, sdmx, sdyrs, sdysd, sdsd, _ = aggregate(sdr_r) sdip,sdimn,sdiav,sdimx,sdiyrs,sdiysd,sdisd,_ = aggregate(sdi_r) scp, scmn, scav, scmx, scyrs, scysd, scsd, _ = aggregate(sca_r) pp, pmn, pav, pmx, pyrs, pysd, psd, _ = aggregate(imp_r) ip, imn, iav, imx, iyrs, iysd, isd, _ = aggregate(inv_r) dp, dmn, dav, dmx, dyrs, dysd, dsd, _ = aggregate(drp_r) rp, rmn, rav, rmx, ryrs, rysd, rsd, _ = aggregate(drt_r) ap, amn, aav, amx, ayrs, aysd, asd, _ = aggregate(all_r) import csv, os, sys csv_path = os.path.join(os.path.dirname(os.path.abspath(sys.argv[0])), f"rio_grande_model_v10_{start_year}.csv") scen_data = [ (wor_name, wp, wsd, wmn, wav, wmx, wyrs, wysd, wor_drift*100, "Worsening"), (win_name, wip, wisd, wimn, wiav, wimx, wiyrs, wiysd, win_drift*100, "Worsening+Invasive"), (wdr_name, wdp, wdsd, wdmn, wdav, wdmx, wdyrs, wdysd, wdr_drift*100, "Worsening+Drip"), (wdi_name, wdip,wdisd,wdimn,wdiav,wdimx,wdiyrs,wdiysd,wdi_drift*100, "Worsening+Canal"), (wca_name, wcp, wcsd, wcmn, wcav, wcmx, wcyrs, wcysd, wca_drift*100, "Worsening+AllThree"), (bas_name, bp, bsd, bmn, bav, bmx, byrs, bysd, bas_drift*100, "Stable"), (sin_name, sip, sisd, simn, siav, simx, siyrs, siysd, sin_drift*100, "Stable+Invasive"), (sdr_name, sdp, sdsd, sdmn, sdav, sdmx, sdyrs, sdysd, sdr_drift*100, "Stable+Drip"), (sdi_name, sdip,sdisd,sdimn,sdiav,sdimx,sdiyrs,sdiysd,sdi_drift*100, "Stable+Canal"), (sca_name, scp, scsd, scmn, scav, scmx, scyrs, scysd, sca_drift*100, "Stable+AllThree"), (imp_name, pp, psd, pmn, pav, pmx, pyrs, pysd, imp_drift*100, "Improving"), (inv_name, ip, isd, imn, iav, imx, iyrs, iysd, inv_drift*100, "Improving+Invasive"), (drip_name,dp, dsd, dmn, dav, dmx, dyrs, dysd, drip_drift*100, "Improving+Drip"), (dirt_name,rp, rsd, rmn, rav, rmx, ryrs, rysd, dirt_drift*100, "Improving+Canal"), (all_name, ap, asd, amn, aav, amx, ayrs, aysd, all_drift*100, "Improving+AllThree"), ] all_years = list(range(start_year, start_year + years_to_predict)) with open(csv_path, "w", newline="") as f: w = csv.writer(f) w.writerow(["=== SHEET 1: OVERALL SCENARIO SUMMARY ==="]) w.writerow(["Scenario", "Flow_Added_CFS", "Net_Drift_pct_per_yr", "P_dry_overall_pct", "SD_overall_pct", "Stage_min_ft", "Stage_avg_ft", "Stage_max_ft", "vs_Stable_pp"]) tc = STRAT_INVASIVE_CFS + STRAT_DRIP_MAX_CFS + STRAT_DIRT_CANAL_CFS flow_map = { wor_name: "0", win_name: f"+{STRAT_INVASIVE_CFS:.0f}", wdr_name: f"+{STRAT_DRIP_MAX_CFS:.0f} (phased 10yr)", wdi_name: f"+{STRAT_DIRT_CANAL_CFS:.0f}", wca_name: f"+{tc:.0f}", bas_name: "0", sin_name: f"+{STRAT_INVASIVE_CFS:.0f}", sdr_name: f"+{STRAT_DRIP_MAX_CFS:.0f} (phased 10yr)", sdi_name: f"+{STRAT_DIRT_CANAL_CFS:.0f}", sca_name: f"+{tc:.0f}", imp_name: "0 (policy only)", inv_name: f"+{STRAT_INVASIVE_CFS:.0f}", drip_name: f"+{STRAT_DRIP_MAX_CFS:.0f} (phased 10yr)", dirt_name: f"+{STRAT_DIRT_CANAL_CFS:.0f}", all_name: f"+{tc:.0f}", } for (sn, prob, sd, smn, sav, smx, ycts, ysd, drift_pct, _) in scen_data: delta = prob - bp if sn != bas_name else 0.0 w.writerow([sn, flow_map[sn], f"{drift_pct:+.1f}", f"{prob:.2f}", f"{sd:.2f}", f"{smn:.3f}", f"{sav:.3f}", f"{smx:.3f}", f"{delta:+.2f}" if sn != bas_name else "baseline"]) w.writerow([]) w.writerow(["=== SHEET 2: PER-YEAR P(DRY) BY SCENARIO ==="]) w.writerow(["Year"] + [col for sn,*_ in scen_data for col in [f"{sn}_P_dry_pct", f"{sn}_SD_pct"]]) for yr in all_years: row = [yr] for (sn, prob, sd, smn, sav, smx, ycts, ysd, drift_pct, _) in scen_data: p = ycts.get(yr, 0) / num_simulations * 100 s = ysd.get(yr, 0.0) row += [f"{p:.2f}", f"{s:.2f}"] w.writerow(row) w.writerow([]) w.writerow(["=== SHEET 3: STRATEGY EFFECT AT KEY YEARS ==="]) yr1 = start_year yr5 = start_year + 4 yr10 = start_year + years_to_predict - 1 w.writerow(["Scenario", "Flow_Added_CFS", f"P_dry_{yr1}_pct", f"SD_{yr1}_pct", f"P_dry_{yr5}_pct", f"SD_{yr5}_pct", f"P_dry_{yr10}_pct", f"SD_{yr10}_pct", f"vs_Stable_at_{yr10}_pp"]) for (sn, prob, sd, smn, sav, smx, ycts, ysd, drift_pct, _) in scen_data: p1 = ycts.get(yr1, 0) / num_simulations * 100 s1 = ysd.get(yr1, 0.0) p5 = ycts.get(yr5, 0) / num_simulations * 100 s5 = ysd.get(yr5, 0.0) p10 = ycts.get(yr10, 0) / num_simulations * 100 s10 = ysd.get(yr10, 0.0) bp10 = byrs.get(yr10, 0) / num_simulations * 100 delta10 = p10 - bp10 if sn != bas_name else 0.0 w.writerow([sn, flow_map[sn], f"{p1:.2f}", f"{s1:.2f}", f"{p5:.2f}", f"{s5:.2f}", f"{p10:.2f}", f"{s10:.2f}", f"{delta10:+.2f}" if sn != bas_name else "baseline"]) w.writerow([]) w.writerow(["=== SHEET 4: MODEL PARAMETERS ==="]) w.writerow(["Parameter", "Value", "Description"]) params = [ ("start_year", start_year, "First year of projection"), ("years_to_predict", years_to_predict, "Projection horizon"), ("num_simulations", num_simulations, "Monte Carlo runs per scenario"), ("dry_threshold_CFS", DRY_THRESHOLD, "CFS below which river is 'dry'"), ("baseline_CU_CFS", BASELINE_CU_CFS, "Annual avg consumptive use, Cochiti-ABQ reach [DRP-045]"), ("AR_order", 12, "Autoregressive model order (monthly)"), ("AR_R2", f"{ar_r2:.4f}", "AR(12) goodness of fit on 1991-2020 gauge data"), ("AR_residual_sigma", f"{AR_RESID_STD:.4f}", "AR(12) residual std dev (log space)"), ("AR_scenario_std", f"{AR_SCENARIO_STD:.4f}", "Monthly shock std used in simulation"), ("AR_shock_clip", AR_SHOCK_CLIP, "Shock clip in sigma units"), ("annual_state_phi", ANNUAL_STATE_PHI, "Inter-annual AR(1) autocorrelation"), ("annual_state_sigma",ANNUAL_STATE_STD, "Inter-annual state std dev (log space)"), ("seed_frac", CURRENT_STATE_FRAC, "Seed deflation vs 2010-2020 mean (current megadrought)"), ("severe_drought_prob",SEVERE_PROB, "Annual prob of D3-D4 severe drought [NM PDSI]"), ("moderate_drought_prob",DROUGHT_PROB, "Annual prob of D2 moderate drought"), ("wet_year_prob", WET_PROB, "Annual prob of above-normal year"), ("severe_mult", SEVERE_MULT, "Flow multiplier in severe drought peak months"), ("drought_mult", DROUGHT_MULT, "Flow multiplier in moderate drought peak months"), ("wet_mult", WET_MULT, "Flow multiplier in wet year peak months"), ("invasive_CFS", STRAT_INVASIVE_CFS, "CFS returned by invasive removal + xeriscape [MRGCD/ABQ Water 2023]"), ("drip_max_CFS", STRAT_DRIP_MAX_CFS, "Max CFS saved by drip irrigation, phased over 10 yr [DRP-045; NRCS]"), ("earthen_canal_CFS", STRAT_DIRT_CANAL_CFS, "CFS net baseflow gain from earthen canals [USBR]"), ] for p, v, d in params: w.writerow([p, v, d]) print(f" {GREEN}CSV saved: {csv_path}{RESET}\n") print(f" {BOLD}Probability river goes dry at least once, " f"{start_year}\u2013{start_year+years_to_predict-1}:{RESET}") print(f" {DIM}(\u00b1SD = binomial SD across {num_simulations:,} simulations){RESET}\n") print_scenario(wor_name, wp, wsd, wmn, wav, wmx, RED) print_scenario(win_name, wip, wisd, wimn, wiav, wimx, YELLOW) print_scenario(wdr_name, wdp, wdsd, wdmn, wdav, wdmx, YELLOW) print_scenario(wdi_name, wdip,wdisd,wdimn,wdiav,wdimx,YELLOW) print_scenario(wca_name, wcp, wcsd, wcmn, wcav, wcmx, YELLOW) print_scenario(bas_name, bp, bsd, bmn, bav, bmx, RED) print_scenario(sin_name, sip, sisd, simn, siav, simx, YELLOW) print_scenario(sdr_name, sdp, sdsd, sdmn, sdav, sdmx, YELLOW) print_scenario(sdi_name, sdip,sdisd,sdimn,sdiav,sdimx,YELLOW) print_scenario(sca_name, scp, scsd, scmn, scav, scmx, YELLOW) print_scenario(imp_name, pp, psd, pmn, pav, pmx, YELLOW) print_scenario(inv_name, ip, isd, imn, iav, imx, YELLOW) print_scenario(drip_name, dp, dsd, dmn, dav, dmx, YELLOW) print_scenario(dirt_name, rp, rsd, rmn, rav, rmx, YELLOW) print_scenario(all_name, ap, asd, amn, aav, amx, GREEN) print(f"{DIM}{LINE}{RESET}") print(f"\n {BOLD}{CYAN}STRATEGY EFFECT TABLE{RESET}" f" {DIM}(strategies built on Improving base \u00b7 vs Stable baseline){RESET}\n") NW = 42 print(f" {BOLD}{DIM}{'Scenario':<{NW}} {'Flow added':>11} {'P(dry)':>7} {'\u00b1SD':>5} " f"{'vs Stable':>10} {'Avg stage':>9}{RESET}") print(f" {DIM}{'-'*88}{RESET}") tc = STRAT_INVASIVE_CFS+STRAT_DRIP_MAX_CFS+STRAT_DIRT_CANAL_CFS trows = [ (wor_name, "—", wp, wsd, wp-bp, wav), (win_name, f"+{STRAT_INVASIVE_CFS:.0f} CFS", wip, wisd, wip-bp, wiav), (wdr_name, f"+{STRAT_DRIP_MAX_CFS:.0f} CFS", wdp, wdsd, wdp-bp, wdav), (wdi_name, f"+{STRAT_DIRT_CANAL_CFS:.0f} CFS", wdip,wdisd,wdip-bp,wdiav), (wca_name, f"+{tc:.0f} CFS", wcp, wcsd, wcp-bp, wcav), (bas_name, "—", bp, bsd, 0.0, bav), (sin_name, f"+{STRAT_INVASIVE_CFS:.0f} CFS", sip, sisd, sip-bp, siav), (sdr_name, f"+{STRAT_DRIP_MAX_CFS:.0f} CFS", sdp, sdsd, sdp-bp, sdav), (sdi_name, f"+{STRAT_DIRT_CANAL_CFS:.0f} CFS", sdip,sdisd,sdip-bp,sdiav), (sca_name, f"+{tc:.0f} CFS", scp, scsd, scp-bp, scav), (imp_name, "policy only", pp, psd, pp-bp, pav), (inv_name, f"+{STRAT_INVASIVE_CFS:.0f} CFS", ip, isd, ip-bp, iav), (drip_name, f"+{STRAT_DRIP_MAX_CFS:.0f} CFS", dp, dsd, dp-bp, dav), (dirt_name, f"+{STRAT_DIRT_CANAL_CFS:.0f} CFS", rp, rsd, rp-bp, rav), (all_name, f"+{tc:.0f} CFS", ap, asd, ap-bp, aav), ] for sn, fl, prob, sd, delta, stage in trows: pc = _col(prob) if sn == bas_name: ds = f"{'baseline':>10}"; dc = DIM else: sign = "+" if delta > 0 else "" ds = f"{sign}{delta:.1f} pp" dc = GREEN if delta < -5 else (YELLOW if delta < 0 else RED) print(f" {BOLD}{sn:<{NW}}{RESET} {DIM}{fl:>11}{RESET}" f" {pc}{prob:>6.1f}%{RESET} {DIM}\u00b1{sd:>4.1f}%{RESET}" f" {dc}{ds:>10}{RESET} {DIM}{stage:>7.2f} ft{RESET}") print(f"\n {DIM}pp = \u00b1 percentage-point change vs Stable \u00b7 SD = binomial SD{RESET}") print(f" {DIM}Invasive/xerisc : tamarisk ~4,000 bosque acres (+19 CFS) + ABQ xeriscape (+4 CFS) [MRGCD/ABQ Water 2023]{RESET}") print(f" {DIM}Drip irrigation : flood 55%% -> drip 90%% efficiency, phased over 10 yr [DRP-045; USDA NRCS]{RESET}") print(f" {DIM}Earthen canals : concrete->dirt seepage recharges aquifer, ~50-60%% resurfaces [USBR]{RESET}") print(f" {DIM}Combined total : +{tc:.0f} CFS \u00b7 Rio Grande ran dry in ABQ 3\u00d7 in 4 yr [MRGCD/USGS 2025]{RESET}") print(f"\n{DIM}{LINE}{RESET}") print(f"\n {BOLD}{CYAN}WHICH YEARS FACE THE HIGHEST RISK?{RESET}" f" {DIM}(% of {num_simulations:,} sims going dry that year \u00b7 \u00b1SD shown){RESET}\n") print_year_breakdown(wor_name, wyrs, wysd, RED) print_year_breakdown(win_name, wiyrs, wiysd, YELLOW) print_year_breakdown(wdr_name, wdyrs, wdysd, YELLOW) print_year_breakdown(wdi_name, wdiyrs,wdiysd,YELLOW) print_year_breakdown(wca_name, wcyrs, wcysd, YELLOW) print_year_breakdown(bas_name, byrs, bysd, RED) print_year_breakdown(sin_name, siyrs, siysd, YELLOW) print_year_breakdown(sdr_name, sdyrs, sdysd, YELLOW) print_year_breakdown(sdi_name, sdiyrs,sdiysd,YELLOW) print_year_breakdown(sca_name, scyrs, scysd, YELLOW) print_year_breakdown(imp_name, pyrs, pysd, YELLOW) print_year_breakdown(inv_name, iyrs, iysd, YELLOW) print_year_breakdown(drip_name, dyrs, dysd, YELLOW) print_year_breakdown(dirt_name, ryrs, rysd, YELLOW) print_year_breakdown(all_name, ayrs, aysd, GREEN) print(f"\n{BLUE}{WAVE}{RESET}\n") except Exception as e: print(f"\n\033[31m{'='*60}\033[0m") print(f"\033[31mERROR: {e}\033[0m") traceback.print_exc() print(f"\033[31m{'='*60}\033[0m\n") try: input(" Press Enter to close...") except (EOFError, KeyboardInterrupt): pass