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] - params.template get<DeathsPerSevere<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})] = (params.template get<DeathsPerSevere<FP>>()[i] + deathsPerSevereAdjusted) /
394  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
395  // InfectedCritical
396  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalNaive, InfectionState::DeadNaive>(
397  {i})] = params.template get<DeathsPerCritical<FP>>()[i] /
398  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrNi];
399 
400  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalNaive,
401  InfectionState::TemporaryImmunePartialImmunity>({i})] =
402  (1 - params.template get<DeathsPerCritical<FP>>()[i]) /
403  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrNi];
404 
405  // Waning immunity
406  flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
407  InfectionState::SusceptibleNaive>({i})] =
408  1 / params.template get<TimeWaningPartialImmunity<FP>>()[i] * y[SPIi];
409  flows[this->template get_flat_flow_index<InfectionState::TemporaryImmunePartialImmunity,
410  InfectionState::SusceptiblePartialImmunity>({i})] =
411  1 / params.template get<TimeTemporaryImmunityPI<FP>>()[i] * y[TImm1];
412 
413  // /**** path of partially immune ***/
414 
415  // Exposed
416  flows[this->template get_flat_flow_index<InfectionState::ExposedPartialImmunity,
417  InfectionState::InfectedNoSymptomsPartialImmunity>({i})] +=
418  y[EPIi] / params.template get<TimeExposed<FP>>()[i];
419 
420  // InfectedNoSymptoms
421  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunity,
422  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
423  (1 - (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
424  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
425  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPIi];
426  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunity,
427  InfectionState::InfectedSymptomsPartialImmunity>({i})] =
428  (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
429  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
430  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPIi];
431  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed,
432  InfectionState::InfectedSymptomsPartialImmunityConfirmed>({i})] =
433  (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
434  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
435  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPICi];
436  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed,
437  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
438  (1 - (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
439  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
440  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPICi];
441 
442  // InfectedSymptoms
443  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunity,
444  InfectionState::InfectedSeverePartialImmunity>({i})] =
445  reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity *
446  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
447  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPIi];
448 
449  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunity,
450  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
451  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity) *
452  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
453  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPIi];
454 
455  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunityConfirmed,
456  InfectionState::InfectedSeverePartialImmunity>({i})] =
457  reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity *
458  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
459  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPICi];
460 
461  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunityConfirmed,
462  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
463  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity) *
464  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
465  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPICi];
466 
467  // InfectedSevere
468  flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
469  InfectionState::InfectedCriticalPartialImmunity>({i})] =
470  reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity *
471  criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
472 
473  flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
474  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
475  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
476  (params.template get<CriticalPerSevere<FP>>()[i] +
477  params.template get<DeathsPerSevere<FP>>()[i])) /
478  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
479 
480  flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
481  InfectionState::DeadPartialImmunity>({i})] =
482  (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
483  (params.template get<DeathsPerSevere<FP>>()[i] + deathsPerSevereAdjusted) /
484  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
485  // InfectedCritical
486  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalPartialImmunity,
487  InfectionState::DeadPartialImmunity>({i})] =
488  (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
489  params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
490  y[ICrPIi];
491 
492  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalPartialImmunity,
493  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
494  (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
495  params.template get<DeathsPerCritical<FP>>()[i]) /
496  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrPIi];
497 
498  // /**** path of improved immunity ***/
499  // Exposed
500  flows[this->template get_flat_flow_index<InfectionState::ExposedImprovedImmunity,
501  InfectionState::InfectedNoSymptomsImprovedImmunity>({i})] +=
502  y[EIIi] / params.template get<TimeExposed<FP>>()[i];
503 
504  // InfectedNoSymptoms
505  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunity,
506  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
507  (1 - (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
508  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
509  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIIi];
510  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunity,
511  InfectionState::InfectedSymptomsImprovedImmunity>({i})] =
512  (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
513  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
514  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIIi];
515  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed,
516  InfectionState::InfectedSymptomsImprovedImmunityConfirmed>({i})] =
517  (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
518  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
519  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIICi];
520  flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed,
521  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
522  (1 - (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
523  (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
524  (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIICi];
525 
526  // InfectedSymptoms
527  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunity,
528  InfectionState::InfectedSevereImprovedImmunity>({i})] =
529  reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity *
530  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
531  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIIi];
532 
533  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunity,
534  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
535  (1 - (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity) *
536  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
537  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIIi];
538 
539  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunityConfirmed,
540  InfectionState::InfectedSevereImprovedImmunity>({i})] =
541  reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity *
542  params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
543  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIICi];
544 
545  flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunityConfirmed,
546  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
547  (1 - (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity) *
548  params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
549  (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIICi];
550 
551  // InfectedSevere
552  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
553  InfectionState::InfectedCriticalImprovedImmunity>({i})] =
554  reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity *
555  criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
556 
557  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
558  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
559  (1 -
560  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
561  (params.template get<CriticalPerSevere<FP>>()[i] +
562  params.template get<DeathsPerSevere<FP>>()[i])) /
563  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
564 
565  flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
566  InfectionState::DeadImprovedImmunity>({i})] =
567  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
568  (params.template get<DeathsPerSevere<FP>>()[i] + deathsPerSevereAdjusted) /
569  params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
570 
571  // InfectedCritical
572  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalImprovedImmunity,
573  InfectionState::DeadImprovedImmunity>({i})] =
574  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
575  params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
576  y[ICrIIi];
577 
578  flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalImprovedImmunity,
579  InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
580  (1 -
581  (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
582  params.template get<DeathsPerCritical<FP>>()[i]) /
583  params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrIIi];
584 
585  // Waning immunity
586  flows[this->template get_flat_flow_index<InfectionState::TemporaryImmuneImprovedImmunity,
587  InfectionState::SusceptibleImprovedImmunity>({i})] =
588  1 / params.template get<TimeTemporaryImmunityII<FP>>()[i] * y[TImm2];
589 
590  flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
591  InfectionState::SusceptiblePartialImmunity>({i})] =
592  1 / params.template get<TimeWaningImprovedImmunity<FP>>()[i] * y[SIIi];
593  }
594  }
595 
607  Eigen::VectorX<FP> vaccinations_at(const FP t,
608  const CustomIndexArray<FP, AgeGroup, SimulationDay>& daily_vaccinations,
609  const FP eps = 0.15) const
610  {
611  using std::floor;
612 
613  auto const& params = this->parameters;
614  const FP ub = floor(t) + 1.0;
615  const FP lb = ub - eps;
616 
617  const auto max_time = static_cast<size_t>(daily_vaccinations.template size<SimulationDay>()) - 1;
618 
619  Eigen::VectorX<FP> smoothed_vaccinations((size_t)params.get_num_groups());
620  smoothed_vaccinations.setZero();
621 
622  // if daily_vaccinations is not available for the current time point, we return zero vaccinations.
623  if (max_time <= floor(t)) {
624  mio::log_warning("Vaccination data not available for time point ", t, ". Returning zero vaccinations.");
625  return smoothed_vaccinations;
626  }
627  if (t >= lb) {
628  for (AgeGroup age = AgeGroup(0); age < params.get_num_groups(); age++) {
629  // if ub + 1 is out of bounds, we use the value at ub
630  FP ubp1 = floor(ub + 1.0);
631  if (max_time < ubp1) {
632  ubp1 = floor(ub);
633  }
634  const auto num_vaccinations_ub = daily_vaccinations[{age, SimulationDay(size_t(floor(ubp1)))}] -
635  daily_vaccinations[{age, SimulationDay(size_t(floor(ub)))}];
636  const auto num_vaccinations_lb = daily_vaccinations[{age, SimulationDay(size_t(floor(lb + 1.0)))}] -
637  daily_vaccinations[{age, SimulationDay(size_t(floor(lb)))}];
638  smoothed_vaccinations[(size_t)age] =
639  smoother_cosine<FP>(t, lb, ub, num_vaccinations_lb, num_vaccinations_ub);
640  }
641  }
642  else {
643  for (auto age = AgeGroup(0); age < params.get_num_groups(); age++) {
644  smoothed_vaccinations[(size_t)age] = daily_vaccinations[{age, SimulationDay(size_t(floor(t + 1)))}] -
645  daily_vaccinations[{age, SimulationDay(size_t(floor(t)))}];
646  }
647  }
648  return smoothed_vaccinations;
649  }
650 
655  template <class IOContext>
656  void serialize(IOContext& io) const
657  {
658  auto obj = io.create_object("Model");
659  obj.add_element("Parameters", this->parameters);
660  obj.add_element("Populations", this->populations);
661  }
662 
667  template <class IOContext>
668  static IOResult<Model> deserialize(IOContext& io)
669  {
670  auto obj = io.expect_object("Model");
671  auto par = obj.expect_element("Parameters", Tag<ParameterSet>{});
672  auto pop = obj.expect_element("Populations", Tag<Populations>{});
673  return apply(
674  io,
675  [](auto&& par_, auto&& pop_) {
676  return Model{pop_, par_};
677  },
678  par, pop);
679  }
680 }; // namespace osecirts
681 
682 //forward declaration, see below.
683 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
684 class Simulation;
685 
693 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
694 FP get_infections_relative(const Simulation<FP, Base>& model, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y);
695 
701 template <typename FP, class BaseT>
702 class Simulation : public BaseT
703 {
704 public:
711  Simulation(mio::osecirts::Model<FP> const& model, FP t0 = 0., FP dt = 0.1)
712  : BaseT(model, t0, dt)
713  {
714  }
715 
729  void apply_variant(const FP t, const CustomIndexArray<UncertainValue<FP>, AgeGroup> base_infectiousness)
730  {
731  using std::min;
732 
733  auto start_day = this->get_model().parameters.template get<StartDay<FP>>();
734  auto start_day_new_variant = this->get_model().parameters.template get<StartDayNewVariant<FP>>();
735 
736  if (start_day + t >= start_day_new_variant - 1e-10) {
737  const FP days_variant = t - (start_day_new_variant - start_day);
738  const FP share_new_variant = min<FP>(1.0, 0.01 * pow(2, (1. / 7) * days_variant));
739  const auto num_groups = this->get_model().parameters.get_num_groups();
740  for (auto i = AgeGroup(0); i < num_groups; ++i) {
741  FP new_transmission = (1 - share_new_variant) * base_infectiousness[i] +
742  share_new_variant * base_infectiousness[i] *
743  this->get_model().parameters.template get<InfectiousnessNewVariant<FP>>()[i];
744  this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>()[i] = new_transmission;
745  }
746  }
747  }
748 
756  Eigen::Ref<Eigen::VectorX<FP>> advance(FP tmax)
757  {
758  using std::floor;
759  using std::min;
760 
761  auto& dyn_npis = this->get_model().parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
762  auto& contact_patterns = this->get_model().parameters.template get<ContactPatterns<FP>>();
763  // const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups();
764 
765  // in the apply_variant function, we adjust the TransmissionProbabilityOnContact parameter. We need to store
766  // the base value to use it in the apply_variant function and also to reset the parameter after the simulation.
767  auto base_infectiousness = this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>();
768 
769  FP delay_npi_implementation;
770  FP t = BaseT::get_result().get_last_time();
771  while (t < tmax) {
772 
773  if (t > 0) {
774  delay_npi_implementation = FP(dyn_npis.get_implementation_delay());
775  }
776  else {
777  // DynamicNPIs for t=0 are treated as 'from-start NPIs'. I.e., do not enforce delay.
778  delay_npi_implementation = 0;
779  }
780 
781  if (dyn_npis.get_thresholds().size() > 0) {
782  FP direc_begin = FP(dyn_npis.get_directive_begin());
783  FP direc_end = FP(dyn_npis.get_directive_end());
784  if (floating_point_greater_equal(t, direc_begin, 1e-10) && t < direc_end) {
785  auto inf_rel = get_infections_relative<FP>(*this, t, this->get_result().get_last_value()) *
786  dyn_npis.get_base_value();
787  auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel);
788  const bool npi_expired = floating_point_greater_equal(t, FP(m_dynamic_npi.second), 1e-10);
789  if (exceeded_threshold != dyn_npis.get_thresholds().end() &&
790  (exceeded_threshold->first > m_dynamic_npi.first || npi_expired)) {
791  //old npi was weaker or is expired
792 
793  // Keep-alive: if the NPI expired but the threshold is still exceeded at the
794  // same level, renew immediately without delay to avoid a gap.
795  // Apply implementation delay only if stronger NPI needed.
796  bool is_expiry_renewal =
797  npi_expired && !(exceeded_threshold->first > m_dynamic_npi.first);
798  FP effective_delay = is_expiry_renewal ? FP(0) : delay_npi_implementation;
799  if (t + effective_delay < direc_end) {
800  auto t_start = SimulationTime<FP>(t + effective_delay);
801  // set the end to the minimum of start+duration and the end of the directive
802  auto t_end = SimulationTime<FP>(min<FP>(direc_end, FP(t_start + 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  // For new NPIs (t_start > 0): shift t_start_damping by +1 so the smooth
807  // transition window [t_start, t_start+1] lies in the future, consistent with
808  // predefined dampings. For t_start = 0 (global t0): the window [-1, 0]
809  // is in the past and no shift is needed.
810  // For keep-alive renewals (is_expiry_renewal): do not shift t_start_damping.
811  // Since t_start == t_end_damping_old, the new start damping has the same
812  // (time, level, type) as the previous entry.
813  // Therefore, the contact matrix stays constant at the NPI level with no dip.
814  auto t_start_damping = (!is_expiry_renewal && FP(t_start) > FP(0))
815  ? SimulationTime<FP>(FP(t_start) + FP(1))
816  : t_start;
817  auto t_end_damping = SimulationTime<FP>(FP(t_end) + FP(1));
818  implement_dynamic_npis<FP>(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
819  t_start_damping, t_end_damping, [](auto& g) {
820  return make_contact_damping_matrix(g);
821  });
822  }
823  }
824  }
825  }
826 
827  auto dt_eff = min<FP>(1.0, tmax - t);
828  // Clamp step size at the NPI end time to avoid stepping over the lifting phase.
829  // This ensures the keep-alive check activates exactly at t_end, preventing a brief dip.
830  FP npi_end = FP(m_dynamic_npi.second);
831  if (t < npi_end) {
832  dt_eff = min<FP>(dt_eff, npi_end - t);
833  }
834  BaseT::advance(t + dt_eff);
835  if (t + 0.5 + dt_eff - floor(t + 0.5) >= 1) {
836  this->apply_variant(t, base_infectiousness);
837  }
838  t = t + dt_eff;
839  }
840  // reset TransmissionProbabilityOnContact. This is important for the graph simulation where the advance
841  // function is called multiple times for the same model.
842  this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>() = base_infectiousness;
843 
844  return this->get_result().get_last_value();
845  }
846 
847 private:
848  std::pair<FP, SimulationTime<FP>> m_dynamic_npi = {-std::numeric_limits<FP>::max(), SimulationTime<FP>(0)};
849 };
850 
862 template <typename FP>
863 inline auto simulate(FP t0, FP tmax, FP dt, const Model<FP>& model,
864  std::unique_ptr<OdeIntegratorCore<FP>>&& integrator_core = nullptr)
865 {
866  return mio::simulate<FP, Model<FP>, Simulation<FP>>(t0, tmax, dt, model, std::move(integrator_core));
867 }
868 
880 template <typename FP>
881 inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model<FP>& model,
882  std::unique_ptr<OdeIntegratorCore<FP>>&& integrator_core = nullptr)
883 {
884  return mio::simulate_flows<FP, Model<FP>, Simulation<FP, mio::FlowSimulation<FP, Model<FP>>>>(
885  t0, tmax, dt, model, std::move(integrator_core));
886 }
887 
888 //see declaration above.
889 template <typename FP, class Base>
890 FP get_infections_relative(const Simulation<FP, Base>& sim, FP /*t*/, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
891 {
892  FP sum_inf = 0;
893  for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) {
894  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsNaive});
895  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsNaiveConfirmed});
896  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsPartialImmunity});
897  sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsImprovedImmunity});
898  sum_inf +=
899  sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
900  sum_inf +=
901  sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
902  }
903  FP inf_rel = sum_inf / sim.get_model().populations.get_total();
904 
905  return inf_rel;
906 }
907 
918 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
919 auto get_migration_factors(const Simulation<Base>& sim, FP /*t*/, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
920 {
921  auto& params = sim.get_model().parameters;
922  //parameters as arrays
923  auto& p_asymp = params.template get<RecoveredPerInfectedNoSymptoms<FP>>().array().template cast<FP>();
924  auto& p_inf = params.template get<RiskOfInfectionFromSymptomatic<FP>>().array().template cast<FP>();
925  auto& p_inf_max = params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>().array().template cast<FP>();
926  //slice of InfectedNoSymptoms
927  auto y_INS = slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsNaive),
928  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)}) +
929  slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsPartialImmunity),
930  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)}) +
931  slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsImprovedImmunity),
932  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)});
933 
934  //compute isolation, same as infection risk from main model
935  auto test_and_trace_required =
936  ((1 - p_asymp) / params.template get<TimeInfectedNoSymptoms<FP>>().array().template cast<FP>() * y_INS.array())
937  .sum();
938  auto riskFromInfectedSymptomatic =
939  smoother_cosine<FP>(test_and_trace_required, FP(params.template get<TestAndTraceCapacity<FP>>()),
940  params.template get<TestAndTraceCapacity<FP>>() * 5, p_inf.matrix(), p_inf_max.matrix());
941 
942  //set factor for infected
943  auto factors = Eigen::VectorX<FP>::Ones(y.rows()).eval();
944  slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsNaive), Eigen::Index(size_t(params.get_num_groups())),
945  Eigen::Index(InfectionState::Count)})
946  .array() = riskFromInfectedSymptomatic;
947  slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsPartialImmunity),
948  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)})
949  .array() = riskFromInfectedSymptomatic;
950  slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsImprovedImmunity),
951  Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)})
952  .array() = riskFromInfectedSymptomatic;
953  return factors;
954 }
955 
965 template <typename FP, class Base = mio::Simulation<FP, Model<FP>>>
966 auto test_commuters(Simulation<FP, Base>& sim, Eigen::Ref<Eigen::VectorX<FP>> migrated, FP time)
967 {
968  auto& model = sim.get_model();
969  FP nondetection = 1.0;
970  if (time >= model.parameters.get_start_commuter_detection() &&
971  time < model.parameters.get_end_commuter_detection()) {
972  nondetection = (FP)model.parameters.get_commuter_nondetection();
973  }
974  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); ++i) {
975  auto ISyNi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsNaive});
976  auto ISyNCi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsNaiveConfirmed});
977  auto INSNi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaive});
978  auto INSNCi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaiveConfirmed});
979 
980  auto ISPIi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunity});
981  auto ISPICi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
982  auto INSPIi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunity});
983  auto INSPICi =
984  model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
985 
986  auto ISyIIi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunity});
987  auto ISyIICi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
988  auto INSIIi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunity});
989  auto INSIICi =
990  model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
991 
992  //put detected commuters in their own compartment so they don't contribute to infections in their home node
993  sim.get_result().get_last_value()[ISyNi] -= migrated[ISyNi] * (1 - nondetection);
994  sim.get_result().get_last_value()[ISyNCi] += migrated[ISyNi] * (1 - nondetection);
995  sim.get_result().get_last_value()[INSNi] -= migrated[INSNi] * (1 - nondetection);
996  sim.get_result().get_last_value()[INSNCi] += migrated[INSNi] * (1 - nondetection);
997 
998  sim.get_result().get_last_value()[ISPIi] -= migrated[ISPIi] * (1 - nondetection);
999  sim.get_result().get_last_value()[ISPICi] += migrated[ISPIi] * (1 - nondetection);
1000  sim.get_result().get_last_value()[INSPIi] -= migrated[INSPIi] * (1 - nondetection);
1001  sim.get_result().get_last_value()[INSPICi] += migrated[INSPIi] * (1 - nondetection);
1002 
1003  sim.get_result().get_last_value()[ISyIIi] -= migrated[ISyIIi] * (1 - nondetection);
1004  sim.get_result().get_last_value()[ISyIICi] += migrated[ISyIIi] * (1 - nondetection);
1005  sim.get_result().get_last_value()[INSIIi] -= migrated[INSIIi] * (1 - nondetection);
1006  sim.get_result().get_last_value()[INSIICi] += migrated[INSIIi] * (1 - nondetection);
1007 
1008  //reduce the number of commuters
1009  migrated[ISyNi] *= nondetection;
1010  migrated[INSNi] *= nondetection;
1011 
1012  migrated[ISPIi] *= nondetection;
1013  migrated[INSPIi] *= nondetection;
1014 
1015  migrated[ISyIIi] *= nondetection;
1016  migrated[INSIIi] *= nondetection;
1017  }
1018 }
1019 
1020 } // namespace osecirts
1021 } // namespace mio
1022 
1023 #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
Typesafe wrapper for a floating-point simulation time value (in days).
Definition: damping.h:78
mio::Simulation
A class for simulating a CompartmentalModel.
Definition: memilio/compartments/simulation.h:43
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:607
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:656
mio::osecirts::Model::deserialize
static IOResult< Model > deserialize(IOContext &io)
deserialize an object of this class.
Definition: ode_secirts/model.h:668
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:771
mio::osecirts::Simulation
specialization of compartment model simulation for the SECIRTS model.
Definition: ode_secirts/model.h:703
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:711
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:729
mio::osecirts::Simulation::advance
Eigen::Ref< Eigen::VectorX< FP > > advance(FP tmax)
advance simulation to tmax.
Definition: ode_secirts/model.h:756
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:966
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:919
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:863
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:881
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:890
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:408
mio::i
auto i
Definition: io.h:810
mio::log_warning
void log_warning(spdlog::string_view_t fmt, const Args &... args)
Definition: logging.h:126
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:482
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::floating_point_greater_equal
bool floating_point_greater_equal(FP v1, FP v2, FP abs_tol=0, FP rel_tol=std::numeric_limits< FP >::min())
compare two floating point values with tolerances.
Definition: floating_point.h:136
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:304
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:354
infection_state.h
parameters.h
mio::AgeGroup
Typesafe index representing an age group.
Definition: age_group.h:40
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:617
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