Merge pull request #405 from ClimateImpactLab/quantreg_sccs

Add quantreg full cloud of points across batches

Author: kemccusker

src/dscim/menu/main_recipe.py

@@ -26,7 +26,7 @@ class MainRecipe(StackedDamages, ABC):
     Parameters
     ----------
     discounting_type : str
-        Choice of discounting: ``euler_gwr``, ``euler_ramsey``, ``constant``, ``naive_ramsey``,
+        Choice of discounting: ``euler_gwr``, ``euler_ramsey``, ``constant``, ``constant_gwr``, ``naive_ramsey``,
         ``naive_gwr``, ``gwr_gwr``.
     discrete_discounting: boolean
         Discounting is discrete if ``True``, else continuous (default is ``False``).
@@ -52,6 +52,7 @@ class MainRecipe(StackedDamages, ABC):
     DISCOUNT_TYPES = [
         "constant",
         "constant_model_collapsed",
+        "constant_gwr",
         "naive_ramsey",
         "euler_ramsey",
         "naive_gwr",
@@ -253,7 +254,11 @@ def __init__(
         # 'constant_model_collapsed' should be here except that we allow
         # for a collapsed-model Ramsey rate to be calculated (for labour
         # and energy purposes)
-        if self.discounting_type in ["constant", "constant_model_collapsed"]:
+        if self.discounting_type in [
+            "constant",
+            "constant_model_collapsed",
+            "constant_gwr",
+        ]:
             self.stream_discount_factors = None
 
         # assert formulas for which clip_gmsl is implemented
@@ -318,8 +323,10 @@ def scc():
             self.logger.info("Processing SCC calculation ...")
             if self.fit_type == "quantreg":
                 self.full_uncertainty_iqr
-                self.calculate_scc
-                self.stat_uncertainty_iqr
+                # stat_uncertainty_iqr function expects collapsed SCCs, so a fair aggregation is required
+                if len(self.fair_aggregation) > 0:
+                    self.calculate_scc
+                    self.stat_uncertainty_iqr
             else:
                 if len(self.fair_aggregation) > 0:
                     self.stream_discount_factors
@@ -1093,7 +1100,7 @@ def discounted_damages(self, damages, discrate):
             xr.Dataset
         """
 
-        if discrate in ["constant", "constant_model_collapsed"]:
+        if discrate in ["constant", "constant_model_collapsed", "constant_gwr"]:
             if self.discrete_discounting:
                 discrate_damages = [
                     damages * (1 / (1 + r)) ** (damages.year - self.climate.pulse_year)

src/dscim/preprocessing/preprocessing.py

@@ -21,6 +21,7 @@ def ce_from_chunk(
     zero,
     socioec,
     ce_batch_coords,
+    quantreg=False,
 ):
     year = chunk.year.values
     ssp = chunk.ssp.values
@@ -44,13 +45,19 @@ def ce_from_chunk(
         raise NotImplementedError("Pass 'cc' or 'no_cc' to reduction.")
 
     if recipe == "adding_up":
-        result = mean_func(
-            np.maximum(
+        if not quantreg:
+            result = mean_func(
+                np.maximum(
+                    calculation,
+                    bottom_code,
+                ),
+                "batch",
+            )
+        else:
+            result = np.maximum(
                 calculation,
                 bottom_code,
-            ),
-            "batch",
-        )
+            )
     elif recipe == "risk_aversion":
         result = ce_func(
             np.maximum(
@@ -73,13 +80,17 @@ def reduce_damages(
     socioec,
     bottom_coding_gdppc=39.39265060424805,
     zero=False,
+    quantreg=False,
 ):
     if recipe == "adding_up":
         assert (
             eta is None
         ), "Adding up does not take an eta argument. Please set to None."
     # client = Client(n_workers=35, memory_limit="9G", threads_per_worker=1)
 
+    if recipe == "risk_aversion":
+        assert not quantreg, "Quantile regression is not compatible with risk aversion. Please set quantreg to False."
+
     with open(config) as stream:
         c = yaml.safe_load(stream)
         params = c["sectors"][sector]
@@ -112,10 +123,15 @@ def reduce_damages(
                     "model": 1,
                     "ssp": 1,
                 }
-
-            ce_batch_dims = [i for i in gdppc.dims] + [
-                i for i in ds.dims if i not in gdppc.dims and i != "batch"
-            ]
+            if quantreg:
+                chunkies["batch"] = 1
+                ce_batch_dims = [i for i in gdppc.dims] + [
+                    i for i in ds.dims if i not in gdppc.dims
+                ]
+            else:
+                ce_batch_dims = [i for i in gdppc.dims] + [
+                    i for i in ds.dims if i not in gdppc.dims and i != "batch"
+                ]
             ce_batch_coords = {c: ds[c].values for c in ce_batch_dims}
             ce_batch_coords["region"] = [
                 i for i in gdppc.region.values if i in ce_batch_coords["region"]
@@ -143,6 +159,7 @@ def reduce_damages(
             zero=zero,
             socioec=socioec,
             ce_batch_coords=ce_batch_coords,
+            quantreg=quantreg,
         ),
         template=template,
     )