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