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 : #include "ide_secir/model.h"
21 : #include "ide_secir/parameters.h"
22 : #include "infection_state.h"
23 : #include "memilio/config.h"
24 : #include "memilio/epidemiology/age_group.h"
25 : #include "memilio/utils/custom_index_array.h"
26 : #include "memilio/utils/logging.h"
27 : #include "memilio/math/eigen.h"
28 :
29 : #include "vector"
30 : #include <algorithm>
31 : #include <cstddef>
32 :
33 : namespace mio
34 : {
35 : namespace isecir
36 : {
37 :
38 23 : Model::Model(TimeSeries<ScalarType>&& transitions_init, CustomIndexArray<ScalarType, AgeGroup> N_init,
39 : CustomIndexArray<ScalarType, AgeGroup> deaths_init, size_t num_agegroups,
40 23 : CustomIndexArray<ScalarType, AgeGroup> total_confirmed_cases_init)
41 23 : : parameters{Parameters(AgeGroup(num_agegroups))}
42 23 : , transitions{std::move(transitions_init)}
43 23 : , populations{TimeSeries<ScalarType>(Eigen::Index(InfectionState::Count) * num_agegroups)}
44 23 : , total_confirmed_cases{total_confirmed_cases_init}
45 23 : , m_N{N_init}
46 46 : , m_num_agegroups{num_agegroups}
47 :
48 : {
49 : // Assert that input arguments for the total population have the correct size regarding
50 : // age groups.
51 23 : assert((size_t)m_N.size() == m_num_agegroups);
52 :
53 23 : if (transitions.get_num_time_points() > 0) {
54 : // Add first time point in m_populations according to last time point in m_transitions which is where we start
55 : // the simulation.
56 63 : populations.add_time_point<Eigen::VectorX<ScalarType>>(
57 21 : transitions.get_last_time(),
58 42 : TimeSeries<ScalarType>::Vector::Constant((size_t)InfectionState::Count * m_num_agegroups, 0.));
59 : }
60 : else {
61 : // Initialize populations with zero as the first point of time if no data is provided for the transitions.
62 : // This can happen for example in the case of initialization with real data.
63 4 : populations.add_time_point<Eigen::VectorX<ScalarType>>(
64 4 : 0, TimeSeries<ScalarType>::Vector::Constant((size_t)InfectionState::Count * m_num_agegroups, 0.));
65 : }
66 :
67 : // Set deaths at simulation start time t0.
68 58 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); group++) {
69 35 : int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
70 35 : populations[Eigen::Index(0)][Di] = deaths_init[group];
71 : }
72 23 : }
73 :
74 12 : bool Model::check_constraints(ScalarType dt) const
75 : {
76 :
77 12 : if (!((size_t)transitions.get_num_elements() == (size_t)InfectionTransition::Count * m_num_agegroups)) {
78 1 : log_error("A variable given for model construction is not valid. Number of elements in vector of "
79 : "transitions does not match the required number.");
80 1 : return true;
81 : }
82 :
83 23 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
84 :
85 116 : for (int i = 0; i < (int)InfectionState::Count; i++) {
86 104 : int index = get_state_flat_index(i, group);
87 104 : if (populations[0][index] < 0) {
88 1 : log_error("Initialization failed. Initial values for populations are less than zero.");
89 1 : return true;
90 : }
91 : }
92 : }
93 :
94 : // It may be possible to run the simulation with fewer time points, but this number ensures that it is possible.
95 10 : if (transitions.get_num_time_points() < (Eigen::Index)std::ceil(get_global_support_max(dt) / dt)) {
96 1 : log_error("Initialization failed. Not enough time points for the transitions given before start of "
97 : "simulation.");
98 1 : return true;
99 : }
100 :
101 19 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
102 :
103 97 : for (int i = 0; i < transitions.get_num_time_points(); i++) {
104 948 : for (int j = 0; j < (int)InfectionTransition::Count; j++) {
105 862 : int index = get_transition_flat_index(j, group);
106 862 : if (transitions[i][index] < 0) {
107 1 : log_error(
108 : "Initialization failed. One or more initial value for the transitions is less than zero.");
109 1 : return true;
110 : }
111 : }
112 : }
113 : }
114 8 : if (transitions.get_last_time() != populations.get_last_time()) {
115 1 : log_error("Last time point of TimeSeries for the transitions does not match last time point of "
116 : "TimeSeries for "
117 : "compartments. Both of these time points have to agree for a sensible simulation.");
118 1 : return true;
119 : }
120 :
121 7 : if (populations.get_num_time_points() != 1) {
122 1 : log_error("The TimeSeries for the compartments contains more than one time point. It is unclear how to "
123 : "initialize.");
124 1 : return true;
125 : }
126 :
127 6 : if ((size_t)total_confirmed_cases.size() > 0 && (size_t)total_confirmed_cases.size() != m_num_agegroups) {
128 1 : log_error("Initialization failed. Number of elements in total_confirmed_cases does not match the number "
129 : "of age groups.");
130 1 : return true;
131 : }
132 :
133 5 : if ((size_t)total_confirmed_cases.size() > 0) {
134 3 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
135 2 : if (total_confirmed_cases[group] < 0) {
136 1 : log_error("Initialization failed. One or more value of total_confirmed_cases is less than zero.");
137 1 : return true;
138 : }
139 : }
140 : }
141 :
142 4 : return parameters.check_constraints();
143 : }
144 :
145 : // Note that this function computes the global_support_max via the get_support_max() method and does not make use
146 : // of the vector m_transitiondistributions_support_max. This is because the global_support_max is already used in
147 : // check_constraints and we cannot ensure that the vector has already been computed when checking for constraints
148 : // (which usually happens before setting the initial transitions and simulating).
149 18 : ScalarType Model::get_global_support_max(ScalarType dt) const
150 : {
151 18 : ScalarType global_support_max = 0.;
152 18 : ScalarType global_support_max_new = 0.;
153 58 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
154 80 : global_support_max_new = std::max(
155 40 : {parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::ExposedToInfectedNoSymptoms]
156 40 : .get_support_max(dt, m_tol),
157 : parameters
158 40 : .get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
159 40 : .get_support_max(dt, m_tol),
160 40 : parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedNoSymptomsToRecovered]
161 40 : .get_support_max(dt, m_tol),
162 : parameters
163 40 : .get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSymptomsToInfectedSevere]
164 40 : .get_support_max(dt, m_tol),
165 40 : parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSymptomsToRecovered]
166 40 : .get_support_max(dt, m_tol),
167 : parameters
168 40 : .get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSevereToInfectedCritical]
169 40 : .get_support_max(dt, m_tol),
170 40 : parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedSevereToRecovered]
171 40 : .get_support_max(dt, m_tol),
172 40 : parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedCriticalToDead]
173 40 : .get_support_max(dt, m_tol),
174 40 : parameters.get<TransitionDistributions>()[group][(int)InfectionTransition::InfectedCriticalToRecovered]
175 40 : .get_support_max(dt, m_tol)});
176 40 : if (global_support_max_new > global_support_max) {
177 18 : global_support_max = global_support_max_new;
178 : }
179 : }
180 18 : return global_support_max;
181 : }
182 :
183 : // ---- Functionality to calculate the sizes of the compartments for time t0. ----
184 85 : void Model::compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState, AgeGroup group,
185 : Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1,
186 : int idx_TransitionDistribution2)
187 : {
188 :
189 85 : ScalarType sum = 0;
190 85 : ScalarType calc_time = 0;
191 : // Determine relevant calculation area and corresponding index.
192 85 : if ((1 - parameters.get<TransitionProbabilities>()[group][idx_TransitionDistribution1]) > 0) {
193 59 : calc_time = std::max(m_transitiondistributions_support_max[group][idx_TransitionDistribution1],
194 59 : m_transitiondistributions_support_max[group][idx_TransitionDistribution2]);
195 : }
196 : else {
197 26 : calc_time = m_transitiondistributions_support_max[group][idx_TransitionDistribution1];
198 : }
199 :
200 85 : Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;
201 :
202 85 : Eigen::Index num_time_points = transitions.get_num_time_points();
203 :
204 : // Index referring to transitions.
205 85 : int transition_index = get_transition_flat_index(idx_IncomingFlow, (size_t)group);
206 : // Index referring to populations.
207 85 : int state_index = get_state_flat_index(idx_InfectionState, (size_t)group);
208 :
209 208 : for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) {
210 :
211 123 : ScalarType state_age = (num_time_points - 1 - i) * dt;
212 :
213 123 : sum += (parameters.get<TransitionProbabilities>()[group][idx_TransitionDistribution1] *
214 123 : parameters.get<TransitionDistributions>()[group][idx_TransitionDistribution1].eval(state_age) +
215 123 : (1 - parameters.get<TransitionProbabilities>()[group][idx_TransitionDistribution1]) *
216 123 : parameters.get<TransitionDistributions>()[group][idx_TransitionDistribution2].eval(state_age)) *
217 123 : transitions[i + 1][transition_index];
218 : }
219 :
220 85 : populations.get_last_value()[state_index] = sum;
221 85 : }
222 :
223 15 : void Model::initial_compute_compartments_infection(ScalarType dt)
224 : {
225 : // We need to compute the initial compartments for every AgeGroup.
226 32 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
227 : // Exposed
228 17 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::Exposed), group,
229 : Eigen::Index(InfectionTransition::SusceptibleToExposed),
230 : (int)InfectionTransition::ExposedToInfectedNoSymptoms);
231 : // InfectedNoSymptoms
232 17 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedNoSymptoms), group,
233 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms),
234 : (int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
235 : (int)InfectionTransition::InfectedNoSymptomsToRecovered);
236 : // InfectedSymptoms
237 17 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedSymptoms), group,
238 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms),
239 : (int)InfectionTransition::InfectedSymptomsToInfectedSevere,
240 : (int)InfectionTransition::InfectedSymptomsToRecovered);
241 : // InfectedSevere
242 17 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedSevere), group,
243 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
244 : (int)InfectionTransition::InfectedSevereToInfectedCritical,
245 : (int)InfectionTransition::InfectedSevereToRecovered);
246 : // InfectedCritical
247 17 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedCritical), group,
248 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
249 : (int)InfectionTransition::InfectedCriticalToDead,
250 : (int)InfectionTransition::InfectedCriticalToRecovered);
251 : }
252 15 : }
253 :
254 15 : void Model::initial_compute_compartments(ScalarType dt)
255 : {
256 : // The initialization method only affects the Susceptible and Recovered compartments.
257 : // It is possible to calculate the sizes of the other compartments in advance because only the initial values of
258 : // the transitions are used.
259 15 : initial_compute_compartments_infection(dt);
260 :
261 : // We store in two Booleans if there are Susceptibles or Recovered given for every age group.
262 15 : bool susceptibles_given = true;
263 16 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
264 15 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
265 15 : if (populations[Eigen::Index(0)][Si] < 1e-12) {
266 14 : susceptibles_given = false;
267 14 : break;
268 : }
269 : }
270 15 : bool recovered_given = true;
271 16 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
272 15 : int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
273 15 : if (populations[Eigen::Index(0)][Ri] < 1e-12) {
274 14 : recovered_given = false;
275 14 : break;
276 : }
277 : }
278 :
279 : // We check which Initialization method we want to use.
280 21 : if (!(total_confirmed_cases == CustomIndexArray<ScalarType, AgeGroup>()) &&
281 6 : std::all_of(total_confirmed_cases.begin(), total_confirmed_cases.end(), [](ScalarType x) {
282 6 : return x > 1e-12;
283 : })) {
284 1 : m_initialization_method = 1;
285 2 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
286 1 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
287 1 : int Ei = get_state_flat_index(Eigen::Index(InfectionState::Exposed), group);
288 1 : int INSi = get_state_flat_index(Eigen::Index(InfectionState::InfectedNoSymptoms), group);
289 1 : int ISyi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSymptoms), group);
290 1 : int ISevi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSevere), group);
291 1 : int ICri = get_state_flat_index(Eigen::Index(InfectionState::InfectedCritical), group);
292 1 : int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
293 1 : int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
294 : // The scheme of the ODE model for initialization is applied here.
295 4 : populations[Eigen::Index(0)][Ri] = total_confirmed_cases[group] - populations[Eigen::Index(0)][ISyi] -
296 3 : populations[Eigen::Index(0)][ISevi] -
297 3 : populations[Eigen::Index(0)][ICri] -
298 2 : populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
299 :
300 4 : populations[Eigen::Index(0)][Si] = m_N[group] - populations[Eigen::Index(0)][Ei] -
301 4 : populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
302 3 : populations[Eigen::Index(0)][ISevi] -
303 4 : populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Ri] -
304 2 : populations[Eigen::Index(0)][Di];
305 : }
306 : }
307 :
308 14 : else if (susceptibles_given) {
309 : // Take initialized value for Susceptibles if value can't be calculated via the standard formula.
310 1 : m_initialization_method = 2;
311 : // R; need an initial value for R, therefore do not calculate via compute_recovered()
312 2 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
313 1 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
314 1 : int Ei = get_state_flat_index(Eigen::Index(InfectionState::Exposed), group);
315 1 : int INSi = get_state_flat_index(Eigen::Index(InfectionState::InfectedNoSymptoms), group);
316 1 : int ISyi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSymptoms), group);
317 1 : int ISevi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSevere), group);
318 1 : int ICri = get_state_flat_index(Eigen::Index(InfectionState::InfectedCritical), group);
319 1 : int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
320 1 : int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
321 4 : populations[Eigen::Index(0)][Ri] = m_N[group] - populations[Eigen::Index(0)][Si] -
322 4 : populations[Eigen::Index(0)][Ei] - populations[Eigen::Index(0)][INSi] -
323 3 : populations[Eigen::Index(0)][ISyi] -
324 3 : populations[Eigen::Index(0)][ISevi] -
325 3 : populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Di];
326 : }
327 : }
328 13 : else if (recovered_given) {
329 : // If value for Recovered is initialized and standard method is not applicable (i.e., no value for total infections
330 : // or Susceptibles given directly), calculate Susceptibles via other compartments.
331 1 : m_initialization_method = 3;
332 2 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
333 1 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
334 1 : int Ei = get_state_flat_index(Eigen::Index(InfectionState::Exposed), group);
335 1 : int INSi = get_state_flat_index(Eigen::Index(InfectionState::InfectedNoSymptoms), group);
336 1 : int ISyi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSymptoms), group);
337 1 : int ISevi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSevere), group);
338 1 : int ICri = get_state_flat_index(Eigen::Index(InfectionState::InfectedCritical), group);
339 1 : int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
340 1 : int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
341 4 : populations[Eigen::Index(0)][Si] = m_N[group] - populations[Eigen::Index(0)][Ei] -
342 4 : populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
343 3 : populations[Eigen::Index(0)][ISevi] -
344 4 : populations[Eigen::Index(0)][ICri] - populations[Eigen::Index(0)][Ri] -
345 2 : populations[Eigen::Index(0)][Di];
346 : }
347 : }
348 : else {
349 : // Compute Susceptibles at t0 and m_forceofinfection at time t0-dt as initial values for discretization scheme.
350 : // Use m_forceofinfection at t0-dt to be consistent with further calculations of S (see compute_susceptibles()),
351 : // where also the value of m_forceofinfection for the previous timestep is used.
352 12 : compute_forceofinfection(dt, true);
353 12 : if (std::all_of(m_forceofinfection.begin(), m_forceofinfection.end(), [](ScalarType x) {
354 14 : return x > 1e-12;
355 : })) {
356 10 : m_initialization_method = 4;
357 22 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
358 12 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
359 12 : int Ei = get_state_flat_index(Eigen::Index(InfectionState::Exposed), group);
360 12 : int INSi = get_state_flat_index(Eigen::Index(InfectionState::InfectedNoSymptoms), group);
361 12 : int ISyi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSymptoms), group);
362 12 : int ISevi = get_state_flat_index(Eigen::Index(InfectionState::InfectedSevere), group);
363 12 : int ICri = get_state_flat_index(Eigen::Index(InfectionState::InfectedCritical), group);
364 12 : int Ri = get_state_flat_index(Eigen::Index(InfectionState::Recovered), group);
365 12 : int Di = get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
366 :
367 12 : int StEi = get_transition_flat_index(Eigen::Index(InfectionTransition::SusceptibleToExposed), group);
368 : /* Attention: With an inappropriate combination of parameters and initial conditions, it can happen that S
369 : becomes greater than N when this formula is used. In this case, at least one compartment is initialized
370 : with a number less than zero, so that a log_error is thrown.
371 : However, this initialization method is consistent with the numerical solver of the model equations,
372 : so it may sometimes make sense to use this method. */
373 24 : populations[Eigen::Index(0)][Si] =
374 24 : transitions.get_last_value()[StEi] / (dt * m_forceofinfection[group]);
375 :
376 : // Recovered; calculated as the difference between the total population and the sum of the other
377 : // compartment sizes.
378 12 : populations[Eigen::Index(0)][Ri] =
379 36 : m_N[group] - populations[Eigen::Index(0)][Si] - populations[Eigen::Index(0)][Ei] -
380 48 : populations[Eigen::Index(0)][INSi] - populations[Eigen::Index(0)][ISyi] -
381 48 : populations[Eigen::Index(0)][ISevi] - populations[Eigen::Index(0)][ICri] -
382 24 : populations[Eigen::Index(0)][Di];
383 : }
384 : }
385 : else {
386 2 : m_initialization_method = -1;
387 2 : log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither "
388 : "Susceptibles nor Recovered or total confirmed cases for time 0 were set. One of them should be "
389 : "larger 0.");
390 : }
391 : }
392 : // Verify that the initialization produces an appropriate result.
393 : // Another check would be if the sum of the compartments is equal to N, but in all initialization methods one of
394 : // the compartments is initialized via N - the others.
395 : // This also means that if a compartment is greater than N, we will always have one or more compartments less than
396 : // zero.
397 : // Check if all compartments are non negative.
398 32 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
399 153 : for (Eigen::Index i = 0; i < (Eigen::Index)InfectionState::Count; i++) {
400 136 : int idx = get_state_flat_index(i, group);
401 136 : if (populations[0][idx] < 0) {
402 1 : m_initialization_method = -2;
403 1 : log_error(
404 : "Initialization failed. One or more initial values for populations are less than zero. It may "
405 : "help to use a different method for initialization.");
406 : }
407 : }
408 : }
409 :
410 : // Compute m_forceofinfection at time t0 needed for further simulation.
411 15 : compute_forceofinfection(dt);
412 15 : }
413 :
414 : // ---- Functionality for the iterations of a simulation. ----
415 108 : void Model::compute_susceptibles(ScalarType dt)
416 : {
417 : // We need to compute to compute the Susceptibles in every AgeGroup.
418 226 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
419 118 : Eigen::Index num_time_points = populations.get_num_time_points();
420 : // Using number of Susceptibles from previous time step and force of infection from previous time step:
421 : // Compute current number of Susceptibles and store Susceptibles in populations.
422 118 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
423 118 : populations.get_last_value()[Si] = populations[num_time_points - 2][Si] / (1 + dt * m_forceofinfection[group]);
424 : }
425 108 : }
426 :
427 2049 : void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt,
428 : Eigen::Index current_time_index, AgeGroup group)
429 : {
430 2049 : ScalarType sum = 0;
431 : /* In order to satisfy TransitionDistribution(dt*i) = 0 for all i >= k, k is determined by the maximum support of
432 : the distribution.
433 : Since we are using a backwards difference scheme to compute the derivative, we have that the
434 : derivative of TransitionDistribution(dt*i) = 0 for all i >= k+1.
435 :
436 : Hence calc_time_index goes until std::ceil(support_max/dt) since for std::ceil(support_max/dt)+1 all terms are
437 : already zero.
438 : This needs to be adjusted if we are changing the finite difference scheme */
439 : Eigen::Index calc_time_index =
440 2049 : (Eigen::Index)std::ceil(m_transitiondistributions_support_max[group][idx_InfectionTransitions] / dt);
441 :
442 2049 : int transition_idx = get_transition_flat_index(idx_InfectionTransitions, size_t(group));
443 9451 : for (Eigen::Index i = current_time_index - calc_time_index; i < current_time_index; i++) {
444 : // (current_time_index - i) is the index corresponding to time the individuals have already spent in this state.
445 7402 : sum += m_transitiondistributions_derivative[group][idx_InfectionTransitions][current_time_index - i] *
446 7402 : transitions[i + 1][idx_IncomingFlow];
447 : }
448 :
449 2049 : transitions.get_value(current_time_index)[transition_idx] =
450 2049 : (-dt) * parameters.get<TransitionProbabilities>()[group][idx_InfectionTransitions] * sum;
451 2049 : }
452 :
453 1062 : void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt,
454 : AgeGroup group)
455 : {
456 1062 : Eigen::Index current_time_index = transitions.get_num_time_points() - 1;
457 1062 : compute_flow(idx_InfectionTransitions, idx_IncomingFlow, dt, current_time_index, group);
458 1062 : }
459 :
460 108 : void Model::flows_current_timestep(ScalarType dt)
461 : {
462 : // Compute flows for every AgeGroup.
463 226 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
464 118 : int StEi = get_transition_flat_index(Eigen::Index(InfectionTransition::SusceptibleToExposed), group);
465 118 : int Si = get_state_flat_index(Eigen::Index(InfectionState::Susceptible), group);
466 :
467 : // Calculate flow SusceptibleToExposed with force of infection from previous time step and Susceptibles from current time step.
468 118 : transitions.get_last_value()[StEi] = dt * m_forceofinfection[group] * populations.get_last_value()[Si];
469 :
470 : // Calculate all other flows with compute_flow.
471 : // Flow from Exposed to InfectedNoSymptoms
472 118 : compute_flow(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms),
473 : Eigen::Index(InfectionTransition::SusceptibleToExposed), dt, group);
474 : // Flow from InfectedNoSymptoms to InfectedSymptoms
475 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms),
476 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, group);
477 : // Flow from InfectedNoSymptoms to Recovered
478 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered),
479 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, group);
480 : // Flow from InfectedSymptoms to InfectedSevere
481 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
482 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt, group);
483 : // Flow from InfectedSymptoms to Recovered
484 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered),
485 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt, group);
486 : // Flow from InfectedSevere to InfectedCritical
487 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
488 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, group);
489 : // Flow from InfectedSevere to Recovered
490 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered),
491 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, group);
492 : // Flow from InfectedCritical to Dead
493 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead),
494 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, group);
495 : // Flow from InfectedCritical to Recovered
496 118 : compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered),
497 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, group);
498 : }
499 108 : }
500 :
501 826 : void Model::update_compartment_from_flow(InfectionState infectionState,
502 : std::vector<InfectionTransition> const& IncomingFlows,
503 : std::vector<InfectionTransition> const& OutgoingFlows, AgeGroup group)
504 : {
505 826 : int state_idx = get_state_flat_index(Eigen::Index(infectionState), group);
506 :
507 826 : Eigen::Index num_time_points = populations.get_num_time_points();
508 826 : ScalarType updated_compartment = populations[num_time_points - 2][state_idx];
509 2006 : for (const InfectionTransition& inflow : IncomingFlows) {
510 1180 : int inflow_idx = get_transition_flat_index(Eigen::Index(inflow), group);
511 1180 : updated_compartment += transitions.get_last_value()[inflow_idx];
512 : }
513 1888 : for (const InfectionTransition& outflow : OutgoingFlows) {
514 1062 : int outflow_idx = get_transition_flat_index(Eigen::Index(outflow), group);
515 1062 : updated_compartment -= transitions.get_last_value()[outflow_idx];
516 : }
517 826 : populations.get_last_value()[state_idx] = updated_compartment;
518 826 : }
519 :
520 108 : void Model::update_compartments()
521 : {
522 : // Update compartments for every AgeGroup.
523 226 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
524 : // Exposed
525 118 : update_compartment_from_flow(InfectionState::Exposed, {InfectionTransition::SusceptibleToExposed},
526 : {InfectionTransition::ExposedToInfectedNoSymptoms}, group);
527 :
528 : // InfectedNoSymptoms
529 118 : update_compartment_from_flow(InfectionState::InfectedNoSymptoms,
530 : {InfectionTransition::ExposedToInfectedNoSymptoms},
531 : {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
532 : InfectionTransition::InfectedNoSymptomsToRecovered},
533 : group);
534 :
535 : // InfectedSymptoms
536 118 : update_compartment_from_flow(
537 : InfectionState::InfectedSymptoms, {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms},
538 : {InfectionTransition::InfectedSymptomsToInfectedSevere, InfectionTransition::InfectedSymptomsToRecovered},
539 : group);
540 :
541 : // InfectedSevere
542 118 : update_compartment_from_flow(
543 : InfectionState::InfectedSevere, {InfectionTransition::InfectedSymptomsToInfectedSevere},
544 : {InfectionTransition::InfectedSevereToInfectedCritical, InfectionTransition::InfectedSevereToRecovered},
545 : group);
546 :
547 : // InfectedCritical
548 118 : update_compartment_from_flow(
549 : InfectionState::InfectedCritical, {InfectionTransition::InfectedSevereToInfectedCritical},
550 : {InfectionTransition::InfectedCriticalToDead, InfectionTransition::InfectedCriticalToRecovered}, group);
551 :
552 : // Recovered
553 354 : update_compartment_from_flow(InfectionState::Recovered,
554 : {
555 : InfectionTransition::InfectedNoSymptomsToRecovered,
556 : InfectionTransition::InfectedSymptomsToRecovered,
557 : InfectionTransition::InfectedSevereToRecovered,
558 : InfectionTransition::InfectedCriticalToRecovered,
559 : },
560 236 : std::vector<InfectionTransition>(), group);
561 :
562 : // Dead
563 354 : update_compartment_from_flow(InfectionState::Dead, {InfectionTransition::InfectedCriticalToDead},
564 236 : std::vector<InfectionTransition>(), group);
565 : }
566 108 : }
567 :
568 135 : void Model::compute_forceofinfection(ScalarType dt, bool initialization)
569 : {
570 :
571 135 : m_forceofinfection = CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(m_num_agegroups), 0.);
572 : // We compute the force of infection for every AgeGroup.
573 :
574 284 : for (AgeGroup i = AgeGroup(0); i < AgeGroup(m_num_agegroups); ++i) {
575 :
576 : Eigen::Index num_time_points;
577 : ScalarType current_time;
578 :
579 149 : if (initialization) {
580 : // Determine m_forceofinfection at time t0-dt which is the penultimate timepoint in transition.
581 14 : num_time_points = transitions.get_num_time_points() - 1;
582 : // Get time of penultimate timepoint in transitions.
583 14 : current_time = transitions.get_time(num_time_points - 1);
584 : }
585 : else {
586 : // Determine m_forceofinfection for current last time in transitions.
587 135 : num_time_points = transitions.get_num_time_points();
588 135 : current_time = transitions.get_last_time();
589 : }
590 : //We compute the Season Value.
591 : ScalarType season_val =
592 : 1 +
593 149 : parameters.get<Seasonality>() *
594 149 : sin(3.141592653589793 * (std::fmod((parameters.get<StartDay>() + current_time), 365.0) / 182.5 + 0.5));
595 : // To include contacts between all age groups we sum over all age groups.
596 340 : for (AgeGroup j = AgeGroup(0); j < AgeGroup(m_num_agegroups); ++j) {
597 : // Determine the relevant calculation area = union of the supports of the relevant transition distributions.
598 955 : ScalarType calc_time = std::max(
599 191 : {m_transitiondistributions_support_max[j]
600 191 : [(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms],
601 191 : m_transitiondistributions_support_max[j][(int)InfectionTransition::InfectedNoSymptomsToRecovered],
602 191 : m_transitiondistributions_support_max[j][(int)InfectionTransition::InfectedSymptomsToInfectedSevere],
603 191 : m_transitiondistributions_support_max[j][(int)InfectionTransition::InfectedSymptomsToRecovered]});
604 : // Corresponding index.
605 : // Need calc_time_index timesteps in sum,
606 : // subtract 1 because in the last summand all TransitionDistributions evaluate to 0 (by definition of support_max).
607 191 : Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;
608 :
609 191 : int Dj = get_state_flat_index(Eigen::Index(InfectionState::Dead), j);
610 191 : int ICrtDj = get_transition_flat_index(Eigen::Index(InfectionTransition::InfectedCriticalToDead), j);
611 :
612 191 : int EtINSj = get_transition_flat_index(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), j);
613 : int INStISyj =
614 191 : get_transition_flat_index(Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), j);
615 :
616 : // We store the number of deaths for every AgeGroup.
617 : ScalarType deaths_j;
618 191 : if (initialization) {
619 : // Determine the number of individuals in Dead compartment at time t0-dt.
620 20 : deaths_j = populations[Eigen::Index(0)][Dj] - transitions.get_last_value()[ICrtDj];
621 : }
622 : else {
623 171 : deaths_j = populations.get_last_value()[Dj];
624 : }
625 :
626 191 : ScalarType sum = 0;
627 :
628 : // Sum over all relevant time points.
629 573 : for (Eigen::Index l = num_time_points - 1 - calc_time_index; l < num_time_points - 1; l++) {
630 382 : Eigen::Index state_age_index = num_time_points - 1 - l;
631 382 : ScalarType state_age = state_age_index * dt;
632 382 : sum += season_val * parameters.get<TransmissionProbabilityOnContact>()[i].eval(state_age) *
633 1146 : parameters.get<ContactPatterns>().get_cont_freq_mat().get_matrix_at(current_time)(
634 764 : static_cast<Eigen::Index>((size_t)i), static_cast<Eigen::Index>((size_t)j)) *
635 382 : (m_transitiondistributions_in_forceofinfection[j][0][num_time_points - l - 1] *
636 764 : transitions[l + 1][EtINSj] *
637 382 : parameters.get<RelativeTransmissionNoSymptoms>()[j].eval(state_age) +
638 382 : m_transitiondistributions_in_forceofinfection[j][1][num_time_points - l - 1] *
639 382 : transitions[l + 1][INStISyj] *
640 382 : parameters.get<RiskOfInfectionFromSymptomatic>()[j].eval(state_age));
641 : }
642 : const double divNj =
643 191 : (m_N[j] - deaths_j < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / (m_N[j] - deaths_j);
644 191 : m_forceofinfection[i] += divNj * sum;
645 : }
646 : }
647 135 : }
648 :
649 : // ---- Functionality to set vectors with necessary information regarding TransitionDistributions. ----
650 21 : void Model::set_transitiondistributions_support_max(ScalarType dt)
651 : {
652 84 : m_transitiondistributions_support_max = CustomIndexArray<std::vector<ScalarType>, AgeGroup>(
653 84 : AgeGroup(m_num_agegroups), std::vector<ScalarType>((int)InfectionTransition::Count, 0.));
654 : // We do not consider the transition SusceptibleToExposed as it is not needed in the computations.
655 59 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
656 380 : for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) {
657 342 : m_transitiondistributions_support_max[group][transition] =
658 684 : parameters.get<TransitionDistributions>()[AgeGroup(group)][transition].get_support_max(dt, m_tol);
659 : }
660 : }
661 21 : }
662 :
663 21 : void Model::set_transitiondistributions_derivative(ScalarType dt)
664 : {
665 84 : m_transitiondistributions_derivative = CustomIndexArray<std::vector<std::vector<ScalarType>>, AgeGroup>(
666 42 : AgeGroup(m_num_agegroups),
667 63 : std::vector<std::vector<ScalarType>>((int)InfectionTransition::Count, std::vector<ScalarType>(1, 0.)));
668 : // We do not consider the transition SusceptibleToExposed as it is not needed in the computations.
669 59 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
670 380 : for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) {
671 : Eigen::Index support_max_index =
672 342 : (Eigen::Index)std::ceil(m_transitiondistributions_support_max[group][transition] / dt);
673 : // Create vec_derivative that contains the value of the approximated derivative for all necessary time points.
674 : // Here, we evaluate the derivative at time points t_0, ..., t_{support_max_index}.
675 684 : std::vector<ScalarType> vec_derivative(support_max_index + 1, 0.);
676 :
677 1807 : for (Eigen::Index i = 0; i <= support_max_index; i++) {
678 : // Compute state_age for considered index.
679 1465 : ScalarType state_age = (ScalarType)i * dt;
680 : // Compute derivative.
681 1465 : vec_derivative[i] =
682 1465 : (parameters.get<TransitionDistributions>()[group][transition].eval(state_age) -
683 1465 : parameters.get<TransitionDistributions>()[group][transition].eval(state_age - dt)) /
684 1465 : dt;
685 : }
686 342 : m_transitiondistributions_derivative[group][transition] = vec_derivative;
687 342 : }
688 : }
689 21 : }
690 :
691 15 : void Model::set_transitiondistributions_in_forceofinfection(ScalarType dt)
692 : {
693 60 : m_transitiondistributions_in_forceofinfection = CustomIndexArray<std::vector<std::vector<ScalarType>>, AgeGroup>(
694 60 : AgeGroup(m_num_agegroups), std::vector<std::vector<ScalarType>>(2, std::vector<ScalarType>(1, 0.)));
695 : // Relevant transitions for force of infection term.
696 15 : std::vector<std::vector<int>> relevant_transitions = {
697 : {(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
698 : (int)InfectionTransition::InfectedNoSymptomsToRecovered},
699 : {(int)InfectionTransition::InfectedSymptomsToInfectedSevere,
700 105 : (int)InfectionTransition::InfectedSymptomsToRecovered}};
701 :
702 : // Determine the relevant calculation area = union of the supports of the relevant transition distributions.
703 :
704 32 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(m_num_agegroups); ++group) {
705 68 : ScalarType calc_time = std::max({m_transitiondistributions_support_max[group][relevant_transitions[0][0]],
706 17 : m_transitiondistributions_support_max[group][relevant_transitions[0][1]],
707 17 : m_transitiondistributions_support_max[group][relevant_transitions[1][0]],
708 17 : m_transitiondistributions_support_max[group][relevant_transitions[1][1]]});
709 : // Corresponding index.
710 : // Need to evaluate survival functions at t_0, ..., t_{calc_time_index} for computation of force of infection,
711 : // subtract 1 because in the last summand all TransitionDistributions evaluate to 0 (by definition of support_max).
712 17 : Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;
713 : // Compute contributions from survival functions and transition probabilities starting in InfectedNoSymptoms and
714 : // InfectedSymptoms, respectively, on force of infection term.
715 51 : for (int contribution = 0; contribution < 2; contribution++) {
716 68 : std::vector<ScalarType> vec_contribution_to_foi(calc_time_index + 1, 0.);
717 116 : for (Eigen::Index i = 0; i <= calc_time_index; i++) {
718 : // Compute state_age for all indices from t_0, ..., t_{calc_time_index}.
719 82 : ScalarType state_age = i * dt;
720 82 : vec_contribution_to_foi[i] =
721 82 : parameters.get<TransitionProbabilities>()[group][relevant_transitions[contribution][0]] *
722 82 : parameters.get<TransitionDistributions>()[group][relevant_transitions[contribution][0]].eval(
723 82 : state_age) +
724 82 : parameters.get<TransitionProbabilities>()[group][relevant_transitions[contribution][1]] *
725 82 : parameters.get<TransitionDistributions>()[group][relevant_transitions[contribution][1]].eval(
726 : state_age);
727 : }
728 :
729 34 : m_transitiondistributions_in_forceofinfection[group][contribution] = vec_contribution_to_foi;
730 34 : }
731 : }
732 30 : }
733 :
734 : } // namespace isecir
735 : } // namespace mio
|
