Add polynomial s shape impact function

climada/entity/impact_funcs/base.py

@@ -273,3 +273,74 @@
         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):
+        """S-shape polynomial impact function hinging on four parameter.
+
+        .. math:: f(I) = \\frac{luk(I)**exponent}{1 + luk(I)**exponent}}\\cdot scale
+
+        .. math:: luk(I) = \\frac{max[I - threshold, 0]}{half_point - threshold}
+
+        This function is inspired by Emanuel et al. (2011)
+        https://doi.org/10.1175/WCAS-D-11-00007.1
+
+        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.
+            Should in general be larger than 0 for computational efficiency
+            of impacts.
+        half_point : float
+            Intensity at which 50% of maxixmum impact is expected.
+            If smaller than threshold, mdd = 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. Must be larger than 0.
+        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()
+
+        Returns
+        -------
+        impf : climada.entity.impact_funcs.ImpactFunc
+            s-shapep 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 = luk**exponent / (1 + luk**exponent) * scale
+        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,50 @@ 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)
+
+        with self.assertRaises(ValueError):
+            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)