Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 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 19 : Model()
64 19 : : Base(Populations(), ParameterSet(num_groups))
65 : {
66 19 : }
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 237 : void get_derivatives(Eigen::Ref<const Vector<ScalarType>> pop, Eigen::Ref<const Vector<ScalarType>> y, ScalarType t,
89 : Eigen::Ref<Vector<ScalarType>> dydt) const override
90 : {
91 : // Vectors are sorted such that we first have all InfectionState%s for AgeGroup 0,
92 : // afterwards all for AgeGroup 1 and so on.
93 237 : dydt.setZero();
94 237 : get_derivatives_impl(pop, y, t, dydt);
95 237 : }
96 :
97 : /**
98 : * @brief Cumulates a simulation result with subcompartments to produce a result that divides the population only
99 : * into the infection states defined in InfectionState.
100 : *
101 : * If the model is used for simulation, we will get a result in form of a TimeSeries with infection states divided
102 : * in subcompartments.
103 : * The function calculates a TimeSeries without subcompartments from another TimeSeries with subcompartments.
104 : * This is done by summing up the numbers in the subcompartments.
105 : * @param[in] subcompartments_ts Result of a simulation with the model.
106 : * @return Result of the simulation divided in infection states without subcompartments.
107 : * Returns TimeSeries with values -1 if calculation is not possible.
108 : */
109 5 : TimeSeries<ScalarType> calculate_compartments(const TimeSeries<ScalarType>& subcompartments_ts) const
110 : {
111 5 : Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count;
112 5 : Eigen::Index num_compartments = count_InfStates * num_groups;
113 5 : TimeSeries<ScalarType> compartments_ts(num_compartments);
114 5 : if (!(this->populations.get_num_compartments() == (size_t)subcompartments_ts.get_num_elements())) {
115 2 : log_error("Result does not match InfectionState of the Model.");
116 1 : Vector<ScalarType> wrong_size = Vector<ScalarType>::Constant(num_compartments, -1);
117 1 : compartments_ts.add_time_point(-1, wrong_size);
118 1 : return compartments_ts;
119 1 : }
120 4 : Vector<ScalarType> compartments(num_compartments);
121 23 : for (Eigen::Index timepoint = 0; timepoint < subcompartments_ts.get_num_time_points(); ++timepoint) {
122 19 : compress_vector(subcompartments_ts[timepoint], compartments);
123 19 : compartments_ts.add_time_point(subcompartments_ts.get_time(timepoint), compartments);
124 : }
125 :
126 4 : return compartments_ts;
127 5 : }
128 :
129 : /**
130 : * @brief Checks that the model satisfies all constraints (e.g. parameter or population constraints).
131 : * @return Returns true if one or more constraints are not satisfied, false otherwise.
132 : */
133 12 : bool check_constraints() const
134 : {
135 :
136 21 : return (check_constraints_impl() || this->parameters.check_constraints() ||
137 21 : this->populations.check_constraints());
138 : }
139 :
140 : private:
141 : /**
142 : * @brief Converts a vector with subcompartments in a vector without subcompartments
143 : * by summing up subcompartment values.
144 : * This is done recursively for each group which corresponds to a slice of the vector.
145 : *
146 : * @tparam group The group specifying the slice of the vector being considered.
147 : * @param[in] subcompartments The vector that should be converted.
148 : * @param[out] compartments Reference to the vector where the output is stored.
149 : */
150 : template <size_t Group = 0>
151 32 : void compress_vector(const Vector<ScalarType>& subcompartments, Vector<ScalarType>& compartments) const
152 : {
153 : static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
154 : using LctStateGroup = type_at_index_t<Group, LctStates...>;
155 :
156 : // Specify first indices of the vector slices considered.
157 : // First index of the group Group in an vector including all subcompartments.
158 32 : Eigen::Index first_index_group_subcomps = this->populations.template get_first_index_of_group<Group>();
159 32 : Eigen::Index count_InfStates = (Eigen::Index)InfectionState::Count;
160 : // First index of the group Group in an vector including all compartments without a resolution
161 : // in subcompartments.
162 32 : Eigen::Index first_index_group_comps = Group * count_InfStates;
163 :
164 : // Use segment of the vector subcompartments of each InfectionState and sum up the values of subcompartments.
165 32 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::Susceptible] =
166 : subcompartments[first_index_group_subcomps];
167 64 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::Exposed] =
168 : subcompartments
169 32 : .segment(first_index_group_subcomps +
170 32 : LctStateGroup::template get_first_index<InfectionState::Exposed>(),
171 32 : LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>())
172 64 : .sum();
173 64 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedNoSymptoms] =
174 : subcompartments
175 32 : .segment(first_index_group_subcomps +
176 32 : LctStateGroup::template get_first_index<InfectionState::InfectedNoSymptoms>(),
177 32 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
178 64 : .sum();
179 64 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedSymptoms] =
180 : subcompartments
181 32 : .segment(first_index_group_subcomps +
182 32 : LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>(),
183 32 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>())
184 64 : .sum();
185 64 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedSevere] =
186 : subcompartments
187 32 : .segment(first_index_group_subcomps +
188 32 : LctStateGroup::template get_first_index<InfectionState::InfectedSevere>(),
189 32 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>())
190 64 : .sum();
191 64 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::InfectedCritical] =
192 : subcompartments
193 32 : .segment(first_index_group_subcomps +
194 32 : LctStateGroup::template get_first_index<InfectionState::InfectedCritical>(),
195 32 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
196 64 : .sum();
197 32 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::Recovered] =
198 64 : subcompartments[first_index_group_subcomps +
199 32 : LctStateGroup::template get_first_index<InfectionState::Recovered>()];
200 32 : compartments[first_index_group_comps + (Eigen::Index)InfectionState::Dead] =
201 64 : subcompartments[first_index_group_subcomps +
202 32 : LctStateGroup::template get_first_index<InfectionState::Dead>()];
203 : // Function call for next group if applicable.
204 : if constexpr (Group + 1 < num_groups) {
205 13 : compress_vector<Group + 1>(subcompartments, compartments);
206 : }
207 32 : }
208 :
209 : /**
210 : * @brief Evaluates the right-hand-side f of the ODE dydt = f(y, t) recursively for each group.
211 : *
212 : * See also the function get_derivative.
213 : * For each group, one slice of the output vector is calculated.
214 : * @tparam group The group specifying the slice of the vector being considered.
215 : * @param[in] pop The current state of the population in the geographic unit we are considering.
216 : * @param[in] y The current state of the model (or a subpopulation) as a flat array.
217 : * @param[in] t The current time.
218 : * @param[out] dydt A reference to the calculated output.
219 : */
220 : template <size_t Group = 0>
221 264 : void get_derivatives_impl(Eigen::Ref<const Vector<ScalarType>> pop, Eigen::Ref<const Vector<ScalarType>> y,
222 : ScalarType t, Eigen::Ref<Vector<ScalarType>> dydt) const
223 : {
224 : static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
225 : using LctStateGroup = type_at_index_t<Group, LctStates...>;
226 :
227 264 : size_t first_index_group = this->populations.template get_first_index_of_group<Group>();
228 264 : auto params = this->parameters;
229 264 : ScalarType flow = 0;
230 :
231 : // Indices of first subcompartment of the InfectionState for the group in the vectors.
232 264 : size_t Ei_first_index = first_index_group + LctStateGroup::template get_first_index<InfectionState::Exposed>();
233 264 : size_t INSi_first_index =
234 264 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedNoSymptoms>();
235 264 : size_t ISyi_first_index =
236 264 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>();
237 264 : size_t ISevi_first_index =
238 264 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSevere>();
239 264 : size_t ICri_first_index =
240 264 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedCritical>();
241 264 : size_t Ri = first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>();
242 264 : size_t Di = first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>();
243 :
244 : // Calculate derivative of the Susceptible compartment.
245 264 : interact<Group, 0>(pop, y, t, dydt);
246 :
247 : // Calculate derivative of the Exposed compartment.
248 264 : dydt[Ei_first_index] = -dydt[first_index_group];
249 653 : for (size_t subcomp = 0; subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>();
250 : subcomp++) {
251 : // Variable flow stores the value of the flow from one subcompartment to the next one.
252 : // Ei_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_first_index
253 : // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms.
254 389 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>() *
255 389 : (1 / params.template get<TimeExposed>()[Group]) * y[Ei_first_index + subcomp];
256 : // Subtract flow from dydt[Ei_first_index + subcomp] and add to next subcompartment.
257 389 : dydt[Ei_first_index + subcomp] -= flow;
258 389 : dydt[Ei_first_index + subcomp + 1] = flow;
259 : }
260 :
261 : // Calculate derivative of the InfectedNoSymptoms compartment.
262 264 : for (size_t subcomp = 0;
263 778 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
264 : subcomp++) {
265 514 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() *
266 514 : (1 / params.template get<TimeInfectedNoSymptoms>()[Group]) * y[INSi_first_index + subcomp];
267 514 : dydt[INSi_first_index + subcomp] -= flow;
268 514 : dydt[INSi_first_index + subcomp + 1] = flow;
269 : }
270 :
271 : // Calculate derivative of the InfectedSymptoms compartment.
272 : // Flow from last (sub-) compartment of C must be split between
273 : // the first subcompartment of InfectedSymptoms and Recovered.
274 264 : dydt[Ri] = dydt[ISyi_first_index] * params.template get<RecoveredPerInfectedNoSymptoms>()[Group];
275 264 : dydt[ISyi_first_index] =
276 264 : dydt[ISyi_first_index] * (1 - params.template get<RecoveredPerInfectedNoSymptoms>()[Group]);
277 264 : for (size_t subcomp = 0;
278 533 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>(); subcomp++) {
279 269 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() *
280 269 : (1 / params.template get<TimeInfectedSymptoms>()[Group]) * y[ISyi_first_index + subcomp];
281 269 : dydt[ISyi_first_index + subcomp] -= flow;
282 269 : dydt[ISyi_first_index + subcomp + 1] = flow;
283 : }
284 :
285 : // Calculate derivative of the InfectedSevere compartment.
286 264 : dydt[Ri] += dydt[ISevi_first_index] * (1 - params.template get<SeverePerInfectedSymptoms>()[Group]);
287 264 : dydt[ISevi_first_index] = dydt[ISevi_first_index] * params.template get<SeverePerInfectedSymptoms>()[Group];
288 264 : for (size_t subcomp = 0;
289 533 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>(); subcomp++) {
290 269 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() *
291 269 : (1 / params.template get<TimeInfectedSevere>()[Group]) * y[ISevi_first_index + subcomp];
292 269 : dydt[ISevi_first_index + subcomp] -= flow;
293 269 : dydt[ISevi_first_index + subcomp + 1] = flow;
294 : }
295 :
296 : // Calculate derivative of the InfectedCritical compartment.
297 264 : dydt[Ri] += dydt[ICri_first_index] * (1 - params.template get<CriticalPerSevere>()[Group]);
298 264 : dydt[ICri_first_index] = dydt[ICri_first_index] * params.template get<CriticalPerSevere>()[Group];
299 264 : for (size_t subcomp = 0;
300 749 : subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() - 1;
301 : subcomp++) {
302 485 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() *
303 485 : (1 / params.template get<TimeInfectedCritical>()[Group]) * y[ICri_first_index + subcomp];
304 485 : dydt[ICri_first_index + subcomp] -= flow;
305 485 : dydt[ICri_first_index + subcomp + 1] = flow;
306 : }
307 : // Last flow from InfectedCritical has to be divided between Recovered and Dead.
308 : // Must be calculated separately in order not to overwrite the already calculated values for Recovered.
309 264 : flow = (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() *
310 264 : (1 / params.template get<TimeInfectedCritical>()[Group]) * y[Ri - 1];
311 264 : dydt[Ri - 1] -= flow;
312 264 : dydt[Ri] = dydt[Ri] + (1 - params.template get<DeathsPerCritical>()[Group]) * flow;
313 264 : dydt[Di] = params.template get<DeathsPerCritical>()[Group] * flow;
314 :
315 : // Function call for next group if applicable.
316 : if constexpr (Group + 1 < num_groups) {
317 27 : get_derivatives_impl<Group + 1>(pop, y, t, dydt);
318 : }
319 528 : }
320 :
321 : /**
322 : * @brief Calculates the derivative of the Susceptible compartment for Group1.
323 : *
324 : * This is done recursively by calculating the interaction terms with each group.
325 : * @tparam Group1 The group for which the derivative of the Susceptible compartment should be calculated.
326 : * @tparam Group2 The group that Group1 interacts with.
327 : * @param[in] pop The current state of the population in the geographic unit we are considering.
328 : * @param[in] y The current state of the model (or a subpopulation) as a flat array.
329 : * @param[in] t The current time.
330 : * @param[out] dydt A reference to the calculated output.
331 : */
332 : template <size_t Group1, size_t Group2 = 0>
333 344 : void interact(Eigen::Ref<const Vector<ScalarType>> pop, Eigen::Ref<const Vector<ScalarType>> y, ScalarType t,
334 : Eigen::Ref<Vector<ScalarType>> dydt) const
335 : {
336 : static_assert((Group1 < num_groups) && (Group1 >= 0) && (Group2 < num_groups) && (Group2 >= 0),
337 : "The template parameters Group1 & Group2 should be valid.");
338 : using LctStateGroup2 = type_at_index_t<Group2, LctStates...>;
339 344 : size_t Si_1 = this->populations.template get_first_index_of_group<Group1>();
340 344 : ScalarType infectedNoSymptoms_2 = 0;
341 344 : ScalarType infectedSymptoms_2 = 0;
342 344 : auto params = this->parameters;
343 :
344 344 : size_t first_index_group2 = this->populations.template get_first_index_of_group<Group2>();
345 :
346 : // Calculate sum of all subcompartments for InfectedNoSymptoms of Group2.
347 344 : infectedNoSymptoms_2 =
348 344 : pop.segment(first_index_group2 +
349 344 : LctStateGroup2::template get_first_index<InfectionState::InfectedNoSymptoms>(),
350 344 : LctStateGroup2::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
351 344 : .sum();
352 : // Calculate sum of all subcompartments for InfectedSymptoms of Group2.
353 344 : infectedSymptoms_2 =
354 344 : pop.segment(first_index_group2 +
355 344 : LctStateGroup2::template get_first_index<InfectionState::InfectedSymptoms>(),
356 344 : LctStateGroup2::template get_num_subcompartments<InfectionState::InfectedSymptoms>())
357 344 : .sum();
358 : // Size of the Subpopulation Group2 without dead people.
359 344 : double N_2 = pop.segment(first_index_group2, LctStateGroup2::Count - 1).sum();
360 344 : ScalarType season_val = 1 + params.template get<Seasonality>() *
361 344 : sin(3.141592653589793 * ((params.template get<StartDay>() + t) / 182.5 + 0.5));
362 688 : dydt[Si_1] += -y[Si_1] / N_2 * season_val * params.template get<TransmissionProbabilityOnContact>()[Group1] *
363 688 : params.template get<ContactPatterns>().get_cont_freq_mat().get_matrix_at(t)(
364 344 : static_cast<Eigen::Index>(Group1), static_cast<Eigen::Index>(Group2)) *
365 344 : (params.template get<RelativeTransmissionNoSymptoms>()[Group2] * infectedNoSymptoms_2 +
366 344 : params.template get<RiskOfInfectionFromSymptomatic>()[Group2] * infectedSymptoms_2);
367 : // Function call for next interacting group if applicable.
368 : if constexpr (Group2 + 1 < num_groups) {
369 80 : interact<Group1, Group2 + 1>(pop, y, t, dydt);
370 : }
371 688 : }
372 :
373 : /**
374 : * @brief Checks whether LctState of a group satisfies all constraints.
375 : * Recursively, it checks that all groups satisfy the constraints.
376 : *
377 : * @tparam group The group for which the constraints should be checked.
378 : * @return Returns true if one or more constraints are not satisfied, false otherwise.
379 : */
380 : template <size_t Group = 0>
381 22 : bool check_constraints_impl() const
382 : {
383 : static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
384 : using LctStateGroup = type_at_index_t<Group, LctStates...>;
385 :
386 22 : if (LctStateGroup::template get_num_subcompartments<InfectionState::Susceptible>() != 1) {
387 2 : log_warning("Constraint check: The number of subcompartments for Susceptibles of group {} should be one!",
388 : Group);
389 1 : return true;
390 : }
391 21 : if (LctStateGroup::template get_num_subcompartments<InfectionState::Recovered>() != 1) {
392 2 : log_warning("Constraint check: The number of subcompartments for Recovered of group {} should be one!",
393 : Group);
394 1 : return true;
395 : }
396 20 : if (LctStateGroup::template get_num_subcompartments<InfectionState::Dead>() != 1) {
397 2 : log_warning("Constraint check: The number of subcompartments for Dead of group {} should be one!", Group);
398 1 : return true;
399 : }
400 :
401 : if constexpr (Group == num_groups - 1) {
402 9 : return false;
403 : }
404 : else {
405 10 : return check_constraints_impl<Group + 1>();
406 : }
407 : }
408 : };
409 :
410 : } // namespace lsecir
411 : } // namespace mio
412 :
413 : #endif // LCTSECIR_MODEL_H
|
