parameters_io.h Source File

CPP API: parameters_io.h Source File
models/ode_secir/parameters_io.h
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2020-2026 MEmilio
3 *
4 * Authors: Wadim Koslow, Daniel Abele, Martin J. Kuehn, 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 #ifndef ODESECIR_PARAMETERS_IO_H
21 #define ODESECIR_PARAMETERS_IO_H
22 
23 #include "memilio/config.h"
24 #include <cassert>
25 
26 #ifdef MEMILIO_HAS_JSONCPP
27 
28 #include "ode_secir/model.h"
29 #include "memilio/io/epi_data.h"
30 #include "memilio/io/parameters_io.h"
31 #include "memilio/io/result_io.h"
32 #include "memilio/math/math_utils.h"
33 
34 namespace mio
35 {
36 
37 namespace osecir
38 {
39 
40 namespace details
41 {
52 IOResult<void> read_confirmed_cases_data(
53  std::vector<ConfirmedCasesDataEntry>& rki_data, std::vector<int> const& vregion, Date date,
54  std::vector<std::vector<double>>& vnum_Exposed, std::vector<std::vector<double>>& vnum_InfectedNoSymptoms,
55  std::vector<std::vector<double>>& vnum_InfectedSymptoms, std::vector<std::vector<double>>& vnum_InfectedSevere,
56  std::vector<std::vector<double>>& vnum_icu, std::vector<std::vector<double>>& vnum_death,
57  std::vector<std::vector<double>>& vnum_rec, const std::vector<std::vector<int>>& vt_Exposed,
58  const std::vector<std::vector<int>>& vt_InfectedNoSymptoms,
59  const std::vector<std::vector<int>>& vt_InfectedSymptoms, const std::vector<std::vector<int>>& vt_InfectedSevere,
60  const std::vector<std::vector<int>>& vt_InfectedCritical, const std::vector<std::vector<double>>& vmu_C_R,
61  const std::vector<std::vector<double>>& vmu_I_H, const std::vector<std::vector<double>>& vmu_H_U,
62  const std::vector<double>& scaling_factor_inf);
63 
73 template <typename FP>
74 IOResult<void> set_confirmed_cases_data(std::vector<Model<FP>>& model, std::vector<ConfirmedCasesDataEntry>& case_data,
75  const std::vector<int>& region, Date date,
76  const std::vector<double>& scaling_factor_inf)
77 {
78  const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size();
79  // allow single scalar scaling that is broadcast to all age groups
80  assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups);
81 
82  // Set scaling factors to match num age groups
83  std::vector<double> scaling_factor_inf_full;
84  if (scaling_factor_inf.size() == 1) {
85  scaling_factor_inf_full.assign(num_age_groups, scaling_factor_inf[0]);
86  }
87  else {
88  scaling_factor_inf_full = scaling_factor_inf;
89  }
90 
91  std::vector<std::vector<int>> t_InfectedNoSymptoms{model.size()};
92  std::vector<std::vector<int>> t_Exposed{model.size()};
93  std::vector<std::vector<int>> t_InfectedSymptoms{model.size()};
94  std::vector<std::vector<int>> t_InfectedSevere{model.size()};
95  std::vector<std::vector<int>> t_InfectedCritical{model.size()};
96 
97  std::vector<std::vector<double>> mu_C_R{model.size()};
98  std::vector<std::vector<double>> mu_I_H{model.size()};
99  std::vector<std::vector<double>> mu_H_U{model.size()};
100  std::vector<std::vector<double>> mu_U_D{model.size()};
101 
102  for (size_t node = 0; node < model.size(); ++node) {
103  const size_t model_groups = (size_t)model[node].parameters.get_num_groups();
104  assert(model_groups == 1 || model_groups == num_age_groups);
105  for (size_t ag = 0; ag < num_age_groups; ag++) {
106  // If the model has fewer groups than casedata entries available,
107  // reuse group 0 parameters for all RKI age groups
108  const size_t group = (model_groups == num_age_groups) ? ag : 0;
109 
110  t_Exposed[node].push_back(
111  static_cast<int>(std::round(model[node].parameters.template get<TimeExposed<FP>>()[(AgeGroup)group])));
112  t_InfectedNoSymptoms[node].push_back(static_cast<int>(
113  std::round(model[node].parameters.template get<TimeInfectedNoSymptoms<FP>>()[(AgeGroup)group])));
114  t_InfectedSymptoms[node].push_back(static_cast<int>(
115  std::round(model[node].parameters.template get<TimeInfectedSymptoms<FP>>()[(AgeGroup)group])));
116  t_InfectedSevere[node].push_back(static_cast<int>(
117  std::round(model[node].parameters.template get<TimeInfectedSevere<FP>>()[(AgeGroup)group])));
118  t_InfectedCritical[node].push_back(static_cast<int>(
119  std::round(model[node].parameters.template get<TimeInfectedCritical<FP>>()[(AgeGroup)group])));
120 
121  mu_C_R[node].push_back(
122  model[node].parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[(AgeGroup)group]);
123  mu_I_H[node].push_back(
124  model[node].parameters.template get<SeverePerInfectedSymptoms<FP>>()[(AgeGroup)group]);
125  mu_H_U[node].push_back(model[node].parameters.template get<CriticalPerSevere<FP>>()[(AgeGroup)group]);
126  mu_U_D[node].push_back(model[node].parameters.template get<DeathsPerCritical<FP>>()[(AgeGroup)group]);
127  }
128  }
129  std::vector<std::vector<double>> num_InfectedSymptoms(model.size(), std::vector<double>(num_age_groups, 0.0));
130  std::vector<std::vector<double>> num_death(model.size(), std::vector<double>(num_age_groups, 0.0));
131  std::vector<std::vector<double>> num_rec(model.size(), std::vector<double>(num_age_groups, 0.0));
132  std::vector<std::vector<double>> num_Exposed(model.size(), std::vector<double>(num_age_groups, 0.0));
133  std::vector<std::vector<double>> num_InfectedNoSymptoms(model.size(), std::vector<double>(num_age_groups, 0.0));
134  std::vector<std::vector<double>> num_InfectedSevere(model.size(), std::vector<double>(num_age_groups, 0.0));
135  std::vector<std::vector<double>> num_icu(model.size(), std::vector<double>(num_age_groups, 0.0));
136 
137  BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
138  num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
139  t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
140  t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf_full));
141 
142  for (size_t node = 0; node < model.size(); node++) {
143 
144  if (std::accumulate(num_InfectedSymptoms[node].begin(), num_InfectedSymptoms[node].end(), FP(0.0),
145  [](const FP& a, const FP& b) {
146  return evaluate_intermediate<FP>(a + b);
147  }) > 0.0) {
148  size_t num_groups = (size_t)model[node].parameters.get_num_groups();
149  if (num_groups == num_age_groups) {
150  for (size_t i = 0; i < num_groups; i++) {
151  model[node].populations[{AgeGroup(i), InfectionState::Exposed}] = num_Exposed[node][i];
152  model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptoms}] =
153  num_InfectedNoSymptoms[node][i];
154  model[node].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsConfirmed}] = 0;
155  model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptoms}] =
156  num_InfectedSymptoms[node][i];
157  model[node].populations[{AgeGroup(i), InfectionState::InfectedSymptomsConfirmed}] = 0;
158  model[node].populations[{AgeGroup(i), InfectionState::InfectedSevere}] =
159  num_InfectedSevere[node][i];
160  // Only set the number of ICU patients here, if the date is not available in the data.
161  if (!is_divi_data_available(date)) {
162  model[node].populations[{AgeGroup(i), InfectionState::InfectedCritical}] = num_icu[node][i];
163  }
164  model[node].populations[{AgeGroup(i), InfectionState::Dead}] = num_death[node][i];
165  model[node].populations[{AgeGroup(i), InfectionState::Recovered}] = num_rec[node][i];
166  }
167  }
168  else {
169  const auto sum_vec = [](const std::vector<FP>& v) {
170  return std::accumulate(v.begin(), v.end(), FP(0.0), [](const FP& a, const FP& b) {
171  return evaluate_intermediate<FP>(a + b);
172  });
173  };
174  const size_t i0 = 0;
175  model[node].populations[{AgeGroup(i0), InfectionState::Exposed}] = sum_vec(num_Exposed[node]);
176  model[node].populations[{AgeGroup(i0), InfectionState::InfectedNoSymptoms}] =
177  sum_vec(num_InfectedNoSymptoms[node]);
178  model[node].populations[{AgeGroup(i0), InfectionState::InfectedNoSymptomsConfirmed}] = 0;
179  model[node].populations[{AgeGroup(i0), InfectionState::InfectedSymptoms}] =
180  sum_vec(num_InfectedSymptoms[node]);
181  model[node].populations[{AgeGroup(i0), InfectionState::InfectedSymptomsConfirmed}] = 0;
182  model[node].populations[{AgeGroup(i0), InfectionState::InfectedSevere}] =
183  sum_vec(num_InfectedSevere[node]);
184  if (!is_divi_data_available(date)) {
185  model[node].populations[{AgeGroup(i0), InfectionState::InfectedCritical}] = sum_vec(num_icu[node]);
186  }
187  model[node].populations[{AgeGroup(i0), InfectionState::Dead}] = sum_vec(num_death[node]);
188  model[node].populations[{AgeGroup(i0), InfectionState::Recovered}] = sum_vec(num_rec[node]);
189  }
190  }
191  else {
192  log_warning("No infections reported on date {} for region {}. Population data has not been set.", date,
193  region[node]);
194  }
195  }
196  return success();
197 }
198 
209 template <typename FP>
210 IOResult<void> set_confirmed_cases_data(std::vector<Model<FP>>& model, const std::string& path,
211  std::vector<int> const& region, Date date,
212  const std::vector<double>& scaling_factor_inf)
213 {
214  BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path));
215  BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf));
216  return success();
217 }
218 
228 template <typename FP>
229 IOResult<void> set_divi_data(std::vector<Model<FP>>& model, const std::string& path, const std::vector<int>& vregion,
230  Date date, double scaling_factor_icu)
231 {
232  // DIVI dataset will no longer be updated from CW29 2024 on.
233  if (!is_divi_data_available(date)) {
234  log_warning("No DIVI data available for date: {}. "
235  "ICU compartment will be set based on Case data.",
236  date);
237  return success();
238  }
239  std::vector<double> sum_mu_I_U(vregion.size(), 0);
240  std::vector<std::vector<double>> mu_I_U{model.size()};
241  for (size_t region = 0; region < vregion.size(); region++) {
242  auto num_groups = model[region].parameters.get_num_groups();
243  for (auto i = AgeGroup(0); i < num_groups; i++) {
244  sum_mu_I_U[region] += model[region].parameters.template get<CriticalPerSevere<FP>>()[i] *
245  model[region].parameters.template get<SeverePerInfectedSymptoms<FP>>()[i];
246  mu_I_U[region].push_back(model[region].parameters.template get<CriticalPerSevere<FP>>()[i] *
247  model[region].parameters.template get<SeverePerInfectedSymptoms<FP>>()[i]);
248  }
249  }
250  std::vector<double> num_icu(model.size(), 0.0);
251  BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu));
252 
253  for (size_t region = 0; region < vregion.size(); region++) {
254  auto num_groups = model[region].parameters.get_num_groups();
255  for (auto i = AgeGroup(0); i < num_groups; i++) {
256  model[region].populations[{i, InfectionState::InfectedCritical}] =
257  scaling_factor_icu * num_icu[region] * mu_I_U[region][(size_t)i] / sum_mu_I_U[region];
258  }
259  }
260 
261  return success();
262 }
263 
271 template <typename FP>
272 IOResult<void> set_population_data(std::vector<Model<FP>>& model, const std::vector<std::vector<FP>>& num_population,
273  const std::vector<int>& vregion)
274 {
275  assert(num_population.size() == vregion.size());
276  assert(model.size() == vregion.size());
277  for (size_t region = 0; region < vregion.size(); region++) {
278  const auto model_groups = (size_t)model[region].parameters.get_num_groups();
279  const auto data_groups = num_population[region].size();
280  assert(data_groups == model_groups || (model_groups == 1 && data_groups >= 1));
281 
282  if (data_groups == model_groups) {
283  for (auto i = AgeGroup(0); i < model[region].parameters.get_num_groups(); i++) {
284  model[region].populations.template set_difference_from_group_total<AgeGroup>(
285  {i, InfectionState::Susceptible}, num_population[region][(size_t)i]);
286  }
287  }
288  else if (model_groups == 1 && data_groups >= 1) {
289  const FP total = std::accumulate(num_population[region].begin(), num_population[region].end(), FP(0.0),
290  [](const FP& a, const FP& b) {
291  return evaluate_intermediate<FP>(a + b);
292  });
293  model[region].populations.template set_difference_from_group_total<AgeGroup>(
294  {AgeGroup(0), InfectionState::Susceptible}, total);
295  }
296  }
297  return success();
298 }
299 
307 template <typename FP>
308 IOResult<void> set_population_data(std::vector<Model<FP>>& model, const std::string& path,
309  const std::vector<int>& vregion)
310 {
311  BOOST_OUTCOME_TRY(const auto&& num_population, read_population_data(path, vregion));
312  BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion));
313  return success();
314 }
315 
316 } //namespace details
317 
318 #ifdef MEMILIO_HAS_HDF5
319 
336 template <class Model>
337 IOResult<void> export_input_data_county_timeseries(
338  std::vector<Model> models, const std::string& results_dir, std::vector<int> const& region, Date date,
339  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu, int num_days,
340  const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path)
341 {
342  const auto num_age_groups = (size_t)models[0].parameters.get_num_groups();
343  // allow scalar scaling factor as convenience for 1-group models
344  assert(scaling_factor_inf.size() == 1 || scaling_factor_inf.size() == num_age_groups);
345  assert(models.size() == region.size());
346  std::vector<TimeSeries<double>> extrapolated_data(
347  region.size(), TimeSeries<double>::zero(num_days + 1, (size_t)InfectionState::Count * num_age_groups));
348 
349  BOOST_OUTCOME_TRY(auto&& num_population, mio::read_population_data(population_data_path, region));
350  BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(confirmed_cases_path));
351 
352  for (int t = 0; t <= num_days; ++t) {
353  auto offset_day = offset_date_by_days(date, t);
354 
355  BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, region, offset_day, scaling_factor_icu));
356  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, region, offset_day, scaling_factor_inf));
357  BOOST_OUTCOME_TRY(details::set_population_data(models, num_population, region));
358  for (size_t r = 0; r < region.size(); r++) {
359  extrapolated_data[r][t] = models[r].get_initial_values();
360  }
361  }
362 
363  BOOST_OUTCOME_TRY(
364  save_result(extrapolated_data, region, (int)num_age_groups, path_join(results_dir, "Results_rki.h5")));
365 
366  auto rki_data_sum = sum_nodes(std::vector<std::vector<TimeSeries<double>>>{extrapolated_data});
367  BOOST_OUTCOME_TRY(
368  save_result({rki_data_sum[0][0]}, {0}, (int)num_age_groups, path_join(results_dir, "Results_rki_sum.h5")));
369 
370  return success();
371 }
372 #else
373 template <class Model>
374 IOResult<void> export_input_data_county_timeseries(
375  std::vector<Model> models, const std::string& results_dir, std::vector<int> const& region, Date date,
376  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu, int num_days,
377  const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path)
378 {
379  mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
380  return success();
381 }
382 #endif // MEMILIO_HAS_HDF5
383 
392 template <class Model>
393 IOResult<void> read_input_data_germany(std::vector<Model>& model, Date date,
394  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
395  const std::string& pydata_dir)
396 {
397  BOOST_OUTCOME_TRY(
398  details::set_divi_data(model, path_join(pydata_dir, "germany_divi.json"), {0}, date, scaling_factor_icu));
399  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_age_ma7.json"), {0},
400  date, scaling_factor_inf));
401  BOOST_OUTCOME_TRY(
402  details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), {0}));
403  return success();
404 }
405 
415 template <class Model>
416 IOResult<void> read_input_data_state(std::vector<Model>& model, Date date, std::vector<int>& state,
417  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
418  const std::string& pydata_dir)
419 {
420 
421  BOOST_OUTCOME_TRY(
422  details::set_divi_data(model, path_join(pydata_dir, "state_divi.json"), state, date, scaling_factor_icu));
423  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_state_age_ma7.json"),
424  state, date, scaling_factor_inf));
425  BOOST_OUTCOME_TRY(
426  details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), state));
427  return success();
428 }
429 
441 template <class Model>
442 IOResult<void> read_input_data_county(std::vector<Model>& model, Date date, const std::vector<int>& county,
443  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
444  const std::string& pydata_dir, int num_days = 0, bool export_time_series = false)
445 {
446  BOOST_OUTCOME_TRY(
447  details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, scaling_factor_icu));
448  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_county_age_ma7.json"),
449  county, date, scaling_factor_inf));
450  BOOST_OUTCOME_TRY(
451  details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"), county));
452 
453  if (export_time_series) {
454  // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
455  // Run time equals run time of the previous functions times the num_days !
456  // (This only represents the vectorization of the previous function over all simulation days...)
457  log_warning("Exporting time series of extrapolated real data. This may take some minutes. "
458  "For simulation runs over the same time period, deactivate it.");
459  BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
460  model, pydata_dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days,
461  path_join(pydata_dir, "county_divi_ma7.json"), path_join(pydata_dir, "cases_all_county_age_ma7.json"),
462  path_join(pydata_dir, "county_current_population.json")));
463  }
464  return success();
465 }
466 
477 template <class Model>
478 IOResult<void> read_input_data(std::vector<Model>& model, Date date, const std::vector<int>& node_ids,
479  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
480  const std::string& pydata_dir, int num_days = 0, bool export_time_series = false)
481 {
482  BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date,
483  scaling_factor_icu));
484  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), node_ids,
485  date, scaling_factor_inf));
486  BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"), node_ids));
487 
488  if (export_time_series) {
489  // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
490  // Run time equals run time of the previous functions times the num_days !
491  // (This only represents the vectorization of the previous function over all simulation days...)
492  log_warning("Exporting time series of extrapolated real data. This may take some minutes. "
493  "For simulation runs over the same time period, deactivate it.");
494  BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
495  model, pydata_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days,
496  path_join(pydata_dir, "critical_cases.json"), path_join(pydata_dir, "confirmed_cases.json"),
497  path_join(pydata_dir, "population_data.json")));
498  }
499  return success();
500 }
501 
502 } // namespace osecir
503 } // namespace mio
504 
505 #endif // MEMILIO_HAS_JSONCPP
506 
507 #endif // ODESECIR_PARAMETERS_IO_H
mio::TimeSeries::zero
static TimeSeries zero(Eigen::Index num_time_points, Eigen::Index num_elements)
constructs TimeSeries instance and initializes it with zeros
Definition: time_series.h:152
epi_data.h
math_utils.h
parameters_io.h
mio
A collection of classes to simplify handling of matrix shapes in meta programming.
Definition: models/abm/analyze_result.h:30
mio::end
requires details::IsElementReference< M > RowMajorIterator< M, false > end(M &m)
create a non-const end iterator for the matrix m.
Definition: eigen_util.h:449
mio::i
auto i
Definition: io.h:809
mio::log_warning
void log_warning(spdlog::string_view_t fmt, const Args &... args)
Definition: logging.h:112
mio::success
auto success()
Create an object that is implicitly convertible to a succesful IOResult<void>.
Definition: io.h:359
mio::begin
requires details::IsElementReference< M > RowMajorIterator< M, false > begin(M &m)
create a non-const iterator to first element of the matrix m.
Definition: eigen_util.h:421
mio::offset_date_by_days
Date offset_date_by_days(Date date, int offset_days)
Computes the new date corresponding to a given date and a offset in days.
Definition: date.h:242
mio::path_join
std::string path_join(String &&base, Strings &&... app)
join one ore more strings with path separators.
Definition: stl_util.h:268
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
mio::sum_nodes
std::vector< std::vector< TimeSeries< FP > > > sum_nodes(const std::vector< std::vector< TimeSeries< FP >>> &ensemble_result)
Definition: memilio/data/analyze_result.h:268
model.h
result_io.h
Generated by doxygen 1.9.1