model.h Source File

CPP API: model.h Source File
ode_secirts/model.h
Go to the documentation of this file.
1 /*
2 * * Copyright (C) 2020-2026 MEmilio
3 *
4 * Authors: Henrik Zunker, Wadim Koslow, Daniel Abele, Martin J. Kühn
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 #ifndef MIO_ODE_SECIRTS_MODEL_H
21 #define MIO_ODE_SECIRTS_MODEL_H
22 
23 #include "memilio/compartments/flow_model.h"
24 #include "memilio/compartments/simulation.h"
25 #include "memilio/compartments/flow_simulation.h"
26 #include "memilio/epidemiology/populations.h"
27 #include "ode_secirts/infection_state.h"
28 #include "ode_secirts/parameters.h"
29 #include "memilio/math/smoother.h"
30 #include "memilio/math/eigen_util.h"
31 
32 #include <numbers>
33 
34 namespace mio
35 {
36 namespace osecirts
37 {
38 // clang-format off
39 using Flows = TypeList<
40  //naive
41  Flow<InfectionState::SusceptibleNaive, InfectionState::ExposedNaive>,
42  Flow<InfectionState::SusceptibleNaive, InfectionState::TemporaryImmunePartialImmunity>,
43  Flow<InfectionState::ExposedNaive, InfectionState::InfectedNoSymptomsNaive>,
44  Flow<InfectionState::InfectedNoSymptomsNaive, InfectionState::InfectedSymptomsNaive>,
45  Flow<InfectionState::InfectedNoSymptomsNaive, InfectionState::TemporaryImmunePartialImmunity>,
46  Flow<InfectionState::InfectedNoSymptomsNaiveConfirmed, InfectionState::InfectedSymptomsNaiveConfirmed>,
47  Flow<InfectionState::InfectedNoSymptomsNaiveConfirmed, InfectionState::TemporaryImmunePartialImmunity>,
48  Flow<InfectionState::InfectedSymptomsNaive, InfectionState::InfectedSevereNaive>,
49  Flow<InfectionState::InfectedSymptomsNaive, InfectionState::TemporaryImmunePartialImmunity>,
50  Flow<InfectionState::InfectedSymptomsNaiveConfirmed, InfectionState::InfectedSevereNaive>,
51  Flow<InfectionState::InfectedSymptomsNaiveConfirmed, InfectionState::TemporaryImmunePartialImmunity>,
52  Flow<InfectionState::InfectedSevereNaive, InfectionState::InfectedCriticalNaive>,
53  Flow<InfectionState::InfectedSevereNaive, InfectionState::TemporaryImmunePartialImmunity>,
54  Flow<InfectionState::InfectedSevereNaive, InfectionState::DeadNaive>,
55  Flow<InfectionState::InfectedCriticalNaive, InfectionState::DeadNaive>,
56  Flow<InfectionState::InfectedCriticalNaive, InfectionState::TemporaryImmunePartialImmunity>,
57  //partial immunity
58  Flow<InfectionState::SusceptiblePartialImmunity, InfectionState::ExposedPartialImmunity>,
59  Flow<InfectionState::SusceptiblePartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
60  Flow<InfectionState::ExposedPartialImmunity, InfectionState::InfectedNoSymptomsPartialImmunity>,
61  Flow<InfectionState::InfectedNoSymptomsPartialImmunity, InfectionState::InfectedSymptomsPartialImmunity>,
62  Flow<InfectionState::InfectedNoSymptomsPartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
63  Flow<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, InfectionState::InfectedSymptomsPartialImmunityConfirmed>,
64  Flow<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
65  Flow<InfectionState::InfectedSymptomsPartialImmunity, InfectionState::InfectedSeverePartialImmunity>,
66  Flow<InfectionState::InfectedSymptomsPartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
67  Flow<InfectionState::InfectedSymptomsPartialImmunityConfirmed, InfectionState::InfectedSeverePartialImmunity>,
68  Flow<InfectionState::InfectedSymptomsPartialImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
69  Flow<InfectionState::InfectedSeverePartialImmunity, InfectionState::InfectedCriticalPartialImmunity>,
70  Flow<InfectionState::InfectedSeverePartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
71  Flow<InfectionState::InfectedSeverePartialImmunity, InfectionState::DeadPartialImmunity>,
72  Flow<InfectionState::InfectedCriticalPartialImmunity, InfectionState::DeadPartialImmunity>,
73  Flow<InfectionState::InfectedCriticalPartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
74  //improved immunity
75  Flow<InfectionState::SusceptibleImprovedImmunity, InfectionState::ExposedImprovedImmunity>,
76  Flow<InfectionState::SusceptibleImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
77  Flow<InfectionState::ExposedImprovedImmunity, InfectionState::InfectedNoSymptomsImprovedImmunity>,
78  Flow<InfectionState::InfectedNoSymptomsImprovedImmunity, InfectionState::InfectedSymptomsImprovedImmunity>,
79  Flow<InfectionState::InfectedNoSymptomsImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
80  Flow<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed, InfectionState::InfectedSymptomsImprovedImmunityConfirmed>,
81  Flow<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
82  Flow<InfectionState::InfectedSymptomsImprovedImmunity, InfectionState::InfectedSevereImprovedImmunity>,
83  Flow<InfectionState::InfectedSymptomsImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
84  Flow<InfectionState::InfectedSymptomsImprovedImmunityConfirmed, InfectionState::InfectedSevereImprovedImmunity>,
85  Flow<InfectionState::InfectedSymptomsImprovedImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
86  Flow<InfectionState::InfectedSevereImprovedImmunity, InfectionState::InfectedCriticalImprovedImmunity>,
87  Flow<InfectionState::InfectedSevereImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
88  Flow<InfectionState::InfectedSevereImprovedImmunity, InfectionState::DeadImprovedImmunity>,
89  Flow<InfectionState::InfectedCriticalImprovedImmunity, InfectionState::DeadImprovedImmunity>,
90  Flow<InfectionState::InfectedCriticalImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
91 
92  // waning
93  Flow<InfectionState::TemporaryImmunePartialImmunity, InfectionState::SusceptiblePartialImmunity>,
94  Flow<InfectionState::TemporaryImmuneImprovedImmunity, InfectionState::SusceptibleImprovedImmunity>,
95  Flow<InfectionState::SusceptibleImprovedImmunity, InfectionState::SusceptiblePartialImmunity>,
96  Flow<InfectionState::SusceptiblePartialImmunity, InfectionState::SusceptibleNaive>>;
97 // clang-format on
98 
99 template <typename FP>
100 class Model
101  : public FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>
102 {
103  using Base = FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>;
104 
105 public:
106  using typename Base::ParameterSet;
107  using typename Base::Populations;
108 
109  Model(const Populations& pop, const ParameterSet& params)
110  : Base(pop, params)
111  {
112  }
113 
114  Model(int num_agegroups)
115  : Model(Populations({AgeGroup(num_agegroups), InfectionState::Count}), ParameterSet(AgeGroup(num_agegroups)))
116  {
117  }
118 
119  void get_flows(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
120  Eigen::Ref<Eigen::VectorX<FP>> flows) const override
121  {
122  using std::floor;
123  using std::min;
124 
125  auto const& params = this->parameters;
126  AgeGroup n_agegroups = params.get_num_groups();
127 
128  ContactMatrixGroup<FP> const& contact_matrix = params.template get<ContactPatterns<FP>>();
129 
130  FP icu_occupancy = 0.0;
131  FP test_and_trace_required = 0.0;
132  for (auto i = AgeGroup(0); i < n_agegroups; ++i) {
133  // naive flow to symptomatic in unit time
134  test_and_trace_required +=
135  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
136  params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
137  (this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsNaive}) +
138  this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsNaiveConfirmed}));
139  // partial immunity flow to symptomatic in unit time
140  test_and_trace_required +=
141  (params.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[i] /
142  params.template get<ReducExposedPartialImmunity<FP>>()[i]) *
143  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
144  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
145  params.template get<ReducTimeInfectedMild<FP>>()[i]) *
146  (this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsPartialImmunity}) +
147  this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}));
148  // improved immunity flow to symptomatic in unit time
149  test_and_trace_required +=
150  (params.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[i] /
151  params.template get<ReducExposedImprovedImmunity<FP>>()[i]) *
152  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
153  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
154  params.template get<ReducTimeInfectedMild<FP>>()[i]) *
155  (this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsImprovedImmunity}) +
156  this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}));
157  icu_occupancy += this->populations.get_from(pop, {i, InfectionState::InfectedCriticalNaive}) +
158  this->populations.get_from(pop, {i, InfectionState::InfectedCriticalPartialImmunity}) +
159  this->populations.get_from(pop, {i, InfectionState::InfectedCriticalImprovedImmunity});
160  }
161 
162  // get vaccinations
163  auto const partial_vaccination = vaccinations_at(t, params.template get<DailyPartialVaccinations<FP>>());
164  auto const full_vaccination = vaccinations_at(t, params.template get<DailyFullVaccinations<FP>>());
165  auto const booster_vaccination = vaccinations_at(t, params.template get<DailyBoosterVaccinations<FP>>());
166 
167  for (auto i = AgeGroup(0); i < n_agegroups; i++) {
168 
169  size_t SNi = this->populations.get_flat_index({i, InfectionState::SusceptibleNaive});
170  size_t ENi = this->populations.get_flat_index({i, InfectionState::ExposedNaive});
171  size_t INSNi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaive});
172  size_t ISyNi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsNaive});
173  size_t ISevNi = this->populations.get_flat_index({i, InfectionState::InfectedSevereNaive});
174  size_t ICrNi = this->populations.get_flat_index({i, InfectionState::InfectedCriticalNaive});
175 
176  size_t INSNCi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaiveConfirmed});
177  size_t ISyNCi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsNaiveConfirmed});
178 
179  size_t SPIi = this->populations.get_flat_index({i, InfectionState::SusceptiblePartialImmunity});
180  size_t EPIi = this->populations.get_flat_index({i, InfectionState::ExposedPartialImmunity});
181  size_t INSPIi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunity});
182  size_t ISyPIi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunity});
183  size_t ISevPIi = this->populations.get_flat_index({i, InfectionState::InfectedSeverePartialImmunity});
184  size_t ICrPIi = this->populations.get_flat_index({i, InfectionState::InfectedCriticalPartialImmunity});
185 
186  size_t INSPICi =
187  this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
188  size_t ISyPICi =
189  this->populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
190 
191  size_t EIIi = this->populations.get_flat_index({i, InfectionState::ExposedImprovedImmunity});
192  size_t INSIIi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunity});
193  size_t ISyIIi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunity});
194  size_t ISevIIi = this->populations.get_flat_index({i, InfectionState::InfectedSevereImprovedImmunity});
195  size_t ICrIIi = this->populations.get_flat_index({i, InfectionState::InfectedCriticalImprovedImmunity});
196 
197  size_t INSIICi =
198  this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
199  size_t ISyIICi =
200  this->populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
201  size_t TImm1 = this->populations.get_flat_index({i, InfectionState::TemporaryImmunePartialImmunity});
202  size_t TImm2 = this->populations.get_flat_index({i, InfectionState::TemporaryImmuneImprovedImmunity});
203 
204  size_t SIIi = this->populations.get_flat_index({i, InfectionState::SusceptibleImprovedImmunity});
205 
206  FP reducExposedPartialImmunity = params.template get<ReducExposedPartialImmunity<FP>>()[i];
207  FP reducExposedImprovedImmunity = params.template get<ReducExposedImprovedImmunity<FP>>()[i];
208  FP reducInfectedSymptomsPartialImmunity =
209  params.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[i];
210  FP reducInfectedSymptomsImprovedImmunity =
211  params.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[i];
212  FP reducInfectedSevereCriticalDeadPartialImmunity =
213  params.template get<ReducInfectedSevereCriticalDeadPartialImmunity<FP>>()[i];
214  FP reducInfectedSevereCriticalDeadImprovedImmunity =
215  params.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<FP>>()[i];
216  FP reducTimeInfectedMild = params.template get<ReducTimeInfectedMild<FP>>()[i];
217 
218  //symptomatic are less well quarantined when testing and tracing is overwhelmed so they infect more people
219  auto riskFromInfectedSymptomatic =
220  smoother_cosine<FP>(test_and_trace_required, params.template get<TestAndTraceCapacity<FP>>(),
221  params.template get<TestAndTraceCapacity<FP>>() *
222  params.template get<TestAndTraceCapacityMaxRiskSymptoms<FP>>(),
223  params.template get<RiskOfInfectionFromSymptomatic<FP>>()[i],
224  params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i]);
225 
226  auto riskFromInfectedNoSymptoms =
227  smoother_cosine<FP>(test_and_trace_required, params.template get<TestAndTraceCapacity<FP>>(),
228  params.template get<TestAndTraceCapacity<FP>>() *
229  params.template get<TestAndTraceCapacityMaxRiskNoSymptoms<FP>>(),
230  params.template get<RelativeTransmissionNoSymptoms<FP>>()[i], 1.0);
231 
232  for (auto j = AgeGroup(0); j < n_agegroups; j++) {
233  size_t SNj = this->populations.get_flat_index({j, InfectionState::SusceptibleNaive});
234  size_t ENj = this->populations.get_flat_index({j, InfectionState::ExposedNaive});
235  size_t INSNj = this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsNaive});
236  size_t ISyNj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsNaive});
237  size_t ISevNj = this->populations.get_flat_index({j, InfectionState::InfectedSevereNaive});
238  size_t ICrNj = this->populations.get_flat_index({j, InfectionState::InfectedCriticalNaive});
239  size_t SIIj = this->populations.get_flat_index({j, InfectionState::SusceptibleImprovedImmunity});
240 
241  size_t INSNCj = this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsNaiveConfirmed});
242  size_t ISyNCj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsNaiveConfirmed});
243 
244  size_t SPIj = this->populations.get_flat_index({j, InfectionState::SusceptiblePartialImmunity});
245  size_t EPIj = this->populations.get_flat_index({j, InfectionState::ExposedPartialImmunity});
246  size_t INSPIj =
247  this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsPartialImmunity});
248  size_t ISyPIj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsPartialImmunity});
249  size_t ISevPIj = this->populations.get_flat_index({j, InfectionState::InfectedSeverePartialImmunity});
250  size_t ICrPIj = this->populations.get_flat_index({j, InfectionState::InfectedCriticalPartialImmunity});
251 
252  size_t INSPICj =
253  this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
254  size_t ISyPICj =
255  this->populations.get_flat_index({j, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
256 
257  size_t EIIj = this->populations.get_flat_index({j, InfectionState::ExposedImprovedImmunity});
258  size_t INSIIj =
259  this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsImprovedImmunity});
260  size_t ISyIIj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsImprovedImmunity});
261  size_t ISevIIj = this->populations.get_flat_index({j, InfectionState::InfectedSevereImprovedImmunity});
262  size_t ICrIIj = this->populations.get_flat_index({j, InfectionState::InfectedCriticalImprovedImmunity});
263 
264  size_t INSIICj =
265  this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
266  size_t ISyIICj =
267  this->populations.get_flat_index({j, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
268 
269  // effective contact rate by contact rate between groups i and j and damping j
270  // std::fmod('time', 365.0) is non differentiable. Use std::floor instead to normalize 'time'.
271  FP normalized_time = (params.template get<StartDay<FP>>() + t) -
272  365.0 * floor((params.template get<StartDay<FP>>() + t) / 365.0);
273  FP season_val = (1 + params.template get<Seasonality<FP>>() *
274  sin(std::numbers::pi_v<ScalarType> * (normalized_time / 182.5 + 0.5)));
275 
276  FP cont_freq_eff =
277  season_val * contact_matrix.get_matrix_at(SimulationTime<FP>(t))(
278  static_cast<Eigen::Index>((size_t)i), static_cast<Eigen::Index>((size_t)j));
279  // without died people
280  FP Nj = pop[SNj] + pop[ENj] + pop[INSNj] + pop[ISyNj] + pop[ISevNj] + pop[ICrNj] + pop[INSNCj] +
281  pop[ISyNCj] + pop[SPIj] + pop[EPIj] + pop[INSPIj] + pop[ISyPIj] + pop[ISevPIj] + pop[ICrPIj] +
282  pop[INSPICj] + pop[ISyPICj] + pop[SIIj] + pop[EIIj] + pop[INSIIj] + pop[ISyIIj] + pop[ISevIIj] +
283  pop[ICrIIj] + pop[INSIICj] + pop[ISyIICj];
284 
285  const FP divNj = (Nj < Limits<FP>::zero_tolerance()) ? FP(0.0) : FP(1.0 / Nj);
286 
287  FP ext_inf_force_dummy = cont_freq_eff * divNj *
288  params.template get<TransmissionProbabilityOnContact<FP>>()[(AgeGroup)i] *
289  (riskFromInfectedNoSymptoms * (pop[INSNj] + pop[INSPIj] + pop[INSIIj]) +
290  riskFromInfectedSymptomatic * (pop[ISyNj] + pop[ISyPIj] + pop[ISyIIj]));
291 
292  FP dummy_SN = y[SNi] * ext_inf_force_dummy;
293 
294  FP dummy_SPI = y[SPIi] * reducExposedPartialImmunity * ext_inf_force_dummy;
295 
296  FP dummy_SII = y[SIIi] * reducExposedImprovedImmunity * ext_inf_force_dummy;
297 
298  flows[this->template get_flat_flow_index<InfectionState::SusceptibleNaive,
299  InfectionState::ExposedNaive>({i})] += dummy_SN;
300  flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
301  InfectionState::ExposedPartialImmunity>({i})] += dummy_SPI;
302  flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
303  InfectionState::ExposedImprovedImmunity>({i})] += dummy_SII;
304  }
305 
306  // vaccinations
307  flows[this->template get_flat_flow_index<InfectionState::SusceptibleNaive,
308  InfectionState::TemporaryImmunePartialImmunity>({i})] =
309  min<FP>(y[SNi] - flows[this->template get_flat_flow_index<InfectionState::SusceptibleNaive,
310  InfectionState::ExposedNaive>({i})],
311  partial_vaccination[static_cast<size_t>(i)]);
312 
313  flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
314  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
315  min<FP>(y[SPIi] -
316  flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
317  InfectionState::ExposedPartialImmunity>({i})],
318  full_vaccination[static_cast<size_t>(i)]);
319 
320  flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
321  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
322  min<FP>(y[SIIi] -
323  flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
324  InfectionState::ExposedImprovedImmunity>({i})],
325  booster_vaccination[static_cast<size_t>(i)]);
326 
327  // ICU capacity shortage is close
328  // TODO: if this is used with vaccination model, it has to be adapted if CriticalPerSevere
329  // is different for different vaccination status. This is not the case here and in addition, ICUCapacity
330  // is set to infinity and this functionality is deactivated, so this is OK for the moment.
331  FP criticalPerSevereAdjusted = smoother_cosine<FP>(
332  icu_occupancy, 0.90 * params.template get<ICUCapacity<FP>>(), params.template get<ICUCapacity<FP>>(),
333  params.template get<CriticalPerSevere<FP>>()[i], 0);
334 
335  FP deathsPerSevereAdjusted = params.template get<CriticalPerSevere<FP>>()[i] - criticalPerSevereAdjusted;
336 
337  /**** path of immune-naive ***/
338  // Exposed
339  flows[this->template get_flat_flow_index<InfectionState::ExposedNaive,
340  InfectionState::InfectedNoSymptomsNaive>({i})] +=
341  y[ENi] / params.template get<TimeExposed<FP>>()[i];
342 
343  // InfectedNoSymptoms
344  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaive,
345  InfectionState::TemporaryImmunePartialImmunity>({i})] =
346  params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i] *
347  (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSNi];
348  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaive,
349  InfectionState::InfectedSymptomsNaive>({i})] =
350  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
351  params.template get<TimeInfectedNoSymptoms<FP>>()[i] * y[INSNi];
352  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaiveConfirmed,
353  InfectionState::InfectedSymptomsNaiveConfirmed>({i})] =
354  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
355  params.template get<TimeInfectedNoSymptoms<FP>>()[i] * y[INSNCi];
356  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaiveConfirmed,
357  InfectionState::TemporaryImmunePartialImmunity>({i})] =
358  params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i] *
359  (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSNCi];
360 
361  // // InfectedSymptoms
362  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaive,
363  InfectionState::InfectedSevereNaive>({i})] =
364  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
365  params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNi];
366 
367  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaive,
368  InfectionState::TemporaryImmunePartialImmunity>({i})] =
369  (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
370  params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNi];
371 
372  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaiveConfirmed,
373  InfectionState::InfectedSevereNaive>({i})] =
374  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
375  params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNCi];
376 
377  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaiveConfirmed,
378  InfectionState::TemporaryImmunePartialImmunity>({i})] =
379  (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
380  params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNCi];
381 
382  // InfectedSevere
383  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereNaive,
384  InfectionState::InfectedCriticalNaive>({i})] =
385  criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
386 
387  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereNaive,
388  InfectionState::TemporaryImmunePartialImmunity>({i})] =
389  (1 - params.template get<CriticalPerSevere<FP>>()[i]) /
390  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
391 
392  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereNaive, InfectionState::DeadNaive>(
393  {i})] = deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
394  // InfectedCritical
395  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalNaive, InfectionState::DeadNaive>(
396  {i})] = params.template get<DeathsPerCritical<FP>>()[i] /
397  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrNi];
398 
399  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalNaive,
400  InfectionState::TemporaryImmunePartialImmunity>({i})] =
401  (1 - params.template get<DeathsPerCritical<FP>>()[i]) /
402  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrNi];
403 
404  // Waning immunity
405  flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
406  InfectionState::SusceptibleNaive>({i})] =
407  1 / params.template get<TimeWaningPartialImmunity<FP>>()[i] * y[SPIi];
408  flows[this->template get_flat_flow_index<InfectionState::TemporaryImmunePartialImmunity,
409  InfectionState::SusceptiblePartialImmunity>({i})] =
410  1 / params.template get<TimeTemporaryImmunityPI<FP>>()[i] * y[TImm1];
411 
412  // /**** path of partially immune ***/
413 
414  // Exposed
415  flows[this->template get_flat_flow_index<InfectionState::ExposedPartialImmunity,
416  InfectionState::InfectedNoSymptomsPartialImmunity>({i})] +=
417  y[EPIi] / params.template get<TimeExposed<FP>>()[i];
418 
419  // InfectedNoSymptoms
420  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunity,
421  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
422  (1 - (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
423  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
424  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPIi];
425  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunity,
426  InfectionState::InfectedSymptomsPartialImmunity>({i})] =
427  (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
428  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
429  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPIi];
430  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed,
431  InfectionState::InfectedSymptomsPartialImmunityConfirmed>({i})] =
432  (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
433  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
434  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPICi];
435  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed,
436  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
437  (1 - (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
438  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
439  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPICi];
440 
441  // InfectedSymptoms
442  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunity,
443  InfectionState::InfectedSeverePartialImmunity>({i})] =
444  reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity *
445  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
446  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPIi];
447 
448  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunity,
449  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
450  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity) *
451  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
452  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPIi];
453 
454  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunityConfirmed,
455  InfectionState::InfectedSeverePartialImmunity>({i})] =
456  reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity *
457  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
458  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPICi];
459 
460  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunityConfirmed,
461  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
462  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity) *
463  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
464  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPICi];
465 
466  // InfectedSevere
467  flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
468  InfectionState::InfectedCriticalPartialImmunity>({i})] =
469  reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity *
470  criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
471 
472  flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
473  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
474  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
475  params.template get<CriticalPerSevere<FP>>()[i]) /
476  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
477 
478  flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
479  InfectionState::DeadPartialImmunity>({i})] =
480  (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
481  deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
482  // InfectedCritical
483  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalPartialImmunity,
484  InfectionState::DeadPartialImmunity>({i})] =
485  (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
486  params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
487  y[ICrPIi];
488 
489  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalPartialImmunity,
490  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
491  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
492  params.template get<DeathsPerCritical<FP>>()[i]) /
493  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrPIi];
494 
495  // /**** path of improved immunity ***/
496  // Exposed
497  flows[this->template get_flat_flow_index<InfectionState::ExposedImprovedImmunity,
498  InfectionState::InfectedNoSymptomsImprovedImmunity>({i})] +=
499  y[EIIi] / params.template get<TimeExposed<FP>>()[i];
500 
501  // InfectedNoSymptoms
502  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunity,
503  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
504  (1 - (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
505  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
506  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIIi];
507  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunity,
508  InfectionState::InfectedSymptomsImprovedImmunity>({i})] =
509  (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
510  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
511  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIIi];
512  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed,
513  InfectionState::InfectedSymptomsImprovedImmunityConfirmed>({i})] =
514  (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
515  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
516  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIICi];
517  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed,
518  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
519  (1 - (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
520  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
521  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIICi];
522 
523  // InfectedSymptoms
524  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunity,
525  InfectionState::InfectedSevereImprovedImmunity>({i})] =
526  reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity *
527  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
528  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIIi];
529 
530  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunity,
531  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
532  (1 - (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity) *
533  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
534  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIIi];
535 
536  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunityConfirmed,
537  InfectionState::InfectedSevereImprovedImmunity>({i})] =
538  reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity *
539  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
540  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIICi];
541 
542  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunityConfirmed,
543  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
544  (1 - (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity) *
545  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
546  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIICi];
547 
548  // InfectedSevere
549  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
550  InfectionState::InfectedCriticalImprovedImmunity>({i})] =
551  reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity *
552  criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
553 
554  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
555  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
556  (1 -
557  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
558  params.template get<CriticalPerSevere<FP>>()[i]) /
559  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
560 
561  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
562  InfectionState::DeadImprovedImmunity>({i})] =
563  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
564  deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
565 
566  // InfectedCritical
567  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalImprovedImmunity,
568  InfectionState::DeadImprovedImmunity>({i})] =
569  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
570  params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
571  y[ICrIIi];
572 
573  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalImprovedImmunity,
574  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
575  (1 -
576  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
577  params.template get<DeathsPerCritical<FP>>()[i]) /
578  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrIIi];
579 
580  // Waning immunity
581  flows[this->template get_flat_flow_index<InfectionState::TemporaryImmuneImprovedImmunity,
582  InfectionState::SusceptibleImprovedImmunity>({i})] =
583  1 / params.template get<TimeTemporaryImmunityII<FP>>()[i] * y[TImm2];
584 
585  flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
586  InfectionState::SusceptiblePartialImmunity>({i})] =
587  1 / params.template get<TimeWaningImprovedImmunity<FP>>()[i] * y[SIIi];
588  }
589  }
590 
602  Eigen::VectorX<FP> vaccinations_at(const FP t,
603  const CustomIndexArray<FP, AgeGroup, SimulationDay>& daily_vaccinations,
604  const FP eps = 0.15) const
605  {
606  using std::floor;
607 
608  auto const& params = this->parameters;
609  const FP ub = floor(t) + 1.0;
610  const FP lb = ub - eps;
611 
612  const auto max_time = static_cast<size_t>(daily_vaccinations.template size<SimulationDay>()) - 1;
613 
614  Eigen::VectorX<FP> smoothed_vaccinations((size_t)params.get_num_groups());
615  smoothed_vaccinations.setZero();
616 
617  // if daily_vaccinations is not available for the current time point, we return zero vaccinations.
618  if (max_time <= floor(t)) {
619  mio::log_warning("Vaccination data not available for time point ", t, ". Returning zero vaccinations.");
620  return smoothed_vaccinations;
621  }
622  if (t >= lb) {
623  for (AgeGroup age = AgeGroup(0); age < params.get_num_groups(); age++) {
624  // if ub + 1 is out of bounds, we use the value at ub
625  FP ubp1 = floor(ub + 1.0);
626  if (max_time < ubp1) {
627  ubp1 = floor(ub);
628  }
629  const auto num_vaccinations_ub = daily_vaccinations[{age, SimulationDay(size_t(floor(ubp1)))}] -
630  daily_vaccinations[{age, SimulationDay(size_t(floor(ub)))}];
631  const auto num_vaccinations_lb = daily_vaccinations[{age, SimulationDay(size_t(floor(lb + 1.0)))}] -
632  daily_vaccinations[{age, SimulationDay(size_t(floor(lb)))}];
633  smoothed_vaccinations[(size_t)age] =
634  smoother_cosine<FP>(t, lb, ub, num_vaccinations_lb, num_vaccinations_ub);
635  }
636  }
637  else {
638  for (auto age = AgeGroup(0); age < params.get_num_groups(); age++) {
639  smoothed_vaccinations[(size_t)age] = daily_vaccinations[{age, SimulationDay(size_t(floor(t + 1)))}] -
640  daily_vaccinations[{age, SimulationDay(size_t(floor(t)))}];
641  }
642  }
643  return smoothed_vaccinations;
644  }
645 
650  template <class IOContext>
651  void serialize(IOContext& io) const
652  {
653  auto obj = io.create_object("Model");
654  obj.add_element("Parameters", this->parameters);
655  obj.add_element("Populations", this->populations);
656  }
657 
662  template <class IOContext>
663  static IOResult<Model> deserialize(IOContext& io)
664  {
665  auto obj = io.expect_object("Model");
666  auto par = obj.expect_element("Parameters", Tag<ParameterSet>{});
667  auto pop = obj.expect_element("Populations", Tag<Populations>{});
668  return apply(
669  io,
670  [](auto&& par_, auto&& pop_) {
671  return Model{pop_, par_};
672  },
673  par, pop);
674  }
675 }; // namespace osecirts
676 
677 //forward declaration, see below.
678 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
679 class Simulation;
680 
688 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
689 FP get_infections_relative(const Simulation<FP, Base>& model, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y);
690 
696 template <typename FP, class BaseT>
697 class Simulation : public BaseT
698 {
699 public:
706  Simulation(mio::osecirts::Model<FP> const& model, FP t0 = 0., FP dt = 0.1)
707  : BaseT(model, t0, dt)
708  , m_t_last_npi_check(t0)
709  {
710  }
711 
725  void apply_variant(const FP t, const CustomIndexArray<UncertainValue<FP>, AgeGroup> base_infectiousness)
726  {
727  using std::min;
728 
729  auto start_day = this->get_model().parameters.template get<StartDay<FP>>();
730  auto start_day_new_variant = this->get_model().parameters.template get<StartDayNewVariant<FP>>();
731 
732  if (start_day + t >= start_day_new_variant - 1e-10) {
733  const FP days_variant = t - (start_day_new_variant - start_day);
734  const FP share_new_variant = min<FP>(1.0, 0.01 * pow(2, (1. / 7) * days_variant));
735  const auto num_groups = this->get_model().parameters.get_num_groups();
736  for (auto i = AgeGroup(0); i < num_groups; ++i) {
737  FP new_transmission = (1 - share_new_variant) * base_infectiousness[i] +
738  share_new_variant * base_infectiousness[i] *
739  this->get_model().parameters.template get<InfectiousnessNewVariant<FP>>()[i];
740  this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>()[i] = new_transmission;
741  }
742  }
743  }
744 
752  Eigen::Ref<Eigen::VectorX<FP>> advance(FP tmax)
753  {
754  using std::floor;
755  using std::min;
756 
757  auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis();
758  auto& dyn_npis = this->get_model().parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
759  auto& contact_patterns = this->get_model().parameters.template get<ContactPatterns<FP>>();
760  // const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups();
761 
762  // in the apply_variant function, we adjust the TransmissionProbabilityOnContact parameter. We need to store
763  // the base value to use it in the apply_variant function and also to reset the parameter after the simulation.
764  auto base_infectiousness = this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>();
765 
766  FP delay_npi_implementation;
767  auto t = BaseT::get_result().get_last_time();
768  const auto dt = dyn_npis.get_interval().get();
769  while (t < tmax) {
770 
771  auto dt_eff = min<FP>({dt, tmax - t, m_t_last_npi_check + dt - t});
772  if (dt_eff >= 1.0) {
773  dt_eff = 1.0;
774  }
775 
776  BaseT::advance(t + dt_eff);
777  if (t + 0.5 + dt_eff - floor(t + 0.5) >= 1) {
778  this->apply_variant(t, base_infectiousness);
779  }
780 
781  if (t > 0) {
782  delay_npi_implementation =
783  this->get_model().parameters.template get<DynamicNPIsImplementationDelay<FP>>();
784  }
785  else {
786  // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay.
787  delay_npi_implementation = 0;
788  }
789  t = t + dt_eff;
790 
791  if (dyn_npis.get_thresholds().size() > 0) {
792  if (floating_point_greater_equal<FP>(t, m_t_last_npi_check + dt)) {
793  if (t < t_end_dyn_npis) {
794  auto inf_rel = get_infections_relative<FP>(*this, t, this->get_result().get_last_value()) *
795  dyn_npis.get_base_value();
796  auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel);
797  if (exceeded_threshold != dyn_npis.get_thresholds().end() &&
798  (exceeded_threshold->first > m_dynamic_npi.first ||
799  t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired
800 
801  auto t_start = SimulationTime<FP>(t + delay_npi_implementation);
802  auto t_end = t_start + SimulationTime<FP>(dyn_npis.get_duration());
803  this->get_model().parameters.get_start_commuter_detection() = (FP)t_start;
804  this->get_model().parameters.get_end_commuter_detection() = (FP)t_end;
805  m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
806  implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
807  t_start, t_end, [](auto& g) {
808  return make_contact_damping_matrix(g);
809  });
810  }
811  }
812  m_t_last_npi_check = t;
813  }
814  }
815  else {
816  m_t_last_npi_check = t;
817  }
818  }
819  // reset TransmissionProbabilityOnContact. This is important for the graph simulation where the advance
820  // function is called multiple times for the same model.
821  this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>() = base_infectiousness;
822 
823  return this->get_result().get_last_value();
824  }
825 
826 private:
827  FP m_t_last_npi_check;
828  std::pair<FP, SimulationTime<FP>> m_dynamic_npi = {-std::numeric_limits<FP>::max(), SimulationTime<FP>(0)};
829 };
830 
842 template <typename FP>
843 inline auto simulate(FP t0, FP tmax, FP dt, const Model<FP>& model,
844  std::unique_ptr<OdeIntegratorCore<FP>>&& integrator_core = nullptr)
845 {
846  return mio::simulate<FP, Model<FP>, Simulation<FP>>(t0, tmax, dt, model, std::move(integrator_core));
847 }
848 
860 template <typename FP>
861 inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model<FP>& model,
862  std::unique_ptr<OdeIntegratorCore<FP>>&& integrator_core = nullptr)
863 {
864  return mio::simulate_flows<FP, Model<FP>, Simulation<FP, mio::FlowSimulation<FP, Model<FP>>>>(
865  t0, tmax, dt, model, std::move(integrator_core));
866 }
867 
868 //see declaration above.
869 template <typename FP, class Base>
870 FP get_infections_relative(const Simulation<FP, Base>& sim, FP /*t*/, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
871 {
872  FP sum_inf = 0;
873  for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) {
874  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsNaive});
875  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsNaiveConfirmed});
876  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsPartialImmunity});
877  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsImprovedImmunity});
878  sum_inf +=
879  sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
880  sum_inf +=
881  sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
882  }
883  FP inf_rel = sum_inf / sim.get_model().populations.get_total();
884 
885  return inf_rel;
886 }
887 
898 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
899 auto get_migration_factors(const Simulation<Base>& sim, FP /*t*/, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
900 {
901  auto& params = sim.get_model().parameters;
902  //parameters as arrays
903  auto& p_asymp = params.template get<RecoveredPerInfectedNoSymptoms<FP>>().array().template cast<FP>();
904  auto& p_inf = params.template get<RiskOfInfectionFromSymptomatic<FP>>().array().template cast<FP>();
905  auto& p_inf_max = params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>().array().template cast<FP>();
906  //slice of InfectedNoSymptoms
907  auto y_INS = slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsNaive),
908  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)}) +
909  slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsPartialImmunity),
910  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)}) +
911  slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsImprovedImmunity),
912  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)});
913 
914  //compute isolation, same as infection risk from main model
915  auto test_and_trace_required =
916  ((1 - p_asymp) / params.template get<TimeInfectedNoSymptoms<FP>>().array().template cast<FP>() * y_INS.array())
917  .sum();
918  auto riskFromInfectedSymptomatic =
919  smoother_cosine<FP>(test_and_trace_required, FP(params.template get<TestAndTraceCapacity<FP>>()),
920  params.template get<TestAndTraceCapacity<FP>>() * 5, p_inf.matrix(), p_inf_max.matrix());
921 
922  //set factor for infected
923  auto factors = Eigen::VectorX<FP>::Ones(y.rows()).eval();
924  slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsNaive), Eigen::Index(size_t(params.get_num_groups())),
925  Eigen::Index(InfectionState::Count)})
926  .array() = riskFromInfectedSymptomatic;
927  slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsPartialImmunity),
928  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)})
929  .array() = riskFromInfectedSymptomatic;
930  slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsImprovedImmunity),
931  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)})
932  .array() = riskFromInfectedSymptomatic;
933  return factors;
934 }
935 
945 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
946 auto test_commuters(Simulation<FP, Base>& sim, Eigen::Ref<Eigen::VectorX<FP>> migrated, FP time)
947 {
948  auto& model = sim.get_model();
949  FP nondetection = 1.0;
950  if (time >= model.parameters.get_start_commuter_detection() &&
951  time < model.parameters.get_end_commuter_detection()) {
952  nondetection = (FP)model.parameters.get_commuter_nondetection();
953  }
954  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); ++i) {
955  auto ISyNi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsNaive});
956  auto ISyNCi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsNaiveConfirmed});
957  auto INSNi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaive});
958  auto INSNCi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaiveConfirmed});
959 
960  auto ISPIi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunity});
961  auto ISPICi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
962  auto INSPIi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunity});
963  auto INSPICi =
964  model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
965 
966  auto ISyIIi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunity});
967  auto ISyIICi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
968  auto INSIIi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunity});
969  auto INSIICi =
970  model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
971 
972  //put detected commuters in their own compartment so they don't contribute to infections in their home node
973  sim.get_result().get_last_value()[ISyNi] -= migrated[ISyNi] * (1 - nondetection);
974  sim.get_result().get_last_value()[ISyNCi] += migrated[ISyNi] * (1 - nondetection);
975  sim.get_result().get_last_value()[INSNi] -= migrated[INSNi] * (1 - nondetection);
976  sim.get_result().get_last_value()[INSNCi] += migrated[INSNi] * (1 - nondetection);
977 
978  sim.get_result().get_last_value()[ISPIi] -= migrated[ISPIi] * (1 - nondetection);
979  sim.get_result().get_last_value()[ISPICi] += migrated[ISPIi] * (1 - nondetection);
980  sim.get_result().get_last_value()[INSPIi] -= migrated[INSPIi] * (1 - nondetection);
981  sim.get_result().get_last_value()[INSPICi] += migrated[INSPIi] * (1 - nondetection);
982 
983  sim.get_result().get_last_value()[ISyIIi] -= migrated[ISyIIi] * (1 - nondetection);
984  sim.get_result().get_last_value()[ISyIICi] += migrated[ISyIIi] * (1 - nondetection);
985  sim.get_result().get_last_value()[INSIIi] -= migrated[INSIIi] * (1 - nondetection);
986  sim.get_result().get_last_value()[INSIICi] += migrated[INSIIi] * (1 - nondetection);
987 
988  //reduce the number of commuters
989  migrated[ISyNi] *= nondetection;
990  migrated[INSNi] *= nondetection;
991 
992  migrated[ISPIi] *= nondetection;
993  migrated[INSPIi] *= nondetection;
994 
995  migrated[ISyIIi] *= nondetection;
996  migrated[INSIIi] *= nondetection;
997  }
998 }
999 
1000 } // namespace osecirts
1001 } // namespace mio
1002 
1003 #endif //MIO_ODE_SECIRTS_MODEL_H
mio::ContactMatrixGroup
represents a collection of contact frequency matrices that whose sum is the total number of contacts.
Definition: contact_matrix.h:536
mio::CustomIndexArray
A class template for an array with custom indices.
Definition: custom_index_array.h:136
mio::FlowModel
A FlowModel is a CompartmentalModel defined by the flows between compartments.
Definition: flow_model.h:74
mio::IntegratorCore
Interface class defining the integration step used in a SystemIntegrator.
Definition: integrator.h:48
mio::Populations
A class template for compartment populations.
Definition: populations.h:55
mio::Populations::get_from
decltype(auto) get_from(Arr &&y, Index const &cats) const
get_from returns the value of a flat container at the flat index corresponding to set of enum values.
Definition: populations.h:123
mio::SimulationDay
Represents the simulation time as an integer index.
Definition: simulation_day.h:32
mio::SimulationTime
double simulation time.
Definition: damping.h:58
mio::abm::Parameters
Parameters of the simulation that are the same everywhere within the Model.
Definition: abm/parameters.h:764
mio::abm::Simulation
Run the Simulation in discrete steps, evolve the Model and report results.
Definition: models/abm/simulation.h:37
mio::osecirts::Model
Definition: ode_secirts/model.h:102
mio::osecirts::Model::vaccinations_at
Eigen::VectorX< FP > vaccinations_at(const FP t, const CustomIndexArray< FP, AgeGroup, SimulationDay > &daily_vaccinations, const FP eps=0.15) const
Calculates smoothed vaccinations for a given time point.
Definition: ode_secirts/model.h:602
mio::osecirts::Model::Model
Model(const Populations &pop, const ParameterSet &params)
Definition: ode_secirts/model.h:109
mio::osecirts::Model::serialize
void serialize(IOContext &io) const
serialize this.
Definition: ode_secirts/model.h:651
mio::osecirts::Model::deserialize
static IOResult< Model > deserialize(IOContext &io)
deserialize an object of this class.
Definition: ode_secirts/model.h:663
mio::osecirts::Model::Model
Model(int num_agegroups)
Definition: ode_secirts/model.h:114
mio::osecirts::Model::get_flows
void get_flows(Eigen::Ref< const Eigen::VectorX< FP >> pop, Eigen::Ref< const Eigen::VectorX< FP >> y, FP t, Eigen::Ref< Eigen::VectorX< FP >> flows) const override
Definition: ode_secirts/model.h:119
mio::osecirts::Parameters
Parameters of the age-resolved SECIRS-type model with high temporary immunity upon immunization and w...
Definition: ode_secirts/parameters.h:769
mio::osecirts::Simulation
specialization of compartment model simulation for the SECIRTS model.
Definition: ode_secirts/model.h:698
mio::osecirts::Simulation::Simulation
Simulation(mio::osecirts::Model< FP > const &model, FP t0=0., FP dt=0.1)
construct a simulation.
Definition: ode_secirts/model.h:706
mio::osecirts::Simulation::apply_variant
void apply_variant(const FP t, const CustomIndexArray< UncertainValue< FP >, AgeGroup > base_infectiousness)
Applies the effect of a new variant of a disease to the transmission probability of the model.
Definition: ode_secirts/model.h:725
mio::osecirts::Simulation::m_t_last_npi_check
FP m_t_last_npi_check
Definition: ode_secirts/model.h:827
mio::osecirts::Simulation::advance
Eigen::Ref< Eigen::VectorX< FP > > advance(FP tmax)
advance simulation to tmax.
Definition: ode_secirts/model.h:752
eigen_util.h
flow_model.h
flow_simulation.h
ad::internal::sin
ad::internal::unary_intermediate< AD_TAPE_REAL, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::operations::ad_sin< AD_TAPE_REAL > > sin(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &x1)
Definition: ad.hpp:913
ad::internal::floor
static double floor(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &x)
Definition: ad.hpp:2451
ad::internal::min
static min_max_return_type< ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > >::type min(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &a, const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &b)
Definition: ad.hpp:2599
ad::internal::pow
ad::internal::binary_intermediate_aa< AD_TAPE_REAL, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::operations::ad_pow_aa< AD_TAPE_REAL > > pow(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &x1, const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &x2)
Definition: ad.hpp:1610
ad::internal::max
static min_max_return_type< ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > >::type max(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &a, const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &b)
Definition: ad.hpp:2596
mio::osecirts::test_commuters
auto test_commuters(Simulation< FP, Base > &sim, Eigen::Ref< Eigen::VectorX< FP >> migrated, FP time)
Adjusts the state of commuters in a model, accounting for detection and mobility effects.
Definition: ode_secirts/model.h:946
mio::osecirts::get_migration_factors
auto get_migration_factors(const Simulation< Base > &sim, FP, const Eigen::Ref< const Eigen::VectorX< FP >> &y)
Get migration factors.
Definition: ode_secirts/model.h:899
mio::osecirts::simulate
auto simulate(FP t0, FP tmax, FP dt, const Model< FP > &model, std::unique_ptr< OdeIntegratorCore< FP >> &&integrator_core=nullptr)
Specialization of simulate for SECIRS-type models using Simulation.
Definition: ode_secirts/model.h:843
mio::osecirts::simulate_flows
auto simulate_flows(FP t0, FP tmax, FP dt, const Model< FP > &model, std::unique_ptr< OdeIntegratorCore< FP >> &&integrator_core=nullptr)
Specialization of simulate for SECIRS-type models using the FlowSimulation.
Definition: ode_secirts/model.h:861
mio::osecirts::get_infections_relative
FP get_infections_relative(const Simulation< FP, Base > &model, FP t, const Eigen::Ref< const Eigen::VectorX< FP >> &y)
get percentage of infections per total population.
Definition: ode_secirts/model.h:870
mio
A collection of classes to simplify handling of matrix shapes in meta programming.
Definition: models/abm/analyze_result.h:30
mio::Tag
boost::outcome_v2::in_place_type_t< T > Tag
Type that is used for overload resolution.
Definition: io.h:407
mio::i
auto i
Definition: io.h:809
mio::log_warning
void log_warning(spdlog::string_view_t fmt, const Args &... args)
Definition: logging.h:112
mio::apply
details::ApplyResultT< F, T... > apply(IOContext &io, F f, const IOResult< T > &... rs)
Evaluate a function with zero or more unpacked IOResults as arguments.
Definition: io.h:481
mio::slice
auto slice(V &&v, Seq< Eigen::Index > elems)
take a regular slice of a row or column vector.
Definition: eigen_util.h:107
mio::implement_dynamic_npis
void implement_dynamic_npis(DampingExprGroup &damping_expr_group, const std::vector< DampingSampling< FP >> &npis, SimulationTime< FP > begin, SimulationTime< FP > end, MakeMatrix &&make_matrix)
implement dynamic NPIs for a time span.
Definition: dynamic_npis.h:279
mio::get
constexpr std::tuple_element< I, std::tuple< Index< CategoryTags >... > >::type & get(Index< CategoryTags... > &i) noexcept
Retrieves the Index (by reference) at the Ith position of a MultiIndex.
Definition: index.h:294
mio::IOResult
boost::outcome_v2::unchecked< T, IOStatus > IOResult
Value-or-error type for operations that return a value but can fail.
Definition: io.h:353
infection_state.h
parameters.h
mio::AgeGroup
The AgeGroup struct is used as a dynamically sized tag for all age dependent categories.
Definition: age_group.h:32
mio::CompartmentalModel< FP, InfectionState, mio::Populations< FP, AgeGroup, InfectionState >, Parameters< FP > >::Populations
mio::Populations< FP, AgeGroup, InfectionState > Populations
Definition: compartmental_model.h:62
mio::Flow
A Flow defines a possible transition between two Compartments in a FlowModel.
Definition: flow.h:33
mio::TypeList
Collection of types. Each type is mapped to an index of type size_t.
Definition: type_list.h:32
mio::osecirts::RecoveredPerInfectedNoSymptoms
The percentage of asymptomatic cases in the SECIRTS model.
Definition: ode_secirts/parameters.h:394
mio::osecirts::ReducExposedPartialImmunity
Factor to reduce infection risk for persons with partial immunity.
Definition: ode_secirts/parameters.h:598
mio::osecirts::TestAndTraceCapacity
The Capacity to test and trace contacts of infected for quarantine per day.
Definition: ode_secirts/parameters.h:119
mio::osecirts::TimeInfectedNoSymptoms
The (mean) time in day unit for asymptomatic cases that are infected but have not yet developed sympt...
Definition: ode_secirts/parameters.h:220
Generated by doxygen 1.9.1