"""Compare the v7 hyperparameter sweep (median member ens010, shared ens148
precip from 1850) with round 4, for the talk of 2026-10-28.

Metrics (comparison rules of 2026-10-23, see CLAUDE.md), on the glaciers
with output in every run (11-01 only):
- LIA area: sum of the per-glacier max of area_min_h over 1800-1875
  (Reinthaler & Paul 2025: 4211 ± 215 km²);
- LIA volume: per glacier, the volume in the year of its max area;
- dV LIA -> 2015 (Reinthaler: −180 ± 39 km³), area ratio LIA / 2015
  (4211 / 1806 = 2.33);
- "anchored" reference LIA volume for the plots: V_model(2003) + 180 −
  |observed dV 2003 -> 2015| (Hugonnet rate x 12 years);
- present day: regional 2000–2020 mismatch against Hugonnet (kg m-2 yr-1),
  V2003 / V_inv, A2003 / A_RGI, dA 2003 -> 2015 (−13.5% observed);
- % area within the calibration tolerance (dyn_calib_reg11.csv).

    python sweep_compare.py  ->  figs/sweep_compare.csv
"""
import glob
import os
import numpy as np
import pandas as pd
import xarray as xr

R = '/home/www/fmaussion/runs/lmr_2025'
HERE = os.path.dirname(os.path.abspath(__file__))
REF_AREA, REF_DV, REF_A2015 = 4211., -180., 1806.

runs = {}
for d in sorted(glob.glob(f'{HERE}/output/sw_*_ys1500')):
    exp = os.path.basename(d)[:-len('_ys1500')]
    f = f'{d}/hist_dyn_calib_reg11.nc'
    if os.path.exists(f):
        runs[exp] = (f, f'{d}/dyn_calib_reg11.csv',
                     f'{R}/lmr_gdirs_v7_202610/exps/{exp}/RGI62/b_240/L3/summary/glacier_statistics_11.csv')
runs['r4 lam8_jjas_dyn (ens183)'] = (
    f'{R}/round4/sensitivity/output/lam8_jjas_dyn/lmr_mira_hist_reg11.nc',
    f'{R}/round4/sensitivity/output/lam8_jjas_dyn/dyn_calib_reg11.csv',
    f'{R}/lmr_gdirs_v5_202610/exps/lam8_jjas/RGI62/b_160/L3/summary/glacier_statistics_11.csv')
print(f'{len(runs)} runs')

D = {}
common = None
for k, (f, c, s) in runs.items():
    with xr.open_dataset(f) as ds:
        ds = ds[['volume', 'area_min_h']].sel(time=slice(1800, 2020)).load()
    ok = ds.volume.sel(time=[1800, 2000, 2003, 2015, 2020]).notnull().all('time').to_pandas()
    ok = ok[ok].index
    ok = ok[ok.str.startswith('RGI60-11.0')]
    common = ok if common is None else common.intersection(ok)
    D[k] = (ds, pd.read_csv(c, index_col=0), pd.read_csv(s, index_col=0, low_memory=False))

st0 = D[next(iter(D))][2]
A_rgi = st0.rgi_area_km2.reindex(common) * 1e6
all_area = st0.loc[st0.index.str.startswith('RGI60-11.0'), 'rgi_area_km2'].sum()
print(f'common set: {len(common)} glaciers, {100 * A_rgi.sum() / 1e6 / all_area:.1f}% of the 11-01 area')

rows = []
for k, (ds, cal, st) in D.items():
    ds = ds.sel(rgi_id=common)
    lia = ds.sel(time=slice(1800, 1875))
    imax = lia.area_min_h.argmax('time')
    a_lia = float(lia.area_min_h.isel(time=imax).sum()) / 1e6
    v_lia = float(lia.volume.isel(time=imax).sum()) / 1e9
    v = ds.volume.sum('rgi_id').to_series() / 1e9
    a = ds.area_min_h.sum('rgi_id').to_series() / 1e6
    ref = cal.dmdtda_mismatch_dynamic_calibration_reference.reindex(common)
    tol = cal.dmdtda_dynamic_calibration_given_error.reindex(common)
    m = 900 * ds.volume
    dm = ((m.sel(time=2020) - m.sel(time=2000)) / A_rgi.to_xarray() / 20).to_pandas()
    mis = dm - ref
    ok = ref.notna()
    target_reg = np.average(ref[ok], weights=A_rgi[ok])
    dv_obs_03_15 = (ref[ok] * A_rgi[ok]).sum() * 12 / 900 / 1e9  # km³ ice, < 0
    inv = st.inv_volume_km3.reindex(common).sum()
    rows.append({
        'run': k,
        'LIA area': a_lia,
        'LIA vol': v_lia,
        'dV LIA-2015': v[2015] - v_lia,
        'area ratio LIA/2015': a_lia / a[2015],
        'V2003': v[2003],
        'anchored ref LIA vol': v[2003] - REF_DV + dv_obs_03_15,
        'mismatch 00-20': np.average(dm[ok], weights=A_rgi[ok]) - target_reg,
        'V2003/Vinv': v[2003] / inv,
        'A2003/Argi': a[2003] / (A_rgi.sum() / 1e6),
        'dA 03-15 %': 100 * (a[2015] / a[2003] - 1),
        '% area within tol': 100 * A_rgi[(mis.abs() <= tol)].sum() / A_rgi.sum(),
    })

df = pd.DataFrame(rows).set_index('run')
df['dLIA area vs ref %'] = 100 * (df['LIA area'] / REF_AREA - 1)
df['dV vs ref %'] = 100 * (df['dV LIA-2015'] / REF_DV - 1)
os.makedirs(f'{HERE}/figs', exist_ok=True)
df.to_csv(f'{HERE}/figs/sweep_compare.csv')
pd.set_option('display.width', 250)
print(f'Reinthaler: LIA area {REF_AREA:.0f} ± 215 km², dV LIA-2015 {REF_DV:.0f} ± 39 km³, '
      f'area ratio {REF_AREA / REF_A2015:.2f}; dA 2003-2015 −13.5%')
print(df.sort_values('LIA area', ascending=False).round(2).to_string())