New sandbox experiments (#38)

* 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

Author: pat-schmitt

agile1d/core/cost_function.py

@@ -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
@@ -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:
@@ -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,
@@ -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)
@@ -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']} "

agile1d/core/data_logging.py

@@ -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
@@ -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 = '''
@@ -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])
@@ -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)

agile1d/core/inversion.py

@@ -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()
 
@@ -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"
@@ -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" \
            "          }" \
            "      }" \

agile1d/core/utils.py

@@ -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]

agile1d/sandbox/calculate_statistics.py

@@ -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
@@ -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] = \
@@ -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:
@@ -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']:
@@ -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()
 
@@ -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']:
@@ -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()
 
@@ -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)
 
@@ -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] = {}
@@ -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']:
@@ -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']:
@@ -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()
 
@@ -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']:
@@ -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()
 
@@ -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)