Skip to content

Commit

Permalink
New sandbox experiments (OGGM#38)
Browse files Browse the repository at this point in the history 
* add three different glacier states per sandbox geometry

* typo

* fix tests

* adapted dashboard for different glacier states

* added initial_flux to data_logger

* added tests for glacier_state mb_models

* small experiment adaptions
  • Loading branch information
@pat-schmitt
pat-schmitt authored Sep 1, 2023
1 parent 672e7bb commit 12a025f
Show file tree
Hide file tree
Showing 11 changed files with 930 additions and 503 deletions.
53 changes: 34 additions & 19 deletions agile1d/core/cost_function.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,6 @@ def get_indices_for_unknown_parameters(data_logger):


def define_scaling_terms(data_logger):

# define scaling terms for observation stuff
if 'scale' in data_logger.obs_scaling_parameters.keys():
observations = data_logger.observations
Expand Down Expand Up @@ -260,10 +259,7 @@ def cost_fct(unknown_parameters, data_logger):

mb_models = initialise_mb_models(unknown_parameters_descaled, data_logger)

if 'smoothed_flux' in data_logger.regularisation_terms.keys():
initial_flux = dynamic_model(flowline).flux_stag[0]
else:
initial_flux = None
initial_flux = dynamic_model(flowline).flux_stag[0]

if data_logger.spinup_type == 'height_shift_spinup':
try:
Expand Down Expand Up @@ -318,7 +314,7 @@ def cost_fct(unknown_parameters, data_logger):
initial_flux, # for regularisation term 'smoothed_flux'
unknown_parameters, # for distance from fg regularisation
data_logger # for reg_parameters and observations
)
)

