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 LCT_SECIR_MODEL_H
22 : #define LCT_SECIR_MODEL_H
23 :
24 : #include "lct_secir/parameters.h"
25 : #include "lct_secir/infection_state.h"
26 : #include "memilio/compartments/compartmentalmodel.h"
27 : #include "memilio/epidemiology/lct_populations.h"
28 : #include "memilio/epidemiology/lct_infection_state.h"
29 : #include "memilio/config.h"
30 : #include "memilio/utils/time_series.h"
31 : #include "memilio/utils/logging.h"
32 : #include "memilio/math/eigen.h"
33 : #include "memilio/utils/type_list.h"
34 : #include "memilio/utils/metaprogramming.h"
35 : namespace mio
36 : {
37 : namespace lsecir
38 : {
39 : /**
40 : * @brief Class that defines an LCT-SECIR model.
41 : *
42 : * @tparam LctStates The LCT model can work with any number of LctStates, where each LctState corresponds to a group,
43 : * e.g. one AgeGroup. The purpose of the LctStates is to define the number of subcompartments for each InfectionState.
44 : * If you do not want to divide the population into groups, just use one LctState.
45 : * If you want to divide the population according to more than one category, e.g. sex and age,
46 : * you have to specify one LctState for each pair of groups, e.g. for (female, A00-A04), (female, A05-A14) etc.
47 : * This is because the number of subcompartments can be different for each group.
48 : * Therefore, the number of LctStates also determines the number of groups.
49 : */
50 : template <class... LctStates>
51 : class Model
52 : : public CompartmentalModel<ScalarType, InfectionState, LctPopulations<ScalarType, LctStates...>, Parameters>
53 : {
54 : public:
55 : using LctStatesGroups = TypeList<LctStates...>;
56 : using Base = CompartmentalModel<ScalarType, InfectionState, LctPopulations<ScalarType, LctStates...>, Parameters>;
57 : using typename Base::ParameterSet;
58 : using typename Base::Populations;
59 : static size_t constexpr num_groups = sizeof...(LctStates);
60 : static_assert(num_groups >= 1, "The number of LctStates provided should be at least one.");
61 :
62 : /// @brief Default constructor.
63 20 : Model()
64 20 : : Base(Populations(), ParameterSet(num_groups))
65 : {
66 20 : }
67 :
68 : /** @brief Constructor using Populations and ParameterSet.
69 : * @param[in] pop An instance of the Populations class.
70 : * @param[in] params Parameters used to be used in the simulation.
71 : */
72 : Model(const Populations& pop, const ParameterSet& params)
73 : : Base(pop, params)
74 : {
75 : }
76 :
77 : /**
78 : * @brief Evaluates the right-hand-side f of the ODE dydt = f(y, t).
79 : *
80 : * The LCT-SECIR model is defined through ordinary differential equations of the form dydt = f(y, t).
81 : * y is a vector containing the number of individuals for each (sub-) compartment.
82 : * This function evaluates the right-hand-side f of the ODE and can be used in an ODE solver.
83 : * @param[in] pop The current state of the population in the geographic unit we are considering.
84 : * @param[in] y The current state of the model (or a subpopulation) as a flat array.
85 : * @param[in] t The current time.
86 : * @param[out] dydt A reference to the calculated output.
87 : */
88 243 : void get_derivatives(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop,
89 : Eigen::Ref<const Eigen::VectorX<ScalarType>> y, ScalarType t,
90 : Eigen::Ref<Eigen::VectorX<ScalarType>> dydt) const override
91 : {
92 : // Vectors are sorted such that we first have all InfectionState%s for AgeGroup 0,
93 : // afterwards all for AgeGroup 1 and so on.
94 243 : dydt.setZero();
95 243 : get_derivatives_impl(pop, y, t, dydt);
96 243 : }
97 :
98 : /**
99 : * @brief Cumulates a simulation result with subcompartments to produce a result that divides the population only
100 : * into the infection states defined in InfectionState.
101 : *
102 : * If the model is used for simulation, we will get a result in form of a TimeSeries with infection states divided
103 : * in subcompartments.
104 : * The function calculates a TimeSeries without subcompartments from another TimeSeries with subcompartments.
105 : * This is done by summing up the numbers in the subcompartments.
106 : * @param[in] subcompartments_ts Result of a simulation with the model.
107 : * @return Result of the simulation divided in infection states without subcompartments.
108 : * Returns TimeSeries with values -1 if calculation is not possible.
109 : */
110 5 : TimeSeries<ScalarType> calculate_compartments(const TimeSeries<ScalarType>& subcompartments_ts) const
111 : {
112 5 : Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count;
113 5 : Eigen::Index num_compartments = count_InfStates * num_groups;
114 5 : TimeSeries<ScalarType> compartments_ts(num_compartments);
115 5 : if (!(this->populations.get_num_compartments() == (size_t)subcompartments_ts.get_num_elements())) {
116 1 : log_error("Result does not match InfectionState of the Model.");
117 1 : Eigen::VectorX<ScalarType> wrong_size = Eigen::VectorX<ScalarType>::Constant(num_compartments, -1);
118 1 : compartments_ts.add_time_point(-1, wrong_size);
119 1 : return compartments_ts;
120 1 : }
121 4 : Eigen::VectorX<ScalarType> compartments(num_compartments);
122 23 : for (Eigen::Index timepoint = 0; timepoint < subcompartments_ts.get_num_time_points(); ++timepoint) {
123 19 : compress_vector(subcompartments_ts[timepoint], compartments);
124 19 : compartments_ts.add_time_point(subcompartments_ts.get_time(timepoint), compartments);
125 : }
126 :
127 4 : return compartments_ts;
128 5 : }
129 :
130 : /**
131 : * @brief Checks that the model satisfies all constraints (e.g. parameter or population constraints).
132 : * @return Returns true if one or more constraints are not satisfied, false otherwise.
133 : */
134 12 : bool check_constraints() const
135 : {
136 :
137 21 : return (check_constraints_impl() || this->parameters.check_constraints() ||
138 21 : this->populations.check_constraints());
139 : }
140 :
141 : private:
142 : /**
143 : * @brief Converts a vector with subcompartments in a vector without subcompartments
144 : * by summing up subcompartment values.
145 : * This is done recursively for each group which corresponds to a slice of the vector.
146 : *
147 : * @tparam group The group specifying the slice of the vector being considered.
148 : * @param[in] subcompartments The vector that should be converted.
149 : * @param[out] compartments Reference to the vector where the output is stored.
150 : */
151 : template <size_t Group = 0>
152 32 : void compress_vector(const Eigen::VectorX<ScalarType>& subcompartments,
153 : Eigen::VectorX<ScalarType>& compartments) const
154 : {
155 : static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
156 : using LctStateGroup = type_at_index_t<Group, LctStates...>;
157 :
158 : // Define first index of the group Group in a vector including all compartments without a resolution
159 : // in subcompartments.
160 32 : Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count;
161 32 : Eigen::Index first_index_group_comps = Group * count_InfStates;
162 :
163 : // Use function from the LctState of the Group to calculate the vector without subcompartments
164 : // using the corresponding vector with subcompartments.
165 64 : compartments.segment(first_index_group_comps, count_InfStates) =
166 96 : LctStateGroup::calculate_compartments(subcompartments.segment(
167 32 : this->populations.template get_first_index_of_group<Group>(), LctStateGroup::Count));
168 :
169 : // Function call for next group if applicable.
170 : if constexpr (Group + 1 < num_groups) {
171 13 : compress_vector<Group + 1>(subcompartments, compartments);
172 : }
173 32 : }
174 :
175 : /**
176 : * @brief Evaluates the right-hand-side f of the ODE dydt = f(y, t) recursively for each group.
177 : *
178 : * See also the function get_derivative.
179 : * For each group, one slice of the output vector is calculated.
180 : * @tparam group The group specifying the slice of the vector being considered.
181 : * @param[in] pop The current state of the population in the geographic unit we are considering.
182 : * @param[in] y The current state of the model (or a subpopulation) as a flat array.
183 : * @param[in] t The current time.
184 : * @param[out] dydt A reference to the calculated output.
185 : */
186 : template <size_t Group = 0>
187 270 : void get_derivatives_impl(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop,
188 : Eigen::Ref<const Eigen::VectorX<ScalarType>> y, ScalarType t,
189 : Eigen::Ref<Eigen::VectorX<ScalarType>> dydt) const
190 : {
191 : static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
192 : using LctStateGroup = type_at_index_t<Group, LctStates...>;
193 :
194 270 : size_t first_index_group = this->populations.template get_first_index_of_group<Group>();
195 270 : auto params = this->parameters;
196 270 : ScalarType flow = 0;
197 :
198 : // Indices of first subcompartment of the InfectionState for the group in the vectors.
199 270 : size_t Ei_first_index = first_index_group + LctStateGroup::template get_first_index<InfectionState::Exposed>();
200 270 : size_t INSi_first_index =
201 270 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedNoSymptoms>();
202 270 : size_t ISyi_first_index =
203 270 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>();
204 270 : size_t ISevi_first_index =
205 270 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSevere>();
206 270 : size_t ICri_first_index =
207 270 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedCritical>();
208 270 : size_t Ri = first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>();
209 270 : size_t Di = first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>();
210 :
211 : // Calculate derivative of the Susceptible compartment.
212 270 : interact<Group, 0>(pop, y, t, dydt);
213 :
214 : // Calculate derivative of the Exposed compartment.
215 270 : dydt[Ei_first_index] = -dydt[first_index_group];
216 665 : for (size_t subcomp = 0; subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>();
217 : subcomp++) {
218 : // Variable flow stores the value of the flow from one subcompartment to the next one.
219 : // Ei_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_first_index
220 : // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms.
221 395 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>() *
222 395 : (1 / params.template get<TimeExposed>()[Group]) * y[Ei_first_index + subcomp];
223 : // Subtract flow from dydt[Ei_first_index + subcomp] and add to next subcompartment.
224 395 : dydt[Ei_first_index + subcomp] -= flow;
225 395 : dydt[Ei_first_index + subcomp + 1] = flow;
226 : }
227 :
228 : // Calculate derivative of the InfectedNoSymptoms compartment.
229 270 : for (size_t subcomp = 0;
230 790 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
231 : subcomp++) {
232 520 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() *
233 520 : (1 / params.template get<TimeInfectedNoSymptoms>()[Group]) * y[INSi_first_index + subcomp];
234 520 : dydt[INSi_first_index + subcomp] -= flow;
235 520 : dydt[INSi_first_index + subcomp + 1] = flow;
236 : }
237 :
238 : // Calculate derivative of the InfectedSymptoms compartment.
239 : // Flow from last (sub-) compartment of C must be split between
240 : // the first subcompartment of InfectedSymptoms and Recovered.
241 270 : dydt[Ri] = dydt[ISyi_first_index] * params.template get<RecoveredPerInfectedNoSymptoms>()[Group];
242 270 : dydt[ISyi_first_index] =
243 270 : dydt[ISyi_first_index] * (1 - params.template get<RecoveredPerInfectedNoSymptoms>()[Group]);
244 270 : for (size_t subcomp = 0;
245 545 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>(); subcomp++) {
246 275 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() *
247 275 : (1 / params.template get<TimeInfectedSymptoms>()[Group]) * y[ISyi_first_index + subcomp];
248 275 : dydt[ISyi_first_index + subcomp] -= flow;
249 275 : dydt[ISyi_first_index + subcomp + 1] = flow;
250 : }
251 :
252 : // Calculate derivative of the InfectedSevere compartment.
253 270 : dydt[Ri] += dydt[ISevi_first_index] * (1 - params.template get<SeverePerInfectedSymptoms>()[Group]);
254 270 : dydt[ISevi_first_index] = dydt[ISevi_first_index] * params.template get<SeverePerInfectedSymptoms>()[Group];
255 270 : for (size_t subcomp = 0;
256 545 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>(); subcomp++) {
257 275 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() *
258 275 : (1 / params.template get<TimeInfectedSevere>()[Group]) * y[ISevi_first_index + subcomp];
259 275 : dydt[ISevi_first_index + subcomp] -= flow;
260 275 : dydt[ISevi_first_index + subcomp + 1] = flow;
261 : }
262 :
263 : // Calculate derivative of the InfectedCritical compartment.
264 270 : dydt[Ri] += dydt[ICri_first_index] * (1 - params.template get<CriticalPerSevere>()[Group]);
265 270 : dydt[ICri_first_index] = dydt[ICri_first_index] * params.template get<CriticalPerSevere>()[Group];
266 270 : for (size_t subcomp = 0;
267 755 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() - 1;
268 : subcomp++) {
269 485 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() *
270 485 : (1 / params.template get<TimeInfectedCritical>()[Group]) * y[ICri_first_index + subcomp];
271 485 : dydt[ICri_first_index + subcomp] -= flow;
272 485 : dydt[ICri_first_index + subcomp + 1] = flow;
273 : }
274 : // Last flow from InfectedCritical has to be divided between Recovered and Dead.
275 : // Must be calculated separately in order not to overwrite the already calculated values for Recovered.
276 270 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() *
277 270 : (1 / params.template get<TimeInfectedCritical>()[Group]) * y[Ri - 1];
278 270 : dydt[Ri - 1] -= flow;
279 270 : dydt[Ri] = dydt[Ri] + (1 - params.template get<DeathsPerCritical>()[Group]) * flow;
280 270 : dydt[Di] = params.template get<DeathsPerCritical>()[Group] * flow;
281 :
282 : // Function call for next group if applicable.
283 : if constexpr (Group + 1 < num_groups) {
284 27 : get_derivatives_impl<Group + 1>(pop, y, t, dydt);
285 : }
286 540 : }
287 :
288 : /**
289 : * @brief Calculates the derivative of the Susceptible compartment for Group1.
290 : *
291 : * This is done recursively by calculating the interaction terms with each group.
292 : * @tparam Group1 The group for which the derivative of the Susceptible compartment should be calculated.
293 : * @tparam Group2 The group that Group1 interacts with.
294 : * @param[in] pop The current state of the population in the geographic unit we are considering.
295 : * @param[in] y The current state of the model (or a subpopulation) as a flat array.
296 : * @param[in] t The current time.
297 : * @param[out] dydt A reference to the calculated output.
298 : */
299 : template <size_t Group1, size_t Group2 = 0>
300 350 : void interact(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop, Eigen::Ref<const Eigen::VectorX<ScalarType>> y,
301 : ScalarType t, Eigen::Ref<Eigen::VectorX<ScalarType>> dydt) const
302 : {
303 : static_assert((Group1 < num_groups) && (Group1 >= 0) && (Group2 < num_groups) && (Group2 >= 0),
304 : "The template parameters Group1 & Group2 should be valid.");
305 : using LctStateGroup2 = type_at_index_t<Group2, LctStates...>;
306 350 : size_t Si_1 = this->populations.template get_first_index_of_group<Group1>();
307 350 : ScalarType infectedNoSymptoms_2 = 0;
308 350 : ScalarType infectedSymptoms_2 = 0;
309 350 : auto params = this->parameters;
310 :
311 350 : size_t first_index_group2 = this->populations.template get_first_index_of_group<Group2>();
312 :
313 : // Calculate sum of all subcompartments for InfectedNoSymptoms of Group2.
314 350 : infectedNoSymptoms_2 =
315 700 : pop.segment(first_index_group2 +
316 350 : LctStateGroup2::template get_first_index<InfectionState::InfectedNoSymptoms>(),
317 : LctStateGroup2::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
318 350 : .sum();
319 : // Calculate sum of all subcompartments for InfectedSymptoms of Group2.
320 350 : infectedSymptoms_2 =
321 700 : pop.segment(first_index_group2 +
322 350 : LctStateGroup2::template get_first_index<InfectionState::InfectedSymptoms>(),
323 : LctStateGroup2::template get_num_subcompartments<InfectionState::InfectedSymptoms>())
324 350 : .sum();
325 : // Size of the Subpopulation Group2 without dead people.
326 350 : double N_2 = pop.segment(first_index_group2, LctStateGroup2::Count - 1).sum();
327 350 : const double divN_2 = (N_2 < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / N_2;
328 350 : ScalarType season_val = 1 + params.template get<Seasonality>() *
329 350 : sin(3.141592653589793 * ((params.template get<StartDay>() + t) / 182.5 + 0.5));
330 700 : dydt[Si_1] += -y[Si_1] * divN_2 * season_val * params.template get<TransmissionProbabilityOnContact>()[Group1] *
331 700 : params.template get<ContactPatterns>().get_cont_freq_mat().get_matrix_at(t)(
332 350 : static_cast<Eigen::Index>(Group1), static_cast<Eigen::Index>(Group2)) *
333 350 : (params.template get<RelativeTransmissionNoSymptoms>()[Group2] * infectedNoSymptoms_2 +
334 350 : params.template get<RiskOfInfectionFromSymptomatic>()[Group2] * infectedSymptoms_2);
335 : // Function call for next interacting group if applicable.
336 : if constexpr (Group2 + 1 < num_groups) {
337 80 : interact<Group1, Group2 + 1>(pop, y, t, dydt);
338 : }
339 700 : }
340 :
341 : /**
342 : * @brief Checks whether LctState of a group satisfies all constraints.
343 : * Recursively, it checks that all groups satisfy the constraints.
344 : *
345 : * @tparam group The group for which the constraints should be checked.
346 : * @return Returns true if one or more constraints are not satisfied, false otherwise.
347 : */
348 : template <size_t Group = 0>
349 22 : bool check_constraints_impl() const
350 : {
351 : static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
352 : using LctStateGroup = type_at_index_t<Group, LctStates...>;
353 :
354 22 : if (LctStateGroup::template get_num_subcompartments<InfectionState::Susceptible>() != 1) {
355 1 : log_warning("Constraint check: The number of subcompartments for Susceptibles of group {} should be one!",
356 : Group);
357 1 : return true;
358 : }
359 21 : if (LctStateGroup::template get_num_subcompartments<InfectionState::Recovered>() != 1) {
360 1 : log_warning("Constraint check: The number of subcompartments for Recovered of group {} should be one!",
361 : Group);
362 1 : return true;
363 : }
364 20 : if (LctStateGroup::template get_num_subcompartments<InfectionState::Dead>() != 1) {
365 1 : log_warning("Constraint check: The number of subcompartments for Dead of group {} should be one!", Group);
366 1 : return true;
367 : }
368 :
369 : if constexpr (Group == num_groups - 1) {
370 9 : return false;
371 : }
372 : else {
373 10 : return check_constraints_impl<Group + 1>();
374 : }
375 : }
376 : };
377 :
378 : } // namespace lsecir
379 : } // namespace mio
380 :
381 : #endif // LCTSECIR_MODEL_H
|
