From 45b2ffd70854cdcd2bff0214a02d8aba93beb9a2 Mon Sep 17 00:00:00 2001
From: Burcin Cakir Erdener <bcakire@kl2.head.cm.kestrel.hpc.nlr.gov>
Date: Fri, 24 Jul 2026 13:16:21 -0600
Subject: [PATCH 01/11] Add CVAR and NCVAR resource adequacy reporting

---
 cases.csv                                 |   3 +-
 cases_test.csv                            |   4 +-
 reeds/resource_adequacy/ra_calcs.py       |  15 +++
 reeds/resource_adequacy/run_pras.jl       |  91 ++++++++++++++----
 reeds/resource_adequacy/stress_periods.py | 110 +++++++++++++++++++---
 runreeds.py                               |  32 ++++---
 6 files changed, 209 insertions(+), 46 deletions(-)

diff --git a/cases.csv b/cases.csv
index 525a78bc7..3b302256d 100644
--- a/cases.csv
+++ b/cases.csv
@@ -273,13 +273,14 @@ GSw_PRM_StressOutages,Whether to apply the availability factor (forced + schedul
 GSw_PRM_StressSeedLoadLevel,Region hierarchy level at which to include peak coincident load days as seeded stress periods (or False to ignore peaks),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,transgrp,
 GSw_PRM_StressSeedMinRElevel,Region hierarchy level at which to include minimum wind and solar capacity factor days as seeded stress periods (or False to ignore min-wind/solar CF days),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,interconnect,
 GSw_PRM_StressStorageCutoff,How to select shoulder stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is 'EUE' or 'capacity' or 'absolute'. If 'EUE' the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if 'cap' it specifies a headspace / storage capacity threshold [fraction]; if 'abs' it specifies absolute number of periods before/after [integer]. Turned off if set to 'off'.,N/A,EUE_0.1,
-GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stress periods; supported metrics (case-insensitive) are: LOLD | LOLE | LOLH | NEUE | Depth | Duration,N/A,NEUE,
+GSw_PRM_StressThresholdMetrics,/-delimited list of resource adequacy metrics. LOLD | LOLE | LOLH | NEUE | Depth | Duration identify stress periods; CVAR and NCVAR are report-only metrics and do not add stress periods or update PRM.,N/A,NEUE,
 GSw_PRM_StressThresholdDepth,Outage depth threshold [MW_EUE/MW_peak_load] (fraction); formulated as HierarchyLevel_Depth where HierarchyLevel is a column in hierarchy.csv; Depth is the max outage magnitude in [MW EUE] / [MW peak demand],N/A,transgrp_0.1,
 GSw_PRM_StressThresholdDuration,Outage duration threshold [hours]; formulated as HierarchyLevel_Duration where HierarchyLevel is a column in hierarchy.csv; Duration is the max outage duration in hours,N/A,transgrp_12,
 GSw_PRM_StressThresholdLOLD,LOLD threshold [event-days/year]; formulated as HierarchyLevel_LOLD where HierarchyLevel is a column in hierarchy.csv; LOLD is loss-of-load event-days per year,N/A,transgrp_0.1,
 GSw_PRM_StressThresholdLOLE,LOLE threshold [events/year]; formulated as HierarchyLevel_LOLE where HierarchyLevel is a column in hierarchy.csv; events is loss-of-load events per year,N/A,transgrp_0.1,
 GSw_PRM_StressThresholdLOLH,LOLH threshold [event-hours/year]; formulated as HierarchyLevel_LOLH where HierarchyLevel is a column in hierarchy.csv; LOLH is loss-of-load event-hours per year,N/A,transgrp_2.4,
 GSw_PRM_StressThresholdNEUE,NEUE threshold [ppm]; formulated as HierarchyLevel_NEUE where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million,N/A,transgrp_1,
+GSw_PRM_CVARAlpha,Alpha value used for CVAR/NCVAR calculations,float,0.95,
 GSw_PRM_UpdateFraction,Fraction to add to the PRM if a region fails RA threshold (only used if GSw_PRM_UpdateMethod = 1),float,0.02,
 GSw_PRM_UpdateMethod,Option to update PRM: (0) no update; (1) static update set by GSw_PRM_UpdateFraction; (2) dynamic update informed by PRAS; (3) dynamic update but only after all new stress periods have been added,0; 1; 2; 3,0,
 GSw_PRMTRADE_level,hierarchy level within which to allow PRM trading,r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region,country,
diff --git a/cases_test.csv b/cases_test.csv
index f3af6b246..a3c8fa38a 100644
--- a/cases_test.csv
+++ b/cases_test.csv
@@ -3,7 +3,7 @@ ignore,1,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
 GSw_Region,cendiv/Pacific,,country/USA,country/USA,country/USA,,st/ID.WY.NE.IA.IL,st/MA,,,,,,,st/WY,interconnect/western,transreg/PJM,st/NY.VT,st/OR,st/NE.NY.PA,st/NE.NY.PA,st/ID.WY.NE.IA.IL,st/KS,country/USA,country/USA,st/NY.NJ,,,,st/MA.RI.CT.NY.NJ.PA.OH,
 endyear,2032,,2050,2050,2050,2029,2060,2026,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,2050,,,,,
 yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,2010..2050..5,2010..2050..5,2010..2060..10,,,2010_2025_2050,2010..2050..5,,,,,
-GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,,,,,,
+GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,z132,,,,,
 GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,,,
 GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,,,
 GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,,
@@ -51,7 +51,7 @@ GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,,,
 GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,
 GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,
 GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,,,
-GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,,
+GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth/CVAR/NCVAR,,,,,
 GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,,
 GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,,
 GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,,
diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py
index cd325d75d..9bc2edff3 100644
--- a/reeds/resource_adequacy/ra_calcs.py
+++ b/reeds/resource_adequacy/ra_calcs.py
@@ -24,6 +24,7 @@ def run_pras(
         write_surplus=False,
         write_energy=False,
         write_shortfall_samples=False,
+        write_shortfall_samples_totals=False,
         write_availability_samples=False,
         **kwargs,
     ):
@@ -68,7 +69,9 @@ def run_pras(
         f"--write_surplus={int(write_surplus)}",
         f"--write_energy={int(write_energy)}",
         f"--write_shortfall_samples={int(write_shortfall_samples)}",
+        f"--write_shortfall_samples_totals={int(write_shortfall_samples_totals)}",
         f"--write_availability_samples={int(write_availability_samples)}",
+        f"--cvar_alpha={float(sw['GSw_PRM_CVARAlpha'])}",
         f"--iteration={iteration}",
         f"--samples={sw['pras_samples']}",
         f"--overwrite={int(overwrite)}",
@@ -153,13 +156,25 @@ def main(t, tnext, casedir, iteration=0):
         2: True,
     }[int(sw['pras'])]
     if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax):
+        stress_metrics = [
+            m.strip().upper()
+            for m in sw.GSw_PRM_StressThresholdMetrics.split('/')
+            if m.strip()
+        ]
+        write_shortfall_samples_totals = any(
+            metric in {'CVAR', 'NCVAR'}
+            for metric in stress_metrics
+        )
+
         result = run_pras(
             casedir, t,
             iteration=iteration,
             write_flow=(True if t == max(solveyears) else False),
             write_energy=True,
             write_shortfall_samples=(True if int(sw.GSw_PRM_UpdateMethod) > 1 else False),
+            write_shortfall_samples_totals=write_shortfall_samples_totals,
         )
