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