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