feat(interpolation): Vastly improves interpolation of metrics

Author: spjuhel

climada/trajectories/interpolation.py

@@ -21,121 +21,149 @@
 """
 
 import logging
-from abc import ABC, abstractmethod
+from abc import ABC
+from typing import Callable
 
 import numpy as np
-from scipy.sparse import csr_matrix, lil_matrix
 
 LOGGER = logging.getLogger(__name__)
 
 
-class InterpolationStrategy(ABC):
-    """Interface for interpolation strategies."""
-
-    @abstractmethod
-    def interpolate(self, imp_E0, imp_E1, time_points: int) -> list: ...
-
-
-class LinearInterpolation(InterpolationStrategy):
-    """Linear interpolation strategy."""
-
-    def interpolate(self, imp_E0, imp_E1, time_points: int):
+def linear_interp_imp_mat(mat_start, mat_end, interpolation_range) -> list:
+    """Linearly interpolates between two impact matrices over an interpolation range.
+
+    Returns a list of `interpolation_range` matrices linearly interpolated between
+    `mat_start` and `mat_end`.
+    """
+    res = []
+    for point in range(interpolation_range):
+        ratio = point / (interpolation_range - 1)
+        mat_interpolated = mat_start + ratio * (mat_end - mat_start)
+        res.append(mat_interpolated)
+    return res
+
+
+def exponential_interp_imp_mat(mat_start, mat_end, interpolation_range, rate) -> list:
+    """Exponentially interpolates between two impact matrices over an interpolation range with a growth rate `rate`.
+
+    Returns a list of `interpolation_range` matrices exponentially (with growth rate `rate`) interpolated between
+    `mat_start` and `mat_end`.
+    """
+    # Convert matrices to logarithmic domain
+    mat_start = mat_start.copy()
+    mat_end = mat_end.copy()
+    mat_start.data = np.log(mat_start.data + np.finfo(float).eps) / np.log(rate)
+    mat_end.data = np.log(mat_end.data + np.finfo(float).eps) / np.log(rate)
+
+    # Perform linear interpolation in the logarithmic domain
+    res = []
+    for point in range(interpolation_range):
+        ratio = point / (interpolation_range - 1)
+        mat_interpolated = mat_start * (1 - ratio) + ratio * mat_end
+        mat_interpolated.data = np.exp(mat_interpolated.data * np.log(rate))
+        res.append(mat_interpolated)
+    return res
+
+
+def linear_interp_arrays(arr_start, arr_end, interpolation_range):
+    """Perform linear interpolation between two arrays (of a scalar metric) over an interpolation range."""
+    prop1 = np.linspace(0, 1, interpolation_range)
+    prop0 = 1 - prop1
+    if arr_start.ndim > 1:
+        prop0, prop1 = prop0.reshape(-1, 1), prop1.reshape(-1, 1)
+
+    return np.multiply(arr_start, prop0) + np.multiply(arr_end, prop1)
+
+
+def exponential_interp_arrays(arr_start, arr_end, interpolation_range, rate):
+    """Perform exponential interpolation between two arrays (of a scalar metric) over an interpolation range with a growth rate `rate`."""
+    prop1 = np.linspace(0, 1, interpolation_range)
+    prop0 = 1 - prop1
+    if arr_start.ndim > 1:
+        prop0, prop1 = prop0.reshape(-1, 1), prop1.reshape(-1, 1)
+
+    return np.exp(
+        (
+            np.multiply(np.log(arr_start) / np.log(rate), prop0)
+            + np.multiply(np.log(arr_end) / np.log(rate), prop1)
+        )
+        * np.log(rate)
+    )
+
+
+def logarithmic_interp_arrays(arr_start, arr_end, interpolation_range):
+    """Perform logarithmic (natural logarithm) interpolation between two arrays (of a scalar metric) over an interpolation range."""
+    prop1 = np.logspace(0, 1, interpolation_range)
+    prop0 = 1 - prop1
+    if arr_start.ndim > 1:
+        prop0, prop1 = prop0.reshape(-1, 1), prop1.reshape(-1, 1)
+
+    return np.multiply(arr_start, prop0) + np.multiply(arr_end, prop1)
+
+
+class InterpolationStrategyBase(ABC):
+    exposure_interp: Callable
+    hazard_interp: Callable
+    vulnerability_interp: Callable
+
+    def interp_exposure_dim(
+        self, imp_E0, imp_E1, interpolation_range: int, **kwargs
+    ) -> list:
+        """Interpolates along the exposure change between two impact matrices.
+
+        Returns a list of `interpolation_range` matrices linearly interpolated between
+        `mat_start` and `mat_end`.
+        """
         try:
-            return self.interpolate_imp_mat(imp_E0, imp_E1, time_points)
-        except ValueError as e:
-            if str(e) == "inconsistent shapes":
+            res = self.exposure_interp(imp_E0, imp_E1, interpolation_range, **kwargs)
+        except ValueError as err:
+            if str(err) == "inconsistent shapes":
                 raise ValueError(
-                    "Interpolation between impact matrices of different shapes"
+                    "Tried to interpolate impact matrices of different shape. A possible reason could be Exposures of different shapes."
                 )
-            else:
-                raise e
-
-    @staticmethod
-    def interpolate_imp_mat(imp0, imp1, time_points):
-        """Interpolate between two impact matrices over a specified time range.
-
-        Parameters
-        ----------
-        imp0 : ImpactCalc
-            The impact calculation for the starting time.
-        imp1 : ImpactCalc
-            The impact calculation for the ending time.
-        time_points:
-            The number of points to interpolate.
-
-        Returns
-        -------
-        list of np.ndarray
-            List of interpolated impact matrices for each time points in the specified range.
-        """
 
-        def interpolate_sm(mat_start, mat_end, time, time_points):
-            """Perform linear interpolation between two matrices for a specified time point."""
-            if time > time_points:
-                raise ValueError("time point must be within the range")
+            raise err
 
-            ratio = time / (time_points - 1)
+        return res
 
-            # Convert the input matrices to a format that allows efficient modification of its elements
-            mat_start = lil_matrix(mat_start)
-            mat_end = lil_matrix(mat_end)
+    def interp_hazard_dim(
+        self, metric_0, metric_1, interpolation_range: int, **kwargs
+    ) -> np.ndarray:
+        return self.hazard_interp(metric_0, metric_1, interpolation_range, **kwargs)
 
-            # Perform the linear interpolation
-            mat_interpolated = mat_start + ratio * (mat_end - mat_start)
+    def interp_vulnerability_dim(
+        self, metric_0, metric_1, interpolation_range: int, **kwargs
+    ) -> np.ndarray:
+        return self.vulnerability_interp(
+            metric_0, metric_1, interpolation_range, **kwargs
+        )
 
-            return csr_matrix(mat_interpolated)
-
-        LOGGER.debug(f"imp0: {imp0.imp_mat.data[0]}, imp1: {imp1.imp_mat.data[0]}")
-        return [
-            interpolate_sm(imp0.imp_mat, imp1.imp_mat, time, time_points)
-            for time in range(time_points)
-        ]
 
+class InterpolationStrategy(InterpolationStrategyBase):
+    """Interface for interpolation strategies."""
 
-class ExponentialInterpolation(InterpolationStrategy):
-    """Exponential interpolation strategy."""
+    def __init__(self, exposure_interp, hazard_interp, vulnerability_interp) -> None:
+        super().__init__()
+        self.exposure_interp = exposure_interp
+        self.hazard_interp = hazard_interp
+        self.vulnerability_interp = vulnerability_interp
 
-    def interpolate(self, imp_E0, imp_E1, time_points: int):
-        return self.interpolate_imp_mat(imp_E0, imp_E1, time_points)
-
-    @staticmethod
-    def interpolate_imp_mat(imp0, imp1, time_points):
-        """Interpolate between two impact matrices over a specified time range.
-
-        Parameters
-        ----------
-        imp0 : ImpactCalc
-            The impact calculation for the starting time.
-        imp1 : ImpactCalc
-            The impact calculation for the ending time.
-        time_points:
-            The number of points to interpolate.
-
-        Returns
-        -------
-        list of np.ndarray
-            List of interpolated impact matrices for each time points in the specified range.
-        """
 
-        def interpolate_sm(mat_start, mat_end, time, time_points):
-            """Perform exponential interpolation between two matrices for a specified time point."""
-            if time > time_points:
-                raise ValueError("time point must be within the range")
-
-            # Convert matrices to logarithmic domain
-            log_mat_start = np.log(mat_start.toarray() + np.finfo(float).eps)
-            log_mat_end = np.log(mat_end.toarray() + np.finfo(float).eps)
+class AllLinearStrategy(InterpolationStrategyBase):
+    """Linear interpolation strategy."""
 
-            # Perform linear interpolation in the logarithmic domain
-            ratio = time / (time_points - 1)
-            log_mat_interpolated = log_mat_start + ratio * (log_mat_end - log_mat_start)
+    def __init__(self) -> None:
+        super().__init__()
+        self.exposure_interp = linear_interp_imp_mat
+        self.hazard_interp = linear_interp_arrays
+        self.vulnerability_interp = linear_interp_arrays
 
-            # Convert back to the original domain using the exponential function
-            mat_interpolated = np.exp(log_mat_interpolated)
 
-            return csr_matrix(mat_interpolated)
+class ExponentialExposureInterpolation(InterpolationStrategyBase):
+    """Exponential interpolation strategy."""
 
-        return [
-            interpolate_sm(imp0.imp_mat, imp1.imp_mat, time, time_points)
-            for time in range(time_points)
-        ]
+    def __init__(self) -> None:
+        super().__init__()
+        self.exposure_interp = exponential_interp_imp_mat
+        self.hazard_interp = linear_interp_arrays
+        self.vulnerability_interp = linear_interp_arrays

climada/trajectories/risk_trajectory.py

@@ -29,12 +29,12 @@
 import pandas as pd
 
 from climada.entity.disc_rates.base import DiscRates
+from climada.trajectories.interpolation import InterpolationStrategyBase
 from climada.trajectories.riskperiod import (
+    AllLinearStrategy,
     CalcRiskPeriod,
     ImpactCalcComputation,
     ImpactComputationStrategy,
-    InterpolationStrategy,
-    LinearInterpolation,
 )
 from climada.trajectories.snapshot import Snapshot
 
@@ -79,7 +79,7 @@ def __init__(
         risk_transf_cover=None,
         risk_transf_attach=None,
         calc_residual: bool = True,
-        interpolation_strategy: InterpolationStrategy | None = None,
+        interpolation_strategy: InterpolationStrategyBase | None = None,
         impact_computation_strategy: ImpactComputationStrategy | None = None,
     ):
         self._reset_metrics()
@@ -94,7 +94,7 @@ def __init__(
         self._risk_transf_cover = risk_transf_cover
         self._risk_transf_attach = risk_transf_attach
         self._calc_residual = calc_residual
-        self._interpolation_strategy = interpolation_strategy or LinearInterpolation()
+        self._interpolation_strategy = interpolation_strategy or AllLinearStrategy()
         self._impact_computation_strategy = (
             impact_computation_strategy or ImpactCalcComputation()
         )
@@ -283,6 +283,9 @@ def _generic_metrics(
                 raise e
         else:
             tmp = tmp.set_index(["date", "group", "measure", "metric"])
+            if "coord_id" in tmp.columns:
+                tmp = tmp.set_index(["coord_id"], append=True)
+
             tmp = tmp[
                 ~tmp.index.duplicated(keep="last")
             ]  # We want to avoid overlap when more than 2 snapshots
@@ -596,6 +599,15 @@ def plot_per_date_waterfall(
         risk_component = self._calc_waterfall_plot_data(
             start_date=start_date, end_date=end_date
         )
+        risk_component = risk_component[
+            [
+                "base risk",
+                "exposure contribution",
+                "hazard contribution",
+                "vulnerability contribution",
+                "interaction contribution",
+            ]
+        ]
         risk_component.plot(ax=ax, kind="bar", stacked=True)
         # Construct y-axis label and title based on parameters
         value_label = "USD"
@@ -654,21 +666,30 @@ def plot_waterfall(
         labels = [
             f"Risk {start_date}",
             f"Exposure {end_date}",
-            f"Hazard {end_date}¹",
+            f"Hazard {end_date}",
+            f"Vulnerability {end_date}",
+            f"Interaction {end_date}",
             f"Total Risk {end_date}",
         ]
         values = [
             risk_component["base risk"],
-            risk_component["delta from exposure"],
-            risk_component["delta from hazard"],
-            risk_component["base risk"]
-            + risk_component["delta from exposure"]
-            + risk_component["delta from hazard"],
+            risk_component["exposure contribution"],
+            risk_component["hazard contribution"],
+            risk_component["vulnerability contribution"],
+            risk_component["interaction contribution"],
+            risk_component.sum(),
         ]
         bottoms = [
             0.0,
             risk_component["base risk"],
-            risk_component["base risk"] + risk_component["delta from exposure"],
+            risk_component["base risk"] + risk_component["exposure contribution"],
+            risk_component["base risk"]
+            + risk_component["exposure contribution"]
+            + risk_component["hazard contribution"],
+            risk_component["base risk"]
+            + risk_component["exposure contribution"]
+            + risk_component["hazard contribution"]
+            + risk_component["vulnerability contribution"],
             0.0,
         ]
 
@@ -677,7 +698,14 @@ def plot_waterfall(
             values,
             bottom=bottoms,
             edgecolor="black",
-            color=["tab:blue", "tab:orange", "tab:green", "tab:red"],
+            color=[
+                "tab:cyan",
+                "tab:orange",
+                "tab:green",
+                "tab:red",
+                "tab:purple",
+                "tab:blue",
+            ],
         )
         for i in range(len(values)):
             ax.text(
@@ -695,16 +723,19 @@ def plot_waterfall(
 
         ax.set_title(title_label)
         ax.set_ylabel(value_label)
-        # ax.tick_params(axis='x', labelrotation=90,)
-        ax.annotate(
-            """¹: The increase in risk due to hazard denotes the difference in risk with future exposure
-and hazard compared to risk with future exposure and present hazard.""",
-            xy=(0.0, -0.15),
-            xycoords="axes fraction",
-            ha="left",
-            va="center",
-            fontsize=8,
+        ax.tick_params(
+            axis="x",
+            labelrotation=90,
         )
+        #         ax.annotate(
+        #             """¹: The increase in risk due to hazard denotes the difference in risk with future exposure
+        # and hazard compared to risk with future exposure and present hazard.""",
+        #             xy=(0.0, -0.15),
+        #             xycoords="axes fraction",
+        #             ha="left",
+        #             va="center",
+        #             fontsize=8,
+        #         )
 
         return ax