Add polynomial s shape impact function

climada/entity/impact_funcs/base.py

@@ -173,12 +173,9 @@
 
         """ Step function type impact function.
 
-        By default, everything is destroyed above the step.
+        By default, the impact is 100% above the step.
         Useful for high resolution modelling.
 
-        This method modifies self (climada.entity.impact_funcs instance)
-        by assigning an id, intensity, mdd and paa to the impact function.
-
         Parameters
         ----------
         intensity: tuple(float, float, float)
@@ -226,12 +223,14 @@
         haz_type: str,
         impf_id: int = 1,
         **kwargs):
-        """Sigmoid type impact function hinging on three parameter.
+        r"""Sigmoid type impact function hinging on three parameter.
 
         This type of impact function is very flexible for any sort of study,
         hazard and resolution. The sigmoid is defined as:
 
-        .. math:: f(x) = \\frac{L}{1+exp^{-k(x-x0)}}
+        .. math::
+
+           f(x) = \frac{L}{1+e^{-k(x-x0)}}
 
         For more information: https://en.wikipedia.org/wiki/Logistic_function
 
@@ -240,7 +239,7 @@
 
         Parameters
         ----------
-        intensity: tuple(float, float, float)
+        intensity : tuple(float, float, float)
             tuple of 3 intensity numbers along np.arange(min, max, step)
         L : float
             "top" of sigmoid
@@ -273,3 +272,88 @@
         LOGGER.warning("The use of ImpactFunc.set_sigmoid_impf is deprecated."
                        "Use ImpactFunc.from_sigmoid_impf instead.")
         self.__dict__ = ImpactFunc.from_sigmoid_impf(*args, **kwargs).__dict__
+
+    @classmethod
+    def from_poly_s_shape(
+        cls,
+        intensity: tuple[float, float, int],
+        threshold: float,
+        half_point: float,
+        scale: float,
+        exponent: float,
+        haz_type: str,
+        impf_id: int = 1,
+        **kwargs):
+        r"""S-shape polynomial impact function hinging on four parameter.
+
+        .. math::
+
+            f(I) = \frac{\textrm{luk}(I)^{\textrm{exponent}}}{
+                1 + \textrm{luk}(I)^{\textrm{exponent}}
+            }
+                \cdot \textrm{scale} \\
+            \textrm{luk}(I) = \frac{\max[I - \textrm{threshold}, 0]}{
+                \textrm{half_point} - \textrm{threshold}
+            }
+
+        This function is inspired by Emanuel et al. (2011)
+        https://doi.org/10.1175/WCAS-D-11-00007.1
+
+        This method only specifies mdd, and paa = 1 for all intensities.
+
+        Parameters
+        ----------
+        intensity : tuple(float, float, float)
+            tuple of 3 intensity numbers along np.linsapce(min, max, num)
+        threshold : float
+            Intensity threshold below which there is no impact.
+            In general choose threshold > 0 for computational efficiency
+            of impacts.
+        half_point : float
+            Intensity at which 50% of maximum impact is expected.
+            If half_point <= threshold, mdd = 0 (and f(I)=0) for all
+            intensities.
+        scale : float
+            Multiplicative factor for the whole function. Typically,
+            this sets the maximum value at large intensities.
+        exponent: float
+            Exponent of the polynomial. Value must be exponent >= 0.
+            Emanuel et al. (2011) uses the value 3.
+        haz_type: str
+            Reference string for the hazard (e.g., 'TC', 'RF', 'WS', ...)
+        impf_id : int, optional, default=1
+            Impact function id
+        kwargs :
+            keyword arguments passed to ImpactFunc()
+
+        Raises
+        ------
+        ValueError : if exponent <= 0
+
+        Returns
+        -------
+        impf : climada.entity.impact_funcs.ImpactFunc
+            s-shaped polynomial impact function
+        """
+        if exponent < 0:
+            raise ValueError('Exponent value must larger than 0')
+
+        inten = np.linspace(*intensity)
+
+        if threshold >= half_point:
+            mdd = np.zeros_like(inten)
+        else:
+            luk = (inten - threshold) / (half_point - threshold)
+            luk[luk < 0] = 0
+            mdd = scale * luk**exponent / (1 + luk**exponent)
+        paa = np.ones_like(inten)
+
+        impf = cls(
+            haz_type=haz_type,
+            id=impf_id,
+            intensity=inten,
+            paa=paa,
+            mdd=mdd,
+            **kwargs
+        )
+        return impf

climada/entity/impact_funcs/test/test_base.py

@@ -47,7 +47,6 @@ def test_from_step(self):
         self.assertEqual(imp_fun.haz_type, 'TC')
         self.assertEqual(imp_fun.id, 2)
 
-
     def test_from_sigmoid(self):
         """Check default impact function: sigmoid function"""
         inten = (0, 100, 5)
@@ -60,6 +59,59 @@ def test_from_sigmoid(self):
         self.assertEqual(imp_fun.haz_type, 'RF')
         self.assertEqual(imp_fun.id, 2)
 
+    def test_from_poly_s_shape(self):
+        """Check default impact function: polynomial s-shape"""
+
+        haz_type = 'RF'
+        threshold = 0.2
+        half_point = 1
+        scale = 0.8
+        exponent = 4
+        impf_id = 2
+        unit = 'm'
+        intensity = (0, 5, 5)
+
+        def test_aux_vars(impf):
+            self.assertTrue(np.array_equal(impf.paa, np.ones(5)))
+            self.assertTrue(np.array_equal(impf.intensity, np.linspace(0, 5, 5)))
+            self.assertEqual(impf.haz_type, haz_type)
+            self.assertEqual(impf.id, impf_id)
+            self.assertEqual(impf.intensity_unit, unit)
+
+        impf = ImpactFunc.from_poly_s_shape(
+            intensity=intensity, threshold=threshold, half_point=half_point, scale=scale,
+            exponent=exponent, haz_type=haz_type, impf_id=impf_id, intensity_unit=unit
+            )
+        # True value can easily be computed with a calculator
+        correct_mdd = np.array([0, 0.59836395, 0.78845941, 0.79794213, 0.79938319])
+        np.testing.assert_array_almost_equal(impf.mdd, correct_mdd)
+        test_aux_vars(impf)
+
+        # If threshold > half_point, mdd should all be 0
+        impf = ImpactFunc.from_poly_s_shape(
+            intensity=intensity, threshold=half_point*2, half_point=half_point, scale=scale,
+            exponent=exponent, haz_type=haz_type, impf_id=impf_id, intensity_unit=unit
+            )
+        np.testing.assert_array_almost_equal(impf.mdd, np.zeros(5))
+        test_aux_vars(impf)
+
+        # If exponent = 0, mdd should be constant
+        impf = ImpactFunc.from_poly_s_shape(
+            intensity=intensity, threshold=threshold, half_point=half_point, scale=scale,
+            exponent=0, haz_type=haz_type, impf_id=impf_id, intensity_unit=unit
+            )
+        np.testing.assert_array_almost_equal(impf.mdd, np.ones(5) * scale / 2)
+        test_aux_vars(impf)
+
+        # If exponent < 0, raise error.
+        with self.assertRaisesRegex(ValueError, "Exponent value"):
+            ImpactFunc.from_poly_s_shape(
+                intensity=intensity, threshold=half_point,
+                half_point=half_point, scale=scale,
+                exponent=-1, haz_type=haz_type,
+                impf_id=impf_id, intensity_unit=unit
+            )
+
 # Execute Tests
 if __name__ == "__main__":
     TESTS = unittest.TestLoader().loadTestsFromTestCase(TestInterpolation)