Added sampling without replacement to RandomMassBalance (#502)

* Added RandomMassBalance sampling without replacement

* added tests and changed parameter name

* forgot important uniqueness test

oggm/core/flowline.py

@@ -1576,6 +1576,7 @@ def run_random_climate(gdir, nyears=1000, y0=None, halfsize=15,
                        climate_input_filesuffix='',
                        output_filesuffix='', init_model_fls=None,
                        zero_initial_glacier=False,
+                       unique_samples=False,
                        **kwargs):
     """Runs the random mass-balance model for a given number of years.
 
@@ -1614,14 +1615,20 @@ def run_random_climate(gdir, nyears=1000, y0=None, halfsize=15,
         present_time_glacier file from the glacier directory)
     zero_initial_glacier : bool
         if true, the ice thickness is set to zero before the simulation
+    unique_samples: bool
+        if true, chosen random mass-balance years will only be available once
+        per random climate period-length
+        if false, every model year will be chosen from the random climate
+        period with the same probability
     kwargs : dict
         kwargs to pass to the FluxBasedModel instance
     """
 
     mb = mbmods.RandomMassBalance(gdir, y0=y0, halfsize=halfsize,
                                   bias=bias, seed=seed,
                                   filename=climate_filename,
-                                  input_filesuffix=climate_input_filesuffix)
+                                  input_filesuffix=climate_input_filesuffix,
+                                  unique_samples=unique_samples)
     if temperature_bias is not None:
         mb.temp_bias = temperature_bias
 

oggm/core/massbalance.py

@@ -557,7 +557,8 @@ class RandomMassBalance(MassBalanceModel):
 
     def __init__(self, gdir, mu_star=None, bias=None, prcp_fac=None,
                  y0=None, halfsize=15, seed=None, filename='climate_monthly',
-                 input_filesuffix='', all_years=False):
+                 input_filesuffix='', all_years=False,
+                 unique_samples=False):
         """Initialize.
 
         Parameters
@@ -590,6 +591,11 @@ def __init__(self, gdir, mu_star=None, bias=None, prcp_fac=None,
         all_years : bool
             if True, overwrites ``y0`` and ``halfsize`` to use all available
             years.
+        unique_samples: bool
+            if true, chosen random mass-balance years will only be available
+            once per random climate period-length
+            if false, every model year will be chosen from the random climate
+            period with the same probability
         """
 
         super(RandomMassBalance, self).__init__()
@@ -613,6 +619,11 @@ def __init__(self, gdir, mu_star=None, bias=None, prcp_fac=None,
         self.rng = np.random.RandomState(seed)
         self._state_yr = dict()
 
+        # Sampling without replacement
+        self.unique_samples = unique_samples
+        if self.unique_samples:
+            self.sampling_years = self.years
+
     @property
     def temp_bias(self):
         """Temperature bias to add to the original series."""
@@ -653,7 +664,22 @@ def get_state_yr(self, year=None):
         """For a given year, get the random year associated to it."""
         year = int(year)
         if year not in self._state_yr:
-            self._state_yr[year] = self.rng.randint(*self.yr_range)
+            if self.unique_samples:
+                # --- Sampling without replacement ---
+                if self.sampling_years.size == 0:
+                    # refill sample pool when all years were picked once
+                    self.sampling_years = self.years
+                # choose one year which was not used in the current period
+                _sample = self.rng.choice(self.sampling_years)
+                # write chosen year to dictionary
+                self._state_yr[year] = _sample
+                # update sample pool: remove the chosen year from it
+                self.sampling_years = np.delete(
+                    self.sampling_years,
+                    np.where(self.sampling_years == _sample))
+            else:
+                # --- Sampling with replacement ---
+                self._state_yr[year] = self.rng.randint(*self.yr_range)
         return self._state_yr[year]
 
     def get_monthly_mb(self, heights, year=None):

oggm/tests/test_models.py

@@ -500,7 +500,7 @@ def test_random_mb(self):
                                                seed=10)
         mb_ts = []
         mb_ts2 = []
