parameter_space.h Source File

CPP API: parameter_space.h Source File
ode_secir/parameter_space.h
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2020-2026 MEmilio
3 *
4 * Authors: Daniel Abele, 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_PARAMETER_SPACE_H
21 #define ODESECIR_PARAMETER_SPACE_H
22 
23 #include "memilio/mobility/metapopulation_mobility_instant.h"
24 #include "memilio/utils/memory.h"
25 #include "memilio/utils/logging.h"
26 #include "memilio/utils/parameter_distributions.h"
27 #include "ode_secir/model.h"
28 
29 #include <assert.h>
30 #include <string>
31 #include <vector>
32 #include <random>
33 #include <memory>
34 
35 namespace mio
36 {
37 namespace osecir
38 {
46 template <typename FP>
47 void set_params_distributions_normal(Model<FP>& model, FP t0, FP tmax, FP dev_rel)
48 {
49  using std::max;
50  using std::min;
51  auto set_distribution = [dev_rel](UncertainValue<FP>& v, FP min_val = 0.001) {
52  auto lower_bound = min<FP>(max<FP>(min_val, (1 - dev_rel * 2.6) * v), 0.1 * std::numeric_limits<FP>::max());
53  auto upper_bound = min<FP>(max<FP>(min_val, (1 + dev_rel * 2.6) * v), 0.5 * std::numeric_limits<FP>::max());
54 
55  if (mio::floating_point_equal<FP>(lower_bound, upper_bound, mio::Limits<FP>::zero_tolerance())) {
56  //MSVC has problems if standard deviation for normal distribution is zero
57  mio::log_debug("Bounded ParameterDistribution has standard deviation close to zero. Therefore constant "
58  "distribution is used.");
59  v.set_distribution(ParameterDistributionConstant(
60  min<FP>(max<FP>(min_val, v.value()), 0.3 * std::numeric_limits<FP>::max())));
61  }
62  else {
63  v.set_distribution(ParameterDistributionNormal(
64  //add add limits for nonsense big values. Also mscv has a problem with a few doubles so this fixes it
65  lower_bound, upper_bound, min<FP>(max<FP>(min_val, v.value()), 0.3 * std::numeric_limits<FP>::max()),
66  min<FP>(max<FP>(min_val, dev_rel * v), std::numeric_limits<FP>::max())));
67  }
68  };
69 
70  set_distribution(model.parameters.template get<Seasonality<FP>>(), 0.0);
71  set_distribution(model.parameters.template get<ICUCapacity<FP>>());
72  set_distribution(model.parameters.template get<TestAndTraceCapacity<FP>>());
73 
74  // populations
75  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
76  for (auto j = Index<InfectionState>(0); j < Index<InfectionState>(InfectionState::Count); j++) {
77 
78  // don't touch S and D
79  if (j == InfectionState::Susceptible || j == InfectionState::Dead) {
80  continue;
81  }
82 
83  // variably sized groups
84  set_distribution(model.populations[{i, j}], 0.0);
85  }
86  }
87 
88  // times
89  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
90 
91  set_distribution(model.parameters.template get<TimeExposed<FP>>()[i]);
92  set_distribution(model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[i]);
93  set_distribution(model.parameters.template get<TimeInfectedSymptoms<FP>>()[i]);
94  set_distribution(model.parameters.template get<TimeInfectedSevere<FP>>()[i]);
95  set_distribution(model.parameters.template get<TimeInfectedCritical<FP>>()[i]);
96 
97  set_distribution(model.parameters.template get<TransmissionProbabilityOnContact<FP>>()[i]);
98  set_distribution(model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[i]);
99  set_distribution(model.parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]);
100  set_distribution(model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[i]);
101  set_distribution(model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i]);
102  set_distribution(model.parameters.template get<DeathsPerCritical<FP>>()[i]);
103  set_distribution(model.parameters.template get<SeverePerInfectedSymptoms<FP>>()[i]);
104  set_distribution(model.parameters.template get<CriticalPerSevere<FP>>()[i]);
105  }
106 
107  // dampings
108  auto matrices = std::vector<size_t>();
109  for (size_t i = 0; i < model.parameters.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_num_matrices();
110  ++i) {
111  matrices.push_back(i);
112  }
113  auto groups = Eigen::VectorX<FP>::Constant(Eigen::Index(model.parameters.get_num_groups().get()), 1.0);
114  model.parameters.template get<ContactPatterns<FP>>().get_dampings().emplace_back(
115  mio::UncertainValue<FP>(0.5), mio::DampingLevel(0), mio::DampingType(0),
116  mio::SimulationTime<FP>(t0 + (tmax - t0) / 2), matrices, groups);
117  set_distribution(model.parameters.template get<ContactPatterns<FP>>().get_dampings()[0].get_value(), 0.0);
118 }
119 
124 template <typename FP>
125 void draw_sample_demographics(Model<FP>& model)
126 {
127  model.parameters.template get<ICUCapacity<FP>>().draw_sample();
128  model.parameters.template get<TestAndTraceCapacity<FP>>().draw_sample();
129 
130  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
131  FP group_total = model.populations.get_group_total(i);
132 
133  model.populations[{i, InfectionState::Exposed}].draw_sample();
134  model.populations[{i, InfectionState::InfectedNoSymptoms}].draw_sample();
135  model.populations[{i, InfectionState::InfectedSymptoms}].draw_sample();
136  model.populations[{i, InfectionState::InfectedSevere}].draw_sample();
137  model.populations[{i, InfectionState::InfectedCritical}].draw_sample();
138  model.populations[{i, InfectionState::Recovered}].draw_sample();
139 
140  // no sampling for dead and total numbers
141  // [...]
142 
143  model.populations.template set_difference_from_group_total<AgeGroup>({i, InfectionState::Susceptible},
144  group_total);
145  model.populations.template set_difference_from_group_total<AgeGroup>({i, InfectionState::Susceptible},
146  model.populations.get_group_total(i));
147  }
148 }
149 
154 template <typename FP>
155 void draw_sample_infection(Model<FP>& model)
156 {
157  model.parameters.template get<Seasonality<FP>>().draw_sample();
158 
159  //not age dependent
160  model.parameters.template get<TimeExposed<FP>>()[AgeGroup(0)].draw_sample();
161  model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[AgeGroup(0)].draw_sample();
162  model.parameters.template get<TimeInfectedSymptoms<FP>>()[AgeGroup(0)].draw_sample();
163  model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[AgeGroup(0)].draw_sample();
164  model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)].draw_sample();
165  model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)].draw_sample();
166 
167  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
168  //not age dependent
169  model.parameters.template get<TimeExposed<FP>>()[i] =
170  model.parameters.template get<TimeExposed<FP>>()[AgeGroup(0)];
171  model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[i] =
172  model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[AgeGroup(0)];
173  model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[i] =
174  model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[AgeGroup(0)];
175  model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[i] =
176  model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)];
177  model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i] =
178  model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)];
179 
180  //age dependent
181  model.parameters.template get<TimeInfectedSymptoms<FP>>()[i].draw_sample();
182  model.parameters.template get<TimeInfectedSevere<FP>>()[i].draw_sample();
183  model.parameters.template get<TimeInfectedCritical<FP>>()[i].draw_sample();
184 
185  model.parameters.template get<TransmissionProbabilityOnContact<FP>>()[i].draw_sample();
186  model.parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i].draw_sample();
187  model.parameters.template get<DeathsPerCritical<FP>>()[i].draw_sample();
188  model.parameters.template get<SeverePerInfectedSymptoms<FP>>()[i].draw_sample();
189  model.parameters.template get<CriticalPerSevere<FP>>()[i].draw_sample();
190  }
191 }
192 
197 template <typename FP>
198 void draw_sample(Model<FP>& model)
199 {
200  draw_sample_infection(model);
201  draw_sample_demographics(model);
202  model.parameters.template get<ContactPatterns<FP>>().draw_sample();
203  model.apply_constraints();
204 }
205 
206 template <typename FP>
207 Graph<Model<FP>, MobilityParameters<FP>> draw_sample(Graph<Model<FP>, MobilityParameters<FP>>& graph)
208 {
209  Graph<Model<FP>, MobilityParameters<FP>> sampled_graph;
210 
211  //sample global parameters
212  auto& shared_params_model = graph.nodes()[0].property;
213  draw_sample_infection(shared_params_model);
214  auto& shared_contacts = shared_params_model.parameters.template get<ContactPatterns<FP>>();
215  shared_contacts.draw_sample_dampings();
216  auto& shared_dynamic_npis = shared_params_model.parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
217  shared_dynamic_npis.draw_sample();
218 
219  for (auto& params_node : graph.nodes()) {
220  auto& node_model = params_node.property;
221 
222  //sample local parameters
223  draw_sample_demographics(params_node.property);
224 
225  //copy global parameters
226  //save demographic parameters so they aren't overwritten
227  auto local_icu_capacity = node_model.parameters.template get<ICUCapacity<FP>>();
228  auto local_tnt_capacity = node_model.parameters.template get<TestAndTraceCapacity<FP>>();
229  auto local_holidays = node_model.parameters.template get<ContactPatterns<FP>>().get_school_holidays();
230  node_model.parameters = shared_params_model.parameters;
231  node_model.parameters.template get<ICUCapacity<FP>>() = local_icu_capacity;
232  node_model.parameters.template get<TestAndTraceCapacity<FP>>() = local_tnt_capacity;
233  node_model.parameters.template get<ContactPatterns<FP>>().get_school_holidays() = local_holidays;
234 
235  node_model.parameters.template get<ContactPatterns<FP>>().make_matrix();
236  node_model.apply_constraints();
237 
238  sampled_graph.add_node(params_node.id, node_model);
239  }
240 
241  for (auto& edge : graph.edges()) {
242  auto edge_params = edge.property;
243  apply_dampings(edge_params.get_coefficients(), shared_contacts.get_dampings(), [&edge_params](auto& v) {
244  return make_mobility_damping_vector(edge_params.get_coefficients().get_shape(), v);
245  });
246  edge_params.set_dynamic_npis_infected(shared_dynamic_npis);
247  sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge_params);
248  }
249 
250  return sampled_graph;
251 }
252 
253 } // namespace osecir
254 } // namespace mio
255 
256 #endif // ODESECIR_PARAMETER_SPACE_H
mio::Graph
generic graph structure
Definition: graph.h:148
mio::Graph::add_edge
requires std::constructible_from< EdgePropertyT, Args... > void add_edge(size_t start_node_idx, size_t end_node_idx, Args &&... args)
add an edge to the graph.
Definition: graph.h:233
mio::Graph::add_node
requires std::constructible_from< NodePropertyT, Args... > void add_node(int id, Args &&... args)
add a node to the graph.
Definition: graph.h:218
mio::Index
An Index with more than one template parameter combines several Index objects.
Definition: index.h:191
mio::MobilityParameters
parameters that influence mobility.
Definition: metapopulation_mobility_instant.h:122
mio::ParameterDistributionConstant
Definition: parameter_distributions.h:813
mio::ParameterDistributionNormal
Definition: parameter_distributions.h:160
mio::SimulationTime
Typesafe wrapper for a floating-point simulation time value (in days).
Definition: damping.h:78
mio::osecir::Model
Definition: ode_secir/model.h:64
mio::UncertainValue::value
FP value() const
Conversion to scalar by returning the scalar contained in UncertainValue.
Definition: uncertain_value.h:118
metapopulation_mobility_instant.h
ad::internal::min
static min_max_return_type< ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > >::type min(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &a, const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &b)
Definition: ad.hpp:2599
ad::internal::max
static min_max_return_type< ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 >, ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > >::type max(const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &a, const ad::internal::active_type< AD_TAPE_REAL, DATA_HANDLER_1 > &b)
Definition: ad.hpp:2596
mio::osecir::draw_sample
void draw_sample(Model< FP > &model)
Draws a sample from Model parameter distributions and stores sample values as SecirParams parameter v...
Definition: ode_secir/parameter_space.h:198
mio::osecir::draw_sample_demographics
void draw_sample_demographics(Model< FP > &model)
draws a sample from the specified distributions for all parameters related to the demographics,...
Definition: ode_secir/parameter_space.h:125
mio::osecir::set_params_distributions_normal
void set_params_distributions_normal(Model< FP > &model, FP t0, FP tmax, FP dev_rel)
Sets alls SecirParams parameters normally distributed, using the current value and a given standard d...
Definition: ode_secir/parameter_space.h:47
mio::osecir::draw_sample_infection
void draw_sample_infection(Model< FP > &model)
draws a sample from the specified distributions for all parameters related to the infection.
Definition: ode_secir/parameter_space.h:155
mio
A collection of classes to simplify handling of matrix shapes in meta programming.
Definition: models/abm/analyze_result.h:30
mio::apply_dampings
void apply_dampings(DampingExpression &damping_expression, const DampingSamplings &dampings, F make_matrix)
add sampled dampings to a damping expression.
Definition: damping_sampling.h:256
mio::i
auto i
Definition: io.h:810
mio::log_debug
void log_debug(spdlog::string_view_t fmt, const Args &... args)
Definition: logging.h:132
mio::get
constexpr std::tuple_element< I, std::tuple< Index< CategoryTags >... > >::type & get(Index< CategoryTags... > &i) noexcept
Retrieves the Index (by reference) at the Ith position of a MultiIndex.
Definition: index.h:304
model.h
parameter_distributions.h
mio::AgeGroup
Typesafe index representing an age group.
Definition: age_group.h:40
mio::CompartmentalModel::apply_constraints
bool apply_constraints()
Checks whether the model satisfies all constraints.
Definition: compartmental_model.h:131
mio::CompartmentalModel::populations
Populations populations
Definition: compartmental_model.h:156
mio::CompartmentalModel::parameters
ParameterSet parameters
Definition: compartmental_model.h:157
mio::Limits
Type specific limits for floating-points.
Definition: memilio/config.h:59
mio::osecir::ContactPatterns
the contact patterns within the society are modelled using an UncertainContactMatrix
Definition: ode_secir/parameters.h:326
mio::osecir::CriticalPerSevere
the percentage of ICU patients per hospitalized patients in the SECIR model
Definition: ode_secir/parameters.h:276
mio::osecir::DeathsPerCritical
the percentage of dead patients per ICU patients in the SECIR model
Definition: ode_secir/parameters.h:310
mio::osecir::ICUCapacity
the icu capacity in the SECIR model
Definition: ode_secir/parameters.h:80
mio::osecir::MaxRiskOfInfectionFromSymptomatic
risk of infection from symptomatic cases increases as test and trace capacity is exceeded.
Definition: ode_secir/parameters.h:244
mio::osecir::RecoveredPerInfectedNoSymptoms
the percentage of asymptomatic cases in the SECIR model
Definition: ode_secir/parameters.h:212
mio::osecir::RelativeTransmissionNoSymptoms
the relative InfectedNoSymptoms infectability
Definition: ode_secir/parameters.h:196
mio::osecir::RiskOfInfectionFromSymptomatic
the risk of infection from symptomatic cases in the SECIR model
Definition: ode_secir/parameters.h:228
mio::osecir::Seasonality
the seasonality in the SECIR model the seasonality is given as (1+k*sin()) where the sine curve is be...
Definition: ode_secir/parameters.h:64
mio::osecir::SeverePerInfectedSymptoms
the percentage of hospitalized patients per infected patients in the SECIR model
Definition: ode_secir/parameters.h:260
mio::osecir::TestAndTraceCapacity
capacity to test and trace contacts of infected for quarantine per day.
Definition: ode_secir/parameters.h:358
mio::osecir::TimeExposed
the (mean) latent time in day unit
Definition: ode_secir/parameters.h:96
mio::osecir::TimeInfectedCritical
the time people are treated by ICU before returning home in the SECIR model in day unit
Definition: ode_secir/parameters.h:164
mio::osecir::TimeInfectedNoSymptoms
the (mean) time in day unit for asymptomatic cases that are infectious but have not yet developed sym...
Definition: ode_secir/parameters.h:113
mio::osecir::TimeInfectedSevere
the time people are 'simply' hospitalized before returning home in the SECIR model in day unit
Definition: ode_secir/parameters.h:147
mio::osecir::TimeInfectedSymptoms
the infectious time for symptomatic cases that are infected but who do not need to be hsopitalized in...
Definition: ode_secir/parameters.h:130
mio::osecir::TransmissionProbabilityOnContact
probability of getting infected from a contact
Definition: ode_secir/parameters.h:180
Generated by doxygen 1.9.1