model.h Source File

CPP API: model.h Source File
glct_secir/model.h
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2020-2026 MEmilio
3 *
4 * Authors: Lena Ploetzke
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 
21 #ifndef MIO_GLCT_SECIR_MODEL_H
22 #define MIO_GLCT_SECIR_MODEL_H
23 
24 #include "glct_secir/parameters.h"
29 #include "memilio/config.h"
30 #include "memilio/utils/logging.h"
32 #include "memilio/math/eigen.h"
33 
34 #include <numbers>
35 
36 namespace mio
37 {
38 namespace glsecir
39 {
40 
50 template <typename FP, size_t NumExposed, size_t NumInfectedNoSymptoms, size_t NumInfectedSymptoms,
51  size_t NumInfectedSevere, size_t NumInfectedCritical>
52 class Model
53  : public CompartmentalModel<
54  FP,
55  LctInfectionState<FP, InfectionState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
56  NumInfectedSevere, NumInfectedCritical, 1, 1>,
57  mio::Populations<FP, LctInfectionState<FP, InfectionState, 1, NumExposed, NumInfectedNoSymptoms,
58  NumInfectedSymptoms, NumInfectedSevere, NumInfectedCritical, 1, 1>>,
59  Parameters<FP>>
60 {
61 
62 public:
63  using LctState = LctInfectionState<FP, InfectionState, 1, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
64  NumInfectedSevere, NumInfectedCritical, 1,
65  1>;
67  using typename Base::ParameterSet;
68  using typename Base::Populations;
69 
73  {
74  }
75 
82  bool check_constraints() const
83  {
84  auto params = this->parameters;
85 
86  // --- Check that the dimensions are consistent. ---
87  if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::Exposed>() !=
88  params.template get<StartingProbabilitiesExposed<FP>>().rows()) {
89  log_error("Constraint check: Dimension of the parameters does not match the number of subcompartments for "
90  "the Exposed "
91  "compartment.");
92  return true;
93  }
94  if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() !=
95  params.template get<StartingProbabilitiesInfectedNoSymptoms<FP>>().rows()) {
96  log_error(
97  "Constraint check: Dimension of the parameters does not match the number of subcompartments for the "
98  "InfectedNoSymptoms compartment.");
99  return true;
100  }
101  if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedSymptoms>() !=
102  params.template get<StartingProbabilitiesInfectedSymptoms<FP>>().rows()) {
103  log_error(
104  "Constraint check: Dimension of the parameters does not match the number of subcompartments for the "
105  "InfectedSymptoms compartment.");
106  return true;
107  }
108  if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedSevere>() !=
109  params.template get<StartingProbabilitiesInfectedSevere<FP>>().rows()) {
110  log_error("Constraint check: Dimension of the parameters does not match the number of subcompartments for "
111  "the InfectedSevere "
112  "compartment.");
113  return true;
114  }
115  if ((Eigen::Index)LctState::template get_num_subcompartments<InfectionState::InfectedCritical>() !=
116  params.template get<StartingProbabilitiesInfectedCritical<FP>>().rows()) {
117  log_error(
118  "Constraint check: Dimension of the parameters does not match the number of subcompartments for the "
119  "InfectedCritical compartment.");
120  return true;
121  }
122 
123  return (params.check_constraints() || this->populations.check_constraints());
124  }
125 
138  void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
139  Eigen::Ref<Eigen::VectorX<FP>> dydt) const override
140  {
141  using std::sin;
142 
143  dydt.setZero();
144 
145  auto params = this->parameters;
146  auto total_population = pop.sum() - pop[LctState::template get_first_index<InfectionState::Dead>()];
147 
148  // Calculate sum of all subcompartments for InfectedNoSymptoms.
149  FP InfectedNoSymptoms_sum =
150  pop.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
151  LctState::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
152  .sum();
153  // Calculate sum of all subcompartments for InfectedSymptoms.
154  FP InfectedSymptoms_sum =
155  pop.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
156  LctState::template get_num_subcompartments<InfectionState::InfectedSymptoms>())
157  .sum();
158 
159  // --- Susceptibles. ---
160  FP season_val =
161  1 + params.template get<Seasonality<FP>>() *
162  sin(std::numbers::pi_v<ScalarType> * ((params.template get<StartDay<FP>>() + t) / 182.5 + 0.5));
163  dydt[0] =
164  -y[0] / total_population * season_val * params.template get<TransmissionProbabilityOnContact<FP>>() *
165  params.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_matrix_at(SimulationTime<FP>(t))(0, 0) *
166  (params.template get<RelativeTransmissionNoSymptoms<FP>>() * InfectedNoSymptoms_sum +
167  params.template get<RiskOfInfectionFromSymptomatic<FP>>() * InfectedSymptoms_sum);
168 
169  // --- Exposed. ---
170  dydt.segment(LctState::template get_first_index<InfectionState::Exposed>(),
171  LctState::template get_num_subcompartments<InfectionState::Exposed>()) -=
172  dydt[0] * params.template get<StartingProbabilitiesExposed<FP>>();
173  dydt.segment(LctState::template get_first_index<InfectionState::Exposed>(),
174  LctState::template get_num_subcompartments<InfectionState::Exposed>()) +=
175  params.template get<TransitionMatrixExposedToInfectedNoSymptoms<FP>>().transpose() *
176  y.segment(LctState::template get_first_index<InfectionState::Exposed>(),
177  LctState::template get_num_subcompartments<InfectionState::Exposed>());
178 
179  // --- InfectedNoSymptoms. ---
180  // Flow from Exposed To InfectedNoSymptoms.
181  dydt.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
182  LctState::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>()) =
184  Eigen::VectorX<FP>::Ones(LctState::template get_num_subcompartments<InfectionState::Exposed>()))
185  .transpose() *
186  y.segment(LctState::template get_first_index<InfectionState::Exposed>(),
187  LctState::template get_num_subcompartments<InfectionState::Exposed>()) *
189  // Flow from InfectedNoSymptoms To InfectedSymptoms.
190  size_t dimensionInfectedNoSymptomsToInfectedSymptoms =
191  params.template get<TransitionMatrixInfectedNoSymptomsToInfectedSymptoms<FP>>().rows();
192  dydt.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
193  dimensionInfectedNoSymptomsToInfectedSymptoms) +=
194  params.template get<TransitionMatrixInfectedNoSymptomsToInfectedSymptoms<FP>>().transpose() *
195  y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
196  dimensionInfectedNoSymptomsToInfectedSymptoms);
197  // Flow from InfectedNoSymptoms To Recovered.
198  size_t dimensionInfectedNoSymptomsToRecovered =
199  params.template get<TransitionMatrixInfectedNoSymptomsToRecovered<FP>>().rows();
200  dydt.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>() +
201  dimensionInfectedNoSymptomsToInfectedSymptoms,
202  dimensionInfectedNoSymptomsToRecovered) +=
203  params.template get<TransitionMatrixInfectedNoSymptomsToRecovered<FP>>().transpose() *
204  y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>() +
205  dimensionInfectedNoSymptomsToInfectedSymptoms,
206  dimensionInfectedNoSymptomsToRecovered);
207  // Add flow directly to Recovered compartment.
208  dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
209  -(params.template get<TransitionMatrixInfectedNoSymptomsToRecovered<FP>>() *
210  Eigen::VectorX<FP>::Ones(dimensionInfectedNoSymptomsToRecovered))
211  .transpose() *
212  y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>() +
213  dimensionInfectedNoSymptomsToInfectedSymptoms,
214  dimensionInfectedNoSymptomsToRecovered);
215 
216  // --- InfectedSymptoms. ---
217  // Flow from InfectedNoSymptoms To InfectedSymptoms.
218  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
219  LctState::template get_num_subcompartments<InfectionState::InfectedSymptoms>()) =
220  -params.template get<StartingProbabilitiesInfectedSymptoms<FP>>() *
221  (params.template get<TransitionMatrixInfectedNoSymptomsToInfectedSymptoms<FP>>() *
222  Eigen::VectorX<FP>::Ones(dimensionInfectedNoSymptomsToInfectedSymptoms))
223  .transpose() *
224  y.segment(LctState::template get_first_index<InfectionState::InfectedNoSymptoms>(),
225  dimensionInfectedNoSymptomsToInfectedSymptoms);
226  // Flow from InfectedSymptoms To InfectedSevere.
227  size_t dimensionInfectedSymptomsToInfectedSevere =
228  params.template get<TransitionMatrixInfectedSymptomsToInfectedSevere<FP>>().rows();
229  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
230  dimensionInfectedSymptomsToInfectedSevere) +=
231  params.template get<TransitionMatrixInfectedSymptomsToInfectedSevere<FP>>().transpose() *
232  y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
233  dimensionInfectedSymptomsToInfectedSevere);
234  // Flow from InfectedSymptoms To Recovered.
235  size_t dimensionInfectedSymptomsToRecovered =
236  params.template get<TransitionMatrixInfectedSymptomsToRecovered<FP>>().rows();
237  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>() +
238  dimensionInfectedSymptomsToInfectedSevere,
239  dimensionInfectedSymptomsToRecovered) +=
240  params.template get<TransitionMatrixInfectedSymptomsToRecovered<FP>>().transpose() *
241  y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>() +
242  dimensionInfectedSymptomsToInfectedSevere,
243  dimensionInfectedSymptomsToRecovered);
244  // Add flow directly to Recovered compartment.
245  dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
246  -(params.template get<TransitionMatrixInfectedSymptomsToRecovered<FP>>() *
247  Eigen::VectorX<FP>::Ones(dimensionInfectedSymptomsToRecovered))
248  .transpose() *
249  y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>() +
250  dimensionInfectedSymptomsToInfectedSevere,
251  dimensionInfectedSymptomsToRecovered);
252 
253  // --- InfectedSevere. ---
254  // Flow from InfectedSymptoms To InfectedSevere.
255  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
256  LctState::template get_num_subcompartments<InfectionState::InfectedSevere>()) =
257  -params.template get<StartingProbabilitiesInfectedSevere<FP>>() *
258  (params.template get<TransitionMatrixInfectedSymptomsToInfectedSevere<FP>>() *
259  Eigen::VectorX<FP>::Ones(dimensionInfectedSymptomsToInfectedSevere))
260  .transpose() *
261  y.segment(LctState::template get_first_index<InfectionState::InfectedSymptoms>(),
262  dimensionInfectedSymptomsToInfectedSevere);
263  // Flow from InfectedSevere To InfectedCritical.
264  size_t dimensionInfectedSevereToInfectedCritical =
265  params.template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>().rows();
266  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
267  dimensionInfectedSevereToInfectedCritical) +=
268  params.template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>().transpose() *
269  y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
270  dimensionInfectedSevereToInfectedCritical);
271  // Flow from InfectedSevere To Dead.
272  size_t dimensionInfectedSevereToDead = params.template get<TransitionMatrixInfectedSevereToDead<FP>>().rows();
273  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
274  dimensionInfectedSevereToInfectedCritical,
275  dimensionInfectedSevereToDead) +=
276  params.template get<TransitionMatrixInfectedSevereToDead<FP>>().transpose() *
277  y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
278  dimensionInfectedSevereToInfectedCritical,
279  dimensionInfectedSevereToDead);
280  // Flow from InfectedSevere To Recovered.
281  size_t dimensionInfectedSevereToRecovered =
282  params.template get<TransitionMatrixInfectedSevereToRecovered<FP>>().rows();
283  dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
284  dimensionInfectedSevereToInfectedCritical + dimensionInfectedSevereToDead,
285  dimensionInfectedSevereToRecovered) +=
286  params.template get<TransitionMatrixInfectedSevereToRecovered<FP>>().transpose() *
287  y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
288  dimensionInfectedSevereToInfectedCritical + dimensionInfectedSevereToDead,
289  dimensionInfectedSevereToRecovered);
290  // Add flow directly to Recovered compartment.
291  dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
292  -(params.template get<TransitionMatrixInfectedSevereToRecovered<FP>>() *
293  Eigen::VectorX<FP>::Ones(dimensionInfectedSevereToRecovered))
294  .transpose() *
295  y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
296  dimensionInfectedSevereToInfectedCritical + dimensionInfectedSevereToDead,
297  dimensionInfectedSevereToRecovered);
298 
299  // --- InfectedCritical. ---
300  // Flow from InfectedSevere To InfectedCritical.
301  dydt.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
302  LctState::template get_num_subcompartments<InfectionState::InfectedCritical>()) =
303  -params.template get<StartingProbabilitiesInfectedCritical<FP>>() *
304  (params.template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>() *
305  Eigen::VectorX<FP>::Ones(dimensionInfectedSevereToInfectedCritical))
306  .transpose() *
307  y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
308  dimensionInfectedSevereToInfectedCritical);
309  // Flow from InfectedCritical To Dead.
310  size_t dimensionInfectedCriticalToDead =
311  params.template get<TransitionMatrixInfectedCriticalToDead<FP>>().rows();
312  dydt.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
313  dimensionInfectedCriticalToDead) +=
314  params.template get<TransitionMatrixInfectedCriticalToDead<FP>>().transpose() *
315  y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
316  dimensionInfectedCriticalToDead);
317  // Flow from InfectedCritical To Recovered.
318  size_t dimensionInfectedCriticalToRecovered =
319  params.template get<TransitionMatrixInfectedCriticalToRecovered<FP>>().rows();
320  dydt.segment(LctState::template get_first_index<InfectionState::InfectedCritical>() +
321  dimensionInfectedCriticalToDead,
322  dimensionInfectedCriticalToRecovered) +=
323  params.template get<TransitionMatrixInfectedCriticalToRecovered<FP>>().transpose() *
324  y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>() +
325  dimensionInfectedCriticalToDead,
326  dimensionInfectedCriticalToRecovered);
327  // Add flow directly to Recovered compartment.
328  dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
329  -(params.template get<TransitionMatrixInfectedCriticalToRecovered<FP>>() *
330  Eigen::VectorX<FP>::Ones(dimensionInfectedCriticalToRecovered))
331  .transpose() *
332  y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>() +
333  dimensionInfectedCriticalToDead,
334  dimensionInfectedCriticalToRecovered);
335 
336  // --- Dead. ---
337  dydt[LctState::template get_first_index<InfectionState::Dead>()] +=
338  -(params.template get<TransitionMatrixInfectedSevereToDead<FP>>() *
339  Eigen::VectorX<FP>::Ones(dimensionInfectedSevereToDead))
340  .transpose() *
341  y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
342  dimensionInfectedSevereToInfectedCritical,
343  dimensionInfectedSevereToDead);
344  dydt[LctState::template get_first_index<InfectionState::Dead>()] +=
345  -(params.template get<TransitionMatrixInfectedCriticalToDead<FP>>() *
346  Eigen::VectorX<FP>::Ones(dimensionInfectedCriticalToDead))
347  .transpose() *
348  y.segment(LctState::template get_first_index<InfectionState::InfectedCritical>(),
349  dimensionInfectedCriticalToDead);
350  }
351 
364  TimeSeries<FP> calculate_compartments(const TimeSeries<FP>& subcompartments_ts) const
365  {
366  TimeSeries<FP> compartments_ts((Eigen::Index)InfectionState::Count);
367  if (!(LctState::Count == subcompartments_ts.get_num_elements())) {
368  log_error("Result does not match InfectionStates of the model.");
369  // Return a TimeSeries with values -1.
370  Eigen::VectorX<FP> error_output = Eigen::VectorX<FP>::Constant((Eigen::Index)InfectionState::Count, -1);
371  compartments_ts.add_time_point(-1, error_output);
372  return compartments_ts;
373  }
374  Eigen::VectorX<FP> compartments((Eigen::Index)InfectionState::Count);
375  for (Eigen::Index i = 0; i < subcompartments_ts.get_num_time_points(); ++i) {
376  compartments_ts.add_time_point(subcompartments_ts.get_time(i),
377  LctState::calculate_compartments(subcompartments_ts[i]));
378  }
379  return compartments_ts;
380  }
381 };
382 
383 } // namespace glsecir
384 } // namespace mio
385 
386 #endif // MIO_GLCT_SECIR_MODEL_H
An Index with more than one template parameter combines several Index objects.
Definition: index.h:191
Provides the functionality to be able to work with subcompartments in an LCT model.
Definition: lct_infection_state.h:42
static constexpr size_t Count
Definition: lct_infection_state.h:114
static Eigen::VectorX< FP > calculate_compartments(const Eigen::VectorX< FP > &subcompartments)
Cumulates a vector with the number of individuals in each subcompartment (with subcompartments accord...
Definition: lct_infection_state.h:93
Typesafe wrapper for a floating-point simulation time value (in days).
Definition: damping.h:78
stores vectors of values at time points (or some other abstract variable) the value at each time poin...
Definition: time_series.h:58
Eigen::Index get_num_elements() const
number of elements of vector at each time point
Definition: time_series.h:205
Eigen::Index get_num_time_points() const
number of time points in the series
Definition: time_series.h:197
FP & get_time(Eigen::Index i)
time of time point at index i
Definition: time_series.h:272
Eigen::Ref< Vector > add_time_point()
add one uninitialized time point
Definition: time_series.h:221
Class that defines an GLCT-SECIR model.
Definition: glct_secir/model.h:60
TimeSeries< FP > calculate_compartments(const TimeSeries< FP > &subcompartments_ts) const
Cumulates a simulation result with subcompartments to produce a result that divides the population on...
Definition: glct_secir/model.h:364
Pop Populations
Definition: compartmental_model.h:62
Params ParameterSet
Definition: compartmental_model.h:63
void get_derivatives(Eigen::Ref< const Eigen::VectorX< FP >> pop, Eigen::Ref< const Eigen::VectorX< FP >> y, FP t, Eigen::Ref< Eigen::VectorX< FP >> dydt) const override
Evaluates the right-hand-side f of the GLCT dydt = f(y, t).
Definition: glct_secir/model.h:138
bool check_constraints() const
Checks that the model satisfies all constraints (e.g.
Definition: glct_secir/model.h:82
Model()
Default constructor.
Definition: glct_secir/model.h:71
Parameters of an GLCT-SECIR model.
Definition: glct_secir/parameters.h:498
ad::internal::unary_intermediate< AD_TAPE_REAL, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::operations::ad_sin< AD_TAPE_REAL > > sin(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &x1)
Definition: ad.hpp:913
InfectionState
The InfectionState enum describes the basic categories for the infection state of persons.
Definition: glct_secir/infection_state.h:31
A collection of classes to simplify handling of matrix shapes in meta programming.
Definition: models/abm/analyze_result.h:30
auto i
Definition: io.h:810
void log_error(spdlog::string_view_t fmt, const Args &... args)
Definition: logging.h:114
constexpr std::tuple_element< I, std::tuple< Index< CategoryTags >... > >::type & get(Index< CategoryTags... > &i) noexcept
Retrieves the Index (by reference) at the Ith position of a MultiIndex.
Definition: index.h:304
CompartmentalModel is a template for a compartmental model for an array of initial populations and a ...
Definition: compartmental_model.h:59
Pop Populations
Definition: compartmental_model.h:62
Params ParameterSet
Definition: compartmental_model.h:63
The contact patterns within the society are modelled using an UncertainContactMatrix.
Definition: glct_secir/parameters.h:402
The risk of infection from symptomatic cases in the GLCT-SECIR model.
Definition: glct_secir/parameters.h:433
The start day in the GLCT-SECIR model.
Definition: glct_secir/parameters.h:452
Vector with the probability to start in any of the subcompartments of the Exposed compartment.
Definition: glct_secir/parameters.h:42
Vector with the probability to start in any of the subcompartments of the InfectedCritical compartmen...
Definition: glct_secir/parameters.h:318
Vector with the probability to start in any of the subcompartments of the InfectedNoSymptoms compartm...
Definition: glct_secir/parameters.h:83
Vector with the probability to start in any of the subcompartments of the InfectedSevere compartment.
Definition: glct_secir/parameters.h:223
Vector with the probability to start in any of the subcompartments of the InfectedSymptoms compartmen...
Definition: glct_secir/parameters.h:153
Transition matrix of the Exposed compartment.
Definition: glct_secir/parameters.h:62