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  set_distribution(model.parameters.template get<DynamicNPIsImplementationDelay<FP>>());
74 
75  // populations
76  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
77  for (auto j = Index<InfectionState>(0); j < Index<InfectionState>(InfectionState::Count); j++) {
78 
79  // don't touch S and D
80  if (j == InfectionState::Susceptible || j == InfectionState::Dead) {
81  continue;
82  }
83 
84  // variably sized groups
85  set_distribution(model.populations[{i, j}], 0.0);
86  }
87  }
88 
89  // times
90  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
91 
92  set_distribution(model.parameters.template get<TimeExposed<FP>>()[i]);
93  set_distribution(model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[i]);
94  set_distribution(model.parameters.template get<TimeInfectedSymptoms<FP>>()[i]);
95  set_distribution(model.parameters.template get<TimeInfectedSevere<FP>>()[i]);
96  set_distribution(model.parameters.template get<TimeInfectedCritical<FP>>()[i]);
97 
98  set_distribution(model.parameters.template get<TransmissionProbabilityOnContact<FP>>()[i]);
99  set_distribution(model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[i]);
100  set_distribution(model.parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]);
101  set_distribution(model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[i]);
102  set_distribution(model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i]);
103  set_distribution(model.parameters.template get<DeathsPerCritical<FP>>()[i]);
104  set_distribution(model.parameters.template get<SeverePerInfectedSymptoms<FP>>()[i]);
105  set_distribution(model.parameters.template get<CriticalPerSevere<FP>>()[i]);
106  }
107 
108  // dampings
109  auto matrices = std::vector<size_t>();
110  for (size_t i = 0; i < model.parameters.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_num_matrices();
111  ++i) {
112  matrices.push_back(i);
113  }
114  auto groups = Eigen::VectorX<FP>::Constant(Eigen::Index(model.parameters.get_num_groups().get()), 1.0);
115  model.parameters.template get<ContactPatterns<FP>>().get_dampings().emplace_back(
116  mio::UncertainValue<FP>(0.5), mio::DampingLevel(0), mio::DampingType(0),
117  mio::SimulationTime<FP>(t0 + (tmax - t0) / 2), matrices, groups);
118  set_distribution(model.parameters.template get<ContactPatterns<FP>>().get_dampings()[0].get_value(), 0.0);
119 }
120 
125 template <typename FP>
126 void draw_sample_demographics(Model<FP>& model)
127 {
128  model.parameters.template get<ICUCapacity<FP>>().draw_sample();
129  model.parameters.template get<TestAndTraceCapacity<FP>>().draw_sample();
130 
131  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
132  FP group_total = model.populations.get_group_total(i);
133 
134  model.populations[{i, InfectionState::Exposed}].draw_sample();
135  model.populations[{i, InfectionState::InfectedNoSymptoms}].draw_sample();
136  model.populations[{i, InfectionState::InfectedSymptoms}].draw_sample();
137  model.populations[{i, InfectionState::InfectedSevere}].draw_sample();
138  model.populations[{i, InfectionState::InfectedCritical}].draw_sample();
139  model.populations[{i, InfectionState::Recovered}].draw_sample();
140 
141  // no sampling for dead and total numbers
142  // [...]
143 
144  model.populations.template set_difference_from_group_total<AgeGroup>({i, InfectionState::Susceptible},
145  group_total);
146  model.populations.template set_difference_from_group_total<AgeGroup>({i, InfectionState::Susceptible},
147  model.populations.get_group_total(i));
148  }
149 }
150 
155 template <typename FP>
156 void draw_sample_infection(Model<FP>& model)
157 {
158  model.parameters.template get<Seasonality<FP>>().draw_sample();
159 
160  //not age dependent
161  model.parameters.template get<TimeExposed<FP>>()[AgeGroup(0)].draw_sample();
162  model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[AgeGroup(0)].draw_sample();
163  model.parameters.template get<TimeInfectedSymptoms<FP>>()[AgeGroup(0)].draw_sample();
164  model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[AgeGroup(0)].draw_sample();
165  model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)].draw_sample();
166  model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)].draw_sample();
167 
168  for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
169  //not age dependent
170  model.parameters.template get<TimeExposed<FP>>()[i] =
171  model.parameters.template get<TimeExposed<FP>>()[AgeGroup(0)];
172  model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[i] =
173  model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[AgeGroup(0)];
174  model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[i] =
175  model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[AgeGroup(0)];
176  model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[i] =
177  model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)];
178  model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i] =
179  model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)];
180 
181  //age dependent
182  model.parameters.template get<TimeInfectedSymptoms<FP>>()[i].draw_sample();
183  model.parameters.template get<TimeInfectedSevere<FP>>()[i].draw_sample();
184  model.parameters.template get<TimeInfectedCritical<FP>>()[i].draw_sample();
185 
186  model.parameters.template get<TransmissionProbabilityOnContact<FP>>()[i].draw_sample();
187  model.parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i].draw_sample();
188  model.parameters.template get<DeathsPerCritical<FP>>()[i].draw_sample();
189  model.parameters.template get<SeverePerInfectedSymptoms<FP>>()[i].draw_sample();
190  model.parameters.template get<CriticalPerSevere<FP>>()[i].draw_sample();
191  }
192 }
193 
198 template <typename FP>
199 void draw_sample(Model<FP>& model)
200 {
201  draw_sample_infection(model);
202  draw_sample_demographics(model);
203  model.parameters.template get<ContactPatterns<FP>>().draw_sample();
204  model.apply_constraints();
205 }
206 
207 template <typename FP>
208 Graph<Model<FP>, MobilityParameters<FP>> draw_sample(Graph<Model<FP>, MobilityParameters<FP>>& graph)
209 {
210  Graph<Model<FP>, MobilityParameters<FP>> sampled_graph;
211 
212  //sample global parameters
213  auto& shared_params_model = graph.nodes()[0].property;
214  draw_sample_infection(shared_params_model);
215  auto& shared_contacts = shared_params_model.parameters.template get<ContactPatterns<FP>>();
216  shared_contacts.draw_sample_dampings();
217  auto& shared_dynamic_npis = shared_params_model.parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
218  shared_dynamic_npis.draw_sample();
219  auto& shared_dynamic_npis_delay = shared_params_model.parameters.template get<DynamicNPIsImplementationDelay<FP>>();
220  shared_dynamic_npis_delay.draw_sample();
221 
222  for (auto& params_node : graph.nodes()) {
223  auto& node_model = params_node.property;
224 
225  //sample local parameters
226  draw_sample_demographics(params_node.property);
227 
228  //copy global parameters
229  //save demographic parameters so they aren't overwritten
230  auto local_icu_capacity = node_model.parameters.template get<ICUCapacity<FP>>();
231  auto local_tnt_capacity = node_model.parameters.template get<TestAndTraceCapacity<FP>>();
232  auto local_holidays = node_model.parameters.template get<ContactPatterns<FP>>().get_school_holidays();
233  node_model.parameters = shared_params_model.parameters;
234  node_model.parameters.template get<ICUCapacity<FP>>() = local_icu_capacity;
235  node_model.parameters.template get<TestAndTraceCapacity<FP>>() = local_tnt_capacity;
236  node_model.parameters.template get<ContactPatterns<FP>>().get_school_holidays() = local_holidays;
237 
238  node_model.parameters.template get<ContactPatterns<FP>>().make_matrix();
239  node_model.apply_constraints();
240 
241  sampled_graph.add_node(params_node.id, node_model);
242  }
243 
244  for (auto& edge : graph.edges()) {
245  auto edge_params = edge.property;
246  apply_dampings(edge_params.get_coefficients(), shared_contacts.get_dampings(), [&edge_params](auto& v) {
247  return make_mobility_damping_vector(edge_params.get_coefficients().get_shape(), v);
248  });
249  edge_params.set_dynamic_npis_infected(shared_dynamic_npis);
250  sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge_params);
251  }
252 
253  return sampled_graph;
254 }
255 
256 } // namespace osecir
257 } // namespace mio
258 
259 #endif // ODESECIR_PARAMETER_SPACE_H
mio::Graph
generic graph structure
Definition: graph.h:153
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:238
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:223
mio::Index
An Index with more than one template parameter combines several Index objects.
Definition: index.h:181
mio::MobilityParameters
parameters that influence mobility.
Definition: metapopulation_mobility_instant.h:123
mio::ParameterDistributionConstant
Definition: parameter_distributions.h:813
mio::ParameterDistributionNormal
Definition: parameter_distributions.h:160
mio::SimulationTime
double simulation time.
Definition: damping.h:58
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:199
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:126
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:156
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:809
mio::log_debug
void log_debug(spdlog::string_view_t fmt, const Args &... args)
Definition: logging.h:118
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:294
model.h
parameter_distributions.h
mio::AgeGroup
The AgeGroup struct is used as a dynamically sized tag for all age dependent categories.
Definition: age_group.h:32
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:308
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:292
mio::osecir::DynamicNPIsImplementationDelay
The delay with which DynamicNPIs are implemented and enforced after exceedance of threshold.
Definition: ode_secir/parameters.h:340
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:356
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