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 : #include "ide_secir/model.h"
21 : #include "ide_secir/parameters.h"
22 : #include "infection_state.h"
23 : #include "memilio/config.h"
24 : #include "memilio/utils/logging.h"
25 : #include "memilio/math/eigen.h"
26 :
27 : #include "vector"
28 :
29 : namespace mio
30 : {
31 : namespace isecir
32 : {
33 :
34 20 : Model::Model(TimeSeries<ScalarType>&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases,
35 20 : const ParameterSet& Parameterset_init)
36 20 : : parameters{Parameterset_init}
37 20 : , m_transitions{std::move(init)}
38 20 : , m_populations{TimeSeries<ScalarType>(Eigen::Index(InfectionState::Count))}
39 20 : , m_total_confirmed_cases{total_confirmed_cases}
40 20 : , m_N{N_init}
41 : {
42 20 : if (m_transitions.get_num_time_points() > 0) {
43 : // Add first time point in m_populations according to last time point in m_transitions which is where we start
44 : // the simulation.
45 57 : m_populations.add_time_point<Eigen::VectorXd>(
46 57 : m_transitions.get_last_time(), TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count, 0.));
47 : }
48 : else {
49 : // Initialize m_populations with zero as the first point of time if no data is provided for the transitions.
50 : // This can happen for example in the case of initialization with real data.
51 2 : m_populations.add_time_point<Eigen::VectorXd>(
52 2 : 0, TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count, 0));
53 : }
54 :
55 : // Set deaths at simulation start time t0.
56 20 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] = deaths;
57 20 : }
58 :
59 : // ---- Functionality to calculate the sizes of the compartments for time t0. ----
60 70 : void Model::compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_InfectionState,
61 : Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1,
62 : int idx_TransitionDistribution2)
63 : {
64 70 : ScalarType sum = 0;
65 70 : ScalarType calc_time = 0;
66 : // Determine relevant calculation area and corresponding index.
67 70 : if ((1 - parameters.get<TransitionProbabilities>()[idx_TransitionDistribution1]) > 0) {
68 47 : calc_time = std::max(m_transitiondistributions_support_max[idx_TransitionDistribution1],
69 47 : m_transitiondistributions_support_max[idx_TransitionDistribution2]);
70 : }
71 : else {
72 23 : calc_time = m_transitiondistributions_support_max[idx_TransitionDistribution1];
73 : }
74 :
75 70 : Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;
76 :
77 70 : Eigen::Index num_time_points = m_transitions.get_num_time_points();
78 :
79 167 : for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) {
80 :
81 97 : ScalarType state_age = (num_time_points - 1 - i) * dt;
82 :
83 97 : sum += (parameters.get<TransitionProbabilities>()[idx_TransitionDistribution1] *
84 97 : parameters.get<TransitionDistributions>()[idx_TransitionDistribution1].eval(state_age) +
85 97 : (1 - parameters.get<TransitionProbabilities>()[idx_TransitionDistribution1]) *
86 97 : parameters.get<TransitionDistributions>()[idx_TransitionDistribution2].eval(state_age)) *
87 97 : m_transitions[i + 1][idx_IncomingFlow];
88 : }
89 :
90 70 : m_populations.get_last_value()[idx_InfectionState] = sum;
91 70 : }
92 :
93 14 : void Model::initial_compute_compartments_infection(ScalarType dt)
94 : {
95 : // Exposed
96 14 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::Exposed),
97 : Eigen::Index(InfectionTransition::SusceptibleToExposed),
98 : (int)InfectionTransition::ExposedToInfectedNoSymptoms);
99 : // InfectedNoSymptoms
100 14 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedNoSymptoms),
101 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms),
102 : (int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
103 : (int)InfectionTransition::InfectedNoSymptomsToRecovered);
104 : // InfectedSymptoms
105 14 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedSymptoms),
106 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms),
107 : (int)InfectionTransition::InfectedSymptomsToInfectedSevere,
108 : (int)InfectionTransition::InfectedSymptomsToRecovered);
109 : // InfectedSevere
110 14 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedSevere),
111 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
112 : (int)InfectionTransition::InfectedSevereToInfectedCritical,
113 : (int)InfectionTransition::InfectedSevereToRecovered);
114 : // InfectedCritical
115 14 : compute_compartment_from_flows(dt, Eigen::Index(InfectionState::InfectedCritical),
116 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
117 : (int)InfectionTransition::InfectedCriticalToDead,
118 : (int)InfectionTransition::InfectedCriticalToRecovered);
119 14 : }
120 :
121 14 : void Model::initial_compute_compartments(ScalarType dt)
122 : {
123 : // The initialization method only affects the Susceptible and Recovered compartments.
124 : // It is possible to calculate the sizes of the other compartments in advance because only the initial values of
125 : // the flows are used.
126 14 : initial_compute_compartments_infection(dt);
127 :
128 14 : if (m_total_confirmed_cases > 1e-12) {
129 1 : m_initialization_method = 1;
130 :
131 : // The scheme of the ODE model for initialization is applied here.
132 1 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
133 3 : m_total_confirmed_cases - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
134 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
135 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
136 2 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
137 :
138 1 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
139 3 : m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
140 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
141 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
142 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
143 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
144 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
145 2 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
146 : }
147 13 : else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) {
148 : // Take initialized value for Susceptibles if value can't be calculated via the standard formula.
149 1 : m_initialization_method = 2;
150 :
151 : // R; need an initial value for R, therefore do not calculate via compute_recovered()
152 1 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
153 3 : m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
154 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
155 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
156 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
157 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
158 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
159 2 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
160 : }
161 12 : else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) {
162 : // If value for Recovered is initialized and standard method is not applicable (i.e., no value for total
163 : // infections or Susceptibles given directly), calculate Susceptibles via other compartments.
164 1 : m_initialization_method = 3;
165 :
166 1 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
167 3 : m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
168 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
169 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
170 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
171 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
172 3 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
173 2 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
174 : }
175 : else {
176 : // Compute Susceptibles at t0 and m_forceofinfection at time t0-dt as initial values for discretization scheme.
177 : // Use m_forceofinfection at t0-dt to be consistent with further calculations of S (see compute_susceptibles()),
178 : // where also the value of m_forceofinfection for the previous timestep is used.
179 11 : compute_forceofinfection(dt, true);
180 11 : if (m_forceofinfection > 1e-12) {
181 9 : m_initialization_method = 4;
182 :
183 : /* Attention: With an inappropriate combination of parameters and initial conditions, it can happen that S
184 : becomes greater than N when this formula is used. In this case, at least one compartment is initialized
185 : with a number less than zero, so that a log_error is thrown.
186 : However, this initialization method is consistent with the numerical solver of the model equations,
187 : so it may sometimes make sense to use this method. */
188 18 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
189 18 : m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] /
190 9 : (dt * m_forceofinfection);
191 :
192 : // Recovered; calculated as the difference between the total population and the sum of the other
193 : // compartment sizes.
194 9 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
195 27 : m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
196 27 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
197 27 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
198 27 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
199 27 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
200 27 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
201 18 : m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
202 : }
203 : else {
204 2 : m_initialization_method = -1;
205 4 : log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither "
206 : "Susceptibles nor Recovered or total confirmed cases for time 0 were set. One of them should be "
207 : "larger 0.");
208 : }
209 : }
210 :
211 : // Verify that the initialization produces an appropriate result.
212 : // Another check would be if the sum of the compartments is equal to N, but in all initialization methods one of
213 : // the compartments is initialized via N - the others.
214 : // This also means that if a compartment is greater than N, we will always have one or more compartments less than
215 : // zero.
216 : // Check if all compartments are non negative.
217 126 : for (int i = 0; i < (int)InfectionState::Count; i++) {
218 112 : if (m_populations[0][i] < 0) {
219 1 : m_initialization_method = -2;
220 2 : log_error("Initialization failed. One or more initial values for populations are less than zero. It may "
221 : "help to use a different method for initialization.");
222 : }
223 : }
224 :
225 : // Compute m_forceofinfection at time t0 needed for further simulation.
226 14 : compute_forceofinfection(dt);
227 14 : }
228 :
229 : // ---- Functionality for the iterations of a simulation. ----
230 103 : void Model::compute_susceptibles(ScalarType dt)
231 : {
232 103 : Eigen::Index num_time_points = m_populations.get_num_time_points();
233 : // Using number of Susceptibles from previous time step and force of infection from previous time step:
234 : // Compute current number of Susceptibles and store Susceptibles in m_populations.
235 206 : m_populations.get_last_value()[Eigen::Index(InfectionState::Susceptible)] =
236 206 : m_populations[num_time_points - 2][Eigen::Index(InfectionState::Susceptible)] / (1 + dt * m_forceofinfection);
237 103 : }
238 :
239 1021 : void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt,
240 : Eigen::Index current_time_index)
241 : {
242 1021 : ScalarType sum = 0;
243 : /* In order to satisfy TransitionDistribution(dt*i) = 0 for all i >= k, k is determined by the maximum support of
244 : the distribution.
245 : Since we are using a backwards difference scheme to compute the derivative, we have that the
246 : derivative of TransitionDistribution(dt*i) = 0 for all i >= k+1.
247 :
248 : Hence calc_time_index goes until std::ceil(support_max/dt) since for std::ceil(support_max/dt)+1 all terms are
249 : already zero.
250 : This needs to be adjusted if we are changing the finite difference scheme */
251 :
252 : Eigen::Index calc_time_index =
253 1021 : (Eigen::Index)std::ceil(m_transitiondistributions_support_max[idx_InfectionTransitions] / dt);
254 :
255 4491 : for (Eigen::Index i = current_time_index - calc_time_index; i < current_time_index; i++) {
256 : // (current_time_index - i) is the index corresponding to time the individuals have already spent in this state.
257 3470 : sum += m_transitiondistributions_derivative[idx_InfectionTransitions][current_time_index - i] *
258 3470 : m_transitions[i + 1][idx_IncomingFlow];
259 : }
260 :
261 1021 : m_transitions.get_value(current_time_index)[Eigen::Index(idx_InfectionTransitions)] =
262 1021 : (-dt) * parameters.get<TransitionProbabilities>()[idx_InfectionTransitions] * sum;
263 1021 : }
264 :
265 927 : void Model::compute_flow(Eigen::Index idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt)
266 : {
267 927 : Eigen::Index current_time_index = m_transitions.get_num_time_points() - 1;
268 927 : compute_flow(idx_InfectionTransitions, idx_IncomingFlow, dt, current_time_index);
269 927 : }
270 :
271 103 : void Model::flows_current_timestep(ScalarType dt)
272 : {
273 : // Calculate flow SusceptibleToExposed with force of infection from previous time step and Susceptibles from
274 : // current time step.
275 103 : m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] =
276 206 : dt * m_forceofinfection * m_populations.get_last_value()[Eigen::Index(InfectionState::Susceptible)];
277 :
278 : // Calculate all other flows with compute_flow.
279 : // Flow from Exposed to InfectedNoSymptoms
280 103 : compute_flow(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms),
281 : Eigen::Index(InfectionTransition::SusceptibleToExposed), dt);
282 : // Flow from InfectedNoSymptoms to InfectedSymptoms
283 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms),
284 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt);
285 : // Flow from InfectedNoSymptoms to Recovered
286 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered),
287 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt);
288 : // Flow from InfectedSymptoms to InfectedSevere
289 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
290 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt);
291 : // Flow from InfectedSymptoms to Recovered
292 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered),
293 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt);
294 : // Flow from InfectedSevere to InfectedCritical
295 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
296 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt);
297 : // Flow from InfectedSevere to Recovered
298 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered),
299 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt);
300 : // Flow from InfectedCritical to Dead
301 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead),
302 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt);
303 : // Flow from InfectedCritical to Recovered
304 103 : compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered),
305 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt);
306 103 : }
307 :
308 103 : void Model::update_compartments()
309 : {
310 : // Exposed
311 103 : update_compartment_from_flow(InfectionState::Exposed, {InfectionTransition::SusceptibleToExposed},
312 : {InfectionTransition::ExposedToInfectedNoSymptoms});
313 :
314 : // InfectedNoSymptoms
315 103 : update_compartment_from_flow(InfectionState::InfectedNoSymptoms, {InfectionTransition::ExposedToInfectedNoSymptoms},
316 : {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
317 : InfectionTransition::InfectedNoSymptomsToRecovered});
318 :
319 : // InfectedSymptoms
320 103 : update_compartment_from_flow(
321 : InfectionState::InfectedSymptoms, {InfectionTransition::InfectedNoSymptomsToInfectedSymptoms},
322 : {InfectionTransition::InfectedSymptomsToInfectedSevere, InfectionTransition::InfectedSymptomsToRecovered});
323 :
324 : // InfectedSevere
325 103 : update_compartment_from_flow(
326 : InfectionState::InfectedSevere, {InfectionTransition::InfectedSymptomsToInfectedSevere},
327 : {InfectionTransition::InfectedSevereToInfectedCritical, InfectionTransition::InfectedSevereToRecovered});
328 :
329 : // InfectedCritical
330 103 : update_compartment_from_flow(
331 : InfectionState::InfectedCritical, {InfectionTransition::InfectedSevereToInfectedCritical},
332 : {InfectionTransition::InfectedCriticalToDead, InfectionTransition::InfectedCriticalToRecovered});
333 :
334 : // Recovered
335 206 : update_compartment_from_flow(
336 : InfectionState::Recovered,
337 : {InfectionTransition::InfectedNoSymptomsToRecovered, InfectionTransition::InfectedSymptomsToRecovered,
338 : InfectionTransition::InfectedSevereToRecovered, InfectionTransition::InfectedCriticalToRecovered},
339 206 : std::vector<InfectionTransition>());
340 :
341 : // Dead
342 206 : update_compartment_from_flow(InfectionState::Dead, {InfectionTransition::InfectedCriticalToDead},
343 206 : std::vector<InfectionTransition>());
344 103 : }
345 :
346 721 : void Model::update_compartment_from_flow(InfectionState infectionState,
347 : std::vector<InfectionTransition> const& IncomingFlows,
348 : std::vector<InfectionTransition> const& OutgoingFlows)
349 : {
350 721 : Eigen::Index num_time_points = m_populations.get_num_time_points();
351 721 : ScalarType updated_compartment = m_populations[num_time_points - 2][Eigen::Index(infectionState)];
352 1751 : for (const InfectionTransition& inflow : IncomingFlows) {
353 1030 : updated_compartment += m_transitions.get_last_value()[Eigen::Index(inflow)];
354 : }
355 1648 : for (const InfectionTransition& outflow : OutgoingFlows) {
356 927 : updated_compartment -= m_transitions.get_last_value()[Eigen::Index(outflow)];
357 : }
358 721 : m_populations.get_last_value()[Eigen::Index(infectionState)] = updated_compartment;
359 721 : }
360 :
361 128 : void Model::compute_forceofinfection(ScalarType dt, bool initialization)
362 : {
363 128 : m_forceofinfection = 0;
364 :
365 : // Determine the relevant calculation area = union of the supports of the relevant transition distributions.
366 : ScalarType calc_time =
367 512 : std::max({m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms],
368 128 : m_transitiondistributions_support_max[(int)InfectionTransition::InfectedNoSymptomsToRecovered],
369 128 : m_transitiondistributions_support_max[(int)InfectionTransition::InfectedSymptomsToInfectedSevere],
370 128 : m_transitiondistributions_support_max[(int)InfectionTransition::InfectedSymptomsToRecovered]});
371 :
372 : // Corresponding index.
373 : // Need calc_time_index timesteps in sum,
374 : // subtract 1 because in the last summand all TransitionDistributions evaluate to 0 (by definition of support_max).
375 128 : Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;
376 :
377 : Eigen::Index num_time_points;
378 : ScalarType current_time;
379 : ScalarType deaths;
380 :
381 128 : if (initialization) {
382 : // Determine m_forceofinfection at time t0-dt which is the penultimate timepoint in m_transitions.
383 11 : num_time_points = m_transitions.get_num_time_points() - 1;
384 : // Get time of penultimate timepoint in m_transitions.
385 11 : current_time = m_transitions.get_time(num_time_points - 1);
386 :
387 : // Determine the number of individuals in Dead compartment at time t0-dt.
388 22 : deaths = m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] -
389 11 : m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)];
390 : }
391 : else {
392 : // Determine m_forceofinfection for current last time in m_transitions.
393 117 : num_time_points = m_transitions.get_num_time_points();
394 117 : current_time = m_transitions.get_last_time();
395 117 : deaths = m_populations.get_last_value()[Eigen::Index(InfectionState::Dead)];
396 : }
397 :
398 384 : for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) {
399 256 : Eigen::Index state_age_index = num_time_points - 1 - i;
400 256 : ScalarType state_age = state_age_index * dt;
401 : ScalarType season_val =
402 256 : 1 + parameters.get<Seasonality>() *
403 256 : sin(3.141592653589793 * ((parameters.get<StartDay>() + current_time) / 182.5 + 0.5));
404 :
405 256 : m_forceofinfection +=
406 256 : season_val * parameters.get<TransmissionProbabilityOnContact>().eval(state_age) *
407 512 : parameters.get<ContactPatterns>().get_cont_freq_mat().get_matrix_at(current_time)(0, 0) *
408 256 : (m_transitiondistributions_in_forceofinfection[0][num_time_points - i - 1] *
409 512 : m_transitions[i + 1][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] *
410 256 : parameters.get<RelativeTransmissionNoSymptoms>().eval(state_age) +
411 256 : m_transitiondistributions_in_forceofinfection[1][num_time_points - i - 1] *
412 512 : m_transitions[i + 1][Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] *
413 256 : parameters.get<RiskOfInfectionFromSymptomatic>().eval(state_age));
414 : }
415 128 : m_forceofinfection = 1 / (m_N - deaths) * m_forceofinfection;
416 128 : }
417 :
418 16 : void Model::set_transitiondistributions_support_max(ScalarType dt)
419 : {
420 : // We do not consider the transition SusceptibleToExposed as it is not needed in the computations.
421 160 : for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) {
422 144 : m_transitiondistributions_support_max[transition] =
423 144 : parameters.get<TransitionDistributions>()[transition].get_support_max(dt, m_tol);
424 : }
425 16 : }
426 :
427 16 : void Model::set_transitiondistributions_derivative(ScalarType dt)
428 : {
429 : // We do not consider the transition SusceptibleToExposed as it is not needed in the computations.
430 160 : for (int transition = 1; transition < (int)InfectionTransition::Count; transition++) {
431 : Eigen::Index support_max_index =
432 144 : (Eigen::Index)std::ceil(m_transitiondistributions_support_max[transition] / dt);
433 : // Create vec_derivative that contains the value of the approximated derivative for all necessary time points.
434 : // Here, we evaluate the derivative at time points t_0, ..., t_{support_max_index}.
435 288 : std::vector<ScalarType> vec_derivative(support_max_index + 1, 0.);
436 :
437 655 : for (Eigen::Index i = 0; i <= support_max_index; i++) {
438 : // Compute state_age for considered index.
439 511 : ScalarType state_age = (ScalarType)i * dt;
440 : // Compute derivative.
441 1022 : vec_derivative[i] = (parameters.get<TransitionDistributions>()[transition].eval(state_age) -
442 511 : parameters.get<TransitionDistributions>()[transition].eval(state_age - dt)) /
443 511 : dt;
444 : }
445 144 : m_transitiondistributions_derivative[transition] = vec_derivative;
446 144 : }
447 16 : }
448 :
449 14 : void Model::set_transitiondistributions_in_forceofinfection(ScalarType dt)
450 : {
451 : // Relevant transitions for force of infection term.
452 14 : std::vector<std::vector<int>> relevant_transitions = {
453 : {(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
454 : (int)InfectionTransition::InfectedNoSymptomsToRecovered},
455 : {(int)InfectionTransition::InfectedSymptomsToInfectedSevere,
456 98 : (int)InfectionTransition::InfectedSymptomsToRecovered}};
457 :
458 : // Determine the relevant calculation area = union of the supports of the relevant transition distributions.
459 56 : ScalarType calc_time = std::max({m_transitiondistributions_support_max[relevant_transitions[0][0]],
460 14 : m_transitiondistributions_support_max[relevant_transitions[0][1]],
461 14 : m_transitiondistributions_support_max[relevant_transitions[1][0]],
462 14 : m_transitiondistributions_support_max[relevant_transitions[1][1]]});
463 :
464 : // Corresponding index.
465 : // Need to evaluate survival functions at t_0, ..., t_{calc_time_index} for computation of force of infection,
466 : // subtract 1 because in the last summand all TransitionDistributions evaluate to 0 (by definition of support_max).
467 14 : Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;
468 :
469 : // Compute contributions from survival functions and transition probabilities starting in InfectedNoSymptoms and
470 : // InfectedSymptoms, respectively, on force of infection term.
471 42 : for (int contribution = 0; contribution < 2; contribution++) {
472 56 : std::vector<ScalarType> vec_contribution_to_foi(calc_time_index + 1, 0.);
473 92 : for (Eigen::Index i = 0; i <= calc_time_index; i++) {
474 : // Compute state_age for all indices from t_0, ..., t_{calc_time_index}.
475 64 : ScalarType state_age = i * dt;
476 64 : vec_contribution_to_foi[i] =
477 64 : parameters.get<TransitionProbabilities>()[relevant_transitions[contribution][0]] *
478 64 : parameters.get<TransitionDistributions>()[relevant_transitions[contribution][0]].eval(state_age) +
479 64 : parameters.get<TransitionProbabilities>()[relevant_transitions[contribution][1]] *
480 64 : parameters.get<TransitionDistributions>()[relevant_transitions[contribution][1]].eval(state_age);
481 : }
482 28 : m_transitiondistributions_in_forceofinfection[contribution] = vec_contribution_to_foi;
483 28 : }
484 28 : }
485 :
486 10 : ScalarType Model::get_global_support_max(ScalarType dt) const
487 : {
488 : // Note that this function computes the global_support_max via the get_support_max() method and does not make use
489 : // of the vector m_transitiondistributions_support_max. This is because the global_support_max is already used in
490 : // check_constraints and we cannot ensure that the vector has already been computed when checking for constraints
491 : // (which usually happens before setting the initial flows and simulating).
492 20 : return std::max(
493 10 : {parameters.get<TransitionDistributions>()[(int)InfectionTransition::ExposedToInfectedNoSymptoms]
494 10 : .get_support_max(dt, m_tol),
495 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
496 10 : .get_support_max(dt, m_tol),
497 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedNoSymptomsToRecovered]
498 10 : .get_support_max(dt, m_tol),
499 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedSymptomsToInfectedSevere]
500 10 : .get_support_max(dt, m_tol),
501 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedSymptomsToRecovered]
502 10 : .get_support_max(dt, m_tol),
503 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedSevereToInfectedCritical]
504 10 : .get_support_max(dt, m_tol),
505 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedSevereToRecovered].get_support_max(
506 10 : dt, m_tol),
507 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedCriticalToDead].get_support_max(
508 10 : dt, m_tol),
509 10 : parameters.get<TransitionDistributions>()[(int)InfectionTransition::InfectedCriticalToRecovered]
510 20 : .get_support_max(dt, m_tol)});
511 : }
512 :
513 : } // namespace isecir
514 : } // namespace mio
|
