Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Anna Wendler, Lena Ploetzke, Martin J. Kuehn
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 IDESECIR_MODEL_H
22 : #define IDESECIR_MODEL_H
23 :
24 : #include "ide_secir/parameters.h"
25 : #include "ide_secir/infection_state.h"
26 : #include "memilio/config.h"
27 : #include "memilio/epidemiology/age_group.h"
28 : #include "memilio/utils/custom_index_array.h"
29 : #include "memilio/utils/time_series.h"
30 :
31 : #include "vector"
32 :
33 : namespace mio
34 : {
35 : namespace isecir
36 : {
37 : class Model
38 : {
39 : using ParameterSet = Parameters;
40 :
41 : public:
42 : /**
43 : * @brief Constructor to create an IDE-SECIR model.
44 : *
45 : * @param[in, out] init TimeSeries with the initial values of the number of individuals,
46 : * which transit within one timestep dt from one compartment to another.
47 : * Possible transitions are specified in #InfectionTransition%s.
48 : * Considered time points should have the distance dt. The last time point determines the start time t0 of the
49 : * simulation.
50 : * The time history must reach a certain point in the past so that the simulation can be performed.
51 : * A warning is displayed if the condition is violated.
52 : * @param[in] N_init A vector, containing the populations of the considered region, for every AgeGroup.
53 : * @param[in] deaths A vector, containing the total number of deaths at time t0, for every AgeGroup.
54 : * @param[in] num_agegroups The number of AgeGroups.
55 : * @param[in] total_confirmed_cases A vector, containing the total confirmed cases at time t0 can be set if it
56 : * should be used for initialization, for every AgeGroup.
57 : */
58 : Model(TimeSeries<ScalarType>&& init, CustomIndexArray<ScalarType, AgeGroup> N_init,
59 : CustomIndexArray<ScalarType, AgeGroup> deaths, size_t num_agegroups,
60 : CustomIndexArray<ScalarType, AgeGroup> total_confirmed_cases = CustomIndexArray<ScalarType, AgeGroup>());
61 :
62 : // ---- Additional functionality such as constraint checking, setters and getters, etc. ----
63 : /**
64 : * @brief Checks constraints on model parameters and initial data.
65 : * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false.
66 : */
67 : bool check_constraints(ScalarType dt) const;
68 :
69 : /**
70 : * @brief Returns a flat index for the InfectionTransition TimeSeries.
71 : *
72 : * In the TimeSeries we store a vector for each time point. In this vector we store the different
73 : * #InfectionTransition%s for every AgeGroup.
74 : * This function is used to get the right index in this vector for a specific AgeGroup and InfectionTransition.
75 : *
76 : * @param[in] transition_idx Index determining which InfectionTransition we want to evaluate.
77 : * @param[in] agegroup The agegroup for which we want to evaluate.
78 : */
79 :
80 5590 : int get_transition_flat_index(Eigen::Index transition_idx, AgeGroup agegroup) const
81 : {
82 5590 : return (static_cast<int>(size_t(agegroup)) * int(InfectionTransition::Count) + int(transition_idx));
83 : }
84 :
85 : /**
86 : * @brief Returns a flat index for the InfectionState TimeSeries.
87 : *
88 : * In the TimeSeries we store a vector for each time point. In this vector we store the different InfectionStates
89 : * for every AgeGroup.
90 : * This function is used to get the right index in this vector for a specific AgeGroup and InfectionState.
91 : *
92 : * @param[in] state_idx Index at which InfectionState we want to evaluate.
93 : * @param[in] agegroup The agegroup for which we want to evaluate.
94 : */
95 1919 : int get_state_flat_index(Eigen::Index state_idx, AgeGroup agegroup) const
96 : {
97 1919 : return (static_cast<int>(size_t(agegroup)) * int(InfectionState::Count) + int(state_idx));
98 : }
99 :
100 : /**
101 : * @brief Getter for the global support_max, i.e. the maximum of support_max over all TransitionDistributions.
102 : *
103 : * This determines how many initial values we need for the flows.
104 : * It may be possible to run the simulation with fewer time points than the value of the global support_max,
105 : * but this number ensures that it is possible.
106 : *
107 : * @param[in] dt Time step size.
108 : *
109 : * @return Global support_max.
110 : */
111 : ScalarType get_global_support_max(ScalarType dt) const;
112 :
113 : /**
114 : * @brief Returns the index of the automatically selected initialization method.
115 : *
116 : * The initialization method is selected automatically based on the different values that need to be set beforehand.
117 : * Infection compartments are always computed through historic flow.
118 : * Initialization methods for Susceptible and Recovered are tested in the following order:
119 : * 1.) If a positive number for the total number of confirmed cases is set, Recovered is set according to that value
120 : * and Susceptible%s are derived.
121 : * 2.) If Susceptible%s are set, Recovered will be derived.
122 : * 3.) If Recovered are set directly, Susceptible%s are derived.
123 : * 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021) to
124 : * set the Susceptible%s.
125 : *
126 : * @return Index representing the initialization method.
127 : */
128 7 : int get_initialization_method_compartments() const
129 : {
130 7 : return m_initialization_method;
131 : }
132 :
133 : /**
134 : * @brief Setter for the tolerance used to calculate the maximum support of the TransitionDistributions.
135 : *
136 : * @param[in] new_tol New tolerance.
137 : */
138 3 : void set_tol_for_support_max(ScalarType new_tol)
139 : {
140 3 : m_tol = new_tol;
141 3 : }
142 :
143 : // ---- Public parameters. ----
144 : ParameterSet parameters{AgeGroup(m_num_agegroups)}; ///< ParameterSet of Model Parameters.
145 : // Attention: m_populations and m_transitions do not necessarily have the same number of time points due to the
146 : // initialization part.
147 : TimeSeries<ScalarType>
148 : m_transitions; ///< TimeSeries containing points of time and the corresponding number of transitions for every
149 : // AgeGroup.
150 : TimeSeries<ScalarType> m_populations; ///< TimeSeries containing points of time and the corresponding number of
151 : // people in defined #InfectionState%s for every AgeGroup.
152 : CustomIndexArray<ScalarType, AgeGroup>
153 : m_total_confirmed_cases; ///< CustomIndexArray that contains the total number of confirmed cases at time t0 for every AgeGroup.
154 :
155 : private:
156 : // ---- Functionality to calculate the sizes of the compartments for time t0. ----
157 : /**
158 : * @brief Compute the compartment specified in idx_InfectionState at the current time -- only using historic flow
159 : * values and disrespecting potential, previous compartment value.
160 : *
161 : * The computation is meaningful for all compartments except Susceptible, Recovered and #Death
162 : * and mostly needed for initialization.
163 : * For Susceptible, Recovered and Dead, use corresponding alternative functions.
164 : *
165 : * @param[in] dt Time discretization step size.
166 : * @param[in] idx_InfectionState Specifies the considered #InfectionState
167 : * @param[in] group The AgeGroup for which we want to compute.
168 : * @param[in] idx_IncomingFlow Specifies the index of the incoming flow to #InfectionState in m_transitions.
169 : * @param[in] idx_TransitionDistribution1 Specifies the index of the first relevant TransitionDistribution,
170 : * related to a flow from the considered #InfectionState to any other #InfectionState.
171 : * This index is also used for related probability.
172 : * @param[in] idx_TransitionDistribution2 Specifies the index of the second relevant TransitionDistribution,
173 : * related to a flow from the considered #InfectionState to any other #InfectionState (in most cases
174 : * to Recovered).
175 : * Related probability is calculated via 1-probability[idx_TransitionDistribution1].
176 : * Sometimes the second index is not needed, e.g., if probability[idx_TransitionDistribution1]=1.
177 : */
178 :
179 : void compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState, AgeGroup group,
180 : Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1,
181 : int idx_TransitionDistribution2 = 0);
182 :
183 : /**
184 : * @brief Computes the values of the infection compartments subset at initialization.
185 : *
186 : * The values for the compartments Exposed, InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and
187 : * InfectedCritical for time t_0 are calculated using the initial data in form of flows.
188 : * Calculated values are stored in m_populations.
189 : *
190 : * @param[in] dt Time discretization step size.
191 : */
192 : void initial_compute_compartments_infection(ScalarType dt);
193 :
194 : /**
195 : * @brief Computes the values of compartments at initialization.
196 : *
197 : * The initialization method is selected automatically based on the different values that need to be set beforehand.
198 : * Infection compartments are always computed through historic flows.
199 : * Initialization methods for Susceptible and Recovered are tested in the following order:
200 : * 1.) If a positive number for the total number of confirmed cases is set, Recovered is set according to that
201 : * value and Susceptible%s are derived.
202 : * 2.) If Susceptible%s are set, Recovered will be derived.
203 : * 3.) If Recovered are set directly, Susceptible%s are derived.
204 : * 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021)
205 : * to set the Susceptible%s.
206 : *
207 : * The function calculate_compartments_initialization() is used in every method for the compartments Exposed,
208 : * InfectedNoSymptoms, InfectedSymptoms, InfectedSevere and InfectedCritical.
209 : * In addition, the force of infection is calculated for start time t_0.
210 : *
211 : * @param[in] dt Time discretization step size.
212 : */
213 : void initial_compute_compartments(ScalarType dt);
214 :
215 : // ---- Functionality for the iterations of a simulation. ----
216 : /**
217 : * @brief Computes number of Susceptible%s for the current last time in m_populations.
218 : *
219 : * Number is computed using previous number of Susceptible%s and the force of infection (also from previous timestep).
220 : * Number is stored at the matching index in m_populations.
221 : * @param[in] dt Time discretization step size.
222 : */
223 : void compute_susceptibles(ScalarType dt);
224 :
225 : /**
226 : * @brief Computes size of a flow.
227 : *
228 : * Computes size of one flow from #InfectionTransition, specified in idx_InfectionTransitions, for the time
229 : * index current_time_index.
230 : *
231 : * @param[in] idx_InfectionTransitions Specifies the considered flow from #InfectionTransition.
232 : * @param[in] idx_IncomingFlow Index of the flow in #InfectionTransition, which goes to the considered starting
233 : * compartment of the flow specified in idx_InfectionTransitions. Size of considered flow is calculated via
234 : * the value of this incoming flow.
235 : * @param[in] dt Time step to compute flow for.
236 : * @param[in] current_time_index The time index the flow should be computed for.
237 : * @param[in] group The Age group for which we want to compute the flow.
238 : */
239 : void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt,
240 : Eigen::Index current_time_index, AgeGroup group);
241 :
242 : /**
243 : * @brief Computes size of a flow for the current last time value in m_transitions.
244 : *
245 : * Computes size of one flow from #InfectionTransition, specified in idx_InfectionTransitions, for the current
246 : * last time value in m_transitions.
247 : *
248 : * @param[in] idx_InfectionTransitions Specifies the considered flow from #InfectionTransition.
249 : * @param[in] idx_IncomingFlow Index of the flow in #InfectionTransition, which goes to the considered starting
250 : * compartment of the flow specified in idx_InfectionTransitions. Size of considered flow is calculated via
251 : * the value of this incoming flow.
252 : * @param[in] dt Time step to compute flow for.
253 : * @param[in] group The Age Group for which we want to compute the flow.
254 : */
255 : void compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt,
256 : AgeGroup group);
257 :
258 : /**
259 : * @brief Sets all required flows for the current last timestep in m_transitions.
260 : *
261 : * New values are stored in m_transitions. Most values are computed via the function compute_flow().
262 : *
263 : * @param[in] dt Time step.
264 : */
265 : void flows_current_timestep(ScalarType dt);
266 :
267 : /**
268 : * @brief Updates the values of one compartment using flows.
269 : *
270 : * New value is stored in m_populations. The value is calculated using the compartment size in the previous
271 : * time step and the related flows of the current time step.
272 : * Therefore the flows of the current time step should be calculated before using this function.
273 : */
274 : void update_compartment_from_flow(InfectionState infectionState,
275 : std::vector<InfectionTransition> const& IncomingFlows,
276 : std::vector<InfectionTransition> const& OutgoingFlows, AgeGroup group);
277 :
278 : /**
279 : * @brief Updates the values of all compartments except Susceptible at initialization.
280 : *
281 : * New values are stored in m_populations. The values are calculated using the compartment size in the previous
282 : * time step and the related flows of the current time step.
283 : * Therefore the flows of the current time step should be calculated before using this function.
284 : *
285 : */
286 : void update_compartments();
287 :
288 : /**
289 : * @brief Computes force of infection for the current last time in m_transitions.
290 : *
291 : * Computed value is stored in m_forceofinfection.
292 : *
293 : * @param[in] dt Time discretization step size.
294 : * @param[in] initialization If true we are in the case of the initialization of the model.
295 : * For this we need forceofinfection at time point t0-dt and not at the current last time
296 : * (given by m_transitions) as in the other time steps.
297 : */
298 : void compute_forceofinfection(ScalarType dt, bool initialization = false);
299 :
300 : // ---- Functionality to set vectors with necessary information regarding TransitionDistributions. ----
301 : /**
302 : * @brief Setter for the vector m_transitiondistributions_support_max that contains the support_max for all
303 : * TransitionDistributions.
304 : *
305 : * This determines how many summands are required when calculating flows, the force of infection or compartments.
306 : *
307 : * @param[in] dt Time step size.
308 : */
309 : void set_transitiondistributions_support_max(ScalarType dt);
310 :
311 : /**
312 : * @brief Setter for the vector m_transitiondistributions_derivative that contains the approximated derivative for
313 : * all TransitionDistributions for all necessary time points.
314 : *
315 : * The derivative is approximated using a backwards difference scheme.
316 : * The number of necessary time points for each TransitionDistribution is determined using
317 : * m_transitiondistributions_support_max.
318 : *
319 : * @param[in] dt Time step size.
320 : */
321 : void set_transitiondistributions_derivative(ScalarType dt);
322 :
323 : /**
324 : * @brief Setter for the vector m_transitiondistributions_in_forceofinfection.
325 : *
326 : * When computing the force of infection, we evaluate the survival functions of the TransitionDistributions
327 : * InfectedNoSymptomsToInfectedSymptoms, InfectedNoSymptomsToRecovered, InfectedSymptomsToInfectedSevere and
328 : * InfectedSymptomsToRecovered, weighted by the corresponding TransitionProbabilities, at the same time points.
329 : * Here, we compute these contributions to the force of infection term and store them in the vector
330 : * m_transitiondistributions_in_forceofinfection so that we can access this vector for all following computations.
331 : *
332 : * @param[in] dt Time step size.
333 : */
334 : void set_transitiondistributions_in_forceofinfection(ScalarType dt);
335 :
336 : // ---- Private parameters. ----
337 : CustomIndexArray<ScalarType, AgeGroup> m_forceofinfection; ///< Force of infection term needed for numerical scheme.
338 : CustomIndexArray<ScalarType, AgeGroup>
339 : m_N; ///< Vector containing the total population size of the considered region for every AgeGroup.
340 : ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the TransitionDistributions.
341 : size_t m_num_agegroups; ///< Number of Age Groups.
342 : int m_initialization_method{0}; ///< Gives the index of the method used for the initialization of the model.
343 : // See also get_initialization_method_compartments() for the number code.
344 : CustomIndexArray<std::vector<ScalarType>, AgeGroup> m_transitiondistributions_support_max; ///< CustomIndexArray
345 : // of Type AgeGroup containing a Vector containing the support_max for all TransitionDistributions.
346 : CustomIndexArray<std::vector<std::vector<ScalarType>>, AgeGroup>
347 : m_transitiondistributions_derivative; ///< CustomIndexArray
348 : //of Type AgeGroup containing a Vector containing the approximated derivative for all TransitionDistributions for
349 : // all necessary time points.
350 : CustomIndexArray<std::vector<std::vector<ScalarType>>, AgeGroup>
351 : m_transitiondistributions_in_forceofinfection; ///< CustomIndexArray
352 : // of Type AgeGroup containing a Vector containing the weighted TransitionDistributions for all necessary time
353 : // points in the force of infection term.
354 :
355 : // ---- Friend classes/functions. ----
356 : // In the Simulation class, the actual simulation is performed which is why it needs access to the here
357 : // defined (and private) functions to solve the model equations.
358 : friend class Simulation;
359 : // In set_initial_flows(), we compute initial flows based on RKI data using the (private) compute_flow() function
360 : // which is why it is defined as a friend function.
361 : friend IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date,
362 : ScalarType scale_confirmed_cases);
363 : };
364 :
365 : } // namespace isecir
366 : } // namespace mio
367 :
368 : #endif // IDESECIR_MODEL_H
|