-        yrs = np.arange(1973, 2003, 1)
+        yrs = np.arange(1973, 2004, 1)
         for yr in yrs:
             mb_ts.append(np.average(mb_ref.get_annual_mb(h, yr) * SEC_IN_YEAR, weights=w))
             mb_ts2.append(np.average(mb_mod.get_annual_mb(h, yr) * SEC_IN_YEAR, weights=w))
@@ -519,6 +519,89 @@ def test_random_mb(self):
         ref_mb = ref_mb / 31
         self.assertTrue(utils.rmsd(ref_mb, my_mb) < 0.1)
 
+    def test_random_mb_unique(self):
+
+        gdir = self.gdir
+        flowline.init_present_time_glacier(gdir)
+
+        ref_mod = massbalance.ConstantMassBalance(gdir,
+                                                  halfsize=15)
+        mb_mod = massbalance.RandomMassBalance(gdir, seed=10,
+                                               unique_samples=True,
+                                               halfsize=15)
+        mb_mod2 = massbalance.RandomMassBalance(gdir, seed=20,
+                                                unique_samples=True,
+                                                halfsize=15)
+        mb_mod3 = massbalance.RandomMassBalance(gdir, seed=20,
+                                                unique_samples=True,
+                                                halfsize=15)
+
+        h, w = gdir.get_inversion_flowline_hw()
+
+        ref_mbh = ref_mod.get_annual_mb(h, None) * SEC_IN_YEAR
+
+        # the same year should be equal
+        r_mbh1 = mb_mod.get_annual_mb(h, 1) * SEC_IN_YEAR
+        r_mbh2 = mb_mod.get_annual_mb(h, 1) * SEC_IN_YEAR
+        np.testing.assert_allclose(r_mbh1, r_mbh2)
+
+        # test 31 years (2*halfsize +1)
+        ny = 31
+        yrs = np.arange(ny)
+        mbts = yrs * 0.
+        r_mbh = 0.
+        r_mbh2 = 0.
+        r_mbh3 = 0.
+        mb_mod3.temp_bias = -0.5
+        annual_previous = -999.
+        for i, yr in enumerate(yrs):
+            # specific mass balance
+            mbts[i] = mb_mod.get_specific_mb(h, w, yr)
+
+            # annual mass balance
+            annual = mb_mod.get_annual_mb(h, yr) * SEC_IN_YEAR
+
+            # annual mass balance must be different than the previous one
+            assert not np.all(np.allclose(annual, annual_previous))
+
+            # sum over all years should be equal to ref_mbh
+            r_mbh += annual
+            r_mbh2 += mb_mod2.get_annual_mb(h, yr) * SEC_IN_YEAR
+
+            # mass balance with temperature bias
+            r_mbh3 += mb_mod3.get_annual_mb(h, yr) * SEC_IN_YEAR
+
+            annual_previous = annual
+
+        r_mbh /= ny
+        r_mbh2 /= ny
+        r_mbh3 /= ny
+
+        # test sums
+        np.testing.assert_allclose(ref_mbh, r_mbh, atol=0.02)
+        np.testing.assert_allclose(r_mbh, r_mbh2, atol=0.02)
+
+        # test uniqueness
+        # size
+        self.assertTrue(len(list(mb_mod._state_yr.values())) ==
+                        np.unique(list(mb_mod._state_yr.values())).size)
+        # size2
+        self.assertTrue(len(list(mb_mod2._state_yr.values())) ==
+                        np.unique(list(mb_mod2._state_yr.values())).size)
+        # state years 1 vs 2
+        self.assertTrue(np.all(np.unique(list(mb_mod._state_yr.values())) ==
+                               np.unique(list(mb_mod2._state_yr.values()))))
+        # state years 1 vs reference model
+        self.assertTrue(np.all(np.unique(list(mb_mod._state_yr.values())) ==
+                               ref_mod.years))
+
+        # test ela vs specific mb
+        elats = mb_mod.get_ela(yrs[:200])
+        assert np.corrcoef(mbts[:200], elats)[0, 1] < -0.95
+
+        # test mass balance with temperature bias
+        self.assertTrue(np.mean(r_mbh) < np.mean(r_mbh3))
+
     def test_uncertain_mb(self):
 
         gdir = self.gdir