parameters_io.h Source File

CPP API: parameters_io.h Source File
models/ode_secirvvs/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. Kühn
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 MIO_ODE_SECIRVVS_PARAMETERS_IO_H
21 #define MIO_ODE_SECIRVVS_PARAMETERS_IO_H
22 
23 #include "memilio/config.h"
24 
25 #ifdef MEMILIO_HAS_JSONCPP
26 
27 #include "ode_secirvvs/model.h"
28 #include "memilio/mobility/graph.h"
29 #include "memilio/io/epi_data.h"
30 #include "memilio/io/parameters_io.h"
31 #include "memilio/io/io.h"
32 #include "memilio/io/result_io.h"
33 #include "memilio/utils/date.h"
34 #include "memilio/utils/logging.h"
35 #include "memilio/math/math_utils.h"
36 
37 namespace mio
38 {
39 namespace osecirvvs
40 {
41 
42 namespace details
43 {
56 IOResult<void> read_confirmed_cases_data(
57  std::string const& path, std::vector<int> const& vregion, Date date, std::vector<std::vector<double>>& num_Exposed,
58  std::vector<std::vector<double>>& num_InfectedNoSymptoms, std::vector<std::vector<double>>& num_InfectedSymptoms,
59  std::vector<std::vector<double>>& num_InfectedSevere, std::vector<std::vector<double>>& num_icu,
60  std::vector<std::vector<double>>& num_death, std::vector<std::vector<double>>& num_rec,
61  const std::vector<std::vector<int>>& t_Exposed, const std::vector<std::vector<int>>& t_InfectedNoSymptoms,
62  const std::vector<std::vector<int>>& t_InfectedSymptoms, const std::vector<std::vector<int>>& t_InfectedSevere,
63  const std::vector<std::vector<int>>& t_InfectedCritical, const std::vector<std::vector<double>>& mu_C_R,
64  const std::vector<std::vector<double>>& mu_I_H, const std::vector<std::vector<double>>& mu_H_U,
65  const std::vector<double>& scaling_factor_inf);
66 
67 IOResult<void> read_confirmed_cases_data(
68  const std::vector<ConfirmedCasesDataEntry>& rki_data, std::vector<int> const& vregion, Date date,
69  std::vector<std::vector<double>>& num_Exposed, std::vector<std::vector<double>>& num_InfectedNoSymptoms,
70  std::vector<std::vector<double>>& num_InfectedSymptoms, std::vector<std::vector<double>>& num_InfectedSevere,
71  std::vector<std::vector<double>>& num_icu, std::vector<std::vector<double>>& num_death,
72  std::vector<std::vector<double>>& num_rec, const std::vector<std::vector<int>>& t_Exposed,
73  const std::vector<std::vector<int>>& t_InfectedNoSymptoms, const std::vector<std::vector<int>>& t_InfectedSymptoms,
74  const std::vector<std::vector<int>>& t_InfectedSevere, const std::vector<std::vector<int>>& t_InfectedCritical,
75  const std::vector<std::vector<double>>& mu_C_R, const std::vector<std::vector<double>>& mu_I_H,
76  const std::vector<std::vector<double>>& mu_H_U, const std::vector<double>& scaling_factor_inf);
89 IOResult<void> read_confirmed_cases_data_fix_recovered(const std::vector<ConfirmedCasesDataEntry>& rki_data,
90  std::vector<int> const& vregion, Date date,
91  std::vector<std::vector<double>>& vnum_rec, double delay = 14.);
92 IOResult<void> read_confirmed_cases_data_fix_recovered(std::string const& path, std::vector<int> const& vregion,
93  Date date, std::vector<std::vector<double>>& vnum_rec,
94  double delay = 14.);
106 template <class Model>
107 IOResult<void> set_confirmed_cases_data(std::vector<Model>& model,
108  const std::vector<ConfirmedCasesDataEntry>& case_data,
109  std::vector<int> const& region, Date date,
110  const std::vector<double>& scaling_factor_inf, bool set_death = false)
111 {
112  auto num_age_groups = (size_t)model[0].parameters.get_num_groups();
113  assert(scaling_factor_inf.size() == num_age_groups); //TODO: allow vector or scalar valued scaling factors
114  assert(ConfirmedCasesDataEntry::age_group_names.size() == num_age_groups);
115 
116  std::vector<std::vector<int>> t_Exposed{model.size()};
117  std::vector<std::vector<int>> t_InfectedNoSymptoms{model.size()};
118  std::vector<std::vector<int>> t_InfectedSymptoms{model.size()};
119  std::vector<std::vector<int>> t_InfectedSevere{model.size()};
120  std::vector<std::vector<int>> t_InfectedCritical{model.size()};
121 
122  std::vector<std::vector<double>> mu_C_R{model.size()};
123  std::vector<std::vector<double>> mu_I_H{model.size()};
124  std::vector<std::vector<double>> mu_H_U{model.size()};
125 
126  std::vector<std::vector<double>> num_InfectedSymptoms(model.size());
127  std::vector<std::vector<double>> num_death(model.size());
128  std::vector<std::vector<double>> num_rec(model.size());
129  std::vector<std::vector<double>> num_Exposed(model.size());
130  std::vector<std::vector<double>> num_InfectedNoSymptoms(model.size());
131  std::vector<std::vector<double>> num_InfectedSevere(model.size());
132  std::vector<std::vector<double>> num_icu(model.size());
133 
134  /*----------- UNVACCINATED -----------*/
135  for (size_t county = 0; county < model.size(); county++) {
136  num_InfectedSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
137  num_death[county] = std::vector<double>(num_age_groups, 0.0);
138  num_rec[county] = std::vector<double>(num_age_groups, 0.0);
139  num_Exposed[county] = std::vector<double>(num_age_groups, 0.0);
140  num_InfectedNoSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
141  num_InfectedSevere[county] = std::vector<double>(num_age_groups, 0.0);
142  num_icu[county] = std::vector<double>(num_age_groups, 0.0);
143  for (size_t group = 0; group < num_age_groups; group++) {
144 
145  t_Exposed[county].push_back(static_cast<int>(
146  std::round(model[county].parameters.template get<TimeExposed<double>>()[(AgeGroup)group])));
147  t_InfectedNoSymptoms[county].push_back(static_cast<int>(
148  std::round(model[county].parameters.template get<TimeInfectedNoSymptoms<double>>()[(AgeGroup)group])));
149  t_InfectedSymptoms[county].push_back(static_cast<int>(
150  std::round(model[county].parameters.template get<TimeInfectedSymptoms<double>>()[(AgeGroup)group])));
151  t_InfectedSevere[county].push_back(static_cast<int>(
152  std::round(model[county].parameters.template get<TimeInfectedSevere<double>>()[(AgeGroup)group])));
153  t_InfectedCritical[county].push_back(static_cast<int>(
154  std::round(model[county].parameters.template get<TimeInfectedCritical<double>>()[(AgeGroup)group])));
155 
156  mu_C_R[county].push_back(
157  model[county].parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[(AgeGroup)group]);
158  mu_I_H[county].push_back(
159  model[county].parameters.template get<SeverePerInfectedSymptoms<double>>()[(AgeGroup)group]);
160  mu_H_U[county].push_back(
161  model[county].parameters.template get<CriticalPerSevere<double>>()[(AgeGroup)group]);
162  }
163  }
164 
165  BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
166  num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
167  t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
168  t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf));
169 
170  for (size_t county = 0; county < model.size(); county++) {
171  // if (std::accumulate(
172  // num_InfectedSymptoms[county].begin(),
173  // num_InfectedSymptoms[county].end(),
174  // double(0.0),
175  // [](const double& a, const double& b) { return evaluate_intermediate<double>(a + b); }
176  // ) > 0)
177  // {
178  size_t num_groups = (size_t)model[county].parameters.get_num_groups();
179  for (size_t i = 0; i < num_groups; i++) {
180  model[county].populations[{AgeGroup(i), InfectionState::ExposedNaive}] = num_Exposed[county][i];
181  model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsNaive}] =
182  num_InfectedNoSymptoms[county][i];
183  model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsNaiveConfirmed}] = 0;
184  model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaive}] =
185  num_InfectedSymptoms[county][i];
186  model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaiveConfirmed}] = 0;
187  model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereNaive}] =
188  num_InfectedSevere[county][i];
189  // Only set the number of ICU patients here, if the date is not available in the data.
190  if (!is_divi_data_available(date)) {
191  model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalNaive}] = num_icu[county][i];
192  }
193  model[county].populations[{AgeGroup(i), InfectionState::SusceptibleImprovedImmunity}] = num_rec[county][i];
194  if (set_death) {
195  // in set_confirmed_cases_data initilization, deaths are now set to 0. In order to visualize
196  // the extrapolated real number of deaths, they have to be set here. In the comparison of data
197  // it has to be paid attention to the fact, the the simulation starts with deaths=0
198  // while this method starts with deaths=number of reported deaths so far...
199  // Additionally, we set the number of reported deaths to DeadNaive since no information on that is
200  // available here.
201  // Do only add deaths after substraction.
202  model[county].populations[{AgeGroup(i), InfectionState::DeadNaive}] = num_death[county][i];
203  }
204  }
205 
206  // }
207  if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), double(0.0),
208  [](const double& a, const double& b) {
209  return evaluate_intermediate<double>(a + b);
210  }) == 0.0) {
211  log_warning(
212  "No infections for unvaccinated reported on date {} for region {}. Population data has not been set.",
213  date, region[county]);
214  }
215  }
216 
217  /*----------- PARTIALLY VACCINATED -----------*/
218  for (size_t county = 0; county < model.size(); county++) {
219  t_InfectedNoSymptoms[county].clear();
220  t_Exposed[county].clear();
221  t_InfectedSymptoms[county].clear();
222  t_InfectedSevere[county].clear();
223  t_InfectedCritical[county].clear();
224 
225  mu_C_R[county].clear();
226  mu_I_H[county].clear();
227  mu_H_U[county].clear();
228 
229  num_InfectedSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
230  num_death[county] = std::vector<double>(num_age_groups, 0.0);
231  num_rec[county] = std::vector<double>(num_age_groups, 0.0);
232  num_Exposed[county] = std::vector<double>(num_age_groups, 0.0);
233  num_InfectedNoSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
234  num_InfectedSevere[county] = std::vector<double>(num_age_groups, 0.0);
235  num_icu[county] = std::vector<double>(num_age_groups, 0.0);
236  for (size_t group = 0; group < num_age_groups; group++) {
237  double reduc_t = model[0].parameters.template get<ReducTimeInfectedMild<double>>()[(AgeGroup)group];
238  t_Exposed[county].push_back(static_cast<int>(
239  std::round(model[county].parameters.template get<TimeExposed<double>>()[(AgeGroup)group])));
240  t_InfectedNoSymptoms[county].push_back(static_cast<int>(std::round(
241  model[county].parameters.template get<TimeInfectedNoSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
242  t_InfectedSymptoms[county].push_back(static_cast<int>(std::round(
243  model[county].parameters.template get<TimeInfectedSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
244  t_InfectedSevere[county].push_back(static_cast<int>(
245  std::round(model[county].parameters.template get<TimeInfectedSevere<double>>()[(AgeGroup)group])));
246  t_InfectedCritical[county].push_back(static_cast<int>(
247  std::round(model[county].parameters.template get<TimeInfectedCritical<double>>()[(AgeGroup)group])));
248 
249  double exp_fac_part_immune =
250  model[county].parameters.template get<ReducExposedPartialImmunity<double>>()[(AgeGroup)group];
251  double inf_fac_part_immune =
252  model[county].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[(AgeGroup)group];
253  double hosp_fac_part_immune =
254  model[county]
255  .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[(AgeGroup)group];
256  double icu_fac_part_immune =
257  model[county]
258  .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[(AgeGroup)group];
259  mu_C_R[county].push_back((
260  1 - inf_fac_part_immune / exp_fac_part_immune *
261  (1 - model[county]
262  .parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[(AgeGroup)group])));
263  mu_I_H[county].push_back(
264  hosp_fac_part_immune / inf_fac_part_immune *
265  model[county].parameters.template get<SeverePerInfectedSymptoms<double>>()[(AgeGroup)group]);
266  // transfer from H to U, D unchanged.
267  mu_H_U[county].push_back(
268  icu_fac_part_immune / hosp_fac_part_immune *
269  model[county].parameters.template get<CriticalPerSevere<double>>()[(AgeGroup)group]);
270  }
271  }
272 
273  BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
274  num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
275  t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
276  t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf));
277 
278  for (size_t county = 0; county < model.size(); county++) {
279  // if (std::accumulate(
280  // num_InfectedSymptoms[county].begin(),
281  // num_InfectedSymptoms[county].end(),
282  // double(0.0),
283  // [](const double& a, const double& b) { return evaluate_intermediate<double>(a + b); }
284  // ) > 0)
285  // {
286  size_t num_groups = (size_t)model[county].parameters.get_num_groups();
287  for (size_t i = 0; i < num_groups; i++) {
288  model[county].populations[{AgeGroup(i), InfectionState::ExposedPartialImmunity}] = num_Exposed[county][i];
289  model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsPartialImmunity}] =
290  num_InfectedNoSymptoms[county][i];
291  model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] = 0;
292  model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunity}] =
293  num_InfectedSymptoms[county][i];
294  model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunityConfirmed}] = 0;
295  model[county].populations[{AgeGroup(i), InfectionState::InfectedSeverePartialImmunity}] =
296  num_InfectedSevere[county][i];
297  // Only set the number of ICU patients here, if the date is not available in the data.
298  if (!is_divi_data_available(date)) {
299  model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalPartialImmunity}] =
300  num_icu[county][i];
301  }
302  }
303  // }
304  if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), double(0.0),
305  [](const double& a, const double& b) {
306  return evaluate_intermediate<double>(a + b);
307  }) == 0.0) {
308  log_warning("No infections for partially vaccinated reported on date {} for region {}. "
309  "Population data has not been set.",
310  date, region[county]);
311  }
312  }
313 
314  /*----------- FULLY VACCINATED -----------*/
315  for (size_t county = 0; county < model.size(); county++) {
316  t_InfectedNoSymptoms[county].clear();
317  t_Exposed[county].clear();
318  t_InfectedSymptoms[county].clear();
319  t_InfectedSevere[county].clear();
320  t_InfectedCritical[county].clear();
321 
322  mu_C_R[county].clear();
323  mu_I_H[county].clear();
324  mu_H_U[county].clear();
325 
326  num_InfectedSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
327  num_death[county] = std::vector<double>(num_age_groups, 0.0);
328  num_rec[county] = std::vector<double>(num_age_groups, 0.0);
329  num_Exposed[county] = std::vector<double>(num_age_groups, 0.0);
330  num_InfectedNoSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
331  num_InfectedSevere[county] = std::vector<double>(num_age_groups, 0.0);
332  num_icu[county] = std::vector<double>(num_age_groups, 0.0);
333  for (size_t group = 0; group < num_age_groups; group++) {
334  double reduc_t = model[0].parameters.template get<ReducTimeInfectedMild<double>>()[(AgeGroup)group];
335  t_Exposed[county].push_back(static_cast<int>(
336  std::round(model[county].parameters.template get<TimeExposed<double>>()[(AgeGroup)group])));
337  t_InfectedNoSymptoms[county].push_back(static_cast<int>(std::round(
338  model[county].parameters.template get<TimeInfectedNoSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
339  t_InfectedSymptoms[county].push_back(static_cast<int>(std::round(
340  model[county].parameters.template get<TimeInfectedSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
341  t_InfectedSevere[county].push_back(static_cast<int>(
342  std::round(model[county].parameters.template get<TimeInfectedSevere<double>>()[(AgeGroup)group])));
343  t_InfectedCritical[county].push_back(static_cast<int>(
344  std::round(model[county].parameters.template get<TimeInfectedCritical<double>>()[(AgeGroup)group])));
345 
346  double reduc_immune_exp =
347  model[county].parameters.template get<ReducExposedImprovedImmunity<double>>()[(AgeGroup)group];
348  double reduc_immune_inf =
349  model[county].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[(AgeGroup)group];
350  double reduc_immune_hosp =
351  model[county].parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[(
352  AgeGroup)group];
353  double reduc_immune_icu =
354  model[county].parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[(
355  AgeGroup)group];
356  mu_C_R[county].push_back((
357  1 - reduc_immune_inf / reduc_immune_exp *
358  (1 - model[county]
359  .parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[(AgeGroup)group])));
360  mu_I_H[county].push_back(
361  reduc_immune_hosp / reduc_immune_inf *
362  model[county].parameters.template get<SeverePerInfectedSymptoms<double>>()[(AgeGroup)group]);
363  // transfer from H to U, D unchanged.
364  mu_H_U[county].push_back(
365  reduc_immune_icu / reduc_immune_hosp *
366  model[county].parameters.template get<CriticalPerSevere<double>>()[(AgeGroup)group]);
367  }
368  }
369 
370  BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
371  num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
372  t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
373  t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf));
374 
375  for (size_t county = 0; county < model.size(); county++) {
376  size_t num_groups = (size_t)model[county].parameters.get_num_groups();
377  for (size_t i = 0; i < num_groups; i++) {
378  model[county].populations[{AgeGroup(i), InfectionState::ExposedImprovedImmunity}] = num_Exposed[county][i];
379  model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsImprovedImmunity}] =
380  num_InfectedNoSymptoms[county][i];
381  model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] = 0;
382  model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunity}] =
383  num_InfectedSymptoms[county][i];
384  model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] = 0;
385  model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereImprovedImmunity}] =
386  num_InfectedSevere[county][i];
387  // Only set the number of ICU patients here, if the date is not available in the data.
388  if (!is_divi_data_available(date)) {
389  model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalImprovedImmunity}] =
390  num_icu[county][i];
391  }
392  }
393 
394  if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), double(0.0),
395  [](const double& a, const double& b) {
396  return evaluate_intermediate<double>(a + b);
397  }) == 0.0) {
398  log_warning("No infections for vaccinated reported on date {} for region {}. "
399  "Population data has not been set.",
400  date, region[county]);
401  }
402  }
403 
404  return success();
405 }
406 
417 template <class Model>
418 IOResult<void> set_confirmed_cases_data(std::vector<Model>& model, const std::string& path,
419  std::vector<int> const& region, Date date,
420  const std::vector<double>& scaling_factor_inf, bool set_death = false)
421 {
422  BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path));
423  BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, set_death));
424  return success();
425 }
426 
435 template <class Model>
436 IOResult<void> set_divi_data(std::vector<Model>& model, const std::string& path, const std::vector<int>& vregion,
437  Date date, double scaling_factor_icu)
438 {
439  // DIVI dataset will no longer be updated from CW29 2024 on.
440  if (!is_divi_data_available(date)) {
441  log_warning("No DIVI data available for date: {}. "
442  "ICU compartment will be set based on Case data.",
443  date);
444  return success();
445  }
446  std::vector<double> sum_mu_I_U(vregion.size(), 0);
447  std::vector<std::vector<double>> mu_I_U{model.size()};
448  for (size_t region = 0; region < vregion.size(); region++) {
449  auto num_groups = model[region].parameters.get_num_groups();
450  for (auto i = AgeGroup(0); i < num_groups; i++) {
451  sum_mu_I_U[region] += model[region].parameters.template get<CriticalPerSevere<double>>()[i] *
452  model[region].parameters.template get<SeverePerInfectedSymptoms<double>>()[i];
453  mu_I_U[region].push_back(model[region].parameters.template get<CriticalPerSevere<double>>()[i] *
454  model[region].parameters.template get<SeverePerInfectedSymptoms<double>>()[i]);
455  }
456  }
457  std::vector<double> num_icu(model.size(), 0.0);
458  BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu));
459 
460  for (size_t region = 0; region < vregion.size(); region++) {
461  auto num_groups = model[region].parameters.get_num_groups();
462  for (auto i = AgeGroup(0); i < num_groups; i++) {
463  model[region].populations[{i, InfectionState::InfectedCriticalNaive}] =
464  scaling_factor_icu * num_icu[region] * mu_I_U[region][(size_t)i] / sum_mu_I_U[region];
465  }
466  }
467 
468  return success();
469 }
470 
479 template <class Model>
480 IOResult<void> set_population_data(std::vector<Model>& model, const std::vector<std::vector<double>>& num_population,
481  const std::vector<ConfirmedCasesDataEntry>& case_data,
482  const std::vector<int>& vregion, Date date)
483 {
484  auto num_age_groups = ConfirmedCasesDataEntry::age_group_names.size();
485  std::vector<std::vector<double>> vnum_rec(model.size(), std::vector<double>(num_age_groups, 0.0));
486 
487  BOOST_OUTCOME_TRY(read_confirmed_cases_data_fix_recovered(case_data, vregion, date, vnum_rec, 14.));
488 
489  for (size_t region = 0; region < vregion.size(); region++) {
490 
491  if (std::accumulate(num_population[region].begin(), num_population[region].end(), double(0.0),
492  [](const double& a, const double& b) {
493  return evaluate_intermediate<double>(a + b);
494  }) > 0) {
495  auto num_groups = model[region].parameters.get_num_groups();
496  for (auto i = AgeGroup(0); i < num_groups; i++) {
497 
498  double S_v = std::min(
499  model[region].parameters.template get<DailyFullVaccinations<double>>()[{i, SimulationDay(0)}] +
500  vnum_rec[region][size_t(i)],
501  num_population[region][size_t(i)]);
502  double S_pv = std::max(
503  model[region].parameters.template get<DailyPartialVaccinations<double>>()[{i, SimulationDay(0)}] -
504  model[region].parameters.template get<DailyFullVaccinations<double>>()[{i, SimulationDay(0)}],
505  0.0); // use std::max with 0
506  double S;
507  if (num_population[region][size_t(i)] - S_pv - S_v < 0.0) {
508  log_warning("Number of vaccinated persons greater than population in county {}, age group {}.",
509  region, size_t(i));
510  S = 0.0;
511  S_v = num_population[region][size_t(i)] - S_pv;
512  }
513  else {
514  S = num_population[region][size_t(i)] - S_pv - S_v;
515  }
516 
517  double denom_E =
518  1 / (S + S_pv * model[region].parameters.template get<ReducExposedPartialImmunity<double>>()[i] +
519  S_v * model[region].parameters.template get<ReducExposedImprovedImmunity<double>>()[i]);
520  double denom_C =
521  1 / (S + S_pv * model[region].parameters.template get<ReducExposedPartialImmunity<double>>()[i] +
522  S_v * model[region].parameters.template get<ReducExposedImprovedImmunity<double>>()[i]);
523  double denom_I =
524  1 /
525  (S +
526  S_pv * model[region].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i] +
527  S_v * model[region].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i]);
528  double denom_HU =
529  1 /
530  (S +
531  S_pv * model[region]
532  .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] +
533  S_v * model[region]
534  .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i]);
535 
536  model[region].populations[{i, InfectionState::ExposedNaive}] =
537  S * model[region].populations[{i, InfectionState::ExposedNaive}] * denom_E;
538  model[region].populations[{i, InfectionState::ExposedPartialImmunity}] =
539  S_pv * model[region].parameters.template get<ReducExposedPartialImmunity<double>>()[i] *
540  model[region].populations[{i, InfectionState::ExposedPartialImmunity}] * denom_E;
541  model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] =
542  S_v * model[region].parameters.template get<ReducExposedImprovedImmunity<double>>()[i] *
543  model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] * denom_E;
544 
545  model[region].populations[{i, InfectionState::InfectedNoSymptomsNaive}] =
546  S * model[region].populations[{i, InfectionState::InfectedNoSymptomsNaive}] * denom_C;
547  model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] =
548  S_pv * model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] * denom_C;
549  model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] =
550  S_v * model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] * denom_C;
551 
552  model[region].populations[{i, InfectionState::InfectedNoSymptomsNaiveConfirmed}] =
553  S * model[region].populations[{i, InfectionState::InfectedNoSymptomsNaiveConfirmed}] * denom_C;
554  model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] =
555  S_pv * model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] *
556  denom_C;
557  model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] =
558  S_v * model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] *
559  denom_C;
560 
561  model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] =
562  S * model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] * denom_I;
563  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] =
564  S_pv * model[region].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i] *
565  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] * denom_I;
566  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] =
567  S_v * model[region].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i] *
568  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] * denom_I;
569 
570  model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] =
571  S * model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] * denom_I;
572  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] =
573  S_pv * model[region].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i] *
574  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] * denom_I;
575  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] =
576  S_v * model[region].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i] *
577  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] * denom_I;
578 
579  model[region].populations[{i, InfectionState::InfectedSevereNaive}] =
580  S * model[region].populations[{i, InfectionState::InfectedSevereNaive}] * denom_HU;
581  model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] =
582  S_pv *
583  model[region].parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] *
584  model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] * denom_HU;
585  model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] =
586  S_v *
587  model[region]
588  .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i] *
589  model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] * denom_HU;
590 
591  model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] =
592  S_pv *
593  model[region].parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] *
594  model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU;
595  model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] =
596  S_v *
597  model[region]
598  .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i] *
599  model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU;
600  model[region].populations[{i, InfectionState::InfectedCriticalNaive}] =
601  S * model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU;
602 
603  model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] =
604  model[region].parameters.template get<DailyFullVaccinations<double>>()[{i, SimulationDay(0)}] +
605  model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] -
606  (model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] +
607  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] +
608  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] +
609  model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] +
610  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] +
611  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] +
612  model[region].populations[{i, InfectionState::InfectedSevereNaive}] +
613  model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] +
614  model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] +
615  model[region].populations[{i, InfectionState::InfectedCriticalNaive}] +
616  model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] +
617  model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] +
618  model[region].populations[{i, InfectionState::DeadNaive}] +
619  model[region].populations[{i, InfectionState::DeadPartialImmunity}] +
620  model[region].populations[{i, InfectionState::DeadImprovedImmunity}]);
621 
622  model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = std::min(
623  S_v - model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] -
624  model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] -
625  model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] -
626  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] -
627  model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] -
628  model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] -
629  model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] -
630  model[region].populations[{i, InfectionState::DeadImprovedImmunity}],
631  std::max(0.0, double(model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}])));
632 
633  model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = std::max(
634  0.0,
635  S_pv - model[region].populations[{i, InfectionState::ExposedPartialImmunity}] -
636  model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] -
637  model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] -
638  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] -
639  model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] -
640  model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] -
641  model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] -
642  model[region].populations[{i, InfectionState::DeadPartialImmunity}]);
643 
644  model[region].populations.template set_difference_from_group_total<AgeGroup>(
645  {i, InfectionState::SusceptibleNaive}, num_population[region][size_t(i)]);
646  }
647 
648  for (auto i = AgeGroup(0); i < AgeGroup(6); i++) {
649  for (auto j = Index<InfectionState>(0); j < InfectionState::Count; ++j) {
650  if (model[region].populations[{i, j}] < 0) {
651  log_warning("Compartment at age group {}, infection state {}, is negative: {}", size_t(i),
652  size_t(j), model[region].populations[{i, j}] / num_population[region][size_t(i)]);
653  }
654  }
655  }
656  }
657  else {
658  log_warning("No population data available for region " + std::to_string(region) +
659  ". Population data has not been set.");
660  }
661  }
662 
663  return success();
664 }
665 
674 template <class Model>
675 IOResult<void> set_population_data(std::vector<Model>& model, const std::string& path, const std::string& path_rki,
676  const std::vector<int>& vregion, Date date)
677 {
678  BOOST_OUTCOME_TRY(auto&& num_population, read_population_data(path, vregion));
679  BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path_rki));
680 
681  BOOST_OUTCOME_TRY(set_population_data(model, num_population, rki_data, vregion, date));
682  return success();
683 }
684 
695 template <typename FP = double>
696 IOResult<void> set_vaccination_data(std::vector<Model<FP>>& model, const std::vector<VaccinationDataEntry>& vacc_data,
697  Date date, const std::vector<int>& vregion, int num_days)
698 {
699  auto num_groups = model[0].parameters.get_num_groups();
700 
701  // type conversion from UncertainValue -> FP -> int
702  auto days_until_effective1 = static_cast<int>(
703  static_cast<FP>(model[0].parameters.template get<DaysUntilEffectivePartialImmunity<FP>>()[AgeGroup(0)]));
704  auto days_until_effective2 = static_cast<int>(
705  static_cast<FP>(model[0].parameters.template get<DaysUntilEffectiveImprovedImmunity<FP>>()[AgeGroup(0)]));
706  auto vaccination_distance =
707  static_cast<int>(static_cast<FP>(model[0].parameters.template get<VaccinationGap<FP>>()[AgeGroup(0)]));
708 
709  // iterate over regions (e.g., counties)
710  for (size_t i = 0; i < model.size(); ++i) {
711  // iterate over age groups in region
712  for (auto g = AgeGroup(0); g < num_groups; ++g) {
713 
714  model[i].parameters.template get<DailyPartialVaccinations<FP>>().resize(SimulationDay(num_days + 1));
715  model[i].parameters.template get<DailyFullVaccinations<FP>>().resize(SimulationDay(num_days + 1));
716  for (auto d = SimulationDay(0); d < SimulationDay(num_days + 1); ++d) {
717  model[i].parameters.template get<DailyPartialVaccinations<FP>>()[{g, d}] = 0.0;
718  model[i].parameters.template get<DailyFullVaccinations<FP>>()[{g, d}] = 0.0;
719  }
720  }
721  }
722 
723  auto max_date_entry = std::max_element(vacc_data.begin(), vacc_data.end(), [](auto&& a, auto&& b) {
724  return a.date < b.date;
725  });
726  if (max_date_entry == vacc_data.end()) {
727  return failure(StatusCode::InvalidFileFormat, "Vaccination data file is empty.");
728  }
729  auto max_date = max_date_entry->date;
730 
731  auto min_date_entry = std::min_element(vacc_data.begin(), vacc_data.end(), [](auto&& a, auto&& b) {
732  return a.date < b.date;
733  });
734  auto min_date = min_date_entry->date;
735  if (min_date > date || max_date < offset_date_by_days(date, num_days)) {
736  log_warning("Vaccination data only available from {} to {}. "
737  "For days before, vaccination data will be set to 0. For days after, "
738  "vaccination data will be set to the last available date.",
739  min_date, max_date);
740  }
741 
742  for (auto&& vacc_data_entry : vacc_data) {
743  auto it = std::find_if(vregion.begin(), vregion.end(), [&vacc_data_entry](auto&& r) {
744  return r == 0 || (vacc_data_entry.county_id && vacc_data_entry.county_id == regions::CountyId(r)) ||
745  (vacc_data_entry.state_id && vacc_data_entry.state_id == regions::StateId(r)) ||
746  (vacc_data_entry.district_id && vacc_data_entry.district_id == regions::DistrictId(r));
747  });
748  auto date_df = vacc_data_entry.date;
749  if (it != vregion.end()) {
750  auto region_idx = size_t(it - vregion.begin());
751  auto age = vacc_data_entry.age_group;
752 
753  for (size_t d = 0; d < (size_t)num_days + 1; ++d) {
754  int days_plus;
755  // In the following, second dose means previous 'full immunization', now 'Grundimmunisierung'.
756  // ---
757  // date: start_date of the simulation (Input from IO call read_input_data_county_vaccmodel())
758  // d: day of simulation, counted from 0 to num_days (for which we need (approximated) vaccination numbers)
759  // root[i]["Vacc_completed"]: accumulated number of total second doses up to day date_df;
760  // taken from input dataframe, single value, per county and age group
761  // ----
762  // An averaged distance between first and second doses (vaccination_distance) is assumed in the following
763  // and the first doses are computed based on the second doses given 'vaccination_distance' days later.
764  // ----
765  // a person whose second dose is reported at start_date + simulation_day - days_until_effective1 + vaccination_distance
766  // had the first dose on start_date + simulation_day - days_until_effective1. Furthermore, he/she has the full protection
767  // of the first dose at day X = start_date + simulation_day
768  // Storing its value in get<DailyPartialVaccinations>() will eventually (in the simulation)
769  // transfer the difference (between get<DailyPartialVaccinations>() at d and d-1) of
770  // N susceptible individuals to 'Susceptible Partially Vaccinated' state at day d; see secir_vaccinated.h
771  auto offset_first_date =
772  offset_date_by_days(date, (int)d - days_until_effective1 + vaccination_distance);
773  if (max_date >= offset_first_date) {
774  // Option 1: considered offset_first_date is available in input data frame
775  if (date_df == offset_first_date) {
776  model[region_idx]
777  .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] =
778  vacc_data_entry.num_vaccinations_completed;
779  }
780  }
781  else { // offset_first_date > max_date
782  // Option 2: considered offset_first_date is NOT available in input data frame
783  // Here, a constant number of first and second doses is assumed, i.e.,
784  // the the number of vaccinationes at day d (N days after max_date) will be:
785  // total number of vaccinations up to day max_date + N * number of vaccinations ON max_date
786  // (where the latter is computed as the difference between the total number at max_date and max_date-1)
787  days_plus = get_offset_in_days(offset_first_date, max_date);
788  if (date_df == offset_date_by_days(max_date, -1)) {
789  model[region_idx]
790  .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] -=
791  days_plus * vacc_data_entry.num_vaccinations_completed;
792  }
793  else if (date_df == max_date) {
794  model[region_idx]
795  .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] +=
796  (days_plus + 1) * vacc_data_entry.num_vaccinations_completed;
797  }
798  }
799 
800  // a person whose second dose is reported at start_date + simulation_day - days_until_effective2
801  // has the full protection of the second dose at day X = start_date + simulation_day
802  // Storing its value in get<DailyFullVaccinations>() will eventually (in the simulation)
803  // transfer the difference (between get<DailyFullVaccinations>() at d and d-1) of
804  // N susceptible, partially vaccinated individuals to 'SusceptibleImprovedImmunity' state at day d; see secir_vaccinated.h
805  auto offset_full_date = offset_date_by_days(date, (int)d - days_until_effective2);
806  if (max_date >= offset_full_date) {
807  // Option 1: considered offset_full_date is available in input data frame
808  if (date_df == offset_full_date) {
809  model[region_idx]
810  .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] =
811  vacc_data_entry.num_vaccinations_completed;
812  }
813  }
814  else { // offset_full_date > max_full_date
815  // Option 2: considered offset_full_date is NOT available in input data frame
816  days_plus = get_offset_in_days(offset_full_date, max_date);
817  if (date_df == offset_date_by_days(max_date, -1)) {
818  model[region_idx]
819  .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] -=
820  days_plus * vacc_data_entry.num_vaccinations_completed;
821  }
822  else if (date_df == max_date) {
823  model[region_idx]
824  .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] +=
825  (days_plus + 1) * vacc_data_entry.num_vaccinations_completed;
826  }
827  }
828  }
829  }
830  }
831  return success();
832 }
833 
844 template <typename FP = double>
845 IOResult<void> set_vaccination_data(std::vector<Model<FP>>& model, const std::string& path, Date date,
846  const std::vector<int>& vregion, int num_days)
847 {
848  // Check if vaccination data is available for the given date range
849  auto end_date = offset_date_by_days(date, num_days);
850  if (!is_vaccination_data_available(date, end_date)) {
851  log_warning("No vaccination data available in range from {} to {}. "
852  "Vaccination data will be set to 0.",
853  date, end_date);
854  // Set vaccination data to 0 for all models
855  for (auto& m : model) {
856  m.parameters.template get<DailyPartialVaccinations<FP>>().resize(SimulationDay(num_days + 1));
857  m.parameters.template get<DailyFullVaccinations<FP>>().resize(SimulationDay(num_days + 1));
858 
859  for (auto d = SimulationDay(0); d < SimulationDay(num_days + 1); ++d) {
860  for (auto a = AgeGroup(0); a < m.parameters.get_num_groups(); ++a) {
861  m.parameters.template get<DailyPartialVaccinations<FP>>()[{a, d}] = 0.0;
862  m.parameters.template get<DailyFullVaccinations<FP>>()[{a, d}] = 0.0;
863  }
864  }
865  }
866  return success();
867  }
868  BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path));
869  BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days));
870  return success();
871 }
872 
873 } // namespace details
874 
875 #ifdef MEMILIO_HAS_HDF5
876 
894 template <class Model>
895 IOResult<void> export_input_data_county_timeseries(
896  std::vector<Model> models, const std::string& results_dir, const std::vector<int>& counties, Date date,
897  const std::vector<double>& scaling_factor_inf, const double scaling_factor_icu, const int num_days,
898  const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path,
899  const std::string& vaccination_data_path = "")
900 {
901  const auto num_groups = (size_t)models[0].parameters.get_num_groups();
902  assert(scaling_factor_inf.size() == num_groups);
903  assert(num_groups == ConfirmedCasesDataEntry::age_group_names.size());
904  assert(models.size() == counties.size());
905  std::vector<TimeSeries<double>> extrapolated_data(
906  models.size(), TimeSeries<double>::zero(num_days + 1, (size_t)InfectionState::Count * num_groups));
907 
908  BOOST_OUTCOME_TRY(auto&& case_data, read_confirmed_cases_data(confirmed_cases_path));
909  BOOST_OUTCOME_TRY(auto&& population_data, read_population_data(population_data_path, counties));
910 
911  // empty vector if set_vaccination_data is not set
912  std::vector<VaccinationDataEntry> vacc_data;
913  if (!vaccination_data_path.empty()) {
914  BOOST_OUTCOME_TRY(vacc_data, read_vaccination_data(vaccination_data_path));
915  }
916 
917  for (int t = 0; t <= num_days; ++t) {
918  auto offset_day = offset_date_by_days(date, t);
919 
920  if (!vaccination_data_path.empty()) {
921  BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days));
922  }
923 
924  // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
925  // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
926  BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu));
927 
928  BOOST_OUTCOME_TRY(
929  details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, true));
930  BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, case_data, counties, offset_day));
931 
932  for (size_t r = 0; r < counties.size(); r++) {
933  extrapolated_data[r][t] = models[r].get_initial_values();
934  // in set_population_data the number of death individuals is subtracted from the SusceptibleImprovedImmunity compartment.
935  // Since we should be independent whether we consider them or not, we add them back here before we save the data.
936  for (size_t age = 0; age < num_groups; age++) {
937  extrapolated_data[r][t][(size_t)InfectionState::SusceptibleImprovedImmunity +
938  age * (size_t)InfectionState::Count] +=
939  extrapolated_data[r][t][(size_t)InfectionState::DeadNaive + age * (size_t)InfectionState::Count];
940  }
941  }
942  }
943  BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast<int>(num_groups),
944  path_join(results_dir, "Results_rki.h5")));
945 
946  auto extrapolated_rki_data_sum = sum_nodes(std::vector<std::vector<TimeSeries<double>>>{extrapolated_data});
947  BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast<int>(num_groups),
948  path_join(results_dir, "Results_rki_sum.h5")));
949 
950  return success();
951 }
952 #else
953 template <class Model>
954 IOResult<void> export_input_data_county_timeseries(std::vector<Model>, const std::string&, const std::vector<int>&,
955  Date, const std::vector<double>&, const double, const int,
956  const std::string&, const std::string&, const std::string&,
957  const std::string&)
958 {
959  mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
960  return success();
961 }
962 
963 #endif //MEMILIO_HAS_HDF5
964 
979 template <class Model>
980 IOResult<void> read_input_data_county(std::vector<Model>& model, Date date, const std::vector<int>& county,
981  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
982  const std::string& pydata_dir, int num_days, bool export_time_series = false)
983 {
984  BOOST_OUTCOME_TRY(details::set_vaccination_data(model, path_join(pydata_dir, "vacc_county_ageinf_ma7.json"), date,
985  county, num_days));
986 
987  // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
988  // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
989  BOOST_OUTCOME_TRY(
990  details::set_divi_data(model, path_join(pydata_dir, "county_divi_ma7.json"), county, date, scaling_factor_icu));
991 
992  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "cases_all_county_age_ma7.json"),
993  county, date, scaling_factor_inf));
994  BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "county_current_population.json"),
995  path_join(pydata_dir, "cases_all_county_age_ma7.json"), county,
996  date));
997 
998  if (export_time_series) {
999  // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
1000  // Run time equals run time of the previous functions times the num_days !
1001  // (This only represents the vectorization of the previous function over all simulation days...)
1002  log_warning("Exporting time series of extrapolated real data. This may take some minutes. "
1003  "For simulation runs over the same time period, deactivate it.");
1004  BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
1005  model, pydata_dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days,
1006  path_join(pydata_dir, "county_divi_ma7.json"), path_join(pydata_dir, "cases_all_county_age_ma7.json"),
1007  path_join(pydata_dir, "county_current_population.json"),
1008  path_join(pydata_dir, "vacc_county_ageinf_ma7.json")));
1009  }
1010 
1011  return success();
1012 }
1013 
1028 template <class Model>
1029 IOResult<void> read_input_data(std::vector<Model>& model, Date date, const std::vector<int>& node_ids,
1030  const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
1031  const std::string& pydata_dir, int num_days, bool export_time_series = false)
1032 {
1033 
1034  BOOST_OUTCOME_TRY(
1035  details::set_vaccination_data(model, path_join(pydata_dir, "vaccination_data.json"), date, node_ids, num_days));
1036 
1037  // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
1038  // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
1039  BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(pydata_dir, "critical_cases.json"), node_ids, date,
1040  scaling_factor_icu));
1041 
1042  BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(pydata_dir, "confirmed_cases.json"), node_ids,
1043  date, scaling_factor_inf));
1044  BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(pydata_dir, "population_data.json"),
1045  path_join(pydata_dir, "confirmed_cases.json"), node_ids, date));
1046 
1047  if (export_time_series) {
1048  // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
1049  // Run time equals run time of the previous functions times the num_days !
1050  // (This only represents the vectorization of the previous function over all simulation days...)
1051  log_warning("Exporting time series of extrapolated real data. This may take some minutes. "
1052  "For simulation runs over the same time period, deactivate it.");
1053  BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
1054  model, pydata_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days,
1055  path_join(pydata_dir, "divi_data.json"), path_join(pydata_dir, "confirmed_cases.json"),
1056  path_join(pydata_dir, "population_data.json"), path_join(pydata_dir, "vaccination_data.json")));
1057  }
1058 
1059  return success();
1060 }
1061 
1071 template <class FP>
1072 IOResult<void> convert_model_data_type(mio::VectorRange<Node<Model<ScalarType>>> model_from,
1073  mio::VectorRange<Node<Model<FP>>> model_to)
1074 {
1075  assert(model_from.size() == model_to.size());
1076  assert((size_t)model_from[0].property.parameters.get_num_groups() ==
1077  (size_t)model_to[0].property.parameters.get_num_groups());
1078  // Todo: Add conversion of ParameterSet and then re-use code from all model parameters io
1079  // OR call set_vacination_data with FP to set correct parameters
1080 
1081  for (size_t region_idx = 0; region_idx < model_from.size(); ++region_idx) {
1082  // convert populations to mio::UncertainValue<FP>
1083  // needs 2 converts as mio::UncertainValue<ScalarType> -> mio::UncertainValue<FP> does not work
1084  model_to[region_idx].property.populations = model_from[region_idx].property.populations.template convert<FP>();
1085  }
1086  return mio::success();
1087 }
1088 
1089 } // namespace osecirvvs
1090 } // namespace mio
1091 
1092 #endif // MEMILIO_HAS_JSONCPP
1093 
1094 #endif // MIO_ODE_SECIRVVS_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
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::mpi::size
int size(Comm comm)
Return the size of the given communicator.
Definition: miompi.cpp:75
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::failure
auto failure(const IOStatus &s)
Create an object that is implicitly convertible to an error IOResult<T>.
Definition: io.h:380
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::VectorRange
std::conditional_t< std::is_const_v< T >, typename mio::Range< typename std::vector< std::remove_const_t< T > >::const_iterator >, typename mio::Range< typename std::vector< std::remove_const_t< T > >::iterator > > VectorRange
Defines generic Range type for IterPair of a vector.
Definition: stl_util.h:341
mio::success
auto success()
Create an object that is implicitly convertible to a succesful IOResult<void>.
Definition: io.h:359
mio::get_offset_in_days
int get_offset_in_days(Date date1, Date date2)
Computes the offset in days given two dates: first date minus second date.
Definition: date.h:310
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