Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 MEmilio
3 : *
4 : * Authors: Daniel Abele, Jan Kleinert, Martin J. Kuehn
5 : *
6 : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
7 : *
8 : * Licensed under the Apache License, Version 2.0 (the "License");
9 : * you may not use this file except in compliance with the License.
10 : * You may obtain a copy of the License at
11 : *
12 : * http://www.apache.org/licenses/LICENSE-2.0
13 : *
14 : * Unless required by applicable law or agreed to in writing, software
15 : * distributed under the License is distributed on an "AS IS" BASIS,
16 : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 : * See the License for the specific language governing permissions and
18 : * limitations under the License.
19 : */
20 : #ifndef ODESECIR_MODEL_H
21 : #define ODESECIR_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_secir/analyze_result.h"
28 : #include "ode_secir/infection_state.h"
29 : #include "ode_secir/parameters.h"
30 : #include "memilio/math/smoother.h"
31 : #include "memilio/math/eigen_util.h"
32 : #include "memilio/math/interpolation.h"
33 :
34 : namespace mio
35 : {
36 : namespace osecir
37 : {
38 :
39 : // Create template specializations for the age resolved
40 : // SECIHURD model
41 : // clang-format off
42 : using Flows = TypeList<Flow<InfectionState::Susceptible, InfectionState::Exposed>,
43 : Flow<InfectionState::Exposed, InfectionState::InfectedNoSymptoms>,
44 : Flow<InfectionState::InfectedNoSymptoms, InfectionState::InfectedSymptoms>,
45 : Flow<InfectionState::InfectedNoSymptoms, InfectionState::Recovered>,
46 : Flow<InfectionState::InfectedNoSymptomsConfirmed, InfectionState::InfectedSymptomsConfirmed>,
47 : Flow<InfectionState::InfectedNoSymptomsConfirmed, InfectionState::Recovered>,
48 : Flow<InfectionState::InfectedSymptoms, InfectionState::InfectedSevere>,
49 : Flow<InfectionState::InfectedSymptoms, InfectionState::Recovered>,
50 : Flow<InfectionState::InfectedSymptomsConfirmed, InfectionState::InfectedSevere>,
51 : Flow<InfectionState::InfectedSymptomsConfirmed, InfectionState::Recovered>,
52 : Flow<InfectionState::InfectedSevere, InfectionState::InfectedCritical>,
53 : Flow<InfectionState::InfectedSevere, InfectionState::Recovered>,
54 : Flow<InfectionState::InfectedSevere, InfectionState::Dead>,
55 : Flow<InfectionState::InfectedCritical, InfectionState::Dead>,
56 : Flow<InfectionState::InfectedCritical, InfectionState::Recovered>>;
57 : // clang-format on
58 :
59 : template <typename FP = ScalarType>
60 : class Model : public FlowModel<FP, InfectionState, Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>
61 : {
62 : using Base = FlowModel<FP, InfectionState, mio::Populations<FP, AgeGroup, InfectionState>, Parameters<FP>, Flows>;
63 :
64 : public:
65 : using typename Base::ParameterSet;
66 : using typename Base::Populations;
67 :
68 371 : Model(const Populations& pop, const ParameterSet& params)
69 371 : : Base(pop, params)
70 : {
71 371 : }
72 :
73 362 : Model(int num_agegroups)
74 724 : : Model(Populations({mio::AgeGroup(num_agegroups), InfectionState::Count}),
75 1086 : ParameterSet(mio::AgeGroup(num_agegroups)))
76 : {
77 362 : }
78 :
79 411453 : void get_flows(Eigen::Ref<const Vector<FP>> pop, Eigen::Ref<const Vector<FP>> y, FP t,
80 : Eigen::Ref<Vector<FP>> flows) const override
81 : {
82 411453 : auto const& params = this->parameters;
83 411453 : AgeGroup n_agegroups = params.get_num_groups();
84 :
85 411453 : ContactMatrixGroup const& contact_matrix = params.template get<ContactPatterns<FP>>();
86 :
87 411453 : auto icu_occupancy = 0.0;
88 411453 : auto test_and_trace_required = 0.0;
89 835098 : for (auto i = AgeGroup(0); i < n_agegroups; ++i) {
90 847290 : test_and_trace_required += (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
91 423645 : params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
92 423645 : this->populations.get_from(pop, {i, InfectionState::InfectedNoSymptoms});
93 423645 : icu_occupancy += this->populations.get_from(pop, {i, InfectionState::InfectedCritical});
94 : }
95 :
96 835098 : for (auto i = AgeGroup(0); i < n_agegroups; i++) {
97 :
98 423645 : size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible});
99 423645 : size_t Ei = this->populations.get_flat_index({i, InfectionState::Exposed});
100 423645 : size_t INSi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptoms});
101 423645 : size_t INSCi = this->populations.get_flat_index({i, InfectionState::InfectedNoSymptomsConfirmed});
102 423645 : size_t ISyi = this->populations.get_flat_index({i, InfectionState::InfectedSymptoms});
103 423645 : size_t ISyCi = this->populations.get_flat_index({i, InfectionState::InfectedSymptomsConfirmed});
104 423645 : size_t ISevi = this->populations.get_flat_index({i, InfectionState::InfectedSevere});
105 423645 : size_t ICri = this->populations.get_flat_index({i, InfectionState::InfectedCritical});
106 :
107 883866 : for (auto j = AgeGroup(0); j < n_agegroups; j++) {
108 460221 : size_t Sj = this->populations.get_flat_index({j, InfectionState::Susceptible});
109 460221 : size_t Ej = this->populations.get_flat_index({j, InfectionState::Exposed});
110 460221 : size_t INSj = this->populations.get_flat_index({j, InfectionState::InfectedNoSymptoms});
111 460221 : size_t ISyj = this->populations.get_flat_index({j, InfectionState::InfectedSymptoms});
112 460221 : size_t ISevj = this->populations.get_flat_index({j, InfectionState::InfectedSevere});
113 460221 : size_t ICrj = this->populations.get_flat_index({j, InfectionState::InfectedCritical});
114 460221 : size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
115 :
116 : //symptomatic are less well quarantined when testing and tracing is overwhelmed so they infect more people
117 : auto riskFromInfectedSymptomatic =
118 1840884 : smoother_cosine(test_and_trace_required, params.template get<TestAndTraceCapacity<FP>>(),
119 920442 : params.template get<TestAndTraceCapacity<FP>>() *
120 460221 : params.template get<TestAndTraceCapacityMaxRisk<FP>>(),
121 460221 : params.template get<RiskOfInfectionFromSymptomatic<FP>>()[j],
122 460221 : params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[j]);
123 :
124 : // effective contact rate by contact rate between groups i and j and damping j
125 460221 : ScalarType season_val =
126 460221 : (1 + params.template get<Seasonality<FP>>() *
127 460221 : sin(3.141592653589793 * ((params.template get<StartDay>() + t) / 182.5 + 0.5)));
128 460221 : ScalarType cont_freq_eff =
129 920442 : season_val * contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
130 460221 : static_cast<Eigen::Index>((size_t)j));
131 460221 : ScalarType Nj =
132 460221 : pop[Sj] + pop[Ej] + pop[INSj] + pop[ISyj] + pop[ISevj] + pop[ICrj] + pop[Rj]; // without died people
133 460221 : const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
134 920442 : ScalarType dummy_S = y[Si] * cont_freq_eff * divNj *
135 460221 : params.template get<TransmissionProbabilityOnContact<FP>>()[i] *
136 460221 : (params.template get<RelativeTransmissionNoSymptoms<FP>>()[j] * pop[INSj] +
137 460221 : riskFromInfectedSymptomatic * pop[ISyj]);
138 :
139 : // Susceptible -> Exposed
140 460221 : flows[this->template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Exposed>({i})] +=
141 : dummy_S;
142 : }
143 :
144 : // ICU capacity shortage is close
145 1270935 : ScalarType criticalPerSevereAdjusted = smoother_cosine(
146 847290 : icu_occupancy, 0.90 * params.template get<ICUCapacity<FP>>(), params.template get<ICUCapacity<FP>>(),
147 423645 : params.template get<CriticalPerSevere<FP>>()[i], 0);
148 :
149 423645 : ScalarType deathsPerSevereAdjusted =
150 423645 : params.template get<CriticalPerSevere<FP>>()[i] - criticalPerSevereAdjusted;
151 :
152 : // Exposed -> InfectedNoSymptoms
153 423645 : flows[this->template get_flat_flow_index<InfectionState::Exposed, InfectionState::InfectedNoSymptoms>(
154 847290 : {i})] = (1 / params.template get<TimeExposed<FP>>()[i]) * y[Ei];
155 :
156 : // InfectedNoSymptoms -> InfectedSymptoms / Recovered
157 847290 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptoms,
158 423645 : InfectionState::InfectedSymptoms>({i})] =
159 423645 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) *
160 423645 : (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSi];
161 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptoms, InfectionState::Recovered>(
162 847290 : {i})] = params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i] *
163 423645 : (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSi];
164 :
165 : // InfectedNoSymptomsConfirmed -> InfectedSymptomsConfirmed / Recovered
166 847290 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsConfirmed,
167 423645 : InfectionState::InfectedSymptomsConfirmed>({i})] =
168 423645 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) *
169 423645 : (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSCi];
170 847290 : flows[this->template get_flat_flow_index<InfectionState::InfectedNoSymptomsConfirmed,
171 423645 : InfectionState::Recovered>({i})] =
172 423645 : params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i] *
173 423645 : (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[i]) * y[INSCi];
174 :
175 : // InfectedSymptoms -> InfectedSevere / Recovered
176 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptoms, InfectionState::InfectedSevere>(
177 1270935 : {i})] = params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
178 847290 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyi];
179 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptoms, InfectionState::Recovered>(
180 1270935 : {i})] = (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
181 847290 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyi];
182 :
183 : // InfectedSymptomsConfirmed -> InfectedSevere / Recovered
184 847290 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsConfirmed,
185 423645 : InfectionState::InfectedSevere>({i})] =
186 847290 : params.template get<SeverePerInfectedSymptoms<FP>>()[i] /
187 847290 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyCi];
188 847290 : flows[this->template get_flat_flow_index<InfectionState::InfectedSymptomsConfirmed,
189 423645 : InfectionState::Recovered>({i})] =
190 847290 : (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[i]) /
191 847290 : params.template get<TimeInfectedSymptoms<FP>>()[i] * y[ISyCi];
192 :
193 : // InfectedSevere -> InfectedCritical / Recovered / Dead
194 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevere, InfectionState::InfectedCritical>(
195 847290 : {i})] = criticalPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevi];
196 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevere, InfectionState::Recovered>({i})] =
197 847290 : (1 - params.template get<CriticalPerSevere<FP>>()[i]) /
198 847290 : params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevi];
199 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedSevere, InfectionState::Dead>({i})] =
200 423645 : deathsPerSevereAdjusted / params.template get<TimeInfectedSevere<FP>>()[i] * y[ISevi];
201 :
202 : // InfectedCritical -> Dead / Recovered
203 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedCritical, InfectionState::Dead>({i})] =
204 423645 : params.template get<DeathsPerCritical<FP>>()[i] / params.template get<TimeInfectedCritical<FP>>()[i] *
205 423645 : y[ICri];
206 423645 : flows[this->template get_flat_flow_index<InfectionState::InfectedCritical, InfectionState::Recovered>(
207 1270935 : {i})] = (1 - params.template get<DeathsPerCritical<FP>>()[i]) /
208 847290 : params.template get<TimeInfectedCritical<FP>>()[i] * y[ICri];
209 : }
210 411453 : }
211 :
212 : /**
213 : * serialize this.
214 : * @see mio::serialize
215 : */
216 : template <class IOContext>
217 23 : void serialize(IOContext& io) const
218 : {
219 69 : auto obj = io.create_object("Model");
220 23 : obj.add_element("Parameters", this->parameters);
221 23 : obj.add_element("Populations", this->populations);
222 46 : }
223 :
224 : /**
225 : * deserialize an object of this class.
226 : * @see mio::deserialize
227 : */
228 : template <class IOContext>
229 9 : static IOResult<Model> deserialize(IOContext& io)
230 : {
231 27 : auto obj = io.expect_object("Model");
232 27 : auto par = obj.expect_element("Parameters", Tag<ParameterSet>{});
233 27 : auto pop = obj.expect_element("Populations", Tag<Populations>{});
234 : return apply(
235 : io,
236 9 : [](auto&& par_, auto&& pop_) {
237 9 : return Model{pop_, par_};
238 : },
239 18 : par, pop);
240 9 : }
241 : };
242 :
243 : //forward declaration, see below.
244 : template <typename FP = ScalarType, class BaseT = mio::Simulation<FP, Model<FP>>>
245 : class Simulation;
246 :
247 : /**
248 : * get percentage of infections per total population.
249 : * @param model the compartment model with initial values.
250 : * @param t current simulation time.
251 : * @param y current value of compartments.
252 : * @tparam FP floating point type, e.g., double.
253 : * @tparam Base simulation type that uses a secir compartment model. see Simulation.
254 : */
255 : template <typename FP = ScalarType, class Base = mio::Simulation<FP, Model<FP>>>
256 : double get_infections_relative(const Simulation<FP, Base>& model, FP t, const Eigen::Ref<const Vector<FP>>& y);
257 :
258 : /**
259 : * specialization of compartment model simulation for secir models.
260 : * @tparam FP floating point type, e.g., double.
261 : * @tparam BaseT simulation type that uses a secir compartment model. default mio::Simulation. For testing purposes only!
262 : */
263 : template <typename FP, class BaseT>
264 : class Simulation : public BaseT
265 : {
266 : public:
267 : /**
268 : * construct a simulation.
269 : * @param model the model to simulate.
270 : * @param t0 start time
271 : * @param dt time steps
272 : */
273 319 : Simulation(mio::osecir::Model<FP> const& model, FP t0 = 0., FP dt = 0.1)
274 : : BaseT(model, t0, dt)
275 319 : , m_t_last_npi_check(t0)
276 : {
277 319 : }
278 :
279 : /**
280 : * @brief advance simulation to tmax.
281 : * Overwrites Simulation::advance and includes a check for dynamic NPIs in regular intervals.
282 : * @see Simulation::advance
283 : * @param tmax next stopping point of simulation
284 : * @return value at tmax
285 : */
286 265 : Eigen::Ref<Vector<FP>> advance(FP tmax)
287 : {
288 265 : auto& t_end_dyn_npis = this->get_model().parameters.get_end_dynamic_npis();
289 265 : auto& dyn_npis = this->get_model().parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
290 265 : auto& contact_patterns = this->get_model().parameters.template get<ContactPatterns<FP>>();
291 :
292 : FP delay_npi_implementation; // delay which is needed to implement a NPI that is criterion-dependent
293 265 : auto t = BaseT::get_result().get_last_time();
294 265 : const auto dt = dyn_npis.get_thresholds().size() > 0 ? dyn_npis.get_interval().get() : tmax;
295 :
296 530 : while (t < tmax) {
297 265 : auto dt_eff = std::min({dt, tmax - t, m_t_last_npi_check + dt - t});
298 :
299 265 : BaseT::advance(t + dt_eff);
300 265 : if (t > 0) {
301 86 : delay_npi_implementation =
302 86 : this->get_model().parameters.template get<DynamicNPIsImplementationDelay<FP>>();
303 : }
304 : else { // DynamicNPIs for t=0 are 'misused' to be from-start NPIs. I.e., do not enforce delay.
305 179 : delay_npi_implementation = 0;
306 : }
307 265 : t = t + dt_eff;
308 :
309 265 : if (dyn_npis.get_thresholds().size() > 0) {
310 35 : if (floating_point_greater_equal(t, m_t_last_npi_check + dt)) {
311 35 : if (t < t_end_dyn_npis) {
312 70 : auto inf_rel = get_infections_relative<FP>(*this, t, this->get_result().get_last_value()) *
313 35 : dyn_npis.get_base_value();
314 35 : auto exceeded_threshold = dyn_npis.get_max_exceeded_threshold(inf_rel);
315 63 : if (exceeded_threshold != dyn_npis.get_thresholds().end() &&
316 28 : (exceeded_threshold->first > m_dynamic_npi.first ||
317 0 : t > ScalarType(m_dynamic_npi.second))) { //old npi was weaker or is expired
318 :
319 28 : auto t_start = SimulationTime(t + delay_npi_implementation);
320 28 : auto t_end = t_start + SimulationTime(ScalarType(dyn_npis.get_duration()));
321 28 : m_dynamic_npi = std::make_pair(exceeded_threshold->first, t_end);
322 28 : implement_dynamic_npis(contact_patterns.get_cont_freq_mat(), exceeded_threshold->second,
323 28 : t_start, t_end, [](auto& g) {
324 28 : return make_contact_damping_matrix(g);
325 : });
326 : }
327 : }
328 35 : m_t_last_npi_check = t;
329 : }
330 : }
331 : else {
332 230 : m_t_last_npi_check = t;
333 : }
334 : }
335 :
336 265 : return this->get_result().get_last_value();
337 : }
338 :
339 : private:
340 : double m_t_last_npi_check;
341 : std::pair<double, SimulationTime> m_dynamic_npi = {-std::numeric_limits<double>::max(), mio::SimulationTime(0)};
342 : };
343 :
344 : /**
345 : * @brief Specialization of simulate for SECIR models using Simulation.
346 : *
347 : * @tparam FP floating point type, e.g., double.
348 : * @param[in] t0 start time.
349 : * @param[in] tmax end time.
350 : * @param[in] dt time step.
351 : * @param[in] model SECIR model to simulate.
352 : * @param[in] integrator optional integrator, uses rk45 if nullptr.
353 : *
354 : * @return Returns the result of the simulation.
355 : */
356 : template <typename FP = ScalarType>
357 137 : inline auto simulate(FP t0, FP tmax, FP dt, const Model<FP>& model,
358 : std::shared_ptr<IntegratorCore<FP>> integrator = nullptr)
359 : {
360 137 : return mio::simulate<FP, Model<FP>, Simulation<>>(t0, tmax, dt, model, integrator);
361 : }
362 :
363 : /**
364 : * @brief Specialization of simulate for SECIR models using the FlowSimulation.
365 : *
366 : * @tparam FP floating point type, e.g., double.
367 : * @param[in] t0 start time.
368 : * @param[in] tmax end time.
369 : * @param[in] dt time step.
370 : * @param[in] model SECIR model to simulate.
371 : * @param[in] integrator optional integrator, uses rk45 if nullptr.
372 : *
373 : * @return Returns the result of the Flowsimulation.
374 : */
375 : template <typename FP = ScalarType>
376 : inline auto simulate_flows(FP t0, FP tmax, FP dt, const Model<FP>& model,
377 : std::shared_ptr<IntegratorCore<FP>> integrator = nullptr)
378 : {
379 : return mio::simulate_flows<FP, Model<FP>, Simulation<FP, mio::FlowSimulation<FP, Model<FP>>>>(t0, tmax, dt, model,
380 : integrator);
381 : }
382 :
383 : //see declaration above.
384 : template <typename FP, class Base>
385 44 : double get_infections_relative(const Simulation<FP, Base>& sim, FP /* t*/, const Eigen::Ref<const Vector<FP>>& y)
386 : {
387 44 : double sum_inf = 0;
388 106 : for (auto i = AgeGroup(0); i < sim.get_model().parameters.get_num_groups(); ++i) {
389 62 : sum_inf += sim.get_model().populations.get_from(y, {i, InfectionState::InfectedSymptoms});
390 : }
391 44 : auto inf_rel = sum_inf / sim.get_model().populations.get_total();
392 :
393 44 : return inf_rel;
394 : }
395 :
396 : /**
397 : *@brief Computes the reproduction number at a given index time of the Model output obtained by the Simulation.
398 : *@param t_idx The index time at which the reproduction number is computed.
399 : *@param sim The simulation holding the SECIR model
400 : *@tparam Base simulation type that uses a SECIR compartment model. see Simulation.
401 : *@returns The computed reproduction number at the provided index time.
402 : */
403 :
404 : template <typename FP, class Base>
405 144 : IOResult<FP> get_reproduction_number(size_t t_idx, const Simulation<FP, Base>& sim)
406 : {
407 :
408 144 : if (!(t_idx < static_cast<size_t>(sim.get_result().get_num_time_points()))) {
409 9 : return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries");
410 : }
411 :
412 135 : auto const& params = sim.get_model().parameters;
413 135 : const size_t num_groups = (size_t)sim.get_model().parameters.get_num_groups();
414 : //The infected compartments in the SECIR Model are: Exposed, Carrier, Infected, Hospitalized, ICU in respective agegroups
415 135 : const size_t num_infected_compartments = 5;
416 135 : const size_t total_infected_compartments = num_infected_compartments * num_groups;
417 135 : const double pi = std::acos(-1);
418 :
419 : //F encodes new Infections and V encodes transition times in the next-generation matrix calculation of R_t
420 135 : Eigen::MatrixXd F(total_infected_compartments, total_infected_compartments);
421 135 : Eigen::MatrixXd V(total_infected_compartments, total_infected_compartments);
422 135 : F = Eigen::MatrixXd::Zero(total_infected_compartments,
423 : total_infected_compartments); //Initialize matrices F and V with zeroes
424 135 : V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments);
425 :
426 135 : auto test_and_trace_required = 0.0;
427 135 : auto icu_occupancy = 0.0;
428 540 : for (auto i = AgeGroup(0); i < (mio::AgeGroup)num_groups; ++i) {
429 405 : test_and_trace_required +=
430 810 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]) /
431 810 : params.template get<TimeInfectedNoSymptoms<FP>>()[i] *
432 405 : sim.get_result().get_value(
433 405 : t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})];
434 405 : icu_occupancy += sim.get_result().get_value(
435 810 : t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedCritical})];
436 : }
437 :
438 135 : double season_val =
439 135 : (1 + params.template get<Seasonality<FP>>() *
440 135 : sin(pi *
441 135 : ((sim.get_model().parameters.template get<StartDay>() + sim.get_result().get_time(t_idx)) / 182.5 +
442 : 0.5)));
443 135 : ContactMatrixGroup const& contact_matrix = sim.get_model().parameters.template get<ContactPatterns<FP>>();
444 :
445 135 : Eigen::MatrixXd cont_freq_eff(num_groups, num_groups);
446 135 : Eigen::MatrixXd riskFromInfectedSymptomatic_derivatives(num_groups, num_groups);
447 135 : Eigen::VectorXd divN(num_groups);
448 135 : Eigen::VectorXd riskFromInfectedSymptomatic(num_groups);
449 :
450 540 : for (mio::AgeGroup k = 0; k < (mio::AgeGroup)num_groups; k++) {
451 405 : double temp = sim.get_result().get_value(
452 810 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Susceptible})] +
453 405 : sim.get_result().get_value(
454 810 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Exposed})] +
455 405 : sim.get_result().get_value(
456 810 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedNoSymptoms})] +
457 405 : sim.get_result().get_value(
458 810 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSymptoms})] +
459 405 : sim.get_result().get_value(
460 810 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSevere})] +
461 405 : sim.get_result().get_value(
462 810 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedCritical})] +
463 405 : sim.get_result().get_value(
464 405 : t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Recovered})];
465 405 : if (temp == 0) {
466 18 : temp = 1;
467 : }
468 405 : divN[(size_t)k] = 1 / temp;
469 :
470 2025 : riskFromInfectedSymptomatic[(size_t)k] = smoother_cosine(
471 405 : test_and_trace_required, params.template get<TestAndTraceCapacity<FP>>(),
472 810 : (params.template get<TestAndTraceCapacity<FP>>()) * params.template get<TestAndTraceCapacityMaxRisk<FP>>(),
473 405 : params.template get<RiskOfInfectionFromSymptomatic<FP>>()[k],
474 405 : params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[k]);
475 :
476 1620 : for (mio::AgeGroup l = 0; l < (mio::AgeGroup)num_groups; l++) {
477 1377 : if (test_and_trace_required < params.template get<TestAndTraceCapacity<FP>>() ||
478 324 : test_and_trace_required > params.template get<TestAndTraceCapacityMaxRisk<FP>>() *
479 162 : params.template get<TestAndTraceCapacity<FP>>()) {
480 1134 : riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) = 0;
481 : }
482 : else {
483 81 : riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) =
484 81 : -0.5 * pi *
485 162 : (params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[k] -
486 81 : params.template get<RiskOfInfectionFromSymptomatic<FP>>()[k]) /
487 81 : (4 * params.template get<TestAndTraceCapacity<FP>>()) *
488 162 : (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[l]) /
489 81 : params.template get<TimeInfectedNoSymptoms<FP>>()[l] *
490 81 : std::sin(pi / (4 * params.template get<TestAndTraceCapacity<FP>>()) *
491 81 : (test_and_trace_required - params.template get<TestAndTraceCapacity<FP>>()));
492 : }
493 : }
494 :
495 1620 : for (Eigen::Index l = 0; l < (Eigen::Index)num_groups; l++) {
496 2430 : cont_freq_eff(l, (size_t)k) =
497 3645 : season_val * contact_matrix.get_matrix_at(static_cast<double>(t_idx))(
498 1215 : static_cast<Eigen::Index>((size_t)l), static_cast<Eigen::Index>((size_t)k));
499 : }
500 : }
501 :
502 : //Check criterion if matrix V will be invertible by checking if subblock J is invertible
503 135 : Eigen::MatrixXd J(num_groups, num_groups);
504 135 : J = Eigen::MatrixXd::Zero(num_groups, num_groups);
505 540 : for (size_t i = 0; i < num_groups; i++) {
506 405 : J(i, i) = 1 / (params.template get<TimeInfectedCritical<FP>>()[(mio::AgeGroup)i]);
507 :
508 432 : if (!(icu_occupancy < 0.9 * params.template get<ICUCapacity<FP>>() ||
509 27 : icu_occupancy > (double)(params.template get<ICUCapacity<FP>>()))) {
510 108 : for (size_t j = 0; j < num_groups; j++) {
511 405 : J(i, j) -= sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
512 243 : {(mio::AgeGroup)i, InfectionState::InfectedSevere})] /
513 243 : params.template get<TimeInfectedSevere<FP>>()[(mio::AgeGroup)i] * 5 *
514 243 : params.template get<CriticalPerSevere<FP>>()[(mio::AgeGroup)i] * pi /
515 162 : (params.template get<ICUCapacity<FP>>()) *
516 81 : std::sin(pi / (0.1 * params.template get<ICUCapacity<FP>>()) *
517 81 : (icu_occupancy - 0.9 * params.template get<ICUCapacity<FP>>()));
518 : }
519 : }
520 : }
521 :
522 : //Check, if J is invertible
523 135 : if (J.determinant() == 0) {
524 9 : return mio::failure(mio::StatusCode::UnknownError, "Matrix V is not invertible");
525 : }
526 :
527 : //Initialize the matrix F
528 504 : for (size_t i = 0; i < num_groups; i++) {
529 1512 : for (size_t j = 0; j < num_groups; j++) {
530 :
531 1134 : double temp = 0;
532 4536 : for (Eigen::Index k = 0; k < (Eigen::Index)num_groups; k++) {
533 3402 : temp += cont_freq_eff(i, k) *
534 13608 : sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
535 6804 : {(mio::AgeGroup)k, InfectionState::InfectedSymptoms})] *
536 6804 : riskFromInfectedSymptomatic_derivatives(k, j) * divN[k];
537 : }
538 :
539 2268 : F(i, j + num_groups) =
540 4536 : sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
541 3402 : {(mio::AgeGroup)i, InfectionState::Susceptible})] *
542 2268 : params.template get<TransmissionProbabilityOnContact<FP>>()[(mio::AgeGroup)i] *
543 2268 : (cont_freq_eff(i, j) * params.template get<RelativeTransmissionNoSymptoms<FP>>()[(mio::AgeGroup)j] *
544 1134 : divN[(size_t)j] +
545 : temp);
546 : }
547 :
548 1512 : for (size_t j = 0; j < num_groups; j++) {
549 5670 : F(i, j + 2 * num_groups) = sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
550 3402 : {(mio::AgeGroup)i, InfectionState::Susceptible})] *
551 2268 : params.template get<TransmissionProbabilityOnContact<FP>>()[(mio::AgeGroup)i] *
552 1134 : cont_freq_eff(i, j) * riskFromInfectedSymptomatic[(size_t)j] * divN[(size_t)j];
553 : }
554 : }
555 :
556 : //Initialize the matrix V
557 504 : for (Eigen::Index i = 0; i < (Eigen::Index)num_groups; i++) {
558 1134 : double criticalPerSevereAdjusted = smoother_cosine(
559 756 : icu_occupancy, 0.90 * params.template get<ICUCapacity<FP>>(), params.template get<ICUCapacity<FP>>(),
560 756 : params.template get<CriticalPerSevere<FP>>()[(mio::AgeGroup)i], 0);
561 :
562 378 : V(i, i) = 1 / params.template get<TimeExposed<FP>>()[(mio::AgeGroup)i];
563 378 : V(i + num_groups, i) = -1 / params.template get<TimeExposed<FP>>()[(mio::AgeGroup)i];
564 378 : V(i + num_groups, i + num_groups) = 1 / params.template get<TimeInfectedNoSymptoms<FP>>()[(mio::AgeGroup)i];
565 378 : V(i + 2 * num_groups, i + num_groups) =
566 1134 : -(1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[(mio::AgeGroup)i]) /
567 756 : params.template get<TimeInfectedNoSymptoms<FP>>()[(mio::AgeGroup)i];
568 378 : V(i + 2 * num_groups, i + 2 * num_groups) =
569 756 : (1 / params.template get<TimeInfectedSymptoms<FP>>()[(mio::AgeGroup)i]);
570 378 : V(i + 3 * num_groups, i + 2 * num_groups) =
571 1134 : -params.template get<SeverePerInfectedSymptoms<FP>>()[(mio::AgeGroup)i] /
572 756 : params.template get<TimeInfectedSymptoms<FP>>()[(mio::AgeGroup)i];
573 378 : V(i + 3 * num_groups, i + 3 * num_groups) =
574 756 : 1 / (params.template get<TimeInfectedSevere<FP>>()[(mio::AgeGroup)i]);
575 378 : V(i + 4 * num_groups, i + 3 * num_groups) =
576 756 : -criticalPerSevereAdjusted / (params.template get<TimeInfectedSevere<FP>>()[(mio::AgeGroup)i]);
577 :
578 1512 : for (size_t j = 0; j < num_groups; j++) {
579 1134 : V(i + 4 * num_groups, j + 4 * num_groups) = J(i, j);
580 : }
581 : }
582 :
583 126 : V = V.inverse();
584 :
585 : //Compute F*V
586 126 : Eigen::MatrixXd NextGenMatrix(num_infected_compartments * num_groups, 5 * num_groups);
587 126 : NextGenMatrix = F * V;
588 :
589 : //Compute the largest eigenvalue in absolute value
590 126 : Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;
591 :
592 126 : ces.compute(NextGenMatrix);
593 126 : const Eigen::VectorXcd eigen_vals = ces.eigenvalues();
594 :
595 126 : Eigen::VectorXd eigen_vals_abs;
596 126 : eigen_vals_abs.resize(eigen_vals.size());
597 :
598 2016 : for (int i = 0; i < eigen_vals.size(); i++) {
599 1890 : eigen_vals_abs[i] = std::abs(eigen_vals[i]);
600 : }
601 126 : return mio::success(eigen_vals_abs.maxCoeff());
602 135 : }
603 : /**
604 : *@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation.
605 : *@param sim The Model Simulation.
606 : *@tparam Base simulation type that uses a SECIR compartment model. see Simulation.
607 : *@returns Eigen::Vector containing all reproduction numbers
608 : */
609 :
610 : template <typename FP, class Base>
611 : Vector<FP> get_reproduction_numbers(const Simulation<FP, Base>& sim)
612 : {
613 : Vector<FP> temp(sim.get_result().get_num_time_points());
614 : for (int i = 0; i < sim.get_result().get_num_time_points(); i++) {
615 : temp[i] = get_reproduction_number((size_t)i, sim).value();
616 : }
617 : return temp;
618 : }
619 :
620 : /**
621 : *@brief @brief Computes the reproduction number at a given time point of the Simulation. If the particular time point is not part of the output, a linearly interpolated value is returned.
622 : *@param t_value The time point at which the reproduction number should be computed.
623 : *@param sim The Model Simulation.
624 : *@tparam Base simulation type that uses a SECIR compartment model. see Simulation.
625 : *@returns The computed reproduction number at the provided time point, potentially using linear interpolation.
626 : */
627 : template <class Base>
628 54 : IOResult<ScalarType> get_reproduction_number(ScalarType t_value, const Simulation<Base>& sim)
629 : {
630 54 : if (t_value < sim.get_result().get_time(0) || t_value > sim.get_result().get_last_time()) {
631 : return mio::failure(mio::StatusCode::OutOfRange,
632 18 : "Cannot interpolate reproduction number outside computed horizon of the TimeSeries");
633 : }
634 :
635 36 : if (t_value == sim.get_result().get_time(0)) {
636 9 : return mio::success(get_reproduction_number((size_t)0, sim).value());
637 : }
638 :
639 27 : const auto& times = sim.get_result().get_times();
640 :
641 27 : auto time_late = std::distance(times.begin(), std::lower_bound(times.begin(), times.end(), t_value));
642 :
643 27 : ScalarType y1 = get_reproduction_number(static_cast<size_t>(time_late - 1), sim).value();
644 27 : ScalarType y2 = get_reproduction_number(static_cast<size_t>(time_late), sim).value();
645 :
646 27 : auto result = linear_interpolation(t_value, sim.get_result().get_time(time_late - 1),
647 27 : sim.get_result().get_time(time_late), y1, y2);
648 27 : return mio::success(static_cast<ScalarType>(result));
649 : }
650 :
651 : /**
652 : * Get mobility factors.
653 : * Used by mobility graph simulation.
654 : * Like infection risk, mobility of infected individuals is reduced if they are well isolated.
655 : * @param model the compartment model with initial values.
656 : * @param t current simulation time.
657 : * @param y current value of compartments.
658 : * @return vector expression, same size as y, with mobility factors per compartment.
659 : * @tparam FP floating point type, e.g., double.
660 : * @tparam Base simulation type that uses a secir compartment model; see Simulation.
661 : */
662 : template <typename FP = ScalarType, class Base = mio::Simulation<Model<FP>, FP>>
663 75 : auto get_mobility_factors(const Simulation<Base>& sim, FP /*t*/, const Eigen::Ref<const Vector<FP>>& y)
664 : {
665 75 : auto& params = sim.get_model().parameters;
666 : //parameters as arrays
667 75 : auto&& p_asymp = params.template get<RecoveredPerInfectedNoSymptoms<FP>>().array().template cast<FP>();
668 75 : auto&& p_inf = params.template get<RiskOfInfectionFromSymptomatic<FP>>().array().template cast<FP>();
669 75 : auto&& p_inf_max = params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>().array().template cast<FP>();
670 : //slice of InfectedNoSymptoms
671 150 : auto y_INS = slice(y, {Eigen::Index(InfectionState::InfectedNoSymptoms),
672 150 : Eigen::Index(size_t(params.get_num_groups())), Eigen::Index(InfectionState::Count)});
673 :
674 : //compute isolation, same as infection risk from main model
675 75 : auto test_and_trace_required =
676 225 : ((1 - p_asymp) / params.template get<TimeInfectedNoSymptoms<FP>>().array().template cast<FP>() * y_INS.array())
677 300 : .sum();
678 75 : auto test_and_trace_capacity = double(params.template get<TestAndTraceCapacity<FP>>());
679 75 : auto test_and_trace_capacity_max_risk = double(params.template get<TestAndTraceCapacityMaxRisk<FP>>());
680 75 : auto riskFromInfectedSymptomatic =
681 : smoother_cosine(test_and_trace_required, test_and_trace_capacity,
682 : test_and_trace_capacity * test_and_trace_capacity_max_risk, p_inf.matrix(), p_inf_max.matrix());
683 :
684 : //set factor for infected
685 75 : auto factors = Eigen::VectorXd::Ones(y.rows()).eval();
686 225 : slice(factors, {Eigen::Index(InfectionState::InfectedSymptoms), Eigen::Index(size_t(params.get_num_groups())),
687 : Eigen::Index(InfectionState::Count)})
688 150 : .array() = riskFromInfectedSymptomatic;
689 150 : return factors;
690 75 : }
691 :
692 : template <typename FP = ScalarType, class Base = mio::Simulation<Model<FP>, FP>>
693 9 : auto test_commuters(Simulation<FP, Base>& sim, Eigen::Ref<Vector<FP>> mobile_population, FP time)
694 : {
695 9 : auto& model = sim.get_model();
696 9 : auto nondetection = 1.0;
697 18 : if (time >= model.parameters.get_start_commuter_detection() &&
698 9 : time < model.parameters.get_end_commuter_detection()) {
699 9 : nondetection = (double)model.parameters.get_commuter_nondetection();
700 : }
701 27 : for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); ++i) {
702 18 : auto INSi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptoms});
703 18 : auto INSCi = model.populations.get_flat_index({i, InfectionState::InfectedNoSymptomsConfirmed});
704 18 : auto ISyi = model.populations.get_flat_index({i, InfectionState::InfectedSymptoms});
705 18 : auto ISyCi = model.populations.get_flat_index({i, InfectionState::InfectedSymptomsConfirmed});
706 :
707 : //put detected commuters in their own compartment so they don't contribute to infections in their home node
708 18 : sim.get_result().get_last_value()[INSi] -= mobile_population[INSi] * (1 - nondetection);
709 18 : sim.get_result().get_last_value()[INSCi] += mobile_population[INSi] * (1 - nondetection);
710 18 : sim.get_result().get_last_value()[ISyi] -= mobile_population[ISyi] * (1 - nondetection);
711 18 : sim.get_result().get_last_value()[ISyCi] += mobile_population[ISyi] * (1 - nondetection);
712 :
713 : //reduce the number of commuters
714 18 : mobile_population[ISyi] *= nondetection;
715 18 : mobile_population[INSi] *= nondetection;
716 : }
717 9 : }
718 :
719 : } // namespace osecir
720 : } // namespace mio
721 :
722 : #endif // ODESECIR_MODEL_H
|
