From fa2ef893f7c23faf1a6eeca5a6d4d51c215f73da Mon Sep 17 00:00:00 2001 From: "brandon s. tober" Date: Fri, 5 Dec 2025 23:18:32 -0500 Subject: [PATCH 1/2] Allow for inversion in mcmc elev_change_1d calib --- pygem/bin/run/run_calibration.py | 113 ++++++++++++++++++------------- pygem/mcmc.py | 20 +++--- 2 files changed, 77 insertions(+), 56 deletions(-) diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py index 49f36dc6..21120692 100755 --- a/pygem/bin/run/run_calibration.py +++ b/pygem/bin/run/run_calibration.py @@ -18,6 +18,7 @@ import time import warnings from datetime import timedelta +from functools import partial import gpytorch import matplotlib.pyplot as plt @@ -38,8 +39,9 @@ # read the config pygem_prms = config_manager.read_config() -from oggm import cfg, tasks, utils +from oggm import cfg, tasks, utils, workflow from oggm.core import flowline +from oggm.core.massbalance import apparent_mb_from_any_mb import pygem.oggm_compat as oggm_compat import pygem.pygem_modelsetup as modelsetup @@ -48,7 +50,7 @@ from pygem.plot import graphics from pygem.utils._funcs import interp1d_fill_gaps from pygem.utils.stats import mcmc_stats - +from pygem.shop import debris # %% FUNCTIONS def getparser(): @@ -256,6 +258,50 @@ def mb_mwea_calc( return mb_mwea +def run_inversion(gdir, modelprms, glacier_rgi_table, inversion_filter=pygem_prms['mb']['include_debris']): + """run inversion steps with a given set of model parameters""" + + # Perform inversion based on PyGEM MB using reference directory + mbmod_inv = PyGEMMassBalance( + gdir, + modelprms, + glacier_rgi_table, + fls=gdir.read_pickle('inversion_flowlines'), + option_areaconstant=True, + inversion_filter=inversion_filter, + ) + # Arbitrariliy shift the MB profile up (or down) until mass balance is zero (equilibrium for inversion) + apparent_mb_from_any_mb(gdir, mb_model=mbmod_inv) + + # Non-tidewater glaciers + if not gdir.is_tidewater or not pygem_prms['setup']['include_frontalablation']: + tasks.prepare_for_inversion(gdir) + tasks.mass_conservation_inversion( + gdir, + glen_a=gdir.get_diagnostics()['inversion_glen_a'], + fs=gdir.get_diagnostics()['inversion_fs'], + ) + + # Tidewater glaciers + else: + cfg.PARAMS['use_kcalving_for_inversion'] = True + cfg.PARAMS['use_kcalving_for_run'] = True + tasks.find_inversion_calving_from_any_mb( + gdir, + mb_model=mbmod_inv, + glen_a=gdir.get_diagnostics()['inversion_glen_a'], + fs=gdir.get_diagnostics()['inversion_fs'], + ) + + tasks.init_present_time_glacier(gdir) # adds bins below + if pygem_prms['mb']['include_debris']: + debris.debris_binned( + gdir, fl_str='model_flowlines' + ) # add debris enhancement factors to flowlines + + return gdir.read_pickle('model_flowlines') + + def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls, debug=False): """run the dynamical evolution model with a given set of model parameters""" @@ -265,11 +311,6 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls, debug=False): # 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) # glacier dynamics model if gdir.is_tidewater and pygem_prms['setup']['include_frontalablation']: ev_model = flowline.FluxBasedModel( @@ -279,7 +320,7 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls, debug=False): glen_a=gdir.get_diagnostics()['inversion_glen_a'], fs=gdir.get_diagnostics()['inversion_fs'], is_tidewater=gdir.is_tidewater, - water_level=water_level, + water_level=gdir.get_diagnostics().get('calving_water_level', None), do_kcalving=pygem_prms['setup']['include_frontalablation'], ) else: @@ -289,8 +330,6 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls, debug=False): mb_model=mbmod, glen_a=gdir.get_diagnostics()['inversion_glen_a'], fs=gdir.get_diagnostics()['inversion_fs'], - is_tidewater=gdir.is_tidewater, - water_level=water_level, ) if debug: fig, ax = plt.subplots(1) @@ -298,29 +337,6 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls, debug=False): try: # run glacier dynamics model forward diag, ds = ev_model.run_until_and_store(y1 + 1, fl_diag_path=True) - with np.errstate(invalid='ignore'): - # record frontal ablation for tidewater glaciers and update total mass balance - if gdir.is_tidewater and pygem_prms['setup']['include_frontalablation']: - # glacier-wide frontal ablation (m3 w.e.) - # - note: diag.calving_m3 is cumulative calving, convert to annual calving - calving_m3we_annual = ( - (diag.calving_m3.values[1:] - diag.calving_m3.values[0:-1]) - * pygem_prms['constants']['density_ice'] - / pygem_prms['constants']['density_water'] - ) - # record each year's frontal ablation in m3 w.e. - if pygem_prms['time']['timestep'] == 'monthly': - for n, year in enumerate(np.arange(y0, y1 + 1)): - tstart, tstop = ev_model.mb_model.get_step_inds(year) - ev_model.mb_model.glac_wide_frontalablation[tstop] = calving_m3we_annual[n] - else: - raise ValueError('Need to add functionality for daily timestep') - - # 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 - ) - ev_model.mb_model.ensure_mass_conservation(diag) if debug: graphics.plot_modeloutput_section(ev_model, ax=ax, srfls='--', lnlabel=f'Glacier year {y1}') plt.show() @@ -410,14 +426,15 @@ def calc_thick_change_1d(gdir, fls, mbmod, ds): def mcmc_model_eval( - gdir, modelprms, + gdir, glacier_rgi_table, fls, mbfxn=None, calib_elev_change_1d=False, calib_snowlines_1d=False, calib_meltextent_1d=False, + do_inversion=False, debug=False, ): """ @@ -428,6 +445,9 @@ def mcmc_model_eval( results = {} mbmod = None + if do_inversion: + fls = run_inversion(gdir, modelprms, glacier_rgi_table) + if calib_elev_change_1d: mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls) # note, the binned thickness change is scaled by modeled density in mcmc.mbPosterior.log_likelihood() to calculate modeled surface elevation change @@ -910,9 +930,10 @@ def run(list_packed_vars): calving_df['O1Region'] = [int(x.split('-')[1].split('.')[0]) for x in calving_df.RGIId.values] calving_df_reg = calving_df.loc[calving_df['O1Region'] == int(gdir.rgi_id[6:8]), :] calving_k = np.median(calving_df_reg.calving_k) - # set calving_k in config - cfg.PARAMS['use_kcalving_for_run'] = True + # set calving_k cfg.PARAMS['calving_k'] = calving_k + cfg.PARAMS['use_kcalving_for_inversion'] = True + cfg.PARAMS['use_kcalving_for_run'] = True # note, some tidewater glaciers require a timestep << OGGM default of 60 seconds cfg.PARAMS['cfl_number'] = pygem_prms['sim']['oggm_dynamics']['cfl_number'] cfg.PARAMS['cfl_min_dt'] = pygem_prms['sim']['oggm_dynamics']['cfl_min_dt'] @@ -2092,14 +2113,14 @@ def rho_constraints(**kwargs): mbfxn = None # define args to pass to fxn2eval in mcmc sampler - fxnargs = ( - gdir, - modelprms, - glacier_rgi_table, - fls, - mbfxn, - args.option_calib_elev_change_1d, - ) + fxn2eval = partial(mcmc_model_eval, + gdir=gdir, + glacier_rgi_table=glacier_rgi_table, + fls=fls, + mbfxn=mbfxn, + calib_elev_change_1d=args.option_calib_elev_change_1d, + do_inversion=args.option_calib_elev_change_1d and not args.spinup + ) # 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) @@ -2107,8 +2128,8 @@ def rho_constraints(**kwargs): mb = mcmc.mbPosterior( obs, priors, - fxn2eval=mcmc_model_eval, - fxnargs=fxnargs, + fxn2eval=fxn2eval, + modelprms_dict=modelprms, 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, diff --git a/pygem/mcmc.py b/pygem/mcmc.py index b10ad61c..d2c30a00 100644 --- a/pygem/mcmc.py +++ b/pygem/mcmc.py @@ -151,13 +151,13 @@ def log_uniform_density(x, **kwargs): # mass balance posterior class class mbPosterior: def __init__( - self, obs, priors, fxn2eval, fxnargs=None, calib_glacierwide_mb_mwea=True, potential_fxns=None, **kwargs + self, obs, priors, fxn2eval, modelprms_dict, 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.modelprms_dict = modelprms_dict self.calib_glacierwide_mb_mwea = calib_glacierwide_mb_mwea self.potential_functions = potential_fxns if potential_fxns is not None else [] self.preds = None @@ -200,13 +200,13 @@ def check_priors(self): # update modelprms for evaluation def update_modelprms(self, m): for i, k in enumerate(['tbias', 'kp', 'ddfsnow']): - self.fxnargs[1][k] = float(m[i]) - self.fxnargs[1]['ddfice'] = self.fxnargs[1]['ddfsnow'] / pygem_prms['sim']['params']['ddfsnow_iceratio'] + self.modelprms_dict[k] = float(m[i]) + self.modelprms_dict['ddfice'] = self.modelprms_dict['ddfsnow'] / pygem_prms['sim']['params']['ddfsnow_iceratio'] # get model predictions def get_model_pred(self, m): self.update_modelprms(m) # update modelprms with current step - self.preds = self.fxn2eval(*self.fxnargs) + self.preds = self.fxn2eval(self.modelprms_dict) # convert all values to torch tensors self.preds = {k: torch.tensor(v, dtype=torch.float) for k, v in self.preds.items()} @@ -256,10 +256,10 @@ def log_likelihood(self, m): # compute the log-potential, summing over all declared potential functions. def log_potential(self, m): # --- Base arguments --- - # kp, tbias, ddfsnow, massbal + # tbias, kp, ddfsnow, massbal kwargs = { - 'kp': m[0], - 'tbias': m[1], + 'tbias': m[0], + 'kp': m[1], 'ddfsnow': m[2], 'massbal': self.preds['glacierwide_mb_mwea'], } @@ -267,8 +267,8 @@ def log_potential(self, m): # --- Optional arguments(if len(m) > 3) --- # rhoabl, rhoacc if len(m) > 3: - kwargs['rhoabl'] = m[-2] - kwargs['rhoacc'] = m[-1] + kwargs['rhoabl'] = m[3] + kwargs['rhoacc'] = m[4] # --- Evaluate all potential functions --- return sum(pf(**kwargs) for pf in self.potential_functions) From a9342072692173ee23cea57e3f7b223cd2166a31 Mon Sep 17 00:00:00 2001 From: "brandon s. tober" Date: Fri, 5 Dec 2025 23:21:47 -0500 Subject: [PATCH 2/2] Ruffed --- pygem/bin/run/run_calibration.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py index 21120692..3884f4e4 100755 --- a/pygem/bin/run/run_calibration.py +++ b/pygem/bin/run/run_calibration.py @@ -39,7 +39,7 @@ # read the config pygem_prms = config_manager.read_config() -from oggm import cfg, tasks, utils, workflow +from oggm import cfg, tasks from oggm.core import flowline from oggm.core.massbalance import apparent_mb_from_any_mb @@ -48,9 +48,10 @@ from pygem import class_climate, mcmc from pygem.massbalance import PyGEMMassBalance from pygem.plot import graphics +from pygem.shop import debris from pygem.utils._funcs import interp1d_fill_gaps from pygem.utils.stats import mcmc_stats -from pygem.shop import debris + # %% FUNCTIONS def getparser(): @@ -295,13 +296,11 @@ def run_inversion(gdir, modelprms, glacier_rgi_table, inversion_filter=pygem_prm tasks.init_present_time_glacier(gdir) # adds bins below if pygem_prms['mb']['include_debris']: - debris.debris_binned( - gdir, fl_str='model_flowlines' - ) # add debris enhancement factors to flowlines - + debris.debris_binned(gdir, fl_str='model_flowlines') # add debris enhancement factors to flowlines + return gdir.read_pickle('model_flowlines') - + def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls, debug=False): """run the dynamical evolution model with a given set of model parameters""" @@ -2113,14 +2112,15 @@ def rho_constraints(**kwargs): mbfxn = None # define args to pass to fxn2eval in mcmc sampler - fxn2eval = partial(mcmc_model_eval, - gdir=gdir, - glacier_rgi_table=glacier_rgi_table, - fls=fls, - mbfxn=mbfxn, - calib_elev_change_1d=args.option_calib_elev_change_1d, - do_inversion=args.option_calib_elev_change_1d and not args.spinup - ) + fxn2eval = partial( + mcmc_model_eval, + gdir=gdir, + glacier_rgi_table=glacier_rgi_table, + fls=fls, + mbfxn=mbfxn, + calib_elev_change_1d=args.option_calib_elev_change_1d, + do_inversion=args.option_calib_elev_change_1d and not args.spinup, + ) # 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)