# sum up cost function terms using cost lambda
cost_lambda = torch.tensor(data_logger.cost_lambda,
Expand All @@ -344,6 +340,9 @@ def cost_fct(unknown_parameters, data_logger):
data_logger.save_data_in_datalogger('sfc_h_start', sfc_h_start)
data_logger.save_data_in_datalogger('section_start', section_start)
data_logger.save_data_in_datalogger('flowlines', detach_flowline(final_fl))
# exclude last grid point because initial flux is defined on staggered grid
data_logger.save_data_in_datalogger('initial_flux',
initial_flux.detach().clone()[:-1])
data_logger.save_data_in_datalogger('costs', cost)
data_logger.save_data_in_datalogger('grads', grad)
data_logger.save_data_in_datalogger('c_terms', c_terms)
Expand Down Expand Up @@ -625,20 +624,36 @@ def initialise_mb_models(unknown_parameters,
# -1 because period defined as [y0 - halfsize, y0 + halfsize + 1]
y0 = (y_start + y_end - 1) / 2
halfsize = (y_end - y_start - 1) / 2
mb_models[mb_mdl_set] = {'mb_model':
ConstantMassBalanceTorch(
gdir, y0=y0, halfsize=halfsize,
torch_type=torch_type, device=device),
'years':
mb_models_settings[mb_mdl_set]['years']}
mb_models[mb_mdl_set] = {
'mb_model': ConstantMassBalanceTorch(
gdir, y0=y0, halfsize=halfsize,
torch_type=torch_type, device=device),
'years':
mb_models_settings[mb_mdl_set]['years']}
elif mb_models_settings[mb_mdl_set]['type'] == 'TIModel':
mb_models[mb_mdl_set] = {'mb_model':
MBModelTorchWrapper(
gdir=gdir,
mb_model=MonthlyTIModel(gdir)),
'years':
mb_models_settings[mb_mdl_set]['years']
}
mb_model_args = mb_models_settings[mb_mdl_set]['model_args']
mb_models[mb_mdl_set] = {
'mb_model': MBModelTorchWrapper(
gdir=gdir,
mb_model=MonthlyTIModel(
gdir,
**mb_model_args)
),
'years':
mb_models_settings[mb_mdl_set]['years']
}
elif mb_models_settings[mb_mdl_set]['type'] == 'ConstantModel':
mb_model_args = mb_models_settings[mb_mdl_set]['model_args']
mb_models[mb_mdl_set] = {
'mb_model': MBModelTorchWrapper(
gdir=gdir,
mb_model=ConstantMassBalance(
gdir,
**mb_model_args)
),
'years':
mb_models_settings[mb_mdl_set]['years']
}
else:
raise NotImplementedError("The MassBalance type "
f"{mb_models_settings[mb_mdl_set]['type']} "
Expand Down
5 changes: 4 additions & 1 deletion agile1d/core/data_logging.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,7 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename,
self.time_needed = None
self.grads = None
self.flowlines = None
self.initial_flux = None
self.sfc_h_start = None
self.section_start = None
self.observations_mdl = None
Expand All @@ -157,7 +158,7 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename,
self.first_guess = None

self.filename = gdir.name + '_' + \
inversion_input['experiment_description']
inversion_input['experiment_description']

# create info Text for callback_fct TODO: think about showing the evolution of the c_terms
self.info_text = '''
Expand Down Expand Up @@ -248,6 +249,7 @@ def filter_data_from_optimization(self):
self.reg_terms = self.reg_terms[index]
self.time_needed = self.squeeze_generic(self.time_needed[index])
self.flowlines = self.squeeze_generic(self.flowlines[index])
self.initial_flux = self.squeeze_generic(self.initial_flux[index])
self.sfc_h_start = self.squeeze_generic(self.sfc_h_start[index])
self.section_start = self.squeeze_generic(self.section_start[index])
self.observations_mdl = self.squeeze_generic(self.observations_mdl[index])
Expand All @@ -269,6 +271,7 @@ def create_and_save_dataset(self):
ds['costs'] = (['iteration'], self.costs)
ds['grads'] = (['iteration', 'nr_unknown_parameters'], self.grads)
ds['flowlines'] = (['iteration'], self.flowlines)
ds['initial_flux'] = (['iteration', 'x'], self.initial_flux)
ds['sfc_h_start'] = (['iteration', 'x'], self.sfc_h_start)
ds['section_start'] = (['iteration', 'x'], self.section_start)
ds['observations_mdl'] = (['iteration'], self.observations_mdl)
Expand Down
14 changes: 8 additions & 6 deletions agile1d/core/inversion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,9 +60,11 @@ def add_setting():
"MassBalanceModel must start at least " \
"one year before first given observation year!" \
"Default: {'MB': {'type': 'TIModel'," \
" 'years': np.array([1980, 2020])}}"
" 'years': np.array([1980, 2020])," \
" 'model_args': {}}}"
_default = {'MB': {'type': 'TIModel',
'years': np.array([1980, 2020])}}
'years': np.array([1980, 2020]),
'model_args': {}}}

add_setting()

Expand Down Expand Up @@ -103,10 +105,10 @@ def add_setting():

_key = "bed_h_bounds"
_doc = "Define how large the boundaries for the bed_h are, in relation of " \
"first guess thickness. (e.g. (0.2, 1.4) means the bed height can " \
"first guess thickness. (e.g. (0.4, 1.6) means the bed height can " \
"be between 1.4*fg_thick and 0.2*fg_thick). " \
"Default: (0.2, 1.4)"
_default = (0.2, 1.4)
"Default: (0.4, 1.6)"
_default = (0.4, 1.6)
add_setting()

_key = "max_deviation_surface_h"
Expand Down Expand Up @@ -245,7 +247,7 @@ def add_setting():
"glacier at the inital state." \
"e.g. {'section':" \
" {'extra_grid_points': 10," \
" 'limits': (0.75, 1.25)," \
" 'limits': (0.6, 1.4)," \
" 'fg_years': 1" \
" }" \
" }" \
Expand Down
4 changes: 3 additions & 1 deletion agile1d/core/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,9 @@ def inversion_settings(self):

def write_inversion_settings(self,
control_vars=['bed_h'],
mb_models_settings={'MB1': {'type': 'constant', 'years': np.array([1950, 2016])}},
mb_models_settings={'MB1': {
'type': 'constant',
'years': np.array([1950, 2016])}},
min_w0_m=10.,
observations=None, # {'Area', {'2010': np.array([23])}}
reg_parameters=None, # [0, 0.1, 10]
Expand Down
58 changes: 38 additions & 20 deletions agile1d/sandbox/calculate_statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,8 @@ def add_0d_stats(x, y):
'abs_diff': float(np.abs(x - y))}


def calculate_result_statistics(gdir, data_logger, print_statistic=False):
def calculate_result_statistics(gdir, glacier_state, data_logger,
print_statistic=False):
"""calculate some statistics of the result for analysis"""

# open the dataset of the run to add our calculated statistics
Expand Down Expand Up @@ -69,7 +70,8 @@ def calculate_result_statistics(gdir, data_logger, print_statistic=False):
ds.unknown_parameters[-1][control_indices[control_var]].values
controls_true = {}
fls_true = gdir.read_pickle('model_flowlines',
filesuffix='_agile_true_init')[0]
filesuffix='_agile_true_init_'
f'{glacier_state}')[0]
for control_var in controls_mdl:
if control_var in ['bed_h', 'area_bed_h']:
controls_true[control_var] = \
Expand All @@ -81,7 +83,8 @@ def calculate_result_statistics(gdir, data_logger, print_statistic=False):
controls_true[control_var] = controls_mdl[control_var]
elif control_var in ['section']:
fls_1980_true = gdir.read_pickle('model_flowlines',
filesuffix='_creation_spinup')[0]
filesuffix='_creation_spinup_'
f'{glacier_state}')[0]
controls_true[control_var] = \
fls_1980_true.section[:len(controls_mdl['section'])]
else:
Expand Down Expand Up @@ -117,7 +120,8 @@ def calculate_result_statistics(gdir, data_logger, print_statistic=False):
sfc_h_start = copy.deepcopy(data_logger.sfc_h_start[-1])
fls_start_mdl.surface_h = sfc_h_start
fls_start_true = gdir.read_pickle('model_flowlines',
filesuffix='_creation_spinup')[0]
filesuffix='_creation_spinup_'
f'{glacier_state}')[0]

past_state_stats = {}
for var in ['thick', 'area_m2', 'volume_m3']:
Expand Down Expand Up @@ -150,7 +154,8 @@ def get_volume(fl):
with xr.open_dataset(fp) as ds_diag:
past_evol_mdl = ds_diag.load()
fp = gdir.get_filepath('model_diagnostics',
filesuffix='_agile_true_total_run')
filesuffix='_agile_true_total_run_'
f'{glacier_state}')
with xr.open_dataset(fp) as ds_diag:
past_evol_true = ds_diag.load()

Expand All @@ -165,7 +170,8 @@ def get_volume(fl):
# how well do we match today's glacier state ------------------------------
fls_end_mdl = copy.deepcopy(data_logger.flowlines[-1])
fls_end_true = gdir.read_pickle('model_flowlines',
filesuffix='_agile_true_end')[0]
filesuffix='_agile_true_end_'
f'{glacier_state}')[0]

today_state_stats = {}
for var in ['thick', 'area_m2', 'volume_m3']:
Expand Down Expand Up @@ -196,7 +202,8 @@ def get_volume(fl):
with xr.open_dataset(fp) as ds_diag:
future_evol_mdl = ds_diag.load()
fp = gdir.get_filepath('model_diagnostics',
filesuffix='_agile_true_future')
filesuffix='_agile_true_future_'
f'{glacier_state}')
with xr.open_dataset(fp) as ds_diag:
future_evol_true = ds_diag.load()

Expand All @@ -218,7 +225,8 @@ def get_volume(fl):
# Here print the statistics in comparision to the default run

# open default statistics
fp_default = os.path.join(gdir.dir, 'default_oggm_statistics.pkl')
fp_default = os.path.join(gdir.dir, f'default_oggm_statistics_'
f'{glacier_state}.pkl')
with open(fp_default, 'rb') as handle:
default_stats = pickle.load(handle)

Expand Down Expand Up @@ -271,28 +279,32 @@ def get_volume(fl):
raise NotImplementedError(f'{stat_clean}')


def calculate_default_oggm_statistics(gdir):
def calculate_default_oggm_statistics(gdir, glacier_state):

default_oggm_statistics = {}

for reali in ['dynamic', 'static']:
# open the run files
with xr.open_dataset(
gdir.get_filepath('model_diagnostics',
filesuffix=f'_oggm_{reali}_past')) as ds:
filesuffix=f'_oggm_{reali}_past_'
f'{glacier_state}')) as ds:
diag_past = ds.load()
f = gdir.get_filepath('fl_diagnostics',
filesuffix=f'_oggm_{reali}_past')
filesuffix=f'_oggm_{reali}_past_'
f'{glacier_state}')
with xr.open_dataset(f, group=f'fl_0') as ds:
fl_diag_past = ds.load()
with xr.open_dataset(
gdir.get_filepath('model_diagnostics',
filesuffix=f'_oggm_{reali}_future')) as ds:
filesuffix=f'_oggm_{reali}_future_'
f'{glacier_state}')) as ds:
diag_future = ds.load()

# how well do we match the observations -------------------------------
obs_given = gdir.read_pickle('inversion_input',
filesuffix='_agile_measurements')
filesuffix='_agile_measurements_'
f'{glacier_state}')
obs_stats = {}
for obs_key in obs_given.keys():
obs_stats[obs_key] = {}
Expand Down Expand Up @@ -374,9 +386,11 @@ def calculate_default_oggm_statistics(gdir):
controls_mdl = {}
controls_true = {}
fls_mdl = gdir.read_pickle('model_flowlines',
filesuffix='_oggm_first_guess')[0]
filesuffix='_oggm_first_guess_'
f'{glacier_state}')[0]
fls_true = gdir.read_pickle('model_flowlines',
filesuffix='_agile_true_init')[0]
filesuffix='_agile_true_init_'
f'{glacier_state}')[0]

for control_var in all_control_vars:
if control_var in ['bed_h', 'area_bed_h']:
Expand Down Expand Up @@ -405,7 +419,8 @@ def calculate_default_oggm_statistics(gdir):
# how well do we match the past glacier state -----------------------------
fls_start_mdl = fl_diag_past.sel(time=fl_diag_past.time[0])
fls_start_true = gdir.read_pickle('model_flowlines',
filesuffix='_creation_spinup')[0]
filesuffix='_creation_spinup_'
f'{glacier_state}')[0]

past_state_stats = {}
for var in ['thick', 'area_m2', 'volume_m3']:
Expand Down Expand Up @@ -435,7 +450,8 @@ def get_volume(fl):
# how well do we match the past glacier evolution ---------------------
past_evol_mdl = diag_past
fp = gdir.get_filepath('model_diagnostics',
filesuffix='_agile_true_total_run')
filesuffix='_agile_true_total_run_'
f'{glacier_state}')
with xr.open_dataset(fp) as ds_diag:
past_evol_true = ds_diag.load()

Expand All @@ -450,7 +466,8 @@ def get_volume(fl):
# how well do we match today's glacier state --------------------------
fls_end_mdl = fl_diag_past.sel(time=fl_diag_past.time[-1])
fls_end_true = gdir.read_pickle('model_flowlines',
filesuffix='_agile_true_end')[0]
filesuffix='_agile_true_end_'
f'{glacier_state}')[0]

today_state_stats = {}
for var in ['thick', 'area_m2', 'volume_m3']:
Expand Down Expand Up @@ -480,7 +497,8 @@ def get_volume(fl):
# how well do we match the future glacier evolution -------------------
future_evol_mdl = diag_future
fp = gdir.get_filepath('model_diagnostics',
filesuffix='_agile_true_future')
filesuffix='_agile_true_future_'
f'{glacier_state}')
with xr.open_dataset(fp) as ds_diag:
future_evol_true = ds_diag.load()

Expand All @@ -494,7 +512,7 @@ def get_volume(fl):

# save final default statistics as pickle
out = os.path.join(gdir.dir,
'default_oggm_statistics.pkl')
f'default_oggm_statistics_{glacier_state}.pkl')
with open(out, 'wb') as handle:
pickle.dump(default_oggm_statistics, handle,
protocol=pickle.HIGHEST_PROTOCOL)
Loading

0 comments on commit 12a025f

Please sign in to comment.