Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Lena Ploetzke
5 : *
6 : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
7 : *
8 : * Licensed under the Apache License, Version 2.0 (the "License");
9 : * you may not use this file except in compliance with the License.
10 : * You may obtain a copy of the License at
11 : *
12 : * http://www.apache.org/licenses/LICENSE-2.0
13 : *
14 : * Unless required by applicable law or agreed to in writing, software
15 : * distributed under the License is distributed on an "AS IS" BASIS,
16 : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 : * See the License for the specific language governing permissions and
18 : * limitations under the License.
19 : */
20 :
21 : #ifndef MIO_GLCT_SECIR_MODEL_H
22 : #define MIO_GLCT_SECIR_MODEL_H
23 :
24 : #include "glct_secir/parameters.h"
25 : #include "glct_secir/infection_state.h"
26 : #include "memilio/epidemiology/lct_infection_state.h"
27 : #include "memilio/compartments/compartmentalmodel.h"
28 : #include "memilio/epidemiology/populations.h"
29 : #include "memilio/config.h"
30 : #include "memilio/utils/logging.h"
31 : #include "memilio/utils/time_series.h"
32 : #include "memilio/math/eigen.h"
33 :
34 : namespace mio
35 : {
36 : namespace glsecir
37 : {
38 :
39 : /**
40 : * @brief Class that defines an GLCT-SECIR model.
41 : *
42 : * @tparam NumExposed The number of subcompartments used for the Exposed compartment.
43 : * @tparam NumInfectedNoSymptoms The number of subcompartments used for the InfectedNoSymptoms compartment.
44 : * @tparam NumInfectedSymptoms The number of subcompartments used for the InfectedSymptoms compartment.
45 : * @tparam NumInfectedSevere The number of subcompartments used for the InfectedSevere compartment.
46 : * @tparam NumInfectedCritical The number of subcompartments used for the InfectedCritical compartment.
47 : */
48 : template <size_t NumExposed, size_t NumInfectedNoSymptoms, size_t NumInfectedSymptoms, size_t NumInfectedSevere,
49 : size_t NumInfectedCritical>
50 : class Model
51 : : public CompartmentalModel<
52 : ScalarType,
53 : LctInfectionState<InfectionState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
54 : NumInfectedSevere, NumInfectedCritical, 1, 1>,
55 : mio::Populations<ScalarType,
56 : LctInfectionState<InfectionState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
57 : NumInfectedSevere, NumInfectedCritical, 1, 1>>,
58 : Parameters>
59 : {
60 :
61 : public:
62 : using LctState =
63 : LctInfectionState<InfectionState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms, NumInfectedSevere,
64 : NumInfectedCritical, 1, 1>; ///< This class specifies the number of subcompartments.
65 : using Base = CompartmentalModel<ScalarType, LctState, mio::Populations<ScalarType, LctState>, Parameters>;
66 : using typename Base::ParameterSet;
67 : using typename Base::Populations;
68 :
69 : /// @brief Default constructor.
70 5 : Model()
71 5 : : Base(Populations({Index<LctState>(LctState::Count)}, 0.), ParameterSet())
72 : {
73 5 : }
74 :
75 : /**
76 : * @brief Checks that the model satisfies all constraints (e.g. parameter or population constraints), and
77 : * logs an error if constraints are not satisfied.
78 : *
79 : * @return Returns true if one or more constraints are not satisfied, false otherwise.
80 : */
81 7 : bool check_constraints() const
82 : {
83 7 : auto params = this->parameters;
84 :
85 : // --- Check that the dimensions are consistent. ---
86 14 : if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::Exposed>() !=
87 7 : params.template get<StartingProbabilitiesExposed>().rows()) {
88 1 : log_error("Constraint check: Dimension of the parameters does not match the number of subcompartments for "
89 : "the Exposed "
90 : "compartment.");
91 1 : return true;
92 : }
93 12 : if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() !=
94 6 : params.template get<StartingProbabilitiesInfectedNoSymptoms>().rows()) {
95 1 : log_error(
96 : "Constraint check: Dimension of the parameters does not match the number of subcompartments for the "
97 : "InfectedNoSymptoms compartment.");
98 1 : return true;
99 : }
100 10 : if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedSymptoms>() !=
101 5 : params.template get<StartingProbabilitiesInfectedSymptoms>().rows()) {
102 1 : log_error(
103 : "Constraint check: Dimension of the parameters does not match the number of subcompartments for the "
104 : "InfectedSymptoms compartment.");
105 1 : return true;
106 : }
107 8 : if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedSevere>() !=
108 4 : params.template get<StartingProbabilitiesInfectedSevere>().rows()) {
109 1 : log_error("Constraint check: Dimension of the parameters does not match the number of subcompartments for "
110 : "the InfectedSevere "
111 : "compartment.");
112 1 : return true;
113 : }
114 6 : if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedCritical>() !=
115 3 : params.template get<StartingProbabilitiesInfectedCritical>().rows()) {
116 1 : log_error(
117 : "Constraint check: Dimension of the parameters does not match the number of subcompartments for the "
118 : "InfectedCritical compartment.");
119 1 : return true;
120 : }
121 :
122 2 : return (params.check_constraints() || this->populations.check_constraints());
123 7 : }
124 :
125 : /**
126 : * @brief Evaluates the right-hand-side f of the GLCT dydt = f(y, t).
127 : *
128 : * The GLCT-SECIR model is defined through ordinary differential equations of the form dydt = f(y, t).
129 : * y is a vector containing number of individuals for each (sub-) compartment.
130 : * This function evaluates the right-hand-side f of the ODE and can be used in an ODE solver.
131 : *
132 : * @param[in] pop The current state of the population in the geographic unit we are considering.
133 : * @param[in] y The current state of the model (or a subpopulation) as a flat array.
134 : * @param[in] t The current time.
135 : * @param[out] dydt A reference to the calculated output.
136 : */
137 85 : void get_derivatives(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop,
138 : Eigen::Ref<const Eigen::VectorX<ScalarType>> y, ScalarType t,
139 : Eigen::Ref<Eigen::VectorX<ScalarType>> dydt) const override
140 : {
141 85 : dydt.setZero();
142 :
143 85 : auto params = this->parameters;
144 85 : auto total_population = pop.sum() - pop[LctState::template get_first_index<InfectionState::Dead>()];
145 :
146 : // Calculate sum of all subcompartments for InfectedNoSymptoms.
147 85 : ScalarType InfectedNoSymptoms_sum =
148 85 : pop.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
149 : LctState::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
150 85 : .sum();
151 : // Calculate sum of all subcompartments for InfectedSymptoms.
152 85 : ScalarType InfectedSymptoms_sum =
153 85 : pop.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
154 : LctState::template get_num_subcompartments<InfectionState::InfectedSymptoms>())
155 85 : .sum();
156 :
157 : // --- Susceptibles. ---
158 85 : ScalarType season_val = 1 + params.template get<Seasonality>() *
159 85 : sin(3.141592653589793 * ((params.template get<StartDay>() + t) / 182.5 + 0.5));
160 170 : dydt[0] = -y[0] / total_population * season_val * params.template get<TransmissionProbabilityOnContact>() *
161 170 : params.template get<ContactPatterns>().get_cont_freq_mat().get_matrix_at(t)(0, 0) *
162 85 : (params.template get<RelativeTransmissionNoSymptoms>() * InfectedNoSymptoms_sum +
163 85 : params.template get<RiskOfInfectionFromSymptomatic>() * InfectedSymptoms_sum);
164 :
165 : // --- Exposed. ---
166 255 : dydt.segment(LctState::template get_first_index<InfectionState::Exposed>(),
167 255 : LctState::template get_num_subcompartments<InfectionState::Exposed>()) -=
168 85 : dydt[0] * params.template get<StartingProbabilitiesExposed>();
169 255 : dydt.segment(LctState::template get_first_index<InfectionState::Exposed>(),
170 85 : LctState::template get_num_subcompartments<InfectionState::Exposed>()) +=
171 255 : params.template get<TransitionMatrixExposedToInfectedNoSymptoms>().transpose() *
172 170 : y.segment(LctState::template get_first_index<InfectionState::Exposed>(),
173 : LctState::template get_num_subcompartments<InfectionState::Exposed>());
174 :
175 : // --- InfectedNoSymptoms. ---
176 : // Flow from Exposed To InfectedNoSymptoms.
177 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
178 510 : LctState::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>()) =
179 85 : -(params.template get<TransitionMatrixExposedToInfectedNoSymptoms>() *
180 85 : Eigen::VectorX<ScalarType>::Ones(LctState::template get_num_subcompartments<InfectionState::Exposed>()))
181 : .transpose() *
182 170 : y.segment(LctState::template get_first_index<InfectionState::Exposed>(),
183 : LctState::template get_num_subcompartments<InfectionState::Exposed>()) *
184 85 : params.template get<StartingProbabilitiesInfectedNoSymptoms>();
185 : // Flow from InfectedNoSymptoms To InfectedSymptoms.
186 85 : size_t dimensionInfectedNoSymptomsToInfectedSymptoms =
187 85 : params.template get<TransitionMatrixInfectedNoSymptomsToInfectedSymptoms>().rows();
188 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
189 85 : dimensionInfectedNoSymptomsToInfectedSymptoms) +=
190 255 : params.template get<TransitionMatrixInfectedNoSymptomsToInfectedSymptoms>().transpose() *
191 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
192 : dimensionInfectedNoSymptomsToInfectedSymptoms);
193 : // Flow from InfectedNoSymptoms To Recovered.
194 85 : size_t dimensionInfectedNoSymptomsToRecovered =
195 85 : params.template get<TransitionMatrixInfectedNoSymptomsToRecovered>().rows();
196 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>() +
197 : dimensionInfectedNoSymptomsToInfectedSymptoms,
198 85 : dimensionInfectedNoSymptomsToRecovered) +=
199 255 : params.template get<TransitionMatrixInfectedNoSymptomsToRecovered>().transpose() *
200 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>() +
201 : dimensionInfectedNoSymptomsToInfectedSymptoms,
202 : dimensionInfectedNoSymptomsToRecovered);
203 : // Add flow directly to Recovered compartment.
204 255 : dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
205 85 : -(params.template get<TransitionMatrixInfectedNoSymptomsToRecovered>() *
206 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedNoSymptomsToRecovered))
207 : .transpose() *
208 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>() +
209 : dimensionInfectedNoSymptomsToInfectedSymptoms,
210 : dimensionInfectedNoSymptomsToRecovered);
211 :
212 : // --- InfectedSymptoms. ---
213 : // Flow from InfectedNoSymptoms To InfectedSymptoms.
214 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
215 510 : LctState::template get_num_subcompartments<InfectionState::InfectedSymptoms>()) =
216 85 : -params.template get<StartingProbabilitiesInfectedSymptoms>() *
217 85 : (params.template get<TransitionMatrixInfectedNoSymptomsToInfectedSymptoms>() *
218 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedNoSymptomsToInfectedSymptoms))
219 : .transpose() *
220 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
221 : dimensionInfectedNoSymptomsToInfectedSymptoms);
222 : // Flow from InfectedSymptoms To InfectedSevere.
223 85 : size_t dimensionInfectedSymptomsToInfectedSevere =
224 85 : params.template get<TransitionMatrixInfectedSymptomsToInfectedSevere>().rows();
225 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
226 85 : dimensionInfectedSymptomsToInfectedSevere) +=
227 255 : params.template get<TransitionMatrixInfectedSymptomsToInfectedSevere>().transpose() *
228 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
229 : dimensionInfectedSymptomsToInfectedSevere);
230 : // Flow from InfectedSymptoms To Recovered.
231 85 : size_t dimensionInfectedSymptomsToRecovered =
232 85 : params.template get<TransitionMatrixInfectedSymptomsToRecovered>().rows();
233 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>() +
234 : dimensionInfectedSymptomsToInfectedSevere,
235 85 : dimensionInfectedSymptomsToRecovered) +=
236 255 : params.template get<TransitionMatrixInfectedSymptomsToRecovered>().transpose() *
237 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>() +
238 : dimensionInfectedSymptomsToInfectedSevere,
239 : dimensionInfectedSymptomsToRecovered);
240 : // Add flow directly to Recovered compartment.
241 255 : dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
242 85 : -(params.template get<TransitionMatrixInfectedSymptomsToRecovered>() *
243 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedSymptomsToRecovered))
244 : .transpose() *
245 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>() +
246 : dimensionInfectedSymptomsToInfectedSevere,
247 : dimensionInfectedSymptomsToRecovered);
248 :
249 : // --- InfectedSevere. ---
250 : // Flow from InfectedSymptoms To InfectedSevere.
251 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
252 510 : LctState::template get_num_subcompartments<InfectionState::InfectedSevere>()) =
253 85 : -params.template get<StartingProbabilitiesInfectedSevere>() *
254 85 : (params.template get<TransitionMatrixInfectedSymptomsToInfectedSevere>() *
255 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedSymptomsToInfectedSevere))
256 : .transpose() *
257 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
258 : dimensionInfectedSymptomsToInfectedSevere);
259 : // Flow from InfectedSevere To InfectedCritical.
260 85 : size_t dimensionInfectedSevereToInfectedCritical =
261 85 : params.template get<TransitionMatrixInfectedSevereToInfectedCritical>().rows();
262 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
263 85 : dimensionInfectedSevereToInfectedCritical) +=
264 255 : params.template get<TransitionMatrixInfectedSevereToInfectedCritical>().transpose() *
265 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
266 : dimensionInfectedSevereToInfectedCritical);
267 : // Flow from InfectedSevere To Recovered.
268 85 : size_t dimensionInfectedSevereToRecovered =
269 85 : params.template get<TransitionMatrixInfectedSevereToRecovered>().rows();
270 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
271 : dimensionInfectedSevereToInfectedCritical,
272 85 : dimensionInfectedSevereToRecovered) +=
273 255 : params.template get<TransitionMatrixInfectedSevereToRecovered>().transpose() *
274 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
275 : dimensionInfectedSevereToInfectedCritical,
276 : dimensionInfectedSevereToRecovered);
277 : // Add flow directly to Recovered compartment.
278 255 : dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
279 85 : -(params.template get<TransitionMatrixInfectedSevereToRecovered>() *
280 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedSevereToRecovered))
281 : .transpose() *
282 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
283 : dimensionInfectedSevereToInfectedCritical,
284 : dimensionInfectedSevereToRecovered);
285 :
286 : // --- InfectedCritical. ---
287 : // Flow from InfectedSevere To InfectedCritical.
288 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
289 510 : LctState::template get_num_subcompartments<InfectionState::InfectedCritical>()) =
290 85 : -params.template get<StartingProbabilitiesInfectedCritical>() *
291 85 : (params.template get<TransitionMatrixInfectedSevereToInfectedCritical>() *
292 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedSevereToInfectedCritical))
293 : .transpose() *
294 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
295 : dimensionInfectedSevereToInfectedCritical);
296 : // Flow from InfectedCritical To Dead.
297 85 : size_t dimensionInfectedCriticalToDead = params.template get<TransitionMatrixInfectedCriticalToDead>().rows();
298 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
299 85 : dimensionInfectedCriticalToDead) +=
300 255 : params.template get<TransitionMatrixInfectedCriticalToDead>().transpose() *
301 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
302 : dimensionInfectedCriticalToDead);
303 : // Flow from InfectedCritical To Recovered.
304 85 : size_t dimensionInfectedCriticalToRecovered =
305 85 : params.template get<TransitionMatrixInfectedCriticalToRecovered>().rows();
306 255 : dydt.segment(LctState::template get_first_index<InfectionState::InfectedCritical>() +
307 : dimensionInfectedCriticalToDead,
308 85 : dimensionInfectedCriticalToRecovered) +=
309 255 : params.template get<TransitionMatrixInfectedCriticalToRecovered>().transpose() *
310 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>() +
311 : dimensionInfectedCriticalToDead,
312 : dimensionInfectedCriticalToRecovered);
313 : // Add flow directly to Recovered compartment.
314 255 : dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
315 85 : -(params.template get<TransitionMatrixInfectedCriticalToRecovered>() *
316 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedCriticalToRecovered))
317 : .transpose() *
318 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>() +
319 : dimensionInfectedCriticalToDead,
320 : dimensionInfectedCriticalToRecovered);
321 :
322 : // --- Dead. ---
323 255 : dydt[LctState::template get_first_index<InfectionState::Dead>()] =
324 85 : -(params.template get<TransitionMatrixInfectedCriticalToDead>() *
325 : Eigen::VectorX<ScalarType>::Ones(dimensionInfectedCriticalToDead))
326 : .transpose() *
327 170 : y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
328 : dimensionInfectedCriticalToDead);
329 170 : }
330 :
331 : /**
332 : * @brief Cumulates a simulation result with subcompartments to produce a result that divides the population only
333 : * into the infection states defined in InfectionState.
334 : *
335 : * If the model is used for simulation, we will get a result in form of a TimeSeries with infection states divided
336 : * in subcompartments.
337 : * The function calculates a TimeSeries without subcompartments from another TimeSeries with subcompartments.
338 : * This is done by summing up the corresponding subcompartments.
339 : * @param[in] subcompartments_ts Result of a simulation with the model.
340 : * @return Result of the simulation divided in infection states without subcompartments.
341 : * Returns TimeSeries with values -1 if calculation is not possible.
342 : */
343 2 : TimeSeries<ScalarType> calculate_compartments(const TimeSeries<ScalarType>& subcompartments_ts) const
344 : {
345 2 : TimeSeries<ScalarType> compartments_ts((Eigen::Index)InfectionState::Count);
346 2 : if (!(LctState::Count == subcompartments_ts.get_num_elements())) {
347 1 : log_error("Result does not match InfectionStates of the model.");
348 : // Return a TimeSeries with values -1.
349 1 : Eigen::VectorX<ScalarType> error_output =
350 : Eigen::VectorX<ScalarType>::Constant((Eigen::Index)InfectionState::Count, -1);
351 1 : compartments_ts.add_time_point(-1, error_output);
352 1 : return compartments_ts;
353 1 : }
354 1 : Eigen::VectorX<ScalarType> compartments((Eigen::Index)InfectionState::Count);
355 14 : for (Eigen::Index i = 0; i < subcompartments_ts.get_num_time_points(); ++i) {
356 13 : compartments_ts.add_time_point(subcompartments_ts.get_time(i),
357 : LctState::calculate_compartments(subcompartments_ts[i]));
358 : }
359 1 : return compartments_ts;
360 2 : }
361 : };
362 :
363 : } // namespace glsecir
364 : } // namespace mio
365 :
366 : #endif // MIO_GLCT_SECIR_MODEL_H
|
