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