From 729f2007bf83e25fb2e94562411836ed480b0653 Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Wed, 29 Oct 2025 19:35:15 -0400
Subject: [PATCH 01/11] set up obs and preds of mcmc as dictionary and match
 keys when calculating log likelihood of each step

---
 pygem/bin/run/run_calibration.py | 201 +++++++++++++++++++------------
 pygem/mcmc.py                    |  80 ++++++------
 pygem/plot/graphics.py           |  52 ++++----
 pygem/setup/config.py            |   1 +
 pygem/setup/config.yaml          |   1 +
 5 files changed, 189 insertions(+), 146 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 60c088da..2ce2c746 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -114,13 +114,13 @@ def getparser():
         help='reference period ending year for calibration (typically 2019)',
     )
     (
-        parser.add_argument(
-            '-rgi_glac_number_fn',
-            action='store',
-            type=str,
-            default=None,
-            help='filepath containing list of rgi_glac_number, helpful for running batches on spc',
-        )
+    parser.add_argument(
+        '-rgi_glac_number_fn',
+        action='store',
+        type=str,
+        default=None,
+        help='filepath containing list of rgi_glac_number, helpful for running batches on spc',
+    )
     )
     parser.add_argument(
         '-rgi_glac_number',
@@ -179,6 +179,12 @@ def getparser():
         action='store_true',
         help='Flag to keep glacier lists ordered (default is false)',
     )