+
         if result.returncode:
             raise Exception(
                 f"run_pras.jl returned code {result.returncode}. Check gamslog.txt for error trace."
diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl
index ec516737b..d39c7964d 100644
--- a/reeds/resource_adequacy/run_pras.jl
+++ b/reeds/resource_adequacy/run_pras.jl
@@ -70,10 +70,20 @@ function parse_commandline()
             default = 0
             required = false
         "--write_shortfall_samples"
-            help = "Write the sample-level shortfall"
+            help = "Write per-sample hourly shortfall by region"
             arg_type = Int
             default = 0
             required = false
+        "--write_shortfall_samples_totals"
+            help = "Write per-sample total shortfall by region over the full PRAS time period"
+            arg_type = Int
+            default = 0
+            required = false
+        "--cvar_alpha"
+            help = "Alpha for CVaR (e.g., 0.95)"
+            arg_type = Float64
+            default = 0.95
+            required = false
         "--write_availability_samples"
             help = "Write the sample-level generator and storage availability"
             arg_type = Int
@@ -182,7 +192,7 @@ function run_pras(pras_system_path::String, args::Dict)
     if args["write_energy"] == 1
         resultspec["energy"] = PRAS.StorageEnergy()
     end
-    if args["write_shortfall_samples"] == 1
+    if args["write_shortfall_samples"] == 1 || args["write_shortfall_samples_totals"] == 1
         resultspec["short_samples"] = PRAS.ShortfallSamples()
     end
     if args["write_availability_samples"] == 1
@@ -208,6 +218,38 @@ function run_pras(pras_system_path::String, args::Dict)
     @info "$(PRAS.EUE(results["short"]))"
     @info "$(PRAS.NEUE(results["short"]))"
 
+    #%% Print CVAR and NCVAR for the entire modeled region
+    if args["write_shortfall_samples_totals"] == 1 && haskey(results, "short_samples")
+        _, _, _, _, energyunit = PRAS.get_params(sys)
+        alpha = Float64(args["cvar_alpha"])
+
+        cvar_result = PRAS.CVAR(
+            energyunit,
+            results["short_samples"],
+            alpha,
+        )
+
+        cvar_value = PRAS.val(cvar_result.cvar)
+        cvar_stderr = PRAS.stderror(cvar_result.cvar)
+        cvar_var = cvar_result.var
+
+        ### Normalize CVaR by total load over the full PRAS time period.
+        total_load = sum(sys.regions.load)
+
+        ncvar_value = cvar_value / total_load * 1e6
+        ncvar_stderr = cvar_stderr / total_load * 1e6
+        ncvar_var = cvar_var / total_load * 1e6
+
+        @info(
+            "CVAR = $(cvar_value)±$(cvar_stderr) MWh; " *
+            "VaR = $(cvar_var) MWh; alpha = $(alpha)"
+        )
+        @info(
+            "NCVAR = $(ncvar_value)±$(ncvar_stderr) ppm; " *
+            "VaR = $(ncvar_var) ppm; alpha = $(alpha)"
+        )
+    end
+
     ## Filter out DC regions used for VSC HVDC transmission
     regions = [r for r in sys.regions.names if !(occursin("|", r))]
 
@@ -285,6 +327,7 @@ function run_pras(pras_system_path::String, args::Dict)
         end
         @info("Wrote PRAS surplus to $(surplusfile)")
     end
+
     ### Storage energy
     if args["write_energy"] == 1
         dfenergy = DF.DataFrame()
@@ -302,31 +345,41 @@ function run_pras(pras_system_path::String, args::Dict)
         @info("Wrote PRAS storage energy to $(energyfile)")
     end
 
-    ### Sample-level shortfall
+    ### Per-sample hourly shortfall by region
     if args["write_shortfall_samples"] == 1
-        dictshort = Dict(s => DF.DataFrame() for s = 1:args["samples"])
-        for s in range(1, args["samples"])
-            dictshort[s] = DF.DataFrame(
-                transpose(getindex.(results["short_samples"][:, :], s)),
-                sys.regions.names
-            )
-            # subset to regions (filter out DC regions) 
-            dictshort[s] = dictshort[s][:,findall(regions .∈ Ref(sys.regions.names))]
-        end
-        ## Write it
+        sf = results["short_samples"]
+        ## Use filtered system regions to avoid DC converter pseudo-regions without load.
+        region_names = sf.regions.names
+        idx = [findfirst(==(r), region_names) for r in regions]
+
         shortfile = replace(outfile, ".h5"=>"-shortfall_samples.h5")
         HDF5.h5open(shortfile, "w") do f
-            ## Create a group for each sample. Within each group, write an array for each region.
-            for s in range(1, args["samples"])
+            for s in 1:args["samples"]
                 HDF5.create_group(f, "$s")
-                for column in DF._names(dictshort[s])
-                    f["$s"]["$column", compress=4] = convert(Array, dictshort[s][!, column])
+                for (r, i_r) in zip(regions, idx)
+                    arr = Float64.(sf.shortfall[i_r, :, s])
+                    f["$s"]["$r", compress=4] = arr
                 end
             end
         end
         @info("Wrote PRAS shortfall by sample to $(shortfile)")
     end
 
+    ### Per-sample total shortfall by region over the full PRAS time period, needed for CVaR
+    if args["write_shortfall_samples_totals"] == 1
+        sf = results["short_samples"]
+        ## Use filtered system regions to avoid DC converter pseudo-regions without load.
+        totalsfile = replace(outfile, ".h5"=>"-shortfall_totals_by_sample.h5")
+        HDF5.h5open(totalsfile, "w") do f
+            f["sample", compress=4] = collect(1:args["samples"])
+            f["USA", compress=4] = Float64.(sf[])
+            for r in regions
+                f["$r", compress=4] = Float64.(sf[r])
+            end
+        end
+        @info("Wrote PRAS shortfall totals by sample to $(totalsfile)")
+    end
+
     ### Sample-level generator and storage availability
     if args["write_availability_samples"] == 1
         dictavail = Dict(s => DF.DataFrame() for s = 1:args["samples"])
@@ -457,6 +510,7 @@ if abspath(PROGRAM_FILE) == @__FILE__
     #     "write_surplus" => 0,
     #     "write_energy" => 0,
     #     "write_shortfall_samples" => 1,
+    #     "write_shortfall_samples_totals" => 1,
     #     "write_availability_samples" => 0,
     #     "overwrite" => 1,
     #     "debug" => 0,
@@ -466,6 +520,7 @@ if abspath(PROGRAM_FILE) == @__FILE__
     #     "pras_existing_unit_size" => 1,
     #     "pras_max_unitsize_prm" => 1,
     #     "pras_seed" => 1,
+    #     "cvar_alpha" => 0.95,
     # )
     # reedscase = args["reedscase"]
     # solve_year = args["solve_year"]
@@ -485,4 +540,4 @@ if abspath(PROGRAM_FILE) == @__FILE__
     main(args)
 
     #%%
-end
+end
\ No newline at end of file
diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index 7f41e0fd0..20150ae7e 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -14,6 +14,8 @@
 # import importlib
 # importlib.reload(functions)
 
+CVAR_METRICS = {'CVAR', 'NCVAR'}
+
 
 #%%### Constants
 RA_SWITCHES = {
@@ -142,6 +144,51 @@ def calc_neue(dfeue_agg, dfload_agg):
     neue = dfeue_agg.sum() / dfload_agg.sum() * 1e6
     return neue
 
+def get_cvar_alpha(sw):
+    alpha = float(sw.GSw_PRM_CVARAlpha)
+    if not (0 <= alpha < 1):
+        raise ValueError(f"GSw_PRM_CVARAlpha must be in [0, 1). Got {alpha}")
+    return alpha
+
+
+def get_shortfall_totals_by_sample(case, t, iteration=0):
+    filepath = os.path.join(case, 'handoff', 'PRAS', f'PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5')
+    if not os.path.isfile(filepath):
+        raise FileNotFoundError(f"{filepath} not found. Re-run PRAS with --write_shortfall_samples_totals 1.")
+    df = reeds.io.read_pras_results(filepath)
+    df.columns = df.columns.astype(str)
+    if 'sample' in df.columns:
+        df = df.set_index('sample')
+    elif df.index.name != 'sample':
+        df.index = pd.RangeIndex(1, len(df) + 1, name='sample')
+    return df.apply(pd.to_numeric, errors='coerce').clip(lower=0)
+
+
+def _sample_cvar(samples, alpha=0.95):
+    x = pd.Series(samples).dropna().astype(float)
+    if x.empty:
+        return np.nan
+    n_tail = max(1, int(np.ceil(round((1 - alpha) * len(x), 12))),)
+    return x.sort_values(ascending=False).iloc[:n_tail].mean()
+
+def calc_cvar(shortfall_samples_agg, alpha=0.95):
+    """
+    CVAR from total shortfall by PRAS sample.
+    """
+    cvar = shortfall_samples_agg.apply(
+        lambda s: _sample_cvar(s, alpha=alpha),
+        axis=0,
+    )
+    return cvar
+
+
+def calc_ncvar(cvar, dfload_agg):
+    """
+    NCVAR = CVAR / total load, in ppm.
+    """
+    ncvar = cvar / dfload_agg.sum().reindex(cvar.index) * 1e6
+    return ncvar
+
 
 def calc_peak_eue(dfeue_agg, dfload_agg, norm:Literal['peak','hourly','absolute']='peak'):
     """
@@ -192,6 +239,19 @@ def calc_ra_metrics(
     sw = reeds.io.get_switches(case)
     numyears = len(sw.resource_adequacy_years_list)
 
+    cvar_metrics = {
+        i.strip().upper()
+        for i in sw.GSw_PRM_StressThresholdMetrics.split('/')
+        if i.strip().upper() in CVAR_METRICS
+    }
+
+    if len(cvar_metrics):
+        shortfall_samples = (
+            get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration)
+            .drop(columns=['USA'], errors='ignore')
+        )
+        cvar_alpha = get_cvar_alpha(sw)
+
     ### Loop over aggregation levels and calculate all metrics
     ra_metrics = {}
     for level in levels:
@@ -200,9 +260,9 @@ def calc_ra_metrics(
         rmap = reeds.io.get_rmap(case=case, hierarchy_level=level)
         ## If multiple zones in one level and hour have LOLE, count that as one event,
         ## so take the max LOLE across the zones
-        dflole_agg = dflole.rename(columns=rmap).T.groupby(level=0).max().T
-        dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T
-        dfload_agg = dfload.rename(columns=rmap).T.groupby(level=0).sum().T
+        dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max()
+        dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum()
+        dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum()
         ## Calculate the full-timeseries metrics for each region
         ra_metrics[level, 'lold_peryear'] = calc_lold(dflole_agg) / numyears
         ra_metrics[level, 'lole_peryear'] = calc_lole(dflole_agg) / numyears
@@ -212,6 +272,22 @@ def calc_ra_metrics(
         ra_metrics[level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak')
         ra_metrics[level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'hourly')
         ra_metrics[level, 'euemax_mw'] = calc_peak_eue(dfeue_agg, dfload_agg, 'absolute')
+        if len(cvar_metrics):
+            regions = [c for c in shortfall_samples.columns if c in rmap.index]
+            if len(regions):
+                shortfall_samples_agg = (
+                    shortfall_samples[regions]
+                    .rename(columns=rmap)
+                    .groupby(axis=1, level=0)
+                    .sum()
+                )
+                cvar = calc_cvar(shortfall_samples_agg, alpha=cvar_alpha)
+
+                if 'CVAR' in cvar_metrics:
+                    ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
+
+                if 'NCVAR' in cvar_metrics:
+                    ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg)
 
     ### Combine it
     dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value')
@@ -232,7 +308,7 @@ def get_eue_events(
     events = {}
     for level in levels:
         rmap = reeds.io.get_rmap(case=case, hierarchy_level=level)
-        dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T
+        dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum()
         events[level] = pd.concat({r: get_events(dfeue_agg[r]) for r in dfeue_agg})
     dfout = pd.concat(events, names=['level','region','number'])
     return dfout
@@ -268,7 +344,7 @@ def get_longest_events(
     dates = []
     for i, row in eue_events.iterrows():
         dates.append(
-            pd.Series(index=pd.date_range(row.start, row.end, freq='h'), data=1)
+            pd.Series(index=pd.date_range(row.start, row.end, freq='H'), data=1)
             .resample('D').count()
         )
     if len(dates):
@@ -486,7 +562,7 @@ def get_stress_periods(case, sw, t, iteration):
     dfenergy = (
         dfenergy_unit
         .rename(columns={c: c.split('|')[1] for c in dfenergy_unit.columns})
-        .T.groupby(level=0).sum().T
+        .groupby(axis=1, level=0).sum()
     )
 
     ### Load this year's stress periods so we don't duplicate
@@ -501,23 +577,26 @@ def get_stress_periods(case, sw, t, iteration):
         stressperiods_this_iteration['start']
         + (
             (pd.Timedelta('5D') if sw.GSw_HourlyType == 'wek' else pd.Timedelta('1D'))
-             - pd.Timedelta('1h')
+             - pd.Timedelta('1H')
         )
     )
     ## Get already-modeled stress hours so we can exclude them from the hourly
     ## EUE and LOLE profiles used to determine new stress periods
     covered_hours = [
-        pd.date_range(row.start, row.end, freq='1h')
+        pd.date_range(row.start, row.end, freq='1H')
         for i,row in stressperiods_this_iteration.iterrows()
     ]
     covered_hours = [i for sublist in covered_hours for i in sublist]
-
-    ### Check all stress criteria; for regions that fail, add new stress periods
     _failed = {}
     _high_stress_periods = {}
     _shoulder_periods = {}
 
-    stress_metrics = [i.lower() for i in sw.GSw_PRM_StressThresholdMetrics.split('/')]
+    stress_metrics = []
+    for metric in sw.GSw_PRM_StressThresholdMetrics.split('/'):
+        metric = str(metric).strip()
+        if metric and metric.upper() not in CVAR_METRICS:
+            stress_metrics.append(metric.lower())
+
     for stress_metric in stress_metrics:
         switch = RA_SWITCHES[stress_metric]
         for criterion in sw[switch].split('/'):
@@ -525,9 +604,9 @@ def get_stress_periods(case, sw, t, iteration):
             ## Example: criterion = 'transgrp_1'
             hierarchy_level, metric_threshold = criterion.split('_')
             rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level)
-            dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T.drop(covered_hours)
-            dflole_agg = dflole.rename(columns=rmap).T.groupby(level=0).max().T.drop(covered_hours)
-            dfenergy_agg = dfenergy.rename(columns=rmap).T.groupby(level=0).sum().T.drop(covered_hours)
+            dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum().drop(covered_hours)
+            dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max().drop(covered_hours)
+            dfenergy_agg = dfenergy.rename(columns=rmap).groupby(axis=1, level=0).sum().drop(covered_hours)
             ## Get the stress periods
             dictout = check_threshold_and_choose_periods(
                 stress_metric,
@@ -865,6 +944,9 @@ def main(sw, t, iteration=0, logging=True):
         os.path.join(sw.casedir, 'inputs_case', newstresspath, 'prm.csv'),
     )
 
+    #%% Done
+    return
+
 
 # #%%### Option to run script directly for debugging
 # if __name__ == '__main__':
diff --git a/runreeds.py b/runreeds.py
index 2688a713c..2c5b175bd 100644
--- a/runreeds.py
+++ b/runreeds.py
@@ -300,21 +300,28 @@ def check_compatibility(sw):
         i.lower(): f'GSw_PRM_StressThreshold{i}'
         for i in ['Depth', 'Duration', 'LOLD', 'LOLE', 'LOLH', 'NEUE']
     }
-    used_metrics = [i.lower() for i in sw['GSw_PRM_StressThresholdMetrics'].split('/')]
+    report_only_metrics = {'cvar', 'ncvar'}
+
+    used_metrics = [i.strip().lower() for i in sw['GSw_PRM_StressThresholdMetrics'].split('/') if i.strip()]
     allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r']
 
     for metric in used_metrics:
+        if metric in report_only_metrics:
+            continue
+
         if metric not in ra_switches:
             raise NotImplementedError(f"GSw_PRM_StressThresholdMetrics = {metric} is not supported")
 
         for threshold in sw[ra_switches[metric]].split('/'):
             ## Example: GSw_PRM_StressThresholdNEUE = 'transgrp_1'
-            (hierarchy_level, stress_value) = threshold.split('_')
+            hierarchy_level, stress_value = threshold.split('_')
+
             if hierarchy_level not in allowed_levels:
                 raise ValueError(
                     f"{ra_switches[metric]}: level={hierarchy_level} but must be in:\n"
                     + '\n'.join(allowed_levels)
                 )
+
             if not (float(stress_value) >= 0):
                 raise ValueError(
                     f"stress value in {ra_switches[metric]} must be a positive number "
@@ -835,19 +842,22 @@ def setupEnvironment(
 
     #%% Check whether the ReEDS conda environment is activated
     if (not skip_checks) and (
-        ('reeds' not in os.environ['CONDA_DEFAULT_ENV'].lower())
-        or (not pd.__version__.startswith('3'))
+        ('reeds2' not in os.environ['CONDA_DEFAULT_ENV'].lower())
+        or (not pd.__version__.startswith('2'))
     ):
-        err = (
+        print(
             f"Your environment is {os.environ['CONDA_DEFAULT_ENV']} and your pandas "
-            f"version is {pd.__version__}.\nThe supported environment is 'reeds', with\n"
-            "pandas version 3.x.\n"
+            f"version is {pd.__version__}.\nThe default environment is 'reeds2', with\n"
+            "pandas version 2.x, so the python parts of ReEDS are unlikely to work.\n"
             "To build the environment for the first time, run:\n"
             "    `conda env create -f environment.yml`\n"
             "To activate the created environment, run:\n"
-            "    `conda activate reeds` (or `activate reeds` on Windows)"
+            "    `conda activate reeds2` (or `activate reeds2` on Windows)\n"
+            "Do you want to continue without activating the environment?"
         )
-        raise ValueError(err)
+        confirm_env = str(input("Continue? y/[n]: ") or 'n')
+        if confirm_env not in ['y','Y','yes','Yes','YES']:
+            quit()
 
     #%% Load specified case file, infer other settings from cases.csv
     if cases_suffix in ['', 'default']:
@@ -1275,7 +1285,7 @@ def write_batch_script(
                 OPATH.writelines("module load conda \n")
                 OPATH.writelines("module load gams \n")
 
-            OPATH.writelines("conda activate reeds \n")
+            OPATH.writelines("conda activate reeds2 \n")
             OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n')
 
         #%% Write the input_processing script calls
@@ -1283,7 +1293,6 @@ def write_batch_script(
         for s in [
             'copy_files',
             'mcs_sampler',
-            'climateprep',            
             'hydcf',
             'h2_storage',
             'calc_financial_inputs',
@@ -1292,6 +1301,7 @@ def write_batch_script(
             'writesupplycurves',
             'writedrshift',
             'plantcostprep',
+            'climateprep',
             'hourly_load',
             'recf',
             'forecast',

From 3b0571081a092338976947ca2a1f8b6c75b29861 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 14:45:55 -0400
Subject: [PATCH 02/11] revert runreeds.py

---
 runreeds.py | 48 ++++++++++++++++++++++--------------------------
 1 file changed, 22 insertions(+), 26 deletions(-)

diff --git a/runreeds.py b/runreeds.py
index 2c5b175bd..7e1db8ba8 100644
--- a/runreeds.py
+++ b/runreeds.py
@@ -216,8 +216,11 @@ def check_cases_format(df_cases):
 
 
 def check_compatibility(sw):
-    if int(sw['startyear']) < 2010:
-        raise ValueError(f"startyear = {sw['startyear']} but must be ≥ 2010")
+    if int(sw['startyear']) != 2010:
+        raise ValueError(f"startyear = {sw['startyear']} but must be = 2010")
+
+    if int(sw['GSw_SkipRAyear']) <= int(sw['startyear']):
+            raise ValueError(f"GSw_SkipRAyear = {sw['GSw_SkipRAyear']} but must be > {sw['startyear']}")
 
     if (sw['GSw_HourlyType'] in ['year']) and int(sw['GSw_InterDayLinkage']):
         raise ValueError(
@@ -300,28 +303,21 @@ def check_compatibility(sw):
         i.lower(): f'GSw_PRM_StressThreshold{i}'
         for i in ['Depth', 'Duration', 'LOLD', 'LOLE', 'LOLH', 'NEUE']
     }
-    report_only_metrics = {'cvar', 'ncvar'}
-
-    used_metrics = [i.strip().lower() for i in sw['GSw_PRM_StressThresholdMetrics'].split('/') if i.strip()]
+    used_metrics = [i.lower() for i in sw['GSw_PRM_StressThresholdMetrics'].split('/')]
     allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r']
 
     for metric in used_metrics:
-        if metric in report_only_metrics:
-            continue
-
         if metric not in ra_switches:
             raise NotImplementedError(f"GSw_PRM_StressThresholdMetrics = {metric} is not supported")
 
         for threshold in sw[ra_switches[metric]].split('/'):
             ## Example: GSw_PRM_StressThresholdNEUE = 'transgrp_1'
-            hierarchy_level, stress_value = threshold.split('_')
-
+            (hierarchy_level, stress_value) = threshold.split('_')
             if hierarchy_level not in allowed_levels:
                 raise ValueError(
                     f"{ra_switches[metric]}: level={hierarchy_level} but must be in:\n"
                     + '\n'.join(allowed_levels)
                 )
-
             if not (float(stress_value) >= 0):
                 raise ValueError(
                     f"stress value in {ra_switches[metric]} must be a positive number "
@@ -842,22 +838,19 @@ def setupEnvironment(
 
     #%% Check whether the ReEDS conda environment is activated
     if (not skip_checks) and (
-        ('reeds2' not in os.environ['CONDA_DEFAULT_ENV'].lower())
-        or (not pd.__version__.startswith('2'))
+        ('reeds' not in os.environ['CONDA_DEFAULT_ENV'].lower())
+        or (not pd.__version__.startswith('3'))
     ):
-        print(
+        err = (
             f"Your environment is {os.environ['CONDA_DEFAULT_ENV']} and your pandas "
-            f"version is {pd.__version__}.\nThe default environment is 'reeds2', with\n"
-            "pandas version 2.x, so the python parts of ReEDS are unlikely to work.\n"
+            f"version is {pd.__version__}.\nThe supported environment is 'reeds', with\n"
+            "pandas version 3.x.\n"
             "To build the environment for the first time, run:\n"
             "    `conda env create -f environment.yml`\n"
             "To activate the created environment, run:\n"
-            "    `conda activate reeds2` (or `activate reeds2` on Windows)\n"
-            "Do you want to continue without activating the environment?"
+            "    `conda activate reeds` (or `activate reeds` on Windows)"
         )
-        confirm_env = str(input("Continue? y/[n]: ") or 'n')
-        if confirm_env not in ['y','Y','yes','Yes','YES']:
-            quit()
+        raise ValueError(err)
 
     #%% Load specified case file, infer other settings from cases.csv
     if cases_suffix in ['', 'default']:
@@ -904,7 +897,7 @@ def setupEnvironment(
                     f'Specified single={single} but available cases are:\n'
                     + '\n> '.join([c for c in df_cases.columns])
                 )
-                raise KeyError(err)
+                raise ValueError(err)
         df_cases = df_cases[single.split(',')].copy()
         casenames = single.split(',')
 
@@ -1285,14 +1278,16 @@ def write_batch_script(
                 OPATH.writelines("module load conda \n")
                 OPATH.writelines("module load gams \n")
 
-            OPATH.writelines("conda activate reeds2 \n")
+            OPATH.writelines("conda activate reeds \n")
             OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n')
 
         #%% Write the input_processing script calls
         big_comment('Input processing', OPATH)
         for s in [
             'copy_files',
+            'process_unitdata',
             'mcs_sampler',
+            'climateprep',            
             'hydcf',
             'h2_storage',
             'calc_financial_inputs',
@@ -1301,13 +1296,12 @@ def write_batch_script(
             'writesupplycurves',
             'writedrshift',
             'plantcostprep',
-            'climateprep',
             'hourly_load',
             'recf',
-            'forecast',
             'WriteHintage',
             'transmission',
             'outage_rates',
+            'forecast',
             'hourly_repperiods',
             'h5_to_gdx',
         ]:
@@ -1399,7 +1393,9 @@ def write_batch_script(
         if not LINUXORMAC:
             OPATH.writelines("endlocal\n")
         OPATH.writelines(f'python {logger}\n')
-        OPATH.writelines(f"python {Path('reeds','core','terminus','report_dump.py')} {casedir} -c\n\n")
+        OPATH.writelines(f"python {Path('reeds','core','terminus','report_dump.py')} {casedir} -c\n")
+        OPATH.writelines(writescripterrorcheck('report_dump.py')+'\n')
+
         if int(caseSwitches['diagnose']):
              OPATH.writelines(
                 "python"

From 7e560cd67d82185cc8516f0260e6b47c160270fd Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 14:47:42 -0400
Subject: [PATCH 03/11] revert cases_test.csv

---
 cases_test.csv | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cases_test.csv b/cases_test.csv
index f8a3e0d07..f3a122adc 100644
--- a/cases_test.csv
+++ b/cases_test.csv
@@ -51,7 +51,7 @@ GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,,,
 GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,
 GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,
 GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,,,
-GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth/CVAR/NCVAR,,,,,
+GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,,
 GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,,
 GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,,
 GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,,

From 3a7626cd213854560353202ac6ea96ed0e136149 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 14:50:41 -0400
Subject: [PATCH 04/11] remove CVAR from GSw_PRM_StressThresholdMetrics

---
 cases.csv | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cases.csv b/cases.csv
index e97c94bff..9b1138cb0 100644
--- a/cases.csv
+++ b/cases.csv
@@ -272,7 +272,7 @@ GSw_PRM_StressOutages,Whether to apply the availability factor (forced + schedul
 GSw_PRM_StressSeedLoadLevel,Region hierarchy level at which to include peak coincident load days as seeded stress periods (or False to ignore peaks),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,transgrp,
 GSw_PRM_StressSeedMinRElevel,Region hierarchy level at which to include minimum wind and solar capacity factor days as seeded stress periods (or False to ignore min-wind/solar CF days),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,interconnect,
 GSw_PRM_StressStorageCutoff,How to select shoulder stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is 'EUE' or 'capacity' or 'absolute'. If 'EUE' the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if 'cap' it specifies a headspace / storage capacity threshold [fraction]; if 'abs' it specifies absolute number of periods before/after [integer]. Turned off if set to 'off'.,N/A,EUE_0.1,
-GSw_PRM_StressThresholdMetrics,/-delimited list of resource adequacy metrics. LOLD | LOLE | LOLH | NEUE | Depth | Duration identify stress periods; CVAR and NCVAR are report-only metrics and do not add stress periods or update PRM.,N/A,NEUE,
+GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stress periods; supported metrics (case-insensitive) are: LOLD | LOLE | LOLH | NEUE | Depth | Duration,N/A,NEUE,
 GSw_PRM_StressThresholdDepth,Outage depth threshold [MW_EUE/MW_peak_load] (fraction); formulated as HierarchyLevel_Depth where HierarchyLevel is a column in hierarchy.csv; Depth is the max outage magnitude in [MW EUE] / [MW peak demand],N/A,transgrp_0.1,
 GSw_PRM_StressThresholdDuration,Outage duration threshold [hours]; formulated as HierarchyLevel_Duration where HierarchyLevel is a column in hierarchy.csv; Duration is the max outage duration in hours,N/A,transgrp_12,
 GSw_PRM_StressThresholdLOLD,LOLD threshold [event-days/year]; formulated as HierarchyLevel_LOLD where HierarchyLevel is a column in hierarchy.csv; LOLD is loss-of-load event-days per year,N/A,transgrp_0.1,

From 48f41eba3308978c7573570f59e1aaaf94ae91e4 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 15:16:39 -0400
Subject: [PATCH 05/11] Always calculate CVAR metric

---
 reeds/resource_adequacy/ra_calcs.py       | 15 +------
 reeds/resource_adequacy/stress_periods.py | 54 ++++++++---------------
 2 files changed, 21 insertions(+), 48 deletions(-)

diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py
index 9bc2edff3..ef120aac1 100644
--- a/reeds/resource_adequacy/ra_calcs.py
+++ b/reeds/resource_adequacy/ra_calcs.py
@@ -24,7 +24,7 @@ def run_pras(
         write_surplus=False,
         write_energy=False,
         write_shortfall_samples=False,
-        write_shortfall_samples_totals=False,
+        write_shortfall_samples_totals=True,
         write_availability_samples=False,
         **kwargs,
     ):
@@ -155,24 +155,13 @@ def main(t, tnext, casedir, iteration=0):
         1: True if t == max(solveyears) else False,
         2: True,
     }[int(sw['pras'])]
-    if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax):
-        stress_metrics = [
-            m.strip().upper()
-            for m in sw.GSw_PRM_StressThresholdMetrics.split('/')
-            if m.strip()
-        ]
-        write_shortfall_samples_totals = any(
-            metric in {'CVAR', 'NCVAR'}
-            for metric in stress_metrics
-        )
-
+    if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax):        
         result = run_pras(
             casedir, t,
             iteration=iteration,
             write_flow=(True if t == max(solveyears) else False),
             write_energy=True,
             write_shortfall_samples=(True if int(sw.GSw_PRM_UpdateMethod) > 1 else False),
-            write_shortfall_samples_totals=write_shortfall_samples_totals,
         )
 
         if result.returncode:
diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index 20150ae7e..bcee6b330 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -14,8 +14,6 @@
 # import importlib
 # importlib.reload(functions)
 
-CVAR_METRICS = {'CVAR', 'NCVAR'}
-
 
 #%%### Constants
 RA_SWITCHES = {
@@ -239,18 +237,12 @@ def calc_ra_metrics(
     sw = reeds.io.get_switches(case)
     numyears = len(sw.resource_adequacy_years_list)
 
-    cvar_metrics = {
-        i.strip().upper()
-        for i in sw.GSw_PRM_StressThresholdMetrics.split('/')
-        if i.strip().upper() in CVAR_METRICS
-    }
-
-    if len(cvar_metrics):
-        shortfall_samples = (
-            get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration)
-            .drop(columns=['USA'], errors='ignore')
-        )
-        cvar_alpha = get_cvar_alpha(sw)
+    ### Get total shortfall samples for CVAR calculation
+    shortfall_samples = (
+        get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration)
+        .drop(columns=['USA'], errors='ignore')
+    )
+    cvar_alpha = get_cvar_alpha(sw)
 
     ### Loop over aggregation levels and calculate all metrics
     ra_metrics = {}
@@ -272,22 +264,18 @@ def calc_ra_metrics(
         ra_metrics[level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak')
         ra_metrics[level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'hourly')
         ra_metrics[level, 'euemax_mw'] = calc_peak_eue(dfeue_agg, dfload_agg, 'absolute')
-        if len(cvar_metrics):
-            regions = [c for c in shortfall_samples.columns if c in rmap.index]
-            if len(regions):
-                shortfall_samples_agg = (
-                    shortfall_samples[regions]
-                    .rename(columns=rmap)
-                    .groupby(axis=1, level=0)
-                    .sum()
-                )
-                cvar = calc_cvar(shortfall_samples_agg, alpha=cvar_alpha)
-
-                if 'CVAR' in cvar_metrics:
-                    ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
-
-                if 'NCVAR' in cvar_metrics:
-                    ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg)
+        ## Calculate tail-based metrics (CVAR and NCVAR)
+        regions = [c for c in shortfall_samples.columns if c in rmap.index]
+        if len(regions):
+            shortfall_samples_agg = (
+                shortfall_samples[regions]
+                .rename(columns=rmap)
+                .groupby(axis=1, level=0)
+                .sum()
+            )
+            cvar = calc_cvar(shortfall_samples_agg, alpha=cvar_alpha)
+            ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
+            ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg)
 
     ### Combine it
     dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value')
@@ -591,11 +579,7 @@ def get_stress_periods(case, sw, t, iteration):
     _high_stress_periods = {}
     _shoulder_periods = {}
 
-    stress_metrics = []
-    for metric in sw.GSw_PRM_StressThresholdMetrics.split('/'):
-        metric = str(metric).strip()
-        if metric and metric.upper() not in CVAR_METRICS:
-            stress_metrics.append(metric.lower())
+    stress_metrics = [i.lower() for i in sw.GSw_PRM_StressThresholdMetrics.split('/')]
 
     for stress_metric in stress_metrics:
         switch = RA_SWITCHES[stress_metric]

From c5995a942f98da44f632b49f1c92ba90f7077061 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 15:22:12 -0400
Subject: [PATCH 06/11] revert back to new groupby syntax

---
 reeds/resource_adequacy/stress_periods.py | 23 ++++++++++-------------
 1 file changed, 10 insertions(+), 13 deletions(-)

diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index bcee6b330..56e6906dc 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -252,9 +252,9 @@ def calc_ra_metrics(
         rmap = reeds.io.get_rmap(case=case, hierarchy_level=level)
         ## If multiple zones in one level and hour have LOLE, count that as one event,
         ## so take the max LOLE across the zones
-        dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max()
-        dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum()
-        dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum()
+        dflole_agg = dflole.rename(columns=rmap).T.groupby(level=0).max().T
+        dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T
+        dfload_agg = dfload.rename(columns=rmap).T.groupby(level=0).sum().T
         ## Calculate the full-timeseries metrics for each region
         ra_metrics[level, 'lold_peryear'] = calc_lold(dflole_agg) / numyears
         ra_metrics[level, 'lole_peryear'] = calc_lole(dflole_agg) / numyears
@@ -270,8 +270,8 @@ def calc_ra_metrics(
             shortfall_samples_agg = (
                 shortfall_samples[regions]
                 .rename(columns=rmap)
-                .groupby(axis=1, level=0)
-                .sum()
+                .T.groupby(level=0)
+                .sum().T
             )
             cvar = calc_cvar(shortfall_samples_agg, alpha=cvar_alpha)
             ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
@@ -296,7 +296,7 @@ def get_eue_events(
     events = {}
     for level in levels:
         rmap = reeds.io.get_rmap(case=case, hierarchy_level=level)
-        dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum()
+        dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T
         events[level] = pd.concat({r: get_events(dfeue_agg[r]) for r in dfeue_agg})
     dfout = pd.concat(events, names=['level','region','number'])
     return dfout
@@ -550,7 +550,7 @@ def get_stress_periods(case, sw, t, iteration):
     dfenergy = (
         dfenergy_unit
         .rename(columns={c: c.split('|')[1] for c in dfenergy_unit.columns})
-        .groupby(axis=1, level=0).sum()
+        .T.groupby(level=0).sum().T
     )
 
     ### Load this year's stress periods so we don't duplicate
@@ -588,9 +588,9 @@ def get_stress_periods(case, sw, t, iteration):
             ## Example: criterion = 'transgrp_1'
             hierarchy_level, metric_threshold = criterion.split('_')
             rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level)
-            dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum().drop(covered_hours)
-            dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max().drop(covered_hours)
-            dfenergy_agg = dfenergy.rename(columns=rmap).groupby(axis=1, level=0).sum().drop(covered_hours)
+            dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T.drop(covered_hours)
+            dflole_agg = dflole.rename(columns=rmap).T.groupby(level=0).max().T.drop(covered_hours)
+            dfenergy_agg = dfenergy.rename(columns=rmap).T.groupby(level=0).sum().T.drop(covered_hours)
             ## Get the stress periods
             dictout = check_threshold_and_choose_periods(
                 stress_metric,
@@ -928,9 +928,6 @@ def main(sw, t, iteration=0, logging=True):
         os.path.join(sw.casedir, 'inputs_case', newstresspath, 'prm.csv'),
     )
 
-    #%% Done
-    return
-
 
 # #%%### Option to run script directly for debugging
 # if __name__ == '__main__':

From 17089a7f80d250495ab253dac5ba0d9bf15719b8 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 15:26:48 -0400
Subject: [PATCH 07/11] move up CVAR switch compatability check

---
 reeds/resource_adequacy/stress_periods.py | 9 +--------
 runreeds.py                               | 5 +++++
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index 56e6906dc..4c48b1108 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -142,12 +142,6 @@ def calc_neue(dfeue_agg, dfload_agg):
     neue = dfeue_agg.sum() / dfload_agg.sum() * 1e6
     return neue
 
-def get_cvar_alpha(sw):
-    alpha = float(sw.GSw_PRM_CVARAlpha)
-    if not (0 <= alpha < 1):
-        raise ValueError(f"GSw_PRM_CVARAlpha must be in [0, 1). Got {alpha}")
-    return alpha
-
 
 def get_shortfall_totals_by_sample(case, t, iteration=0):
     filepath = os.path.join(case, 'handoff', 'PRAS', f'PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5')
@@ -242,7 +236,6 @@ def calc_ra_metrics(
         get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration)
         .drop(columns=['USA'], errors='ignore')
     )
-    cvar_alpha = get_cvar_alpha(sw)
 
     ### Loop over aggregation levels and calculate all metrics
     ra_metrics = {}
@@ -273,7 +266,7 @@ def calc_ra_metrics(
                 .T.groupby(level=0)
                 .sum().T
             )
-            cvar = calc_cvar(shortfall_samples_agg, alpha=cvar_alpha)
+            cvar = calc_cvar(shortfall_samples_agg, alpha=float(sw.GSw_PRM_CVARAlpha))
             ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
             ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg)
 
diff --git a/runreeds.py b/runreeds.py
index 7e1db8ba8..61e815d69 100644
--- a/runreeds.py
+++ b/runreeds.py
@@ -323,6 +323,11 @@ def check_compatibility(sw):
                     f"stress value in {ra_switches[metric]} must be a positive number "
                     f"but '{stress_value}' was provided"
                 )
+    
+    ## CVAR value in [0,1)
+    alpha = float(sw['GSw_PRM_CVARAlpha'])
+    if not (0 <= alpha < 1):
+        raise ValueError(f"GSw_PRM_CVARAlpha must be in [0, 1). Got {alpha}")
 
     ### GSw_PRM_UpdateMethod 1-3 (static or PRAS-informed PRM update) is computed from the
     ### NEUE-based shortfall, so it requires NEUE to be an active stress metric

From 3735a61cb1502934240d2c617823e25d47742080 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 15:28:38 -0400
Subject: [PATCH 08/11] change to GSw_PRM_CVARalpha

---
 reeds/resource_adequacy/ra_calcs.py       | 2 +-
 reeds/resource_adequacy/stress_periods.py | 2 +-
 runreeds.py                               | 4 ++--
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py
index ef120aac1..9c3b474e5 100644
--- a/reeds/resource_adequacy/ra_calcs.py
+++ b/reeds/resource_adequacy/ra_calcs.py
@@ -71,7 +71,7 @@ def run_pras(
         f"--write_shortfall_samples={int(write_shortfall_samples)}",
         f"--write_shortfall_samples_totals={int(write_shortfall_samples_totals)}",
         f"--write_availability_samples={int(write_availability_samples)}",
-        f"--cvar_alpha={float(sw['GSw_PRM_CVARAlpha'])}",
+        f"--cvar_alpha={float(sw['GSw_PRM_CVARalpha'])}",
         f"--iteration={iteration}",
         f"--samples={sw['pras_samples']}",
         f"--overwrite={int(overwrite)}",
diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index 4c48b1108..3751d0304 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -266,7 +266,7 @@ def calc_ra_metrics(
                 .T.groupby(level=0)
                 .sum().T
             )
-            cvar = calc_cvar(shortfall_samples_agg, alpha=float(sw.GSw_PRM_CVARAlpha))
+            cvar = calc_cvar(shortfall_samples_agg, alpha=float(sw.GSw_PRM_CVARalpha))
             ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
             ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg)
 
diff --git a/runreeds.py b/runreeds.py
index 61e815d69..69429300c 100644
--- a/runreeds.py
+++ b/runreeds.py
@@ -325,9 +325,9 @@ def check_compatibility(sw):
                 )
     
     ## CVAR value in [0,1)
-    alpha = float(sw['GSw_PRM_CVARAlpha'])
+    alpha = float(sw['GSw_PRM_CVARalpha'])
     if not (0 <= alpha < 1):
-        raise ValueError(f"GSw_PRM_CVARAlpha must be in [0, 1). Got {alpha}")
+        raise ValueError(f"GSw_PRM_CVARalpha must be in [0, 1). Got {alpha}")
 
     ### GSw_PRM_UpdateMethod 1-3 (static or PRAS-informed PRM update) is computed from the
     ### NEUE-based shortfall, so it requires NEUE to be an active stress metric

From 91730d278535b3103705df7aed1aa6c1ede6f994 Mon Sep 17 00:00:00 2001
From: bsergi <bsergi@nlr.gov>
Date: Mon, 10 Aug 2026 15:30:46 -0400
Subject: [PATCH 09/11] Missed a change to GSw_PRM_CVARalpha

---
 cases.csv | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cases.csv b/cases.csv
index 9b1138cb0..a9bc848c1 100644
--- a/cases.csv
+++ b/cases.csv
@@ -279,7 +279,7 @@ GSw_PRM_StressThresholdLOLD,LOLD threshold [event-days/year]; formulated as Hier
 GSw_PRM_StressThresholdLOLE,LOLE threshold [events/year]; formulated as HierarchyLevel_LOLE where HierarchyLevel is a column in hierarchy.csv; events is loss-of-load events per year,N/A,transgrp_0.1,
 GSw_PRM_StressThresholdLOLH,LOLH threshold [event-hours/year]; formulated as HierarchyLevel_LOLH where HierarchyLevel is a column in hierarchy.csv; LOLH is loss-of-load event-hours per year,N/A,transgrp_2.4,
 GSw_PRM_StressThresholdNEUE,NEUE threshold [ppm]; formulated as HierarchyLevel_NEUE where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million,N/A,transgrp_1,
-GSw_PRM_CVARAlpha,Alpha value used for CVAR/NCVAR calculations,float,0.95,
+GSw_PRM_CVARalpha,Alpha value used for CVAR/NCVAR calculations,float,0.95,
 GSw_PRM_UpdateFraction,Fraction to add to the PRM if a region fails RA threshold (only used if GSw_PRM_UpdateMethod = 1),float,0.02,
 GSw_PRM_UpdateMethod,Option to update PRM: (0) no update; (1) static update set by GSw_PRM_UpdateFraction; (2) dynamic update informed by PRAS; (3) dynamic update but only after all new stress periods have been added,0; 1; 2; 3,0,
 GSw_PRMTRADE_level,hierarchy level within which to allow PRM trading,r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region,country,

From 081086f487781b607ba3f9c25b7dad4eb16dab02 Mon Sep 17 00:00:00 2001
From: Burcin Cakir Erdener <bcakire@kl2.head.cm.kestrel.hpc.nlr.gov>
Date: Thu, 27 Aug 2026 13:39:57 -0600
Subject: [PATCH 10/11] Remove Julia-side CVAR dependency

---
 reeds/resource_adequacy/ra_calcs.py       |  3 +-
 reeds/resource_adequacy/run_pras.jl       | 38 -----------------------
 reeds/resource_adequacy/stress_periods.py | 10 ++++--
 3 files changed, 9 insertions(+), 42 deletions(-)

diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py
index 9c3b474e5..6ed95ff21 100644
--- a/reeds/resource_adequacy/ra_calcs.py
+++ b/reeds/resource_adequacy/ra_calcs.py
@@ -71,7 +71,6 @@ def run_pras(
         f"--write_shortfall_samples={int(write_shortfall_samples)}",
         f"--write_shortfall_samples_totals={int(write_shortfall_samples_totals)}",
         f"--write_availability_samples={int(write_availability_samples)}",
-        f"--cvar_alpha={float(sw['GSw_PRM_CVARalpha'])}",
         f"--iteration={iteration}",
         f"--samples={sw['pras_samples']}",
         f"--overwrite={int(overwrite)}",
@@ -155,7 +154,7 @@ def main(t, tnext, casedir, iteration=0):
         1: True if t == max(solveyears) else False,
         2: True,
     }[int(sw['pras'])]
-    if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax):        
+    if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax):
         result = run_pras(
             casedir, t,
             iteration=iteration,
diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl
index d39c7964d..132db16d1 100644
--- a/reeds/resource_adequacy/run_pras.jl
+++ b/reeds/resource_adequacy/run_pras.jl
@@ -79,11 +79,6 @@ function parse_commandline()
             arg_type = Int
             default = 0
             required = false
-        "--cvar_alpha"
-            help = "Alpha for CVaR (e.g., 0.95)"
-            arg_type = Float64
-            default = 0.95
-            required = false
         "--write_availability_samples"
             help = "Write the sample-level generator and storage availability"
             arg_type = Int
@@ -218,38 +213,6 @@ function run_pras(pras_system_path::String, args::Dict)
     @info "$(PRAS.EUE(results["short"]))"
     @info "$(PRAS.NEUE(results["short"]))"
 
-    #%% Print CVAR and NCVAR for the entire modeled region
-    if args["write_shortfall_samples_totals"] == 1 && haskey(results, "short_samples")
-        _, _, _, _, energyunit = PRAS.get_params(sys)
-        alpha = Float64(args["cvar_alpha"])
-
-        cvar_result = PRAS.CVAR(
-            energyunit,
-            results["short_samples"],
-            alpha,
-        )
-
-        cvar_value = PRAS.val(cvar_result.cvar)
-        cvar_stderr = PRAS.stderror(cvar_result.cvar)
-        cvar_var = cvar_result.var
-
-        ### Normalize CVaR by total load over the full PRAS time period.
-        total_load = sum(sys.regions.load)
-
-        ncvar_value = cvar_value / total_load * 1e6
-        ncvar_stderr = cvar_stderr / total_load * 1e6
-        ncvar_var = cvar_var / total_load * 1e6
-
-        @info(
-            "CVAR = $(cvar_value)±$(cvar_stderr) MWh; " *
-            "VaR = $(cvar_var) MWh; alpha = $(alpha)"
-        )
-        @info(
-            "NCVAR = $(ncvar_value)±$(ncvar_stderr) ppm; " *
-            "VaR = $(ncvar_var) ppm; alpha = $(alpha)"
-        )
-    end
-
     ## Filter out DC regions used for VSC HVDC transmission
     regions = [r for r in sys.regions.names if !(occursin("|", r))]
 
@@ -520,7 +483,6 @@ if abspath(PROGRAM_FILE) == @__FILE__
     #     "pras_existing_unit_size" => 1,
     #     "pras_max_unitsize_prm" => 1,
     #     "pras_seed" => 1,
-    #     "cvar_alpha" => 0.95,
     # )
     # reedscase = args["reedscase"]
     # solve_year = args["solve_year"]
diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index 3751d0304..7e2ae4a99 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -160,8 +160,14 @@ def _sample_cvar(samples, alpha=0.95):
     x = pd.Series(samples).dropna().astype(float)
     if x.empty:
         return np.nan
-    n_tail = max(1, int(np.ceil(round((1 - alpha) * len(x), 12))),)
-    return x.sort_values(ascending=False).iloc[:n_tail].mean()
+
+    # Round before applying ceil to remove floating-point noise.
+    # For example, a mathematically exact tail size of 50 may be
+    # represented as 50.00000000000004, which would otherwise select 51 samples.
+    tail_size = round((1 - alpha) * len(x), 12)
+    n_tail = max(1, int(np.ceil(tail_size)))
+
+    return x.nlargest(n_tail).mean()
 
 def calc_cvar(shortfall_samples_agg, alpha=0.95):
     """

From 6a95ef2e6bf2e9a467c64b49f1444846dee713ff Mon Sep 17 00:00:00 2001
From: Burcin Cakir Erdener <bcakire@kl2.head.cm.kestrel.hpc.nlr.gov>
Date: Fri, 28 Aug 2026 20:46:31 -0600
Subject: [PATCH 11/11] Restore pandas 3 hourly frequency syntax

---
 reeds/resource_adequacy/stress_periods.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py
index 7e2ae4a99..e6cecd1e1 100644
--- a/reeds/resource_adequacy/stress_periods.py
+++ b/reeds/resource_adequacy/stress_periods.py
@@ -564,13 +564,13 @@ def get_stress_periods(case, sw, t, iteration):
         stressperiods_this_iteration['start']
         + (
             (pd.Timedelta('5D') if sw.GSw_HourlyType == 'wek' else pd.Timedelta('1D'))
-             - pd.Timedelta('1H')
+             - pd.Timedelta('1h')
         )
     )
     ## Get already-modeled stress hours so we can exclude them from the hourly
     ## EUE and LOLE profiles used to determine new stress periods
     covered_hours = [
-        pd.date_range(row.start, row.end, freq='1H')
+        pd.date_range(row.start, row.end, freq='1h')
         for i,row in stressperiods_this_iteration.iterrows()
     ]
     covered_hours = [i for sublist in covered_hours for i in sublist]