Line data Source code
1 : /*
2 : * * Copyright (C) 2020-2025 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 : namespace mio
33 : {
34 : namespace osecirts
35 : {
36 : // clang-format off
37 : using Flows = TypeList<
38 : //naive
39 : Flow<InfectionState::SusceptibleNaive, InfectionState::ExposedNaive>,
40 : Flow<InfectionState::SusceptibleNaive, InfectionState::TemporaryImmunePartialImmunity>,
41 : Flow<InfectionState::ExposedNaive, InfectionState::InfectedNoSymptomsNaive>,
42 : Flow<InfectionState::InfectedNoSymptomsNaive, InfectionState::InfectedSymptomsNaive>,
43 : Flow<InfectionState::InfectedNoSymptomsNaive, InfectionState::TemporaryImmunePartialImmunity>,
44 : Flow<InfectionState::InfectedNoSymptomsNaiveConfirmed, InfectionState::InfectedSymptomsNaiveConfirmed>,
45 : Flow<InfectionState::InfectedNoSymptomsNaiveConfirmed, InfectionState::TemporaryImmunePartialImmunity>,
46 : Flow<InfectionState::InfectedSymptomsNaive, InfectionState::InfectedSevereNaive>,
47 : Flow<InfectionState::InfectedSymptomsNaive, InfectionState::TemporaryImmunePartialImmunity>,
48 : Flow<InfectionState::InfectedSymptomsNaiveConfirmed, InfectionState::InfectedSevereNaive>,
49 : Flow<InfectionState::InfectedSymptomsNaiveConfirmed, InfectionState::TemporaryImmunePartialImmunity>,
50 : Flow<InfectionState::InfectedSevereNaive, InfectionState::InfectedCriticalNaive>,
51 : Flow<InfectionState::InfectedSevereNaive, InfectionState::TemporaryImmunePartialImmunity>,
52 : Flow<InfectionState::InfectedSevereNaive, InfectionState::DeadNaive>,
53 : Flow<InfectionState::InfectedCriticalNaive, InfectionState::DeadNaive>,
54 : Flow<InfectionState::InfectedCriticalNaive, InfectionState::TemporaryImmunePartialImmunity>,
55 : Flow<InfectionState::TemporaryImmunePartialImmunity, InfectionState::SusceptiblePartialImmunity>,
56 : //partial immunity
57 : Flow<InfectionState::SusceptiblePartialImmunity, InfectionState::ExposedPartialImmunity>,
58 : Flow<InfectionState::SusceptiblePartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
59 : Flow<InfectionState::ExposedPartialImmunity, InfectionState::InfectedNoSymptomsPartialImmunity>,
60 : Flow<InfectionState::InfectedNoSymptomsPartialImmunity, InfectionState::InfectedSymptomsPartialImmunity>,
61 : Flow<InfectionState::InfectedNoSymptomsPartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
62 : Flow<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, InfectionState::InfectedSymptomsPartialImmunityConfirmed>,
63 : Flow<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
64 : Flow<InfectionState::InfectedSymptomsPartialImmunity, InfectionState::InfectedSeverePartialImmunity>,
65 : Flow<InfectionState::InfectedSymptomsPartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
66 : Flow<InfectionState::InfectedSymptomsPartialImmunityConfirmed, InfectionState::InfectedSeverePartialImmunity>,
67 : Flow<InfectionState::InfectedSymptomsPartialImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
68 : Flow<InfectionState::InfectedSeverePartialImmunity, InfectionState::InfectedCriticalPartialImmunity>,
69 : Flow<InfectionState::InfectedSeverePartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
70 : Flow<InfectionState::InfectedSeverePartialImmunity, InfectionState::DeadPartialImmunity>,
71 : Flow<InfectionState::InfectedCriticalPartialImmunity, InfectionState::DeadPartialImmunity>,
72 : Flow<InfectionState::InfectedCriticalPartialImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
73 : //improved immunity
74 : Flow<InfectionState::SusceptibleImprovedImmunity, InfectionState::ExposedImprovedImmunity>,
75 : Flow<InfectionState::SusceptibleImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
76 : Flow<InfectionState::ExposedImprovedImmunity, InfectionState::InfectedNoSymptomsImprovedImmunity>,
77 : Flow<InfectionState::InfectedNoSymptomsImprovedImmunity, InfectionState::InfectedSymptomsImprovedImmunity>,
78 : Flow<InfectionState::InfectedNoSymptomsImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
79 : Flow<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed, InfectionState::InfectedSymptomsImprovedImmunityConfirmed>,
80 : Flow<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
81 : Flow<InfectionState::InfectedSymptomsImprovedImmunity, InfectionState::InfectedSevereImprovedImmunity>,
82 : Flow<InfectionState::InfectedSymptomsImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
83 : Flow<InfectionState::InfectedSymptomsImprovedImmunityConfirmed, InfectionState::InfectedSevereImprovedImmunity>,
84 : Flow<InfectionState::InfectedSymptomsImprovedImmunityConfirmed, InfectionState::TemporaryImmuneImprovedImmunity>,
85 : Flow<InfectionState::InfectedSevereImprovedImmunity, InfectionState::InfectedCriticalImprovedImmunity>,
86 : Flow<InfectionState::InfectedSevereImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
87 : Flow<InfectionState::InfectedSevereImprovedImmunity, InfectionState::DeadImprovedImmunity>,
88 : Flow<InfectionState::InfectedCriticalImprovedImmunity, InfectionState::DeadImprovedImmunity>,
89 : Flow<InfectionState::InfectedCriticalImprovedImmunity, InfectionState::TemporaryImmuneImprovedImmunity>,
90 :
91 : // waning
92 : Flow<InfectionState::TemporaryImmunePartialImmunity, InfectionState::SusceptiblePartialImmunity>,
93 : Flow<InfectionState::TemporaryImmuneImprovedImmunity, InfectionState::SusceptibleImprovedImmunity>,
94 : Flow<InfectionState::SusceptibleImprovedImmunity, InfectionState::SusceptiblePartialImmunity>,
95 : Flow<InfectionState::SusceptiblePartialImmunity, InfectionState::SusceptibleNaive>>;
96 : // clang-format on
97 :
98 : template <typename FP = ScalarType>
99 : class Model
100 : : public FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>
101 : {
102 : using Base = FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>;
103 :
104 : public:
105 : using typename Base::ParameterSet;
106 : using typename Base::Populations;
107 :
108 208 : Model(const Populations& pop, const ParameterSet& params)
109 208 : : Base(pop, params)
110 : {
111 208 : }
112 :
113 208 : Model(int num_agegroups)
114 208 : : Model(Populations({AgeGroup(num_agegroups), InfectionState::Count}), ParameterSet(AgeGroup(num_agegroups)))
115 : {
116 208 : }
117 :
118 3933 : void get_flows(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
119 : Eigen::Ref<Eigen::VectorX<FP>> flows) const override
120 : {
121 3933 : auto const& params = this->parameters;
122 3933 : AgeGroup n_agegroups = params.get_num_groups();
123 :
124 3933 : ContactMatrixGroup const& contact_matrix = params.template get<ContactPatterns<FP>>();
125 :
126 3933 : auto icu_occupancy = 0.0;
127 3933 : auto test_and_trace_required = 0.0;
128 11376 : for (auto i = AgeGroup(0); i < n_agegroups; ++i) {
129 : // naive flow to symptomatic in unit time
130 7443 : test_and_trace_required +=
131 14886 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
132 7443 : params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
133 14886 : (this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsNaive}) +
134 7443 : this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsNaiveConfirmed}));
135 : // partial immunity flow to symptomatic in unit time
136 7443 : test_and_trace_required +=
137 14886 : (params.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[i] /
138 7443 : params.template get<ReducExposedPartialImmunity<FP>>()[i]) *
139 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
140 14886 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
141 7443 : params.template get<ReducTimeInfectedMild<FP>>()[i]) *
142 14886 : (this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsPartialImmunity}) +
143 7443 : this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}));
144 : // improved immunity flow to symptomatic in unit time
145 7443 : test_and_trace_required +=
146 14886 : (params.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[i] /
147 7443 : params.template get<ReducExposedImprovedImmunity<FP>>()[i]) *
148 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
149 14886 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
150 7443 : params.template get<ReducTimeInfectedMild<FP>>()[i]) *
151 14886 : (this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsImprovedImmunity}) +
152 7443 : this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}));
153 14886 : icu_occupancy += this->populations.get_from(pop, {i, InfectionState::InfectedCriticalNaive}) +
154 14886 : this->populations.get_from(pop, {i, InfectionState::InfectedCriticalPartialImmunity}) +
155 7443 : this->populations.get_from(pop, {i, InfectionState::InfectedCriticalImprovedImmunity});
156 : }
157 :
158 : // get vaccinations
159 3933 : auto const partial_vaccination = vaccinations_at(t, params.template get<DailyPartialVaccinations<FP>>());
160 3933 : auto const full_vaccination = vaccinations_at(t, params.template get<DailyFullVaccinations<FP>>());
161 3933 : auto const booster_vaccination = vaccinations_at(t, params.template get<DailyBoosterVaccinations<FP>>());
162 :
163 11376 : for (auto i = AgeGroup(0); i < n_agegroups; i++) {
164 :
165 7443 : size_t SNi = this->populations.get_flat_index({i, InfectionState::SusceptibleNaive});
166 7443 : size_t ENi = this->populations.get_flat_index({i, InfectionState::ExposedNaive});
167 7443 : size_t INSNi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaive});
168 7443 : size_t ISyNi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsNaive});
169 7443 : size_t ISevNi = this->populations.get_flat_index({i, InfectionState::InfectedSevereNaive});
170 7443 : size_t ICrNi = this->populations.get_flat_index({i, InfectionState::InfectedCriticalNaive});
171 :
172 7443 : size_t INSNCi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaiveConfirmed});
173 7443 : size_t ISyNCi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsNaiveConfirmed});
174 :
175 7443 : size_t SPIi = this->populations.get_flat_index({i, InfectionState::SusceptiblePartialImmunity});
176 7443 : size_t EPIi = this->populations.get_flat_index({i, InfectionState::ExposedPartialImmunity});
177 7443 : size_t INSPIi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunity});
178 7443 : size_t ISyPIi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunity});
179 7443 : size_t ISevPIi = this->populations.get_flat_index({i, InfectionState::InfectedSeverePartialImmunity});
180 7443 : size_t ICrPIi = this->populations.get_flat_index({i, InfectionState::InfectedCriticalPartialImmunity});
181 :
182 : size_t INSPICi =
183 7443 : this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
184 : size_t ISyPICi =
185 7443 : this->populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
186 :
187 7443 : size_t EIIi = this->populations.get_flat_index({i, InfectionState::ExposedImprovedImmunity});
188 7443 : size_t INSIIi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunity});
189 7443 : size_t ISyIIi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunity});
190 7443 : size_t ISevIIi = this->populations.get_flat_index({i, InfectionState::InfectedSevereImprovedImmunity});
191 7443 : size_t ICrIIi = this->populations.get_flat_index({i, InfectionState::InfectedCriticalImprovedImmunity});
192 :
193 : size_t INSIICi =
194 7443 : this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
195 : size_t ISyIICi =
196 7443 : this->populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
197 7443 : size_t TImm1 = this->populations.get_flat_index({i, InfectionState::TemporaryImmunePartialImmunity});
198 7443 : size_t TImm2 = this->populations.get_flat_index({i, InfectionState::TemporaryImmuneImprovedImmunity});
199 :
200 7443 : size_t SIIi = this->populations.get_flat_index({i, InfectionState::SusceptibleImprovedImmunity});
201 :
202 7443 : FP reducExposedPartialImmunity = params.template get<ReducExposedPartialImmunity<FP>>()[i];
203 7443 : FP reducExposedImprovedImmunity = params.template get<ReducExposedImprovedImmunity<FP>>()[i];
204 7443 : FP reducInfectedSymptomsPartialImmunity =
205 7443 : params.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[i];
206 7443 : FP reducInfectedSymptomsImprovedImmunity =
207 7443 : params.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[i];
208 7443 : FP reducInfectedSevereCriticalDeadPartialImmunity =
209 7443 : params.template get<ReducInfectedSevereCriticalDeadPartialImmunity<FP>>()[i];
210 7443 : FP reducInfectedSevereCriticalDeadImprovedImmunity =
211 7443 : params.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<FP>>()[i];
212 7443 : FP reducTimeInfectedMild = params.template get<ReducTimeInfectedMild<FP>>()[i];
213 :
214 : //symptomatic are less well quarantined when testing and tracing is overwhelmed so they infect more people
215 : auto riskFromInfectedSymptomatic =
216 29772 : smoother_cosine(test_and_trace_required, params.template get<TestAndTraceCapacity<FP>>(),
217 14886 : params.template get<TestAndTraceCapacity<FP>>() *
218 7443 : params.template get<TestAndTraceCapacityMaxRiskSymptoms<FP>>(),
219 7443 : params.template get<RiskOfInfectionFromSymptomatic<FP>>()[i],
220 7443 : params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i]);
221 :
222 : auto riskFromInfectedNoSymptoms =
223 22329 : smoother_cosine(test_and_trace_required, params.template get<TestAndTraceCapacity<FP>>(),
224 14886 : params.template get<TestAndTraceCapacity<FP>>() *
225 7443 : params.template get<TestAndTraceCapacityMaxRiskNoSymptoms<FP>>(),
226 7443 : params.template get<RelativeTransmissionNoSymptoms<FP>>()[i], 1.0);
227 :
228 21906 : for (auto j = AgeGroup(0); j < n_agegroups; j++) {
229 14463 : size_t SNj = this->populations.get_flat_index({j, InfectionState::SusceptibleNaive});
230 14463 : size_t ENj = this->populations.get_flat_index({j, InfectionState::ExposedNaive});
231 14463 : size_t INSNj = this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsNaive});
232 14463 : size_t ISyNj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsNaive});
233 14463 : size_t ISevNj = this->populations.get_flat_index({j, InfectionState::InfectedSevereNaive});
234 14463 : size_t ICrNj = this->populations.get_flat_index({j, InfectionState::InfectedCriticalNaive});
235 14463 : size_t SIIj = this->populations.get_flat_index({j, InfectionState::SusceptibleImprovedImmunity});
236 :
237 14463 : size_t INSNCj = this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsNaiveConfirmed});
238 14463 : size_t ISyNCj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsNaiveConfirmed});
239 :
240 14463 : size_t SPIj = this->populations.get_flat_index({j, InfectionState::SusceptiblePartialImmunity});
241 14463 : size_t EPIj = this->populations.get_flat_index({j, InfectionState::ExposedPartialImmunity});
242 : size_t INSPIj =
243 14463 : this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsPartialImmunity});
244 14463 : size_t ISyPIj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsPartialImmunity});
245 14463 : size_t ISevPIj = this->populations.get_flat_index({j, InfectionState::InfectedSeverePartialImmunity});
246 14463 : size_t ICrPIj = this->populations.get_flat_index({j, InfectionState::InfectedCriticalPartialImmunity});
247 :
248 : size_t INSPICj =
249 14463 : this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
250 : size_t ISyPICj =
251 14463 : this->populations.get_flat_index({j, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
252 :
253 14463 : size_t EIIj = this->populations.get_flat_index({j, InfectionState::ExposedImprovedImmunity});
254 : size_t INSIIj =
255 14463 : this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsImprovedImmunity});
256 14463 : size_t ISyIIj = this->populations.get_flat_index({j, InfectionState::InfectedSymptomsImprovedImmunity});
257 14463 : size_t ISevIIj = this->populations.get_flat_index({j, InfectionState::InfectedSevereImprovedImmunity});
258 14463 : size_t ICrIIj = this->populations.get_flat_index({j, InfectionState::InfectedCriticalImprovedImmunity});
259 :
260 : size_t INSIICj =
261 14463 : this->populations.get_flat_index({j, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
262 : size_t ISyIICj =
263 14463 : this->populations.get_flat_index({j, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
264 :
265 : // effective contact rate by contact rate between groups i and j and damping j
266 14463 : FP season_val = (1 + params.template get<Seasonality<FP>>() *
267 14463 : sin(3.141592653589793 *
268 14463 : (std::fmod((params.template get<StartDay>() + t), 365.0) / 182.5 + 0.5)));
269 28926 : FP cont_freq_eff = season_val * contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
270 14463 : static_cast<Eigen::Index>((size_t)j));
271 : // without died people
272 14463 : FP Nj = pop[SNj] + pop[ENj] + pop[INSNj] + pop[ISyNj] + pop[ISevNj] + pop[ICrNj] + pop[INSNCj] +
273 14463 : pop[ISyNCj] + pop[SPIj] + pop[EPIj] + pop[INSPIj] + pop[ISyPIj] + pop[ISevPIj] + pop[ICrPIj] +
274 14463 : pop[INSPICj] + pop[ISyPICj] + pop[SIIj] + pop[EIIj] + pop[INSIIj] + pop[ISyIIj] + pop[ISevIIj] +
275 14463 : pop[ICrIIj] + pop[INSIICj] + pop[ISyIICj];
276 :
277 14463 : const FP divNj = (Nj < Limits<FP>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
278 :
279 14463 : FP ext_inf_force_dummy = cont_freq_eff * divNj *
280 14463 : params.template get<TransmissionProbabilityOnContact<FP>>()[(AgeGroup)i] *
281 14463 : (riskFromInfectedNoSymptoms * (pop[INSNj] + pop[INSPIj] + pop[INSIIj]) +
282 14463 : riskFromInfectedSymptomatic * (pop[ISyNj] + pop[ISyPIj] + pop[ISyIIj]));
283 :
284 14463 : FP dummy_SN = y[SNi] * ext_inf_force_dummy;
285 :
286 14463 : FP dummy_SPI = y[SPIi] * reducExposedPartialImmunity * ext_inf_force_dummy;
287 :
288 14463 : FP dummy_SII = y[SIIi] * reducExposedImprovedImmunity * ext_inf_force_dummy;
289 :
290 28926 : flows[this->template get_flat_flow_index<InfectionState::SusceptibleNaive,
291 14463 : InfectionState::ExposedNaive>({i})] += dummy_SN;
292 28926 : flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
293 14463 : InfectionState::ExposedPartialImmunity>({i})] += dummy_SPI;
294 28926 : flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
295 14463 : InfectionState::ExposedImprovedImmunity>({i})] += dummy_SII;
296 : }
297 :
298 : // vaccinations
299 14886 : flows[this->template get_flat_flow_index<InfectionState::SusceptibleNaive,
300 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
301 22329 : std::min(y[SNi] - flows[this->template get_flat_flow_index<InfectionState::SusceptibleNaive,
302 7443 : InfectionState::ExposedNaive>({i})],
303 7443 : partial_vaccination[static_cast<size_t>(i)]);
304 :
305 14886 : flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
306 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
307 22329 : std::min(y[SPIi] -
308 7443 : flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
309 7443 : InfectionState::ExposedPartialImmunity>({i})],
310 7443 : full_vaccination[static_cast<size_t>(i)]);
311 :
312 14886 : flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
313 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
314 22329 : std::min(y[SIIi] -
315 7443 : flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
316 7443 : InfectionState::ExposedImprovedImmunity>({i})],
317 7443 : booster_vaccination[static_cast<size_t>(i)]);
318 :
319 : // ICU capacity shortage is close
320 : // TODO: if this is used with vaccination model, it has to be adapted if CriticalPerSevere
321 : // is different for different vaccination status. This is not the case here and in addition, ICUCapacity
322 : // is set to infinity and this functionality is deactivated, so this is OK for the moment.
323 14886 : FP criticalPerSevereAdjusted = smoother_cosine(icu_occupancy, 0.90 * params.template get<ICUCapacity<FP>>(),
324 7443 : params.template get<ICUCapacity<FP>>(),
325 7443 : params.template get<CriticalPerSevere<FP>>()[i], 0);
326 :
327 7443 : FP deathsPerSevereAdjusted = params.template get<CriticalPerSevere<FP>>()[i] - criticalPerSevereAdjusted;
328 :
329 : /**** path of immune-naive ***/
330 : // Exposed
331 14886 : flows[this->template get_flat_flow_index<InfectionState::ExposedNaive,
332 7443 : InfectionState::InfectedNoSymptomsNaive>({i})] +=
333 7443 : y[ENi] / params.template get<TimeExposed<FP>>()[i];
334 :
335 : // InfectedNoSymptoms
336 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaive,
337 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
338 7443 : params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i] *
339 7443 : (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSNi];
340 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaive,
341 7443 : InfectionState::InfectedSymptomsNaive>({i})] =
342 14886 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
343 14886 : params.template get<TimeInfectedNoSymptoms<FP>>()[i] * y[INSNi];
344 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaiveConfirmed,
345 7443 : InfectionState::InfectedSymptomsNaiveConfirmed>({i})] =
346 14886 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
347 14886 : params.template get<TimeInfectedNoSymptoms<FP>>()[i] * y[INSNCi];
348 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsNaiveConfirmed,
349 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
350 7443 : params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i] *
351 7443 : (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSNCi];
352 :
353 : // // InfectedSymptoms
354 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaive,
355 7443 : InfectionState::InfectedSevereNaive>({i})] =
356 14886 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
357 14886 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNi];
358 :
359 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaive,
360 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
361 14886 : (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
362 14886 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNi];
363 :
364 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaiveConfirmed,
365 7443 : InfectionState::InfectedSevereNaive>({i})] =
366 14886 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
367 14886 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNCi];
368 :
369 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsNaiveConfirmed,
370 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
371 14886 : (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
372 14886 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyNCi];
373 :
374 : // InfectedSevere
375 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevereNaive,
376 7443 : InfectionState::InfectedCriticalNaive>({i})] =
377 7443 : criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
378 :
379 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevereNaive,
380 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
381 14886 : (1 - params.template get<CriticalPerSevere<FP>>()[i]) /
382 14886 : params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
383 :
384 7443 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevereNaive, InfectionState::DeadNaive>(
385 14886 : {i})] = deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevNi];
386 : // InfectedCritical
387 7443 : flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalNaive, InfectionState::DeadNaive>(
388 22329 : {i})] = params.template get<DeathsPerCritical<FP>>()[i] /
389 14886 : params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrNi];
390 :
391 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalNaive,
392 7443 : InfectionState::TemporaryImmunePartialImmunity>({i})] =
393 14886 : (1 - params.template get<DeathsPerCritical<FP>>()[i]) /
394 14886 : params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrNi];
395 :
396 : // Waning immunity
397 14886 : flows[this->template get_flat_flow_index<InfectionState::SusceptiblePartialImmunity,
398 7443 : InfectionState::SusceptibleNaive>({i})] =
399 7443 : 1 / params.template get<TimeWaningPartialImmunity<FP>>()[i] * y[SPIi];
400 14886 : flows[this->template get_flat_flow_index<InfectionState::TemporaryImmunePartialImmunity,
401 7443 : InfectionState::SusceptiblePartialImmunity>({i})] =
402 7443 : 1 / params.template get<TimeTemporaryImmunityPI<FP>>()[i] * y[TImm1];
403 :
404 : // /**** path of partially immune ***/
405 :
406 : // Exposed
407 14886 : flows[this->template get_flat_flow_index<InfectionState::ExposedPartialImmunity,
408 7443 : InfectionState::InfectedNoSymptomsPartialImmunity>({i})] +=
409 7443 : y[EPIi] / params.template get<TimeExposed<FP>>()[i];
410 :
411 : // InfectedNoSymptoms
412 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunity,
413 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
414 14886 : (1 - (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
415 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
416 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPIi];
417 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunity,
418 7443 : InfectionState::InfectedSymptomsPartialImmunity>({i})] =
419 14886 : (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
420 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
421 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPIi];
422 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed,
423 7443 : InfectionState::InfectedSymptomsPartialImmunityConfirmed>({i})] =
424 14886 : (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
425 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
426 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPICi];
427 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsPartialImmunityConfirmed,
428 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
429 14886 : (1 - (reducInfectedSymptomsPartialImmunity / reducExposedPartialImmunity) *
430 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
431 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSPICi];
432 :
433 : // InfectedSymptoms
434 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunity,
435 7443 : InfectionState::InfectedSeverePartialImmunity>({i})] =
436 7443 : reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity *
437 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
438 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPIi];
439 :
440 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunity,
441 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
442 7443 : (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity) *
443 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
444 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPIi];
445 :
446 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunityConfirmed,
447 7443 : InfectionState::InfectedSeverePartialImmunity>({i})] =
448 7443 : reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity *
449 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
450 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPICi];
451 :
452 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsPartialImmunityConfirmed,
453 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
454 7443 : (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSymptomsPartialImmunity) *
455 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
456 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyPICi];
457 :
458 : // InfectedSevere
459 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
460 7443 : InfectionState::InfectedCriticalPartialImmunity>({i})] =
461 7443 : reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity *
462 7443 : criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
463 :
464 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
465 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
466 7443 : (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
467 14886 : params.template get<CriticalPerSevere<FP>>()[i]) /
468 14886 : params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
469 :
470 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSeverePartialImmunity,
471 7443 : InfectionState::DeadPartialImmunity>({i})] =
472 7443 : (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
473 7443 : deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevPIi];
474 : // InfectedCritical
475 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalPartialImmunity,
476 7443 : InfectionState::DeadPartialImmunity>({i})] =
477 7443 : (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
478 14886 : params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
479 7443 : y[ICrPIi];
480 :
481 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalPartialImmunity,
482 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
483 7443 : (1 - (reducInfectedSevereCriticalDeadPartialImmunity / reducInfectedSevereCriticalDeadPartialImmunity) *
484 14886 : params.template get<DeathsPerCritical<FP>>()[i]) /
485 14886 : params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrPIi];
486 :
487 : // /**** path of improved immunity ***/
488 : // Exposed
489 14886 : flows[this->template get_flat_flow_index<InfectionState::ExposedImprovedImmunity,
490 7443 : InfectionState::InfectedNoSymptomsImprovedImmunity>({i})] +=
491 7443 : y[EIIi] / params.template get<TimeExposed<FP>>()[i];
492 :
493 : // InfectedNoSymptoms
494 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunity,
495 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
496 14886 : (1 - (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
497 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
498 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIIi];
499 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunity,
500 7443 : InfectionState::InfectedSymptomsImprovedImmunity>({i})] =
501 14886 : (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
502 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
503 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIIi];
504 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed,
505 7443 : InfectionState::InfectedSymptomsImprovedImmunityConfirmed>({i})] =
506 14886 : (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
507 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
508 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIICi];
509 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed,
510 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
511 14886 : (1 - (reducInfectedSymptomsImprovedImmunity / reducExposedImprovedImmunity) *
512 7443 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i])) /
513 7443 : (params.template get<TimeInfectedNoSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[INSIICi];
514 :
515 : // InfectedSymptoms
516 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunity,
517 7443 : InfectionState::InfectedSevereImprovedImmunity>({i})] =
518 7443 : reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity *
519 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
520 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIIi];
521 :
522 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunity,
523 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
524 7443 : (1 - (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity) *
525 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
526 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIIi];
527 :
528 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunityConfirmed,
529 7443 : InfectionState::InfectedSevereImprovedImmunity>({i})] =
530 7443 : reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity *
531 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
532 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIICi];
533 :
534 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsImprovedImmunityConfirmed,
535 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
536 7443 : (1 - (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSymptomsImprovedImmunity) *
537 7443 : params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
538 7443 : (params.template get<TimeInfectedSymptoms<FP>>()[i] * reducTimeInfectedMild) * y[ISyIICi];
539 :
540 : // InfectedSevere
541 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
542 7443 : InfectionState::InfectedCriticalImprovedImmunity>({i})] =
543 7443 : reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity *
544 7443 : criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
545 :
546 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
547 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
548 7443 : (1 -
549 7443 : (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
550 14886 : params.template get<CriticalPerSevere<FP>>()[i]) /
551 14886 : params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
552 :
553 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevereImprovedImmunity,
554 7443 : InfectionState::DeadImprovedImmunity>({i})] =
555 7443 : (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
556 7443 : deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevIIi];
557 :
558 : // InfectedCritical
559 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalImprovedImmunity,
560 7443 : InfectionState::DeadImprovedImmunity>({i})] =
561 7443 : (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
562 14886 : params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
563 7443 : y[ICrIIi];
564 :
565 14886 : flows[this->template get_flat_flow_index<InfectionState::InfectedCriticalImprovedImmunity,
566 7443 : InfectionState::TemporaryImmuneImprovedImmunity>({i})] =
567 7443 : (1 -
568 7443 : (reducInfectedSevereCriticalDeadImprovedImmunity / reducInfectedSevereCriticalDeadImprovedImmunity) *
569 14886 : params.template get<DeathsPerCritical<FP>>()[i]) /
570 14886 : params.template get<TimeInfectedCritical<FP>>()[i] * y[ICrIIi];
571 :
572 : // Waning immunity
573 14886 : flows[this->template get_flat_flow_index<InfectionState::TemporaryImmuneImprovedImmunity,
574 7443 : InfectionState::SusceptibleImprovedImmunity>({i})] =
575 7443 : 1 / params.template get<TimeTemporaryImmunityII<FP>>()[i] * y[TImm2];
576 :
577 14886 : flows[this->template get_flat_flow_index<InfectionState::SusceptibleImprovedImmunity,
578 7443 : InfectionState::SusceptiblePartialImmunity>({i})] =
579 7443 : 1 / params.template get<TimeWaningImprovedImmunity<FP>>()[i] * y[SIIi];
580 : }
581 7866 : }
582 :
583 : /**
584 : * @brief Calculates smoothed vaccinations for a given time point.
585 : *
586 : * This function calculates the number of vaccinations for each age group at a given time t,
587 : * based on daily vaccinations data. The smoothing is done using a cosine function.
588 : *
589 : * @param t The time in the simulation.
590 : * @param daily_vaccinations The daily vaccinations data, indexed by age group and simulation day.
591 : * @param eps [Default: 0.15] The smoothing factor used in the cosine smoothing function.
592 : * @return A vector containing the number of vaccinations for each age group at time t.
593 : */
594 11898 : Eigen::VectorXd vaccinations_at(const FP t,
595 : const CustomIndexArray<ScalarType, AgeGroup, SimulationDay>& daily_vaccinations,
596 : const FP eps = 0.15) const
597 : {
598 11898 : auto const& params = this->parameters;
599 11898 : const FP ub = (size_t)t + 1.0;
600 11898 : const FP lb = ub - eps;
601 :
602 11898 : const auto max_time = static_cast<size_t>(daily_vaccinations.size<SimulationDay>()) - 1;
603 :
604 11898 : Eigen::VectorXd smoothed_vaccinations((size_t)params.get_num_groups());
605 11898 : smoothed_vaccinations.setZero();
606 :
607 : // if daily_vaccinations is not available for the current time point, we return zero vaccinations.
608 11898 : if (max_time <= (size_t)t) {
609 9 : mio::log_warning("Vaccination data not available for time point ", t, ". Returning zero vaccinations.");
610 9 : return smoothed_vaccinations;
611 : }
612 11889 : if (t >= lb) {
613 10476 : for (AgeGroup age = AgeGroup(0); age < params.get_num_groups(); age++) {
614 : // if ub + 1 is out of bounds, we use the value at ub
615 6912 : auto ubp1 = static_cast<size_t>(ub + 1);
616 6912 : if (max_time < ubp1) {
617 9 : ubp1 = static_cast<size_t>(ub);
618 : }
619 20736 : const auto num_vaccinations_ub = daily_vaccinations[{age, SimulationDay(static_cast<size_t>(ubp1))}] -
620 13824 : daily_vaccinations[{age, SimulationDay(static_cast<size_t>(ub))}];
621 20736 : const auto num_vaccinations_lb = daily_vaccinations[{age, SimulationDay(static_cast<size_t>(lb + 1))}] -
622 13824 : daily_vaccinations[{age, SimulationDay(static_cast<size_t>(lb))}];
623 6912 : smoothed_vaccinations[(size_t)age] =
624 6912 : smoother_cosine(t, lb, ub, num_vaccinations_lb, num_vaccinations_ub);
625 : }
626 : }
627 : else {
628 23832 : for (auto age = AgeGroup(0); age < params.get_num_groups(); age++) {
629 46521 : smoothed_vaccinations[(size_t)age] = daily_vaccinations[{age, SimulationDay((size_t)t + 1)}] -
630 31014 : daily_vaccinations[{age, SimulationDay((size_t)t)}];
631 : }
632 : }
633 11889 : return smoothed_vaccinations;
634 11898 : }
635 :
636 : /**
637 : * serialize this.
638 : * @see mio::serialize
639 : */
640 : template <class IOContext>
641 : void serialize(IOContext& io) const
642 : {
643 : auto obj = io.create_object("Model");
644 : obj.add_element("Parameters", this->parameters);
645 : obj.add_element("Populations", this->populations);
646 : }
647 :
648 : /**
649 : * deserialize an object of this class.
650 : * @see mio::deserialize
651 : */
652 : template <class IOContext>
653 : static IOResult<Model> deserialize(IOContext& io)
654 : {
655 : auto obj = io.expect_object("Model");
656 : auto par = obj.expect_element("Parameters", Tag<ParameterSet>{});
657 : auto pop = obj.expect_element("Populations", Tag<Populations>{});
658 : return apply(
659 : io,
660 : [](auto&& par_, auto&& pop_) {
661 : return Model{pop_, par_};
662 : },
663 : par, pop);
664 : }
665 : }; // namespace osecirts
666 :
667 : //forward declaration, see below.
668 : template <typename FP = ScalarType, class Base = mio::Simulation<FP, Model<FP>>>
669 : class Simulation;
670 :
671 : /**
672 : * get percentage of infections per total population.
673 : * @param model the compartment model with initial values.
674 : * @param t current simulation time.
675 : * @param y current value of compartments.
676 : * @tparam Base simulation type that uses the SECIRS-type compartment model. see Simulation.
677 : */
678 : template <typename FP = ScalarType, class Base = mio::Simulation<FP, Model<FP>>>
679 : FP get_infections_relative(const Simulation<FP, Base>& model, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y);
680 :
681 : /**
682 : * specialization of compartment model simulation for the SECIRTS model.
683 : * @tparam FP floating point type, e.g., double.
684 : * @tparam BaseT simulation type, default mio::Simulation. For testing purposes only!
685 : */
686 : template <typename FP, class BaseT>
687 : class Simulation : public BaseT
688 : {
689 : public:
690 : /**
691 : * construct a simulation.
692 : * @param model the model to simulate.
693 : * @param t0 start time
694 : * @param dt time steps
695 : */
696 63 : Simulation(mio::osecirts::Model<FP> const& model, FP t0 = 0., FP dt = 0.1)
697 : : BaseT(model, t0, dt)
698 63 : , m_t_last_npi_check(t0)
699 : {
700 63 : }
701 :
702 : /**
703 : * @brief Applies the effect of a new variant of a disease to the transmission probability of the model.
704 : *
705 : * This function adjusts the transmission probability of the disease for each age group based on the share of the new variant.
706 : * The share of the new variant is calculated based on the time `t` and the start day of the new variant.
707 : * The transmission probability is then updated for each age group in the model.
708 : *
709 : * Based on Equation (35) and (36) in doi.org/10.1371/journal.pcbi.1010054
710 : *
711 : * @param [in] t The current time.
712 : * @param [in] base_infectiousness The base infectiousness of the old variant for each age group.
713 : */
714 :
715 333 : void apply_variant(const FP t, const CustomIndexArray<UncertainValue<FP>, AgeGroup> base_infectiousness)
716 : {
717 333 : auto start_day = this->get_model().parameters.template get<StartDay>();
718 333 : auto start_day_new_variant = this->get_model().parameters.template get<StartDayNewVariant>();
719 :
720 333 : if (start_day + t >= start_day_new_variant - 1e-10) {
721 27 : const FP days_variant = t - (start_day_new_variant - start_day);
722 27 : const FP share_new_variant = std::min(1.0, 0.01 * pow(2, (1. / 7) * days_variant));
723 27 : const auto num_groups = this->get_model().parameters.get_num_groups();
724 54 : for (auto i = AgeGroup(0); i < num_groups; ++i) {
725 27 : FP new_transmission = (1 - share_new_variant) * base_infectiousness[i] +
726 27 : share_new_variant * base_infectiousness[i] *
727 27 : this->get_model().parameters.template get<InfectiousnessNewVariant<FP>>()[i];
728 27 : this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>()[i] = new_transmission;
729 : }
730 : }
731 333 : }
732 :
733 : /**
734 : * @brief advance simulation to tmax.
735 : * Overwrites Simulation::advance and includes a check for dynamic NPIs in regular intervals.
736 : * @see Simulation::advance
737 : * @param tmax next stopping point of simulation
738 : * @return value at tmax
739 : */
740 36 : Eigen::Ref<Eigen::VectorXd> advance(FP tmax)
741 : {
742 36 : auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis();
743 36 : auto& dyn_npis = this->get_model().parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
744 36 : auto& contact_patterns = this->get_model().parameters.template get<ContactPatterns<FP>>();
745 : // const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups();
746 :
747 : // in the apply_variant function, we adjust the TransmissionProbabilityOnContact parameter. We need to store
748 : // the base value to use it in the apply_variant function and also to reset the parameter after the simulation.
749 36 : auto base_infectiousness = this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>();
750 :
751 : FP delay_npi_implementation;
752 36 : auto t = BaseT::get_result().get_last_time();
753 36 : const auto dt = dyn_npis.get_interval().get();
754 333 : while (t < tmax) {
755 :
756 297 : auto dt_eff = std::min({dt, tmax - t, m_t_last_npi_check + dt - t});
757 297 : if (dt_eff >= 1.0) {
758 288 : dt_eff = 1.0;
759 : }
760 :
761 297 : BaseT::advance(t + dt_eff);
762 297 : if (t + 0.5 + dt_eff - std::floor(t + 0.5) >= 1) {
763 288 : this->apply_variant(t, base_infectiousness);
764 : }
765 :
766 297 : if (t > 0) {
767 261 : delay_npi_implementation =
768 261 : this->get_model().parameters.template get<DynamicNPIsImplementationDelay<FP>>();
769 : }
770 : else {
771 : // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay.
772 36 : delay_npi_implementation = 0;
773 : }
774 297 : t = t + dt_eff;
775 :
776 297 : if (dyn_npis.get_thresholds().size() > 0) {
777 270 : if (floating_point_greater_equal(t, m_t_last_npi_check + dt)) {
778 90 : if (t < t_end_dyn_npis) {
779 54 : auto inf_rel = get_infections_relative<FP>(*this, t, this->get_result().get_last_value()) *
780 27 : dyn_npis.get_base_value();
781 27 : auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel);
782 54 : if (exceeded_threshold != dyn_npis.get_thresholds().end() &&
783 45 : (exceeded_threshold->first > m_dynamic_npi.first ||
784 18 : t > FP(m_dynamic_npi.second))) { //old npi was weaker or is expired
785 :
786 9 : auto t_start = SimulationTime(t + delay_npi_implementation);
787 9 : auto t_end = t_start + SimulationTime(dyn_npis.get_duration());
788 9 : this->get_model().parameters.get_start_commuter_detection() = (FP)t_start;
789 9 : this->get_model().parameters.get_end_commuter_detection() = (FP)t_end;
790 9 : m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
791 9 : implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
792 9 : t_start, t_end, [](auto& g) {
793 9 : return make_contact_damping_matrix(g);
794 : });
795 : }
796 : }
797 90 : m_t_last_npi_check = t;
798 : }
799 : }
800 : else {
801 27 : m_t_last_npi_check = t;
802 : }
803 : }
804 : // reset TransmissionProbabilityOnContact. This is important for the graph simulation where the advance
805 : // function is called multiple times for the same model.
806 36 : this->get_model().parameters.template get<TransmissionProbabilityOnContact<FP>>() = base_infectiousness;
807 :
808 72 : return this->get_result().get_last_value();
809 36 : }
810 :
811 : private:
812 : FP m_t_last_npi_check;
813 : std::pair<FP, SimulationTime> m_dynamic_npi = {-std::numeric_limits<FP>::max(), SimulationTime(0)};
814 : };
815 :
816 : /**
817 : * @brief Specialization of simulate for SECIRS-type models using Simulation.
818 : *
819 : * @param[in] t0 start time.
820 : * @param[in] tmax end time.
821 : * @param[in] dt time step.
822 : * @param[in] model SECIRS-type model to simulate.
823 : * @param[in] integrator optional integrator, uses rk45 if nullptr.
824 : *
825 : * @return Returns the result of the simulation.
826 : */
827 : template <typename FP = ScalarType>
828 9 : inline auto simulate(FP t0, FP tmax, FP dt, const Model<FP>& model,
829 : std::shared_ptr<IntegratorCore<FP>> integrator = nullptr)
830 : {
831 9 : return mio::simulate<FP, Model<FP>, Simulation<FP>>(t0, tmax, dt, model, integrator);
832 : }
833 :
834 : /**
835 : * @brief Specialization of simulate for SECIRS-type models using the FlowSimulation.
836 : *
837 : * @param[in] t0 start time.
838 : * @param[in] tmax end time.
839 : * @param[in] dt time step.
840 : * @param[in] model SECIRS-type model to simulate.
841 : * @param[in] integrator optional integrator, uses rk45 if nullptr.
842 : *
843 : * @return Returns the result of the Flowsimulation.
844 : */
845 : template <typename FP = ScalarType>
846 : inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model<FP>& model,
847 : std::shared_ptr<IntegratorCore<FP>> integrator = nullptr)
848 : {
849 : return mio::simulate_flows<FP, Model<FP>, Simulation<FP, mio::FlowSimulation<FP, Model<FP>>>>(t0, tmax, dt, model,
850 : integrator);
851 : }
852 :
853 : //see declaration above.
854 : template <typename FP, class Base>
855 36 : FP get_infections_relative(const Simulation<FP, Base>& sim, FP /*t*/, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
856 : {
857 36 : FP sum_inf = 0;
858 108 : for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) {
859 72 : sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsNaive});
860 72 : sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsNaiveConfirmed});
861 72 : sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsPartialImmunity});
862 72 : sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsImprovedImmunity});
863 72 : sum_inf +=
864 72 : sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
865 72 : sum_inf +=
866 72 : sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
867 : }
868 36 : auto inf_rel = sum_inf / sim.get_model().populations.get_total();
869 :
870 36 : return inf_rel;
871 : }
872 :
873 : /**
874 : * Get migration factors.
875 : * Used by migration graph simulation.
876 : * Like infection risk, migration of infected individuals is reduced if they are well isolated.
877 : * @param model the compartment model with initial values.
878 : * @param t current simulation time.
879 : * @param y current value of compartments.
880 : * @return vector expression, same size as y, with migration factors per compartment.
881 : * @tparam Base simulation type that uses a SECIRS-type compartment model. see Simulation.
882 : */
883 : template <typename FP = double, class Base = mio::Simulation<Model<FP>, FP>>
884 9 : auto get_migration_factors(const Simulation<Base>& sim, FP /*t*/, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
885 : {
886 9 : auto& params = sim.get_model().parameters;
887 : //parameters as arrays
888 9 : auto& p_asymp = params.template get<RecoveredPerInfectedNoSymptoms<FP>>().array().template cast<double>();
889 9 : auto& p_inf = params.template get<RiskOfInfectionFromSymptomatic<FP>>().array().template cast<double>();
890 9 : auto& p_inf_max = params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>().array().template cast<double>();
891 : //slice of InfectedNoSymptoms
892 54 : auto y_INS = slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsNaive),
893 18 : Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)}) +
894 : slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsPartialImmunity),
895 18 : Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)}) +
896 : slice(y, {Eigen::Index(InfectionState::InfectedNoSymptomsImprovedImmunity),
897 18 : Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)});
898 :
899 : //compute isolation, same as infection risk from main model
900 9 : auto test_and_trace_required =
901 27 : ((1 - p_asymp) / params.template get<TimeInfectedNoSymptoms<FP>>().array().template cast<double>() *
902 : y_INS.array())
903 36 : .sum();
904 27 : auto riskFromInfectedSymptomatic =
905 9 : smoother_cosine(test_and_trace_required, double(params.template get<TestAndTraceCapacity<FP>>()),
906 9 : params.template get<TestAndTraceCapacity<FP>>() * 5, p_inf.matrix(), p_inf_max.matrix());
907 :
908 : //set factor for infected
909 9 : auto factors = Eigen::VectorXd::Ones(y.rows()).eval();
910 27 : slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsNaive), Eigen::Index(size_t(params.get_num_groups())),
911 : Eigen::Index(InfectionState::Count)})
912 18 : .array() = riskFromInfectedSymptomatic;
913 9 : slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsPartialImmunity),
914 18 : Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)})
915 18 : .array() = riskFromInfectedSymptomatic;
916 9 : slice(factors, {Eigen::Index(InfectionState::InfectedSymptomsImprovedImmunity),
917 18 : Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)})
918 18 : .array() = riskFromInfectedSymptomatic;
919 18 : return factors;
920 9 : }
921 :
922 : /**
923 : * @brief Adjusts the state of commuters in a model, accounting for detection and mobility effects.
924 : *
925 : * @tparam FP Floating-point type, e.g., double.
926 : * @tparam Base Simulation type that uses the SECIRTS-type model. see Simulation.
927 : * @param[in,out] sim Simulation object containing the model and result data.
928 : * @param[in,out] migrated Vector representing the number of commuters in each state.
929 : * @param[in] time Current simulation time, used to determine the commuter detection period.
930 : */
931 : template <typename FP = double, class Base = mio::Simulation<Model<FP>, FP>>
932 9 : auto test_commuters(Simulation<FP, Base>& sim, Eigen::Ref<Eigen::VectorX<FP>> migrated, FP time)
933 : {
934 9 : auto& model = sim.get_model();
935 9 : auto nondetection = 1.0;
936 18 : if (time >= model.parameters.get_start_commuter_detection() &&
937 9 : time < model.parameters.get_end_commuter_detection()) {
938 9 : nondetection = (double)model.parameters.get_commuter_nondetection();
939 : }
940 27 : for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); ++i) {
941 18 : auto ISyNi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsNaive});
942 18 : auto ISyNCi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsNaiveConfirmed});
943 18 : auto INSNi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaive});
944 18 : auto INSNCi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsNaiveConfirmed});
945 :
946 18 : auto ISPIi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunity});
947 18 : auto ISPICi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsPartialImmunityConfirmed});
948 18 : auto INSPIi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunity});
949 : auto INSPICi =
950 18 : model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed});
951 :
952 18 : auto ISyIIi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunity});
953 18 : auto ISyIICi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed});
954 18 : auto INSIIi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunity});
955 : auto INSIICi =
956 18 : model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed});
957 :
958 : //put detected commuters in their own compartment so they don't contribute to infections in their home node
959 18 : sim.get_result().get_last_value()[ISyNi] -= migrated[ISyNi] * (1 - nondetection);
960 18 : sim.get_result().get_last_value()[ISyNCi] += migrated[ISyNi] * (1 - nondetection);
961 18 : sim.get_result().get_last_value()[INSNi] -= migrated[INSNi] * (1 - nondetection);
962 18 : sim.get_result().get_last_value()[INSNCi] += migrated[INSNi] * (1 - nondetection);
963 :
964 18 : sim.get_result().get_last_value()[ISPIi] -= migrated[ISPIi] * (1 - nondetection);
965 18 : sim.get_result().get_last_value()[ISPICi] += migrated[ISPIi] * (1 - nondetection);
966 18 : sim.get_result().get_last_value()[INSPIi] -= migrated[INSPIi] * (1 - nondetection);
967 18 : sim.get_result().get_last_value()[INSPICi] += migrated[INSPIi] * (1 - nondetection);
968 :
969 18 : sim.get_result().get_last_value()[ISyIIi] -= migrated[ISyIIi] * (1 - nondetection);
970 18 : sim.get_result().get_last_value()[ISyIICi] += migrated[ISyIIi] * (1 - nondetection);
971 18 : sim.get_result().get_last_value()[INSIIi] -= migrated[INSIIi] * (1 - nondetection);
972 18 : sim.get_result().get_last_value()[INSIICi] += migrated[INSIIi] * (1 - nondetection);
973 :
974 : //reduce the number of commuters
975 18 : migrated[ISyNi] *= nondetection;
976 18 : migrated[INSNi] *= nondetection;
977 :
978 18 : migrated[ISPIi] *= nondetection;
979 18 : migrated[INSPIi] *= nondetection;
980 :
981 18 : migrated[ISyIIi] *= nondetection;
982 18 : migrated[INSIIi] *= nondetection;
983 : }
984 9 : }
985 :
986 : } // namespace osecirts
987 : } // namespace mio
988 :
989 : #endif //MIO_ODE_SECIRTS_MODEL_H
|