+    parser.add_argument(
+        '-option_calib_glacierwide_mb_mwea',
+        action='store_true',
+        default=pygem_prms['calib']['MCMC_params']['option_calib_glacierwide_mb_mwea'],
+        help='Flag to calibrate against average glacierwide mass balance',
+    )
     parser.add_argument(
         '-option_calib_elev_change_1d',
         action='store_true',
@@ -252,27 +258,20 @@ def mb_mwea_calc(
         return mb_mwea
 
 
-def calculate_elev_change_1d(
-    gdir,
-    modelprms,
-    glacier_rgi_table,
-    fls,
-    debug=False,
-):
-    """
-    For a given set of model parameters, run the ice thickness inversion and mass balance model to get binned annual ice thickness change
-    Convert to monthly thickness by assuming that the flux divergence is constant throughout the year
-    """
+def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
+    """run the dynamical evolution model with a given set of model parameters"""
+    
     y0 = gdir.dates_table.year.min()
     y1 = gdir.dates_table.year.max()
 
+    # mass balance model with evolving area
+    mbmod = PyGEMMassBalance(gdir, modelprms, glacier_rgi_table, fls=fls)
+
     # Check that water level is within given bounds
     cls = gdir.read_pickle('inversion_input')[-1]
     th = cls['hgt'][-1]
     vmin, vmax = cfg.PARAMS['free_board_marine_terminating']
     water_level = utils.clip_scalar(0, th - vmax, th - vmin)
-    # mass balance model with evolving area
-    mbmod = PyGEMMassBalance(gdir, modelprms, glacier_rgi_table, fls=fls)
     # glacier dynamics model
     if gdir.is_tidewater and pygem_prms['setup']['include_frontalablation']:
         ev_model = flowline.FluxBasedModel(
@@ -313,22 +312,17 @@ def calculate_elev_change_1d(
                 for n in np.arange(calving_m3we_annual.shape[0]):
                     ev_model.mb_model.glac_wide_frontalablation[12 * n + 11] = calving_m3we_annual[n]
 
-                # Add mass lost from frontal ablation to Glacier-wide total mass balance (m3 w.e.)
+                # add mass lost from frontal ablation to Glacier-wide total mass balance (m3 w.e.)
                 ev_model.mb_model.glac_wide_massbaltotal = (
                     ev_model.mb_model.glac_wide_massbaltotal + ev_model.mb_model.glac_wide_frontalablation
                 )
-
-            mod_glacierwide_mb_mwea = (
-                mbmod.glac_wide_massbaltotal[gdir.mbdata['t1_idx'] : gdir.mbdata['t2_idx'] + 1].sum()
-                / mbmod.glac_wide_area_annual[0]
-                / gdir.mbdata['nyears']
-            )
-
-    # if there is an issue evaluating the dynamics model for a given parameter set in MCMC calibration,
-    # return -inf for mb_mwea and binned_dh, so MCMC calibration won't accept given parameters
     except RuntimeError:
-        return float('-inf'), float('-inf')
+        ds = None
+    
+    return mbmod, ds
+
 
+def calc_elev_change_1d(gdir, mbmod, ds):
     ### get monthly ice thickness
     # grab components of interest
     thickness_m = ds[0].thickness_m.values.T  # glacier thickness [m ice], (nbins, nyears)
@@ -375,7 +369,7 @@ def calculate_elev_change_1d(
     h_monthly_ = np.column_stack([interp1d_fill_gaps(x.copy()) for x in h_monthly.T])
 
     # difference each set of inds in diff_inds_map
-    mod_elev_change_1d = np.column_stack(
+    elev_change_1d = np.column_stack(
         [
             h_monthly_[:, tup[1]] - h_monthly_[:, tup[0]]
             if tup[0] is not None and tup[1] is not None
@@ -383,8 +377,71 @@ def calculate_elev_change_1d(
             for tup in gdir.elev_change_1d['model2obs_inds_map']
         ]
     )
+    return elev_change_1d
+
+
+def evaluate_model_outputs(
+    gdir,
+    modelprms,
+    glacier_rgi_table,
+    fls,
+    mbfxn=None,
+    calib_glacierwide_mb_mwea=True,
+    calib_elev_change_1d=False,
+    calib_snowlines_1d=False,
+    calib_meltextent_1d=False,
+    debug=False,
+):
+    """
+    For a given set of model parameters, evaluate the desired model outputs.
+    Optionally use an emulator function to compute mass balance.
+    Returns a dictionary with only the requested results.
+    """
+    results = {}
+    mbmod = None
+
+    if calib_elev_change_1d:
+        mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls)
+        results["elev_change_1d"] = (
+            calc_elev_change_1d(gdir, mbmod, ds) if ds else float("-inf")
+        )
+
+    if calib_glacierwide_mb_mwea:
+        if mbfxn is not None:
+            # grab current values from modelprms for the emulator
+            mb_args = [
+                modelprms['tbias'],
+                modelprms['kp'],
+                modelprms['ddfsnow']
+            ]
+            glacierwide_mb_mwea = mbfxn(*[mb_args])
+        else:
+            if mbmod is None:
+                glacierwide_mb_mwea = mb_mwea_calc(gdir, modelprms, glacier_rgi_table, fls)
+            else:
+                glacierwide_mb_mwea = (
+                    mbmod.glac_wide_massbaltotal[
+                        gdir.mbdata["t1_idx"]: gdir.mbdata["t2_idx"] + 1
+                    ].sum()
+                    / mbmod.glac_wide_area_annual[0]
+                    / gdir.mbdata["nyears"]
+                )
 
-    return mod_glacierwide_mb_mwea, mod_elev_change_1d
+        results["glacierwide_mb_mwea"] = glacierwide_mb_mwea
+
+    # (add future calibration options here)
+    if calib_snowlines_1d:
+        pass
+        # results["snowlines_1d"] = calc_snowlines_1d(gdir, mbmod)
+
+    if calib_meltextent_1d:
+        pass
+        # results["meltextent_1d"] = calc_meltextent_1d(gdir, mbmod)
+
+    if debug:
+        print("Returned keys:", list(results.keys()))
+
+    return results
 
 
 # class for Gaussian Process model for mass balance emulator
@@ -2010,8 +2067,14 @@ def rho_constraints(**kwargs):
                 # -------------------
                 # --- set up MCMC ---
                 # -------------------
-                # mass balance observation and standard deviation
-                obs = [(torch.tensor([mb_obs_mwea]), torch.tensor([mb_obs_mwea_err]))]
+                # the mcmc.mbPosterior class expects observations to be provided as a dictionary,
+                # where each key corresponds to a variable being calibrated.
+                # each value should be a tuple of the form (observation, variance).
+                obs = {'glacierwide_mb_mwea':(torch.tensor([mb_obs_mwea]), torch.tensor([mb_obs_mwea_err]))}
+                mbfxn = None
+
+                if pygem_prms['calib']['MCMC_params']['option_use_emulator']:
+                    mbfxn = mbEmulator.eval  # returns (mb_mwea)
 
                 # if running full model (no emulator), or calibrating against binned elevation change, several arguments are needed
                 if args.option_calib_elev_change_1d:
@@ -2036,34 +2099,20 @@ def rho_constraints(**kwargs):
                         ),
                     )
 
-                    # calculate inds of data v. model
-                    mbfxn = calculate_elev_change_1d  # returns (mb_mwea, binned_dm)
-                    mbargs = (
-                        gdir,  # arguments for get_binned_dh()
-                        modelprms,
-                        glacier_rgi_table,
-                        fls,
-                    )
-                    # append deltah obs and and sigma obs list
-                    obs.append(
-                        (
-                            torch.tensor(gdir.elev_change_1d['dh']),
-                            torch.tensor(gdir.elev_change_1d['dh_sigma']),
-                        )
-                    )
-                # if there are more observations to calibrate against, simply append a tuple of (obs, variance) to obs list
-                # e.g. obs.append((torch.tensor(dmda_array),torch.tensor(dmda_err_array)))
-                elif pygem_prms['calib']['MCMC_params']['option_use_emulator']:
-                    mbfxn = mbEmulator.eval  # returns (mb_mwea)
-                    mbargs = None  # no additional arguments for mbEmulator.eval()
-                else:
-                    mbfxn = mb_mwea_calc  # returns (mb_mwea)
-                    mbargs = (
-                        gdir,
-                        modelprms,
-                        glacier_rgi_table,
-                        fls,
-                    )  # arguments for mb_mwea_calc()
+                    # add elevation change data to observations dictionary
+                    obs['elev_change_1d'] = (torch.tensor(gdir.elev_change_1d['dh']), torch.tensor(gdir.elev_change_1d['dh_sigma']))
+                # if there are more observations to calibrate against, simply add them as a tuple of (obs, variance) to the obs dictionary
+
+                # define args to pass to fxn2eval in mcmc sampler
+                fxnargs = (
+                    gdir,
+                    modelprms,
+                    glacier_rgi_table,
+                    fls,
+                    mbfxn,
+                    args.option_calib_glacierwide_mb_mwea,
+                    args.option_calib_elev_change_1d,
+                )
 
                 # instantiate mbPosterior given priors, and observed values
                 # note, mbEmulator.eval expects the modelprms to be ordered like so: [tbias, kp, ddfsnow], so priors and initial guesses must also be ordered as such)
@@ -2071,8 +2120,8 @@ def rho_constraints(**kwargs):
                 mb = mcmc.mbPosterior(
                     obs,
                     priors,
-                    mb_func=mbfxn,
-                    mb_args=mbargs,
+                    fxn2eval=evaluate_model_outputs,
+                    fxnargs=fxnargs,
                     potential_fxns=[mb_max, must_melt, rho_constraints],
                     ela=gdir.ela.min() if hasattr(gdir, 'ela') else None,
                     bin_z=gdir.elev_change_1d['bin_centers'] if hasattr(gdir, 'elev_change_1d') else None,
@@ -2089,7 +2138,7 @@ def rho_constraints(**kwargs):
                 if args.option_calib_elev_change_1d:
                     modelprms_export['elev_change_1d'] = {}
                     modelprms_export['elev_change_1d']['bin_edges'] = gdir.elev_change_1d['bin_edges']
-                    modelprms_export['elev_change_1d']['obs'] = [ob.flatten().tolist() for ob in obs[1]]
+                    modelprms_export['elev_change_1d']['obs'] = [ob.flatten().tolist() for ob in obs['elev_change_1d']]
                     modelprms_export['elev_change_1d']['dates'] = [
                         (dt1, dt2) for dt1, dt2 in gdir.elev_change_1d['dates']
                     ]
@@ -2178,9 +2227,9 @@ def rho_constraints(**kwargs):
                             )
 
                         # concatenate mass balance
-                        m_chain = torch.cat((m_chain, torch.tensor(pred_chain[0]).reshape(-1, 1)), dim=1)
+                        m_chain = torch.cat((m_chain, torch.tensor(pred_chain['glacierwide_mb_mwea']).reshape(-1, 1)), dim=1)
                         m_primes = torch.cat(
-                            (m_primes, torch.tensor(pred_primes[0]).reshape(-1, 1)),
+                            (m_primes, torch.tensor(pred_primes['glacierwide_mb_mwea']).reshape(-1, 1)),
                             dim=1,
                         )
 
@@ -2213,23 +2262,23 @@ def rho_constraints(**kwargs):
                                 graphics.plot_mcmc_chain(
                                     m_primes,
                                     m_chain,
-                                    obs[0],
+                                    obs,
                                     ar,
                                     glacier_str,
                                     show=show,
                                     fpath=f'{fp}/{glacier_str}-chain{n_chain}.png',
                                 )
-                                for i in pred_chain.keys():
+                                for k in pred_chain.keys():
                                     graphics.plot_resid_histogram(
-                                        obs[i],
-                                        pred_chain[i],
+                                        obs[k],
+                                        pred_chain[k],
                                         glacier_str,
                                         show=show,
-                                        fpath=f'{fp}/{glacier_str}-chain{n_chain}-residuals-{i}.png',
+                                        fpath=f'{fp}/{glacier_str}-chain{n_chain}-residuals-{k}.png',
                                     )
-                                    if i == 1:
+                                    if k == 'elev_change_1d':
                                         graphics.plot_mcmc_elev_change_1d(
-                                            pred_chain[1],
+                                            pred_chain[k],
                                             fls,
                                             gdir.elev_change_1d,
                                             gdir.ela.min(),
@@ -2253,7 +2302,7 @@ def rho_constraints(**kwargs):
                         modelprms_export['ar'][chain_str] = ar
                         if args.option_calib_elev_change_1d:
                             modelprms_export['elev_change_1d'][chain_str] = [
-                                preds.flatten().tolist() for preds in pred_chain[1]
+                                preds.flatten().tolist() for preds in pred_chain['elev_change_1d']
                             ]
                             modelprms_export['rhoabl'][chain_str] = m_chain[:, 3].tolist()
                             modelprms_export['rhoacc'][chain_str] = m_chain[:, 4].tolist()
diff --git a/pygem/mcmc.py b/pygem/mcmc.py
index f149b9a2..647e222f 100644
--- a/pygem/mcmc.py
+++ b/pygem/mcmc.py
@@ -150,12 +150,12 @@ def log_uniform_density(x, **kwargs):
 
 # mass balance posterior class
 class mbPosterior:
-    def __init__(self, obs, priors, mb_func, mb_args=None, potential_fxns=None, **kwargs):
+    def __init__(self, obs, priors, fxn2eval, fxnargs=None, potential_fxns=None, **kwargs):
         # obs will be passed as a list, where each item is a tuple with the first element being the mean observation, and the second being the variance
         self.obs = obs
         self.priors = copy.deepcopy(priors)
-        self.mb_func = mb_func
-        self.mb_args = mb_args
+        self.fxn2eval = fxn2eval
+        self.fxnargs = fxnargs
         self.potential_functions = potential_fxns if potential_fxns is not None else []
         self.preds = None
         self.check_priors()
@@ -197,19 +197,15 @@ def check_priors(self):
     # update modelprms for evaluation
     def update_modelprms(self, m):
         for i, k in enumerate(['tbias', 'kp', 'ddfsnow']):
-            self.mb_args[1][k] = float(m[i])
-        self.mb_args[1]['ddfice'] = self.mb_args[1]['ddfsnow'] / pygem_prms['sim']['params']['ddfsnow_iceratio']
+            self.fxnargs[1][k] = float(m[i])
+        self.fxnargs[1]['ddfice'] = self.fxnargs[1]['ddfsnow'] / pygem_prms['sim']['params']['ddfsnow_iceratio']
 
-    # get mb_pred
+    # get model predictions
     def get_model_pred(self, m):
-        if self.mb_args:
-            self.update_modelprms(m)
-            self.preds = self.mb_func(*self.mb_args)
-        else:
-            self.preds = self.mb_func([*m])
-        if not isinstance(self.preds, tuple):
-            self.preds = [self.preds]
-        self.preds = [torch.tensor(item) for item in self.preds]  # make all preds torch.tensor() objects
+        self.update_modelprms(m)           # update modelprms with current step
+        self.preds = self.fxn2eval(*self.fxnargs)
+        # convert all values to torch tensors
+        self.preds = {k: torch.tensor(v, dtype=torch.float) for k, v in self.preds.items()}
 
     # get total log prior density
     def log_prior(self, m):
@@ -225,36 +221,30 @@ def log_prior(self, m):
     # get log likelihood
     def log_likelihood(self, m):
         log_likehood = 0
-        for i, pred in enumerate(self.preds):
-            # --- Check for invalid predictions early ---
+        for k, pred in self.preds.items():
+            # --- Check for invalid predictions  ---
             if torch.all(pred == float('-inf')):
                 # Invalid model output -> assign -inf likelihood
                 return torch.tensor([-float('inf')])
 
-            if i == 0:
-                # --- Base case: mass balance likelihood ---
-                log_likehood += log_normal_density(
-                    self.obs[i][0],  # observed values
-                    mu=pred,  # predicted values
-                    sigma=self.obs[i][1],  # observation uncertainty
-                )
-
-            elif i == 1 and len(m) > 3:
-                # --- Extended case: apply density scaling to get binned elevation change ---
+            # if key is `elev_change_1d` scale by density to predict binned surface elevation change
+            if k == 'elev_change_1d':
                 # Create density field, separate values for ablation/accumulation zones
                 rho = np.ones_like(self.bin_z)
                 rho[self.abl_mask] = m[3]  # rhoabl
                 rho[~self.abl_mask] = m[4]  # rhoacc
                 rho = torch.tensor(rho)
-                self.preds[i] = pred = (
-                    self.preds[i] * (pygem_prms['constants']['density_ice'] / rho[:, np.newaxis])
-                )  # scale prediction by model density values (convert from m ice to m thickness change considering modeled density)
-
-                log_likehood += log_normal_density(
-                    self.obs[i][0],  # observations
-                    mu=pred,  # scaled predictions
-                    sigma=self.obs[i][1],  # uncertainty
-                )
+                # scale prediction by model density values (convert from m ice to m surface elevation change considering modeled density)
+                pred *= (pygem_prms['constants']['density_ice'] / rho[:, np.newaxis])
+                # update values in preds dict
+                self.preds[k] = pred
+
+            log_likehood += log_normal_density(
+                self.obs[k][0],  # observations
+                mu=pred,  # scaled predictions
+                sigma=self.obs[k][1],  # uncertainty
+            )
+
         return log_likehood
 
     # compute the log-potential, summing over all declared potential functions.
@@ -265,7 +255,7 @@ def log_potential(self, m):
             'kp': m[0],
             'tbias': m[1],
             'ddfsnow': m[2],
-            'massbal': self.preds[0],
+            'massbal': self.preds['glacierwide_mb_mwea'],
         }
 
         # --- Optional arguments(if len(m) > 3) ---
@@ -332,9 +322,9 @@ def rm_stuck_samples(self):
         self.m_primes = self.m_primes[self.n_rm :]
         self.steps = self.steps[self.n_rm :]
         self.acceptance = self.acceptance[self.n_rm :]
-        for j in self.preds_primes.keys():
-            self.preds_primes[j] = self.preds_primes[j][self.n_rm :]
-            self.preds_chain[j] = self.preds_chain[j][self.n_rm :]
+        for k in self.preds_primes.keys():
+            self.preds_primes[k] = self.preds_primes[k][self.n_rm :]
+            self.preds_chain[k] = self.preds_chain[k][self.n_rm :]
         return
 
     def sample(
@@ -394,12 +384,12 @@ def sample(
                     self.m_chain.append(m_0)
                     self.m_primes.append(m_prime)
                     self.acceptance.append(self.naccept / (i + (thin_factor * self.n_rm)))
-                    for j in range(len(pred_1)):
-                        if j not in self.preds_chain.keys():
-                            self.preds_chain[j] = []
-                            self.preds_primes[j] = []
-                        self.preds_chain[j].append(pred_0[j])
-                        self.preds_primes[j].append(pred_1[j])
+                    for k, values in pred_1.items():
+                        if k not in self.preds_chain.keys():
+                            self.preds_chain[k] = []
+                            self.preds_primes[k] = []
+                        self.preds_chain[k].append(pred_0[k])
+                        self.preds_primes[k].append(pred_1[k])
 
             # trim off any initial steps that are stagnant
             if (i == (n_samples - 1)) and (trim):
diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py
index 74f0bbba..20b49951 100644
--- a/pygem/plot/graphics.py
+++ b/pygem/plot/graphics.py
@@ -113,7 +113,7 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs):
     ax.set_title(title, loc='left')
 
 
-def plot_mcmc_chain(m_primes, m_chain, mb_obs, ar, title, ms=1, fontsize=8, show=False, fpath=None):
+def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=False, fpath=None):
     # Plot the trace of the parameters
     n = m_primes.shape[1]
     fig, axes = plt.subplots(n + 1, 1, figsize=(6, n * 1.5), sharex=True)
@@ -188,30 +188,32 @@ def plot_mcmc_chain(m_primes, m_chain, mb_obs, ar, title, ms=1, fontsize=8, show
         legs.append(l4)
         axes[4].set_ylabel(r'$\rho_{acc}$', fontsize=fontsize)
 
-    axes[-2].fill_between(
-        np.arange(len(ar)),
-        mb_obs[0] - (2 * mb_obs[1]),
-        mb_obs[0] + (2 * mb_obs[1]),
-        color='grey',
-        alpha=0.3,
-    )
-    axes[-2].fill_between(
-        np.arange(len(ar)),
-        mb_obs[0] - mb_obs[1],
-        mb_obs[0] + mb_obs[1],
-        color='grey',
-        alpha=0.3,
-    )
-    axes[-2].plot(m_primes[:, -1], '.', ms=ms, c='tab:blue')
-    axes[-2].plot(m_chain[:, -1], '.', ms=ms, c='tab:orange')
-    axes[-2].plot(
-        [],
-        [],
-        label=f'median={np.median(m_chain[:, -1]):.3f}\niqr={np.subtract(*np.percentile(m_chain[:, -1], [75, 25])):.3f}',
-    )
-    ln2 = axes[-2].legend(loc='upper right', handlelength=0, borderaxespad=0, fontsize=fontsize)
-    legs.append(ln2)
-    axes[-2].set_ylabel(r'$\dot{{b}}$', fontsize=fontsize)
+    if 'glacierwide_mb_mwea' in obs.keys():
+        mb_obs = obs['glacierwide_mb_mwea']
+        axes[-2].fill_between(
+            np.arange(len(ar)),
+            mb_obs[0] - (2 * mb_obs[1]),
+            mb_obs[0] + (2 * mb_obs[1]),
+            color='grey',
+            alpha=0.3,
+        )
+        axes[-2].fill_between(
+            np.arange(len(ar)),
+            mb_obs[0] - mb_obs[1],
+            mb_obs[0] + mb_obs[1],
+            color='grey',
+            alpha=0.3,
+        )
+        axes[-2].plot(m_primes[:, -1], '.', ms=ms, c='tab:blue')
+        axes[-2].plot(m_chain[:, -1], '.', ms=ms, c='tab:orange')
+        axes[-2].plot(
+            [],
+            [],
+            label=f'median={np.median(m_chain[:, -1]):.3f}\niqr={np.subtract(*np.percentile(m_chain[:, -1], [75, 25])):.3f}',
+        )
+        ln2 = axes[-2].legend(loc='upper right', handlelength=0, borderaxespad=0, fontsize=fontsize)
+        legs.append(ln2)
+        axes[-2].set_ylabel(r'$\dot{{b}}$', fontsize=fontsize)
 
     axes[-1].plot(ar, 'tab:orange', lw=1)
     axes[-1].plot(
diff --git a/pygem/setup/config.py b/pygem/setup/config.py
index c60f520c..fc16cfaf 100644
--- a/pygem/setup/config.py
+++ b/pygem/setup/config.py
@@ -220,6 +220,7 @@ def _validate_config(self, config):
         'calib.emulator_params.ftol_opt': float,
         'calib.emulator_params.eps_opt': float,
         'calib.MCMC_params': dict,
+        'calib.MCMC_params.option_calib_glacierwide_mb_mwea': bool,
         'calib.MCMC_params.option_use_emulator': bool,
         'calib.MCMC_params.emulator_sims': int,
         'calib.MCMC_params.tbias_step': float,
diff --git a/pygem/setup/config.yaml b/pygem/setup/config.yaml
index 83bf12b8..bb44ed8c 100644
--- a/pygem/setup/config.yaml
+++ b/pygem/setup/config.yaml
@@ -150,6 +150,7 @@ calib:
   
   # MCMC params
   MCMC_params:
+    option_calib_glacierwide_mb_mwea: true        # option to calibrate against average glacierwide mass balance data (true or false)
     option_use_emulator: true                     # use emulator or full model (if true, calibration must have first been run with option_calibretion=='emulator')
     emulator_sims: 100
     tbias_step: 0.1

From 5937399984448b6bbf6173023a43a11d400ff844 Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Wed, 29 Oct 2025 19:43:16 -0400
Subject: [PATCH 02/11] Formatting

---
 pygem/bin/run/run_calibration.py | 49 +++++++++++++++-----------------
 pygem/mcmc.py                    |  4 +--
 2 files changed, 25 insertions(+), 28 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 2ce2c746..263a0523 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -114,13 +114,13 @@ def getparser():
         help='reference period ending year for calibration (typically 2019)',
     )
     (
-    parser.add_argument(
-        '-rgi_glac_number_fn',
-        action='store',
-        type=str,
-        default=None,
-        help='filepath containing list of rgi_glac_number, helpful for running batches on spc',
-    )
+        parser.add_argument(
+            '-rgi_glac_number_fn',
+            action='store',
+            type=str,
+            default=None,
+            help='filepath containing list of rgi_glac_number, helpful for running batches on spc',
+        )
     )
     parser.add_argument(
         '-rgi_glac_number',
@@ -260,7 +260,7 @@ def mb_mwea_calc(
 
 def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
     """run the dynamical evolution model with a given set of model parameters"""
-    
+
     y0 = gdir.dates_table.year.min()
     y1 = gdir.dates_table.year.max()
 
@@ -318,7 +318,7 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
                 )
     except RuntimeError:
         ds = None
-    
+
     return mbmod, ds
 
 
@@ -402,32 +402,24 @@ def evaluate_model_outputs(
 
     if calib_elev_change_1d:
         mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls)
-        results["elev_change_1d"] = (
-            calc_elev_change_1d(gdir, mbmod, ds) if ds else float("-inf")
-        )
+        results['elev_change_1d'] = calc_elev_change_1d(gdir, mbmod, ds) if ds else float('-inf')
 
     if calib_glacierwide_mb_mwea:
         if mbfxn is not None:
             # grab current values from modelprms for the emulator
-            mb_args = [
-                modelprms['tbias'],
-                modelprms['kp'],
-                modelprms['ddfsnow']
-            ]
+            mb_args = [modelprms['tbias'], modelprms['kp'], modelprms['ddfsnow']]
             glacierwide_mb_mwea = mbfxn(*[mb_args])
         else:
             if mbmod is None:
                 glacierwide_mb_mwea = mb_mwea_calc(gdir, modelprms, glacier_rgi_table, fls)
             else:
                 glacierwide_mb_mwea = (
-                    mbmod.glac_wide_massbaltotal[
-                        gdir.mbdata["t1_idx"]: gdir.mbdata["t2_idx"] + 1
-                    ].sum()
+                    mbmod.glac_wide_massbaltotal[gdir.mbdata['t1_idx'] : gdir.mbdata['t2_idx'] + 1].sum()
                     / mbmod.glac_wide_area_annual[0]
-                    / gdir.mbdata["nyears"]
+                    / gdir.mbdata['nyears']
                 )
 
-        results["glacierwide_mb_mwea"] = glacierwide_mb_mwea
+        results['glacierwide_mb_mwea'] = glacierwide_mb_mwea
 
     # (add future calibration options here)
     if calib_snowlines_1d:
@@ -439,7 +431,7 @@ def evaluate_model_outputs(
         # results["meltextent_1d"] = calc_meltextent_1d(gdir, mbmod)
 
     if debug:
-        print("Returned keys:", list(results.keys()))
+        print('Returned keys:', list(results.keys()))
 
     return results
 
@@ -2070,7 +2062,7 @@ def rho_constraints(**kwargs):
                 # the mcmc.mbPosterior class expects observations to be provided as a dictionary,
                 # where each key corresponds to a variable being calibrated.
                 # each value should be a tuple of the form (observation, variance).
-                obs = {'glacierwide_mb_mwea':(torch.tensor([mb_obs_mwea]), torch.tensor([mb_obs_mwea_err]))}
+                obs = {'glacierwide_mb_mwea': (torch.tensor([mb_obs_mwea]), torch.tensor([mb_obs_mwea_err]))}
                 mbfxn = None
 
                 if pygem_prms['calib']['MCMC_params']['option_use_emulator']:
@@ -2100,7 +2092,10 @@ def rho_constraints(**kwargs):
                     )
 
                     # add elevation change data to observations dictionary
-                    obs['elev_change_1d'] = (torch.tensor(gdir.elev_change_1d['dh']), torch.tensor(gdir.elev_change_1d['dh_sigma']))
+                    obs['elev_change_1d'] = (
+                        torch.tensor(gdir.elev_change_1d['dh']),
+                        torch.tensor(gdir.elev_change_1d['dh_sigma']),
+                    )
                 # if there are more observations to calibrate against, simply add them as a tuple of (obs, variance) to the obs dictionary
 
                 # define args to pass to fxn2eval in mcmc sampler
@@ -2227,7 +2222,9 @@ def rho_constraints(**kwargs):
                             )
 
                         # concatenate mass balance
-                        m_chain = torch.cat((m_chain, torch.tensor(pred_chain['glacierwide_mb_mwea']).reshape(-1, 1)), dim=1)
+                        m_chain = torch.cat(
+                            (m_chain, torch.tensor(pred_chain['glacierwide_mb_mwea']).reshape(-1, 1)), dim=1
+                        )
                         m_primes = torch.cat(
                             (m_primes, torch.tensor(pred_primes['glacierwide_mb_mwea']).reshape(-1, 1)),
                             dim=1,
diff --git a/pygem/mcmc.py b/pygem/mcmc.py
index 647e222f..f84ee5a4 100644
--- a/pygem/mcmc.py
+++ b/pygem/mcmc.py
@@ -202,7 +202,7 @@ def update_modelprms(self, m):
 
     # get model predictions
     def get_model_pred(self, m):
-        self.update_modelprms(m)           # update modelprms with current step
+        self.update_modelprms(m)  # update modelprms with current step
         self.preds = self.fxn2eval(*self.fxnargs)
         # convert all values to torch tensors
         self.preds = {k: torch.tensor(v, dtype=torch.float) for k, v in self.preds.items()}
@@ -235,7 +235,7 @@ def log_likelihood(self, m):
                 rho[~self.abl_mask] = m[4]  # rhoacc
                 rho = torch.tensor(rho)
                 # scale prediction by model density values (convert from m ice to m surface elevation change considering modeled density)
-                pred *= (pygem_prms['constants']['density_ice'] / rho[:, np.newaxis])
+                pred *= pygem_prms['constants']['density_ice'] / rho[:, np.newaxis]
                 # update values in preds dict
                 self.preds[k] = pred
 

From f32c55c52b38d03aa4c3bdb4ec5775b8b8660ceb Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Wed, 29 Oct 2025 20:12:58 -0400
Subject: [PATCH 03/11] Formatting

---
 pygem/bin/run/run_calibration.py | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 263a0523..8a3b2f32 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -113,14 +113,12 @@ def getparser():
         default=pygem_prms['climate']['ref_endyear'],
         help='reference period ending year for calibration (typically 2019)',
     )
-    (
-        parser.add_argument(
-            '-rgi_glac_number_fn',
-            action='store',
-            type=str,
-            default=None,
-            help='filepath containing list of rgi_glac_number, helpful for running batches on spc',
-        )
+    parser.add_argument(
+        '-rgi_glac_number_fn',
+        action='store',
+        type=str,
+        default=None,
+        help='filepath containing list of rgi_glac_number, helpful for running batches on spc',
     )
     parser.add_argument(
         '-rgi_glac_number',
@@ -380,7 +378,7 @@ def calc_elev_change_1d(gdir, mbmod, ds):
     return elev_change_1d
 
 
-def evaluate_model_outputs(
+def mcmc_model_eval(
     gdir,
     modelprms,
     glacier_rgi_table,
@@ -2115,7 +2113,7 @@ def rho_constraints(**kwargs):
                 mb = mcmc.mbPosterior(
                     obs,
                     priors,
-                    fxn2eval=evaluate_model_outputs,
+                    fxn2eval=mcmc_model_eval,
                     fxnargs=fxnargs,
                     potential_fxns=[mb_max, must_melt, rho_constraints],
                     ela=gdir.ela.min() if hasattr(gdir, 'ela') else None,

From b6cc4ce6960dfe17fbf1225e6b5f78b7cc64c7d2 Mon Sep 17 00:00:00 2001
From: "Brandon S. Tober" <tobers.brandon@gmail.com>
Date: Thu, 30 Oct 2025 09:17:02 -0400
Subject: [PATCH 04/11] Add average along-chain residual to mcmc diagnostics
 plot (#157)

Closes #154
---
 pygem/bin/run/run_calibration.py | 15 ++-----
 pygem/plot/graphics.py           | 72 ++++++++++++++++++++++++--------
 2 files changed, 59 insertions(+), 28 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 8a3b2f32..e6587bdd 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -2219,21 +2219,12 @@ def rho_constraints(**kwargs):
                                 f'Chain {n_chain}: failed to produce an unstuck result after {attempts_per_chain} initial guesses.'
                             )
 
-                        # concatenate mass balance
-                        m_chain = torch.cat(
-                            (m_chain, torch.tensor(pred_chain['glacierwide_mb_mwea']).reshape(-1, 1)), dim=1
-                        )
-                        m_primes = torch.cat(
-                            (m_primes, torch.tensor(pred_primes['glacierwide_mb_mwea']).reshape(-1, 1)),
-                            dim=1,
-                        )
-
                         if debug:
                             print(
                                 'mb_mwea_mean:',
-                                np.round(torch.mean(m_chain[:, -1]).item(), 3),
+                                np.round(torch.mean(torch.stack(pred_chain['glacierwide_mb_mwea'])).item(), 3),
                                 'mb_mwea_std:',
-                                np.round(torch.std(m_chain[:, -1]).item(), 3),
+                                np.round(torch.std(torch.stack(pred_chain['glacierwide_mb_mwea'])).item(), 3),
                                 '\nmb_obs_mean:',
                                 np.round(mb_obs_mwea, 3),
                                 'mb_obs_std:',
@@ -2257,6 +2248,8 @@ def rho_constraints(**kwargs):
                                 graphics.plot_mcmc_chain(
                                     m_primes,
                                     m_chain,
+                                    pred_primes,
+                                    pred_chain,
                                     obs,
                                     ar,
                                     glacier_str,
diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py
index 20b49951..b0d37bb0 100644
--- a/pygem/plot/graphics.py
+++ b/pygem/plot/graphics.py
@@ -13,6 +13,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import torch
 from scipy.stats import binned_statistic
 
 from pygem.utils.stats import effective_n
@@ -113,10 +114,15 @@ def plot_modeloutput_section(model=None, ax=None, title='', **kwargs):
     ax.set_title(title, loc='left')
 
 
-def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=False, fpath=None):
+def plot_mcmc_chain(
+    m_primes, m_chain, pred_primes, pred_chain, obs, ar, title, ms=1, fontsize=8, show=False, fpath=None
+):
     # Plot the trace of the parameters
-    n = m_primes.shape[1]
-    fig, axes = plt.subplots(n + 1, 1, figsize=(6, n * 1.5), sharex=True)
+    nparams = m_primes.shape[1]
+    npreds = len(pred_chain.keys())
+    N = nparams + npreds + 1
+    fig, axes = plt.subplots(N, 1, figsize=(6, N * 1), sharex=True)
+    # convert torch objects to numpy
     m_chain = m_chain.detach().numpy()
     m_primes = m_primes.detach().numpy()
 
@@ -125,6 +131,7 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
     # instantiate list to hold legend objs
     legs = []
 
+    # axes[0] will always be tbias
     axes[0].plot(
         [],
         [],
@@ -139,6 +146,7 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
     # axes[0].add_artist(leg)
     axes[0].set_ylabel(r'$T_{bias}$', fontsize=fontsize)
 
+    # axes[1] will always be kp
     axes[1].plot(m_primes[:, 1], '.', ms=ms, c='tab:blue')
     axes[1].plot(m_chain[:, 1], '.', ms=ms, c='tab:orange')
     axes[1].plot(
@@ -150,6 +158,7 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
     legs.append(l1)
     axes[1].set_ylabel(r'$K_p$', fontsize=fontsize)
 
+    # axes[2] will always be ddfsnow
     axes[2].plot(m_primes[:, 2], '.', ms=ms, c='tab:blue')
     axes[2].plot(m_chain[:, 2], '.', ms=ms, c='tab:orange')
     axes[2].plot(
@@ -161,7 +170,8 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
     legs.append(l2)
     axes[2].set_ylabel(r'$fsnow$', fontsize=fontsize)
 
-    if n > 4:
+    if nparams > 3:
+        # axes[3] will be rho_ablation if more than 3 model params
         m_chain[:, 3] = m_chain[:, 3]
         m_primes[:, 3] = m_primes[:, 3]
         axes[3].plot(m_primes[:, 3], '.', ms=ms, c='tab:blue')
@@ -175,6 +185,7 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
         legs.append(l3)
         axes[3].set_ylabel(r'$\rho_{abl}$', fontsize=fontsize)
 
+        # axes[4] will be rho_accumulation if more than 3 model params
         m_chain[:, 4] = m_chain[:, 4]
         m_primes[:, 4] = m_primes[:, 4]
         axes[4].plot(m_primes[:, 4], '.', ms=ms, c='tab:blue')
@@ -188,33 +199,60 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
         legs.append(l4)
         axes[4].set_ylabel(r'$\rho_{acc}$', fontsize=fontsize)
 
-    if 'glacierwide_mb_mwea' in obs.keys():
+    # plot predictions
+    if 'glacierwide_mb_mwea' in pred_primes.keys():
         mb_obs = obs['glacierwide_mb_mwea']
-        axes[-2].fill_between(
+        axes[nparams].fill_between(
             np.arange(len(ar)),
             mb_obs[0] - (2 * mb_obs[1]),
             mb_obs[0] + (2 * mb_obs[1]),
             color='grey',
             alpha=0.3,
         )
-        axes[-2].fill_between(
+        axes[nparams].fill_between(
             np.arange(len(ar)),
             mb_obs[0] - mb_obs[1],
             mb_obs[0] + mb_obs[1],
             color='grey',
             alpha=0.3,
         )
-        axes[-2].plot(m_primes[:, -1], '.', ms=ms, c='tab:blue')
-        axes[-2].plot(m_chain[:, -1], '.', ms=ms, c='tab:orange')
-        axes[-2].plot(
+
+        mb_primes = torch.stack(pred_primes['glacierwide_mb_mwea']).numpy()
+        mb_chain = torch.stack(pred_chain['glacierwide_mb_mwea']).numpy()
+        axes[nparams].plot(mb_primes, '.', ms=ms, c='tab:blue')
+        axes[nparams].plot(mb_chain, '.', ms=ms, c='tab:orange')
+        axes[nparams].plot(
             [],
             [],
-            label=f'median={np.median(m_chain[:, -1]):.3f}\niqr={np.subtract(*np.percentile(m_chain[:, -1], [75, 25])):.3f}',
+            label=f'median={np.median(mb_chain):.3f}\niqr={np.subtract(*np.percentile(mb_chain, [75, 25])):.3f}',
         )
-        ln2 = axes[-2].legend(loc='upper right', handlelength=0, borderaxespad=0, fontsize=fontsize)
+        ln2 = axes[nparams].legend(loc='upper right', handlelength=0, borderaxespad=0, fontsize=fontsize)
         legs.append(ln2)
-        axes[-2].set_ylabel(r'$\dot{{b}}$', fontsize=fontsize)
+        axes[nparams].set_ylabel(r'$\dot{{b}}$', fontsize=fontsize)
+        nparams += 1
+
+    # plot MAE for all other prediction keys
+    for key in pred_primes.keys():
+        if key == 'glacierwide_mb_mwea':
+            continue
+        pred_primes = torch.stack(pred_primes[key]).numpy()
+        pred_chain = torch.stack(pred_chain[key]).numpy()
+        obs_vals = np.array(obs[key][0])
+
+        mae_primes = np.mean(pred_primes - obs_vals, axis=(1, 2))
+        mae_chain = np.mean(pred_chain - obs_vals, axis=(1, 2))
 
+        axes[nparams].plot(mae_primes, '.', ms=ms, c='tab:blue')
+        axes[nparams].plot(mae_chain, '.', ms=ms, c='tab:orange')
+
+        if key == 'elev_change_1d':
+            axes[nparams].set_ylabel(r'$\overline{\hat{dh} - dh}$', fontsize=fontsize)
+        else:
+            axes[nparams].set_ylabel(r'$\overline{\mathrm{pred} - \mathrm{obs}}$', fontsize=fontsize)
+        legs.append(None)
+        nparams += 1
+
+    # axes[-1] will always be acceptance rate
     axes[-1].plot(ar, 'tab:orange', lw=1)
     axes[-1].plot(
         np.convolve(ar, np.ones(100) / 100, mode='valid'),
@@ -230,7 +268,8 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
         ax.xaxis.set_ticks_position('both')
         ax.yaxis.set_ticks_position('both')
         ax.tick_params(axis='both', direction='inout')
-        if i == n:
+        ax.yaxis.set_label_coords(-0.1, 0.5)
+        if i > m_primes.shape[1] - 1:
             continue
         ax.plot([], [], label=f'n_eff={neff[i]}')
         hands, ls = ax.get_legend_handles_labels()
@@ -252,9 +291,8 @@ def plot_mcmc_chain(m_primes, m_chain, obs, ar, title, ms=1, fontsize=8, show=Fa
                 handlelength=0,
                 fontsize=fontsize,
             )
-
-    for i, ax in enumerate(axes):
-        ax.add_artist(legs[i])
+        if legs[i] is not None:
+            ax.add_artist(legs[i])
 
     axes[0].set_xlim([0, m_chain.shape[0]])
     axes[0].set_title(title, fontsize=fontsize)

From 50f1d9a747279404e9365e744d0c3202561d2306 Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Thu, 30 Oct 2025 17:23:44 -0400
Subject: [PATCH 05/11] Plot residual for any pred-obs pairs that aren't
 glacierwide mb, whether 2d or 1d

---
 pygem/plot/graphics.py | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py
index b0d37bb0..0e48d2c5 100644
--- a/pygem/plot/graphics.py
+++ b/pygem/plot/graphics.py
@@ -235,15 +235,24 @@ def plot_mcmc_chain(
     for key in pred_primes.keys():
         if key == 'glacierwide_mb_mwea':
             continue
+
+        # stack predictions first (shape: n_steps x ... x ...) - may end up being 2d or 3d
         pred_primes = torch.stack(pred_primes[key]).numpy()
         pred_chain = torch.stack(pred_chain[key]).numpy()
-        obs_vals = np.array(obs[key][0])
 
-        mae_primes = np.mean(pred_primes - obs_vals, axis=(1, 2))
-        mae_chain = np.mean(pred_chain - obs_vals, axis=(1, 2))
+        # flatten all axes except the first (n_steps) -> 2D array (n_steps, M)
+        pred_primes_flat = pred_primes.reshape(pred_primes.shape[0], -1)
+        pred_chain_flat = pred_chain.reshape(pred_chain.shape[0], -1)
+
+        # make obs array broadcastable (flatten if needed)
+        obs_vals_flat = np.ravel(np.array(obs[key][0]))
+
+        # compute mean residual per step
+        mean_resid_primes = np.mean(pred_primes_flat - obs_vals_flat, axis=1)
+        mean_resid_chain = np.mean(pred_chain_flat - obs_vals_flat, axis=1)
 
-        axes[nparams].plot(mae_primes, '.', ms=ms, c='tab:blue')
-        axes[nparams].plot(mae_chain, '.', ms=ms, c='tab:orange')
+        axes[nparams].plot(mean_resid_primes, '.', ms=ms, c='tab:blue')
+        axes[nparams].plot(mean_resid_chain, '.', ms=ms, c='tab:orange')
 
         if key == 'elev_change_1d':
             axes[nparams].set_ylabel(r'$\overline{\hat{dh} - dh}$', fontsize=fontsize)

From 1938ec3ab5afb54ead73dc6467425723406a85df Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Thu, 30 Oct 2025 18:37:06 -0400
Subject: [PATCH 06/11] use np.nanmean to mask any nans in preds or obs arrays

---
 pygem/bin/run/run_calibration.py | 3 ++-
 pygem/plot/graphics.py           | 4 ++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index e6587bdd..05b170f2 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -314,7 +314,8 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
                 ev_model.mb_model.glac_wide_massbaltotal = (
                     ev_model.mb_model.glac_wide_massbaltotal + ev_model.mb_model.glac_wide_frontalablation
                 )
-    except RuntimeError:
+    # safely catch any errors with dynamical run
+    except Exception:
         ds = None
 
     return mbmod, ds
diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py
index 0e48d2c5..99ed9d84 100644
--- a/pygem/plot/graphics.py
+++ b/pygem/plot/graphics.py
@@ -248,8 +248,8 @@ def plot_mcmc_chain(
         obs_vals_flat = np.ravel(np.array(obs[key][0]))
 
         # compute mean residual per step
-        mean_resid_primes = np.mean(pred_primes_flat - obs_vals_flat, axis=1)
-        mean_resid_chain = np.mean(pred_chain_flat - obs_vals_flat, axis=1)
+        mean_resid_primes = np.nanmean(pred_primes_flat - obs_vals_flat, axis=1)
+        mean_resid_chain = np.nanmean(pred_chain_flat - obs_vals_flat, axis=1)
 
         axes[nparams].plot(mean_resid_primes, '.', ms=ms, c='tab:blue')
         axes[nparams].plot(mean_resid_chain, '.', ms=ms, c='tab:orange')

From d8167317ba502f5040a7522ea4ce3fed49ba7295 Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Thu, 30 Oct 2025 18:39:32 -0400
Subject: [PATCH 07/11] Comment

---
 pygem/plot/graphics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygem/plot/graphics.py b/pygem/plot/graphics.py
index 99ed9d84..d809c5bb 100644
--- a/pygem/plot/graphics.py
+++ b/pygem/plot/graphics.py
@@ -231,7 +231,7 @@ def plot_mcmc_chain(
         axes[nparams].set_ylabel(r'$\dot{{b}}$', fontsize=fontsize)
         nparams += 1
 
-    # plot MAE for all other prediction keys
+    # plot along-chain mean residual for all other prediction keys
     for key in pred_primes.keys():
         if key == 'glacierwide_mb_mwea':
             continue

From a0bc1f4f06a702a275dc78b23bb8d78f4ed7bd04 Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Thu, 30 Oct 2025 18:48:05 -0400
Subject: [PATCH 08/11] Update 01Region -> O1Region in output regional glen_a
 dataset

---
 pygem/bin/run/run_inversion.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/pygem/bin/run/run_inversion.py b/pygem/bin/run/run_inversion.py
index ab8a029c..24ffaa2b 100644
--- a/pygem/bin/run/run_inversion.py
+++ b/pygem/bin/run/run_inversion.py
@@ -55,15 +55,15 @@ def export_regional_results(regions, outpath):
     # sort by the region number
     merged_df = merged_df.sort_values('rnum').drop(columns='rnum')
 
-    # if the file already exists, replace rows with same '01Region'
+    # if the file already exists, replace rows with same 'O1Region'
     if os.path.exists(outpath):
         existing_df = pd.read_csv(outpath)
-        # remove rows with the same '01Region' values as in the new merge
+        # remove rows with the same 'O1Region' values as in the new merge
         merged_df = pd.concat(
-            [existing_df[~existing_df['01Region'].isin(merged_df['01Region'])], merged_df], ignore_index=True
+            [existing_df[~existing_df['O1Region'].isin(merged_df['O1Region'])], merged_df], ignore_index=True
         )
         # re-sort
-        merged_df = merged_df.sort_values('01Region')
+        merged_df = merged_df.sort_values('O1Region')
 
     # export final merged csv
     merged_df.to_csv(outpath, index=False)
@@ -350,7 +350,7 @@ def run(
     # prepare ouptut dataset
     df = pd.Series(
         {
-            '01Region': reg,
+            'O1Region': reg,
             'count': len(glac_no),
             'inversion_glen_a': gdirs[0].get_diagnostics()['inversion_glen_a'],
             'inversion_fs': gdirs[0].get_diagnostics()['inversion_fs'],

From 660574a8484ac976562c35eedd245d97d0dd24dc Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Mon, 3 Nov 2025 09:51:12 -0500
Subject: [PATCH 09/11] Always compute and store glacierwide avg mass balance,
 but optionionally calibrate against

---
 pygem/bin/run/run_calibration.py | 30 ++++++++++++++----------------
 pygem/mcmc.py                    |  8 +++++++-
 2 files changed, 21 insertions(+), 17 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 05b170f2..8b780647 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -385,7 +385,6 @@ def mcmc_model_eval(
     glacier_rgi_table,
     fls,
     mbfxn=None,
-    calib_glacierwide_mb_mwea=True,
     calib_elev_change_1d=False,
     calib_snowlines_1d=False,
     calib_meltextent_1d=False,
@@ -403,22 +402,21 @@ def mcmc_model_eval(
         mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls)
         results['elev_change_1d'] = calc_elev_change_1d(gdir, mbmod, ds) if ds else float('-inf')
 
-    if calib_glacierwide_mb_mwea:
-        if mbfxn is not None:
-            # grab current values from modelprms for the emulator
-            mb_args = [modelprms['tbias'], modelprms['kp'], modelprms['ddfsnow']]
-            glacierwide_mb_mwea = mbfxn(*[mb_args])
+    if mbfxn is not None:
+        # grab current values from modelprms for the emulator
+        mb_args = [modelprms['tbias'], modelprms['kp'], modelprms['ddfsnow']]
+        glacierwide_mb_mwea = mbfxn(*[mb_args])
+    else:
+        if mbmod is None:
+            glacierwide_mb_mwea = mb_mwea_calc(gdir, modelprms, glacier_rgi_table, fls)
         else:
-            if mbmod is None:
-                glacierwide_mb_mwea = mb_mwea_calc(gdir, modelprms, glacier_rgi_table, fls)
-            else:
-                glacierwide_mb_mwea = (
-                    mbmod.glac_wide_massbaltotal[gdir.mbdata['t1_idx'] : gdir.mbdata['t2_idx'] + 1].sum()
-                    / mbmod.glac_wide_area_annual[0]
-                    / gdir.mbdata['nyears']
-                )
+            glacierwide_mb_mwea = (
+                mbmod.glac_wide_massbaltotal[gdir.mbdata['t1_idx'] : gdir.mbdata['t2_idx'] + 1].sum()
+                / mbmod.glac_wide_area_annual[0]
+                / gdir.mbdata['nyears']
+            )
 
-        results['glacierwide_mb_mwea'] = glacierwide_mb_mwea
+    results['glacierwide_mb_mwea'] = glacierwide_mb_mwea
 
     # (add future calibration options here)
     if calib_snowlines_1d:
@@ -2104,7 +2102,6 @@ def rho_constraints(**kwargs):
                     glacier_rgi_table,
                     fls,
                     mbfxn,
-                    args.option_calib_glacierwide_mb_mwea,
                     args.option_calib_elev_change_1d,
                 )
 
@@ -2116,6 +2113,7 @@ def rho_constraints(**kwargs):
                     priors,
                     fxn2eval=mcmc_model_eval,
                     fxnargs=fxnargs,
+                    calib_glacierwide_mb_mwea=args.option_calib_glacierwide_mb_mwea,
                     potential_fxns=[mb_max, must_melt, rho_constraints],
                     ela=gdir.ela.min() if hasattr(gdir, 'ela') else None,
                     bin_z=gdir.elev_change_1d['bin_centers'] if hasattr(gdir, 'elev_change_1d') else None,
diff --git a/pygem/mcmc.py b/pygem/mcmc.py
index f84ee5a4..f57df929 100644
--- a/pygem/mcmc.py
+++ b/pygem/mcmc.py
@@ -150,12 +150,15 @@ def log_uniform_density(x, **kwargs):
 
 # mass balance posterior class
 class mbPosterior:
-    def __init__(self, obs, priors, fxn2eval, fxnargs=None, potential_fxns=None, **kwargs):
+    def __init__(
+        self, obs, priors, fxn2eval, fxnargs=None, calib_glacierwide_mb_mwea=True, potential_fxns=None, **kwargs
+    ):
         # obs will be passed as a list, where each item is a tuple with the first element being the mean observation, and the second being the variance
         self.obs = obs
         self.priors = copy.deepcopy(priors)
         self.fxn2eval = fxn2eval
         self.fxnargs = fxnargs
+        self.calib_glacierwide_mb_mwea = calib_glacierwide_mb_mwea
         self.potential_functions = potential_fxns if potential_fxns is not None else []
         self.preds = None
         self.check_priors()
@@ -239,6 +242,9 @@ def log_likelihood(self, m):
                 # update values in preds dict
                 self.preds[k] = pred
 
+            if k == 'glacierwide_mb_mwea' and not self.calib_glacierwide_mb_mwea:
+                continue  # skip this model output if not calibrating glacierwide mass balance
+
             log_likehood += log_normal_density(
                 self.obs[k][0],  # observations
                 mu=pred,  # scaled predictions

From b886f40eeb932d0effa9cbe99f6370d984f9f97c Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Tue, 4 Nov 2025 09:56:18 -0500
Subject: [PATCH 10/11] Better commenting on calc_thick_change_1d() fxn

---
 pygem/bin/run/run_calibration.py | 55 ++++++++++++++++++++------------
 pygem/setup/config.py            |  1 -
 2 files changed, 35 insertions(+), 21 deletions(-)

diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py
index 74e379b4..e35fe76e 100755
--- a/pygem/bin/run/run_calibration.py
+++ b/pygem/bin/run/run_calibration.py
@@ -321,31 +321,44 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
     return mbmod, ds
 
 
-def calc_elev_change_1d(gdir, mbmod, ds):
-    ### get monthly ice thickness
+def calc_thick_change_1d(gdir, mbmod, ds):
+    """
+    calculate binned change in ice thickness assuming constant annual flux divergence.
+    sub-annual ice thickness is differenced at timesteps coincident with observations.
+    """
+    years_subannual = np.array([d.year for d in gdir.dates_table['date']])
+    yrs = np.unique(years_subannual)
+    nyrs = len(yrs)
     # grab components of interest
-    thickness_m = ds[0].thickness_m.values.T  # glacier thickness [m ice], (nbins, nyears)
+    bin_thick_annual = ds[0].thickness_m.values.T  # glacier thickness [m ice], (nbins, nyears)
 
     # set any < 0 thickness to nan
-    thickness_m[thickness_m <= 0] = np.nan
+    bin_thick_annual[bin_thick_annual <= 0] = np.nan
 
-    # climatic mass balance
-    dotb_monthly = mbmod.glac_bin_massbalclim  # climatic mass balance [m w.e.] per month
+    #  --- Step 1: convert mass balance from m w.e. to m ice  ---
+    bin_massbalclim = mbmod.glac_bin_massbalclim  # climatic mass balance [m w.e.] per step
     # convert to m ice
-    dotb_monthly = dotb_monthly * (pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice'])
+    bin_massbalclim_ice = bin_massbalclim * (
+        pygem_prms['constants']['density_water'] / pygem_prms['constants']['density_ice']
+    )
 
-    ### to get monthly thickness and mass we require monthly flux divergence ###
-    # we'll assume the flux divergence is constant througohut the year
-    # ie. take annual values and divide by 12 - use numpy repeat to repeat values across 12 months
-    flux_div_monthly_mmo = np.repeat(-ds[0].flux_divergence_myr.values.T[:, 1:] / 12, 12, axis=-1)
+    # --- Step 2: expand flux divergence to subannual steps ---
+    # assume flux divergence is constant throughout the year
+    # (divide annual by the number of steps in the binned climatic mass balance to get subannual flux divergence)
+    bin_flux_divergence_annual = -ds[0].flux_divergence_myr.values.T[:, 1:]
+    bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice)
+    for i, year in enumerate(yrs):
+        idx = np.where(years_subannual == year)[0]
+        bin_flux_divergence_subannual[:, idx] = bin_flux_divergence_annual[:, i][:, np.newaxis] / len(idx)
 
-    # get monthly binned change in thickness
-    delta_h_monthly = dotb_monthly - flux_div_monthly_mmo  # [m ice per month]
+    # --- Step 3: compute subannual thickness change ---
+    bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual
 
-    # get binned monthly thickness = running thickness change + initial thickness
-    running_delta_h_monthly = np.cumsum(delta_h_monthly, axis=-1)
-    h_monthly = running_delta_h_monthly + thickness_m[:, 0][:, np.newaxis]
+    # --- Step 4: calculate subannual thickness = running thickness change + initial thickness---
+    running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1)
+    bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, 0][:, np.newaxis]
 
+    # --- Step 5: rebin ---
     # get surface height at the specified reference year
     ref_surface_height = ds[0].bed_h.values + ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values
     # aggregate model bin thicknesses as desired
@@ -365,8 +378,8 @@ def calc_elev_change_1d(gdir, mbmod, ds):
     # interpolate over any empty bins
     bin_thick_subannual = np.column_stack([interp1d_fill_gaps(x.copy()) for x in bin_thick_subannual.T])
 
-    # difference each set of inds in diff_inds_map
-    elev_change_1d = np.column_stack(
+    # --- Step 5: compute binned thickness change ---
+    bin_thick_change = np.column_stack(
         [
             bin_thick_subannual[:, tup[1]] - bin_thick_subannual[:, tup[0]]
             if tup[0] is not None and tup[1] is not None
@@ -374,7 +387,8 @@ def calc_elev_change_1d(gdir, mbmod, ds):
             for tup in gdir.elev_change_1d['model2obs_inds_map']
         ]
     )
-    return elev_change_1d
+
+    return bin_thick_change
 
 
 def mcmc_model_eval(
@@ -398,7 +412,8 @@ def mcmc_model_eval(
 
     if calib_elev_change_1d:
         mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls)
-        results['elev_change_1d'] = calc_elev_change_1d(gdir, mbmod, ds) if ds else float('-inf')
+        # note, the binned thickness change is scaled by modeled density in mcmc.mbPosterior.log_likelihood() to calculate modeled surface elevation change
+        results['elev_change_1d'] = calc_thick_change_1d(gdir, mbmod, ds) if ds else float('-inf')
 
     if mbfxn is not None:
         # grab current values from modelprms for the emulator
diff --git a/pygem/setup/config.py b/pygem/setup/config.py
index e5ddc50e..f0def64a 100644
--- a/pygem/setup/config.py
+++ b/pygem/setup/config.py
@@ -163,7 +163,6 @@ def flatten_dict(d, parent_key=''):
             # Skip patterns
             if any(fnmatch.fnmatch(key, pat) for pat in skip_patterns):
                 continue
-            print(key, value)
 
             path = os.path.join(root, value.strip(os.sep))
 

From 93a0647068e727872015a9c6baeb2dbc6fa344bd Mon Sep 17 00:00:00 2001
From: "brandon s. tober" <tobers.brandon@gmail.com>
Date: Tue, 4 Nov 2025 10:18:07 -0500
Subject: [PATCH 11/11] Re-order

---
 pygem/mcmc.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pygem/mcmc.py b/pygem/mcmc.py
index f57df929..b10ad61c 100644
--- a/pygem/mcmc.py
+++ b/pygem/mcmc.py
@@ -230,6 +230,9 @@ def log_likelihood(self, m):
                 # Invalid model output -> assign -inf likelihood
                 return torch.tensor([-float('inf')])
 
+            if k == 'glacierwide_mb_mwea' and not self.calib_glacierwide_mb_mwea:
+                continue  # skip this model output if not calibrating glacierwide mass balance
+
             # if key is `elev_change_1d` scale by density to predict binned surface elevation change
             if k == 'elev_change_1d':
                 # Create density field, separate values for ablation/accumulation zones
@@ -242,9 +245,6 @@ def log_likelihood(self, m):
                 # update values in preds dict
                 self.preds[k] = pred
 
-            if k == 'glacierwide_mb_mwea' and not self.calib_glacierwide_mb_mwea:
-                continue  # skip this model output if not calibrating glacierwide mass balance
-
             log_likehood += log_normal_density(
                 self.obs[k][0],  # observations
                 mu=pred,  # scaled predictions