refactor bootstrap: first resample, then metric

CHANGELOG.rst

@@ -77,6 +77,9 @@ Internals/Minor Fixes
   climpred.stats (:pr:`354`) `Aaron Spring`_.
 - Require ``cftime v1.1.2``, which modifies their object handling to create 200-400x
   speedups in some basic operations. (:pr:`356`) `Riley X. Brady`_.
+- Resample first and then calculate skill in
+  :py:func:`~climpred.bootstrap.bootstrap_perfect_model` and
+  :py:func:`~climpred.bootstrap.bootstrap_hindcast` (:pr:`355`) `Aaron Spring`_.
 
 
 Documentation

asv_bench/benchmarks/benchmarks_hindcast.py

@@ -1,7 +1,9 @@
 import numpy as np
 import xarray as xr
+from dask.distributed import Client
 
 from climpred.bootstrap import bootstrap_hindcast
+from climpred.metrics import PROBABILISTIC_METRICS
 from climpred.prediction import compute_hindcast
 
 from . import ensure_loaded, parameterized, randn, requires_dask
@@ -11,7 +13,7 @@
 # only take comparisons compatible with probabilistic metrics
 HINDCAST_COMPARISONS = ['m2o']
 
-ITERATIONS = 8
+ITERATIONS = 16
 
 
 class Generate:
@@ -30,16 +32,20 @@ def make_hind_obs(self):
         self.uninit = xr.Dataset()
 
         self.nmember = 3
-        self.nlead = 3
-        self.nx = 64
-        self.ny = 64
+        self.nlead = 5
+        self.nx = 72
+        self.ny = 36
+        self.iterations = ITERATIONS
         self.init_start = 1960
         self.init_end = 2000
         self.ninit = self.init_end - self.init_start
+        self.client = None
 
         FRAC_NAN = 0.0
 
-        inits = np.arange(self.init_start, self.init_end)
+        inits = xr.cftime_range(
+            start=str(self.init_start), end=str(self.init_end - 1), freq='YS'
+        )
         leads = np.arange(1, 1 + self.nlead)
         members = np.arange(1, 1 + self.nmember)
 
@@ -86,6 +92,9 @@ def make_hind_obs(self):
         )
 
         self.hind.attrs = {'history': 'created for xarray benchmarking'}
+        self.hind.lead.attrs['units'] = 'years'
+        self.uninit.time.attrs['units'] = 'years'
+        self.observations.time.attrs['units'] = 'years'
 
 
 class Compute(Generate):
@@ -99,48 +108,60 @@ def setup(self, *args, **kwargs):
     @parameterized(['metric', 'comparison'], (METRICS, HINDCAST_COMPARISONS))
     def time_compute_hindcast(self, metric, comparison):
         """Take time for `compute_hindcast`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else 'init'
         ensure_loaded(
             compute_hindcast(
-                self.hind, self.observations, metric=metric, comparison=comparison,
+                self.hind,
+                self.observations,
+                metric=metric,
+                comparison=comparison,
+                dim=dim,
             )
         )
 
     @parameterized(['metric', 'comparison'], (METRICS, HINDCAST_COMPARISONS))
     def peakmem_compute_hindcast(self, metric, comparison):
         """Take memory peak for `compute_hindcast`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else 'init'
         ensure_loaded(
             compute_hindcast(
-                self.hind, self.observations, metric=metric, comparison=comparison,
+                self.hind,
+                self.observations,
+                metric=metric,
+                comparison=comparison,
+                dim=dim,
             )
         )
 
     @parameterized(['metric', 'comparison'], (METRICS, HINDCAST_COMPARISONS))
     def time_bootstrap_hindcast(self, metric, comparison):
         """Take time for `bootstrap_hindcast`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else 'init'
         ensure_loaded(
             bootstrap_hindcast(
                 self.hind,
                 self.uninit,
                 self.observations,
                 metric=metric,
                 comparison=comparison,
-                iterations=ITERATIONS,
-                dim='member',
+                iterations=self.iterations,
+                dim=dim,
             )
         )
 
     @parameterized(['metric', 'comparison'], (METRICS, HINDCAST_COMPARISONS))
     def peakmem_bootstrap_hindcast(self, metric, comparison):
         """Take memory peak for `bootstrap_hindcast`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else 'init'
         ensure_loaded(
             bootstrap_hindcast(
                 self.hind,
                 self.uninit,
                 self.observations,
                 metric=metric,
                 comparison=comparison,
-                iterations=ITERATIONS,
-                dim='member',
+                iterations=self.iterations,
+                dim=dim,
             )
         )
 
@@ -155,11 +176,24 @@ def setup(self, *args, **kwargs):
         # https://github.com/pydata/xarray/blob/stable/asv_bench/benchmarks/rolling.py
         super().setup(**kwargs)
         # chunk along a spatial dimension to enable embarrasingly parallel computation
-        self.hind = self.hind['var'].chunk({'lon': self.nx // ITERATIONS})
-        self.observations = self.observations['var'].chunk(
-            {'lon': self.nx // ITERATIONS}
-        )
-        self.uninit = self.uninit['var'].chunk({'lon': self.nx // ITERATIONS})
+        self.hind = self.hind['var'].chunk()
+        self.observations = self.observations['var'].chunk()
+        self.uninit = self.uninit['var'].chunk()
+
+
+class ComputeDaskDistributed(ComputeDask):
+    def setup(self, *args, **kwargs):
+        """Benchmark time and peak memory of `compute_hindcast` and
+        `bootstrap_hindcast`. This executes the same tests as `Compute` but
+        on chunked data with dask.distributed.Client."""
+        requires_dask()
+        # magic taken from
+        # https://github.com/pydata/xarray/blob/stable/asv_bench/benchmarks/rolling.py
+        super().setup(**kwargs)
+        self.client = Client()
+
+    def cleanup(self):
+        self.client.shutdown()
 
 
 class ComputeSmall(Compute):
@@ -176,3 +210,4 @@ def setup(self, *args, **kwargs):
         self.hind = self.hind.mean(spatial_dims)
         self.observations = self.observations.mean(spatial_dims)
         self.uninit = self.uninit.mean(spatial_dims)
+        self.iterations = 500

asv_bench/benchmarks/benchmarks_perfect_model.py

@@ -1,7 +1,9 @@
 import numpy as np
 import xarray as xr
+from dask.distributed import Client
 
 from climpred.bootstrap import bootstrap_perfect_model
+from climpred.metrics import PROBABILISTIC_METRICS
 from climpred.prediction import compute_perfect_model
 
 from . import ensure_loaded, parameterized, randn, requires_dask
@@ -11,7 +13,7 @@
 # only take comparisons compatible with probabilistic metrics
 PM_COMPARISONS = ['m2m', 'm2c']
 
-ITERATIONS = 8
+ITERATIONS = 16
 
 
 class Generate:
@@ -26,23 +28,32 @@ def make_initialized_control(self):
         perfect-model experiment."""
         self.ds = xr.Dataset()
         self.control = xr.Dataset()
-        self.nmember = 3
-        self.ninit = 4
-        self.nlead = 3
-        self.nx = 64
-        self.ny = 64
+        self.nmember = 5
+        self.ninit = 6
+        self.nlead = 10
+        self.iterations = ITERATIONS
+        self.nx = 72
+        self.ny = 36
         self.control_start = 3000
         self.control_end = 3300
-        self.ntime = 300
+        self.ntime = self.control_end - self.control_start
+        self.client = None
 
         FRAC_NAN = 0.0
 
-        times = np.arange(self.control_start, self.control_end)
+        times = xr.cftime_range(
+            start=str(self.control_start),
+            periods=self.ntime,
+            freq='YS',
+            calendar='noleap',
+        )
         leads = np.arange(1, 1 + self.nlead)
         members = np.arange(1, 1 + self.nmember)
-        inits = (
-            np.random.choice(self.control_end - self.control_start, self.ninit)
-            + self.control_start
+        inits = xr.cftime_range(
+            start=str(self.control_start),
+            periods=self.ninit,
+            freq='10YS',
+            calendar='noleap',
         )
 
         lons = xr.DataArray(
@@ -80,6 +91,8 @@ def make_initialized_control(self):
         )
 
         self.ds.attrs = {'history': 'created for xarray benchmarking'}
+        self.ds.lead.attrs['units'] = 'years'
+        self.control.time.attrs['units'] = 'years'
 
 
 class Compute(Generate):
@@ -94,44 +107,50 @@ def setup(self, *args, **kwargs):
     @parameterized(['metric', 'comparison'], (METRICS, PM_COMPARISONS))
     def time_compute_perfect_model(self, metric, comparison):
         """Take time for `compute_perfect_model`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else None
         ensure_loaded(
             compute_perfect_model(
-                self.ds, self.control, metric=metric, comparison=comparison
+                self.ds, self.control, metric=metric, comparison=comparison, dim=dim
             )
         )
 
     @parameterized(['metric', 'comparison'], (METRICS, PM_COMPARISONS))
     def peakmem_compute_perfect_model(self, metric, comparison):
         """Take memory peak for `compute_perfect_model`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else None
         ensure_loaded(
             compute_perfect_model(
-                self.ds, self.control, metric=metric, comparison=comparison
+                self.ds, self.control, metric=metric, comparison=comparison, dim=dim
             )
         )
 
     @parameterized(['metric', 'comparison'], (METRICS, PM_COMPARISONS))
     def time_bootstrap_perfect_model(self, metric, comparison):
         """Take time for `bootstrap_perfect_model`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else None
         ensure_loaded(
             bootstrap_perfect_model(
                 self.ds,
                 self.control,
                 metric=metric,
                 comparison=comparison,
-                iterations=ITERATIONS,
+                iterations=self.iterations,
+                dim=dim,
             )
         )
 
     @parameterized(['metric', 'comparison'], (METRICS, PM_COMPARISONS))
     def peakmem_bootstrap_perfect_model(self, metric, comparison):
         """Take memory peak for `bootstrap_perfect_model`."""
+        dim = 'member' if metric in PROBABILISTIC_METRICS else None
         ensure_loaded(
             bootstrap_perfect_model(
                 self.ds,
                 self.control,
                 metric=metric,
                 comparison=comparison,
-                iterations=ITERATIONS,
+                iterations=self.iterations,
+                dim=dim,
             )
         )
 
@@ -146,8 +165,23 @@ def setup(self, *args, **kwargs):
         # https://github.com/pydata/xarray/blob/stable/asv_bench/benchmarks/rolling.py
         super().setup(**kwargs)
         # chunk along a spatial dimension to enable embarrasingly parallel computation
-        self.ds = self.ds['var'].chunk({'lon': self.nx // ITERATIONS})
-        self.control = self.control['var'].chunk({'lon': self.nx // ITERATIONS})
+        self.ds = self.ds['var'].chunk()
+        self.control = self.control['var'].chunk()
+
+
+class ComputeDaskDistributed(ComputeDask):
+    def setup(self, *args, **kwargs):
+        """Benchmark time and peak memory of `compute_perfect_model` and
+        `bootstrap_perfect_model`. This executes the same tests as `Compute` but
+        on chunked data with dask.distributed.Client."""
+        requires_dask()
+        # magic taken from
+        # https://github.com/pydata/xarray/blob/stable/asv_bench/benchmarks/rolling.py
+        super().setup(**kwargs)
+        self.client = Client()
+
+    def cleanup(self):
+        self.client.shutdown()
 
 
 class ComputeSmall(Compute):
@@ -158,7 +192,7 @@ def setup(self, *args, **kwargs):
         # magic taken from
         # https://github.com/pydata/xarray/blob/stable/asv_bench/benchmarks/rolling.py
         super().setup(**kwargs)
-        # chunk along a spatial dimension to enable embarrasingly parallel computation
         spatial_dims = ['lon', 'lat']
         self.ds = self.ds.mean(spatial_dims)
         self.control = self.control.mean(spatial_dims)
+        self.iterations = 500

climpred/alignment.py

@@ -89,11 +89,12 @@ def _maximize_alignment(init_lead_matrix, all_verifs, leads):
     # Probably a way to do this more efficiently since we're doing essentially
     # the same thing at each step.
     verif_dates = {
-        l: lead_dependent_verif_dates.sel(lead=l).dropna('time').to_index()
-        for l in leads
+        lead: lead_dependent_verif_dates.sel(lead=lead).dropna('time').to_index()
+        for lead in leads
     }
     inits = {
-        l: lead_dependent_verif_dates.sel(lead=l).dropna('time')['time'] for l in leads
+        lead: lead_dependent_verif_dates.sel(lead=lead).dropna('time')['time']
+        for lead in leads
     }
     return inits, verif_dates
 
@@ -106,9 +107,10 @@ def _same_inits_alignment(init_lead_matrix, valid_inits, all_verifs, leads, n, f
     """
     verifies_at_all_leads = init_lead_matrix.isin(all_verifs).all('lead')
     inits = valid_inits.where(verifies_at_all_leads, drop=True)
-    inits = {l: inits for l in leads}
+    inits = {lead: inits for lead in leads}
     verif_dates = {
-        l: shift_cftime_index(inits[l], 'time', n, freq) for (l, n) in zip(leads, n)
+        lead: shift_cftime_index(inits[lead], 'time', n, freq)
+        for (lead, n) in zip(leads, n)
     }
     return inits, verif_dates
 
@@ -132,10 +134,12 @@ def _same_verifs_alignment(init_lead_matrix, valid_inits, all_verifs, leads, n,
     verif_dates = xr.concat(common_set_of_verifs, 'time').to_index()
     inits_that_verify_with_verif_dates = init_lead_matrix.isin(verif_dates)
     inits = {
-        l: valid_inits.where(inits_that_verify_with_verif_dates.sel(lead=l), drop=True)
-        for l in leads
+        lead: valid_inits.where(
+            inits_that_verify_with_verif_dates.sel(lead=lead), drop=True
+        )
+        for lead in leads
     }
-    verif_dates = {l: verif_dates for l in leads}
+    verif_dates = {lead: verif_dates for lead in leads}
     return inits, verif_dates