Skip to content

Commit 31ea96d

Browse files
spjuhelspjuhel
committed
feat(trajectory): interpolation and impactcalc files, eai_gdf metric
1 parent 0d6335b commit 31ea96d

File tree

4 files changed

+330
-212
lines changed

4 files changed

+330
-212
lines changed

‎climada/trajectories/impact_calc_strat.py‎

Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
"""
2+
This file is part of CLIMADA.
3+
4+
Copyright (C) 2017 ETH Zurich, CLIMADA contributors listed in AUTHORS.
5+
6+
CLIMADA is free software: you can redistribute it and/or modify it under the
7+
terms of the GNU General Public License as published by the Free
8+
Software Foundation, version 3.
9+
10+
CLIMADA is distributed in the hope that it will be useful, but WITHOUT ANY
11+
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
12+
PARTICULAR PURPOSE. See the GNU General Public License for more details.
13+
14+
You should have received a copy of the GNU General Public License along
15+
with CLIMADA. If not, see <https://www.gnu.org/licenses/>.
16+
17+
---
18+
19+
This modules implements the Snapshot and SnapshotsCollection classes.
20+
21+
"""
22+
23+
import copy
24+
from abc import ABC, abstractmethod
25+
26+
import numpy as np
27+
28+
from climada.engine.impact import Impact
29+
from climada.engine.impact_calc import ImpactCalc
30+
from climada.trajectories.snapshot import Snapshot
31+
32+
33+
class ImpactComputationStrategy(ABC):
34+
"""Interface for impact computation strategies."""
35+
36+
@abstractmethod
37+
def compute_impacts(
38+
self,
39+
snapshot0: Snapshot,
40+
snapshot1: Snapshot,
41+
risk_transf_attach: float | None,
42+
risk_transf_cover: float | None,
43+
calc_residual: bool,
44+
) -> tuple:
45+
pass
46+
47+
48+
class ImpactCalcComputation(ImpactComputationStrategy):
49+
"""Default impact computation strategy."""
50+
51+
def compute_impacts(
52+
self,
53+
snapshot0: Snapshot,
54+
snapshot1: Snapshot,
55+
risk_transf_attach: float | None,
56+
risk_transf_cover: float | None,
57+
calc_residual: bool = False,
58+
):
59+
impacts = self._calculate_impacts_for_snapshots(snapshot0, snapshot1)
60+
self._apply_risk_transfer(
61+
impacts, risk_transf_attach, risk_transf_cover, calc_residual
62+
)
63+
return impacts
64+
65+
def _calculate_impacts_for_snapshots(
66+
self, snapshot0: Snapshot, snapshot1: Snapshot
67+
):
68+
"""Calculate impacts for the given snapshots and impact function set."""
69+
imp_E0H0 = ImpactCalc(
70+
snapshot0.exposure, snapshot0.impfset, snapshot0.hazard
71+
).impact()
72+
imp_E1H0 = ImpactCalc(
73+
snapshot1.exposure, snapshot1.impfset, snapshot0.hazard
74+
).impact()
75+
imp_E0H1 = ImpactCalc(
76+
snapshot0.exposure, snapshot0.impfset, snapshot1.hazard
77+
).impact()
78+
imp_E1H1 = ImpactCalc(
79+
snapshot1.exposure, snapshot1.impfset, snapshot1.hazard
80+
).impact()
81+
return imp_E0H0, imp_E1H0, imp_E0H1, imp_E1H1
82+
83+
def _apply_risk_transfer(
84+
self,
85+
impacts: tuple[Impact, Impact, Impact, Impact],
86+
risk_transf_attach: float | None,
87+
risk_transf_cover: float | None,
88+
calc_residual: bool,
89+
):
90+
"""Apply risk transfer to the calculated impacts."""
91+
if risk_transf_attach is not None and risk_transf_cover is not None:
92+
for imp in impacts:
93+
imp.imp_mat = self.calculate_residual_or_risk_transfer_impact_matrix(
94+
imp.imp_mat, risk_transf_attach, risk_transf_cover, calc_residual
95+
)
96+
97+
def calculate_residual_or_risk_transfer_impact_matrix(
98+
self, imp_mat, risk_transf_attach, risk_transf_cover, calc_residual
99+
):
100+
"""
101+
Calculate either the residual or the risk transfer impact matrix.
102+
103+
The impact matrix is adjusted based on the total impact for each event.
104+
When calculating the residual impact, the result is the total impact minus
105+
the risk layer. The risk layer is defined as the minimum of the cover and
106+
the maximum of the difference between the total impact and the attachment.
107+
If `calc_residual` is False, the function returns the risk layer matrix
108+
instead of the residual.
109+
110+
Parameters
111+
----------
112+
imp_mat : scipy.sparse.csr_matrix
113+
The original impact matrix to be scaled.
114+
attachment : float, optional
115+
The attachment point for the risk layer.
116+
cover : float, optional
117+
The maximum coverage for the risk layer.
118+
calc_residual : bool, default=True
119+
Determines if the function calculates the residual (if True) or the
120+
risk layer (if False).
121+
122+
Returns
123+
-------
124+
scipy.sparse.csr_matrix
125+
The adjusted impact matrix, either residual or risk transfer.
126+
127+
Example
128+
-------
129+
>>> calc_residual_or_risk_transf_imp_mat(imp_mat, attachment=100, cover=500, calc_residual=True)
130+
Residual impact matrix with applied risk layer adjustments.
131+
"""
132+
if risk_transf_attach and risk_transf_cover:
133+
# Make a copy of the impact matrix
134+
imp_mat = copy.deepcopy(imp_mat)
135+
# Calculate the total impact per event
136+
total_at_event = imp_mat.sum(axis=1).A1
137+
# Risk layer at event
138+
transfer_at_event = np.minimum(
139+
np.maximum(total_at_event - risk_transf_attach, 0), risk_transf_cover
140+
)
141+
# Resiudal impact
142+
residual_at_event = np.maximum(total_at_event - transfer_at_event, 0)
143+
144+
# Calculate either the residual or transfer impact matrix
145+
# Choose the denominator to rescale the impact values
146+
if calc_residual:
147+
# Rescale the impact values
148+
numerator = residual_at_event
149+
else:
150+
# Rescale the impact values
151+
numerator = transfer_at_event
152+
153+
# Rescale the impact values
154+
rescale_impact_values = np.divide(
155+
numerator,
156+
total_at_event,
157+
out=np.zeros_like(numerator, dtype=float),
158+
where=total_at_event != 0,
159+
)
160+
161+
# The multiplication is broadcasted across the columns for each row
162+
result_matrix = imp_mat.multiply(rescale_impact_values[:, np.newaxis])
163+
164+
return result_matrix
165+
166+
else:
167+
168+
return imp_mat

‎climada/trajectories/interpolation.py‎

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
import logging
2+
from abc import ABC, abstractmethod
3+
4+
import numpy as np
5+
from scipy.sparse import lil_matrix
6+
7+
LOGGER = logging.getLogger(__name__)
8+
9+
10+
class InterpolationStrategy(ABC):
11+
"""Interface for interpolation strategies."""
12+
13+
@abstractmethod
14+
def interpolate(self, imp_E0, imp_E1, time_points: int) -> list: ...
15+
16+
17+
class LinearInterpolation(InterpolationStrategy):
18+
"""Linear interpolation strategy."""
19+
20+
def interpolate(self, imp_E0, imp_E1, time_points: int):
21+
try:
22+
return self.interpolate_imp_mat(imp_E0, imp_E1, time_points)
23+
except ValueError as e:
24+
if str(e) == "inconsistent shape":
25+
raise ValueError(
26+
"Interpolation between impact matrices of different shapes"
27+
)
28+
else:
29+
raise e
30+
31+
@staticmethod
32+
def interpolate_imp_mat(imp0, imp1, time_points):
33+
"""Interpolate between two impact matrices over a specified time range.
34+
35+
Parameters
36+
----------
37+
imp0 : ImpactCalc
38+
The impact calculation for the starting time.
39+
imp1 : ImpactCalc
40+
The impact calculation for the ending time.
41+
time_points:
42+
The number of points to interpolate.
43+
44+
Returns
45+
-------
46+
list of np.ndarray
47+
List of interpolated impact matrices for each time points in the specified range.
48+
"""
49+
50+
def interpolate_sm(mat_start, mat_end, time, time_points):
51+
"""Perform linear interpolation between two matrices for a specified time point."""
52+
if time > time_points:
53+
raise ValueError("time point must be within the range")
54+
55+
ratio = time / (time_points - 1)
56+
57+
# Convert the input matrices to a format that allows efficient modification of its elements
58+
mat_start = lil_matrix(mat_start)
59+
mat_end = lil_matrix(mat_end)
60+
61+
# Perform the linear interpolation
62+
mat_interpolated = mat_start + ratio * (mat_end - mat_start)
63+
64+
return mat_interpolated
65+
66+
LOGGER.debug(f"imp0: {imp0.imp_mat.data[0]}, imp1: {imp1.imp_mat.data[0]}")
67+
return [
68+
interpolate_sm(imp0.imp_mat, imp1.imp_mat, time, time_points)
69+
for time in range(time_points)
70+
]
71+
72+
73+
class ExponentialInterpolation:
74+
"""Exponential interpolation strategy."""
75+
76+
def interpolate(self, imp_E0, imp_E1, time_points: int):
77+
try:
78+
return self.interpolate_imp_mat(imp_E0, imp_E1, time_points)
79+
except ValueError as e:
80+
if str(e) == "inconsistent shape":
81+
raise ValueError(
82+
"Interpolation between impact matrices of different shapes"
83+
)
84+
else:
85+
raise e
86+
87+
@staticmethod
88+
def interpolate_imp_mat(imp0, imp1, time_points):
89+
"""Interpolate between two impact matrices over a specified time range.
90+
91+
Parameters
92+
----------
93+
imp0 : ImpactCalc
94+
The impact calculation for the starting time.
95+
imp1 : ImpactCalc
96+
The impact calculation for the ending time.
97+
time_points:
98+
The number of points to interpolate.
99+
100+
Returns
101+
-------
102+
list of np.ndarray
103+
List of interpolated impact matrices for each time points in the specified range.
104+
"""
105+
106+
def interpolate_sm(mat_start, mat_end, time, time_points):
107+
"""Perform exponential interpolation between two matrices for a specified time point."""
108+
if time > time_points:
109+
raise ValueError("time point must be within the range")
110+
111+
# Convert matrices to logarithmic domain
112+
log_mat_start = np.log(mat_start.toarray() + np.finfo(float).eps)
113+
log_mat_end = np.log(mat_end.toarray() + np.finfo(float).eps)
114+
115+
# Perform linear interpolation in the logarithmic domain
116+
ratio = time / (time_points - 1)
117+
log_mat_interpolated = log_mat_start + ratio * (log_mat_end - log_mat_start)
118+
119+
# Convert back to the original domain using the exponential function
120+
mat_interpolated = np.exp(log_mat_interpolated)
121+
122+
return lil_matrix(mat_interpolated)
123+
124+
return [
125+
interpolate_sm(imp0.imp_mat, imp1.imp_mat, time, time_points)
126+
for time in range(time_points)
127+
]
128+
129+
130+
# Example usage
131+
# Assuming imp0 and imp1 are instances of ImpactCalc with imp_mat attributes as sparse matrices
132+
# interpolator = ExponentialInterpolation()
133+
# interpolated_matrices = interpolator.interpolate(imp0, imp1, 100)

‎climada/trajectories/risk_trajectory.py‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -78,12 +78,7 @@ def __init__(
7878
interpolation_strategy: InterpolationStrategy | None = None,
7979
impact_computation_strategy: ImpactComputationStrategy | None = None,
8080
):
81-
self._aai_metrics = None
82-
self._return_periods_metrics = None
83-
self._risk_components_metrics = None
84-
self._aai_per_group_metrics = None
85-
self._all_risk_metrics = None
86-
self._metrics_up_to_date = False
81+
self._reset_metrics()
8782
self._risk_period_up_to_date: bool = False
8883
self._snapshots = snapshots_list
8984
self._all_groups_name = all_groups_name
@@ -103,6 +98,7 @@ def __init__(
10398
self._risk_periods_calculators = self._calc_risk_periods(snapshots_list)
10499

105100
def _reset_metrics(self):
101+
self._eai_metrics = None
106102
self._aai_metrics = None
107103
self._return_periods_metrics = None
108104
self._risk_components_metrics = None

0 commit comments

Comments
 (0)