Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 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/io/epi_data.h"
29 : #include "memilio/io/io.h"
30 : #include "memilio/io/result_io.h"
31 : #include "memilio/utils/date.h"
32 : #include "memilio/utils/logging.h"
33 :
34 : namespace mio
35 : {
36 : namespace osecirvvs
37 : {
38 :
39 : namespace details
40 : {
41 : /**
42 : * @brief Reads subpopulations of infection states from transformed RKI cases file.
43 : * @param[in] path Path to transformed RKI cases file.
44 : * @param[in] vregion Vector of keys of the region of interest.
45 : * @param[in] date Date for which the arrays are initialized.
46 : * @param[in, out] num_* Output vector for number of people in the corresponding compartement.
47 : * @param[in] vt_* Average time it takes to get from one compartement to another (vector with one element per age group).
48 : * @param[in] vmu_* Probabilities to get from one compartement to another (vector with one element per age group).
49 : * @param[in] scaling_factor_inf Factor for scaling the confirmed cases to account for estimated undetected cases.
50 : * @see mio::read_confirmed_cases_data
51 : * @{
52 : */
53 : IOResult<void> read_confirmed_cases_data(
54 : std::string const& path, std::vector<int> const& vregion, Date date, std::vector<std::vector<double>>& num_Exposed,
55 : std::vector<std::vector<double>>& num_InfectedNoSymptoms, std::vector<std::vector<double>>& num_InfectedSymptoms,
56 : std::vector<std::vector<double>>& num_InfectedSevere, std::vector<std::vector<double>>& num_icu,
57 : std::vector<std::vector<double>>& num_death, std::vector<std::vector<double>>& num_rec,
58 : const std::vector<std::vector<int>>& t_Exposed, const std::vector<std::vector<int>>& t_InfectedNoSymptoms,
59 : const std::vector<std::vector<int>>& t_InfectedSymptoms, const std::vector<std::vector<int>>& t_InfectedSevere,
60 : const std::vector<std::vector<int>>& t_InfectedCritical, const std::vector<std::vector<double>>& mu_C_R,
61 : const std::vector<std::vector<double>>& mu_I_H, const std::vector<std::vector<double>>& mu_H_U,
62 : const std::vector<double>& scaling_factor_inf);
63 :
64 : IOResult<void> read_confirmed_cases_data(
65 : const std::vector<ConfirmedCasesDataEntry>& rki_data, std::vector<int> const& vregion, Date date,
66 : std::vector<std::vector<double>>& num_Exposed, std::vector<std::vector<double>>& num_InfectedNoSymptoms,
67 : std::vector<std::vector<double>>& num_InfectedSymptoms, std::vector<std::vector<double>>& num_InfectedSevere,
68 : std::vector<std::vector<double>>& num_icu, std::vector<std::vector<double>>& num_death,
69 : std::vector<std::vector<double>>& num_rec, const std::vector<std::vector<int>>& t_Exposed,
70 : const std::vector<std::vector<int>>& t_InfectedNoSymptoms, const std::vector<std::vector<int>>& t_InfectedSymptoms,
71 : const std::vector<std::vector<int>>& t_InfectedSevere, const std::vector<std::vector<int>>& t_InfectedCritical,
72 : const std::vector<std::vector<double>>& mu_C_R, const std::vector<std::vector<double>>& mu_I_H,
73 : const std::vector<std::vector<double>>& mu_H_U, const std::vector<double>& scaling_factor_inf);
74 : /**@}*/
75 :
76 : /**
77 : * @brief Reads confirmed cases data and translates data of day t0-delay to recovered compartment.
78 : * @param[in] path Path to RKI confirmed cases file.
79 : * @param[in] vregion Vector of keys of the region of interest.
80 : * @param[in] date Date for which the arrays are initialized.
81 : * @param[in, out] num_rec Output vector for number of people in the compartement recovered.
82 : * @param[in] delay Number of days in the past the are used to set recovered compartment.
83 : * @see mio::read_confirmed_cases_data
84 : * @{
85 : */
86 : IOResult<void> read_confirmed_cases_data_fix_recovered(const std::vector<ConfirmedCasesDataEntry>& rki_data,
87 : std::vector<int> const& vregion, Date date,
88 : std::vector<std::vector<double>>& vnum_rec, double delay = 14.);
89 : IOResult<void> read_confirmed_cases_data_fix_recovered(std::string const& path, std::vector<int> const& vregion,
90 : Date date, std::vector<std::vector<double>>& vnum_rec,
91 : double delay = 14.);
92 : /**@}*/
93 :
94 : /**
95 : * @brief Sets the confirmed cases data for a vector of models based on input data.
96 : * @param[in, out] model Vector of objects in which the data is set.
97 : * @param[in] case_data Vector of case data. Each inner vector represents a different region.
98 : * @param[in] region Vector of keys of the region of interest.
99 : * @param[in] date Date for which the arrays are initialized.
100 : * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of RKI data.
101 : * @param[in] set_death If true, set the number of deaths.
102 : */
103 : template <class Model>
104 99 : IOResult<void> set_confirmed_cases_data(std::vector<Model>& model,
105 : const std::vector<ConfirmedCasesDataEntry>& case_data,
106 : std::vector<int> const& region, Date date,
107 : const std::vector<double>& scaling_factor_inf, bool set_death = false)
108 : {
109 99 : auto num_age_groups = (size_t)model[0].parameters.get_num_groups();
110 99 : assert(scaling_factor_inf.size() == num_age_groups); //TODO: allow vector or scalar valued scaling factors
111 99 : assert(ConfirmedCasesDataEntry::age_group_names.size() == num_age_groups);
112 :
113 198 : std::vector<std::vector<int>> t_Exposed{model.size()};
114 198 : std::vector<std::vector<int>> t_InfectedNoSymptoms{model.size()};
115 198 : std::vector<std::vector<int>> t_InfectedSymptoms{model.size()};
116 198 : std::vector<std::vector<int>> t_InfectedSevere{model.size()};
117 198 : std::vector<std::vector<int>> t_InfectedCritical{model.size()};
118 :
119 198 : std::vector<std::vector<double>> mu_C_R{model.size()};
120 198 : std::vector<std::vector<double>> mu_I_H{model.size()};
121 198 : std::vector<std::vector<double>> mu_H_U{model.size()};
122 :
123 198 : std::vector<std::vector<double>> num_InfectedSymptoms(model.size());
124 198 : std::vector<std::vector<double>> num_death(model.size());
125 198 : std::vector<std::vector<double>> num_rec(model.size());
126 198 : std::vector<std::vector<double>> num_Exposed(model.size());
127 198 : std::vector<std::vector<double>> num_InfectedNoSymptoms(model.size());
128 198 : std::vector<std::vector<double>> num_InfectedSevere(model.size());
129 198 : std::vector<std::vector<double>> num_icu(model.size());
130 :
131 : /*----------- UNVACCINATED -----------*/
132 198 : for (size_t county = 0; county < model.size(); county++) {
133 99 : num_InfectedSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
134 99 : num_death[county] = std::vector<double>(num_age_groups, 0.0);
135 99 : num_rec[county] = std::vector<double>(num_age_groups, 0.0);
136 99 : num_Exposed[county] = std::vector<double>(num_age_groups, 0.0);
137 99 : num_InfectedNoSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
138 99 : num_InfectedSevere[county] = std::vector<double>(num_age_groups, 0.0);
139 99 : num_icu[county] = std::vector<double>(num_age_groups, 0.0);
140 693 : for (size_t group = 0; group < num_age_groups; group++) {
141 :
142 1188 : t_Exposed[county].push_back(static_cast<int>(
143 1188 : std::round(model[county].parameters.template get<TimeExposed<double>>()[(AgeGroup)group])));
144 1188 : t_InfectedNoSymptoms[county].push_back(static_cast<int>(
145 1188 : std::round(model[county].parameters.template get<TimeInfectedNoSymptoms<double>>()[(AgeGroup)group])));
146 1188 : t_InfectedSymptoms[county].push_back(static_cast<int>(
147 1188 : std::round(model[county].parameters.template get<TimeInfectedSymptoms<double>>()[(AgeGroup)group])));
148 1188 : t_InfectedSevere[county].push_back(static_cast<int>(
149 1188 : std::round(model[county].parameters.template get<TimeInfectedSevere<double>>()[(AgeGroup)group])));
150 1188 : t_InfectedCritical[county].push_back(static_cast<int>(
151 1188 : std::round(model[county].parameters.template get<TimeInfectedCritical<double>>()[(AgeGroup)group])));
152 :
153 1188 : mu_C_R[county].push_back(
154 1188 : model[county].parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[(AgeGroup)group]);
155 1188 : mu_I_H[county].push_back(
156 1188 : model[county].parameters.template get<SeverePerInfectedSymptoms<double>>()[(AgeGroup)group]);
157 1188 : mu_H_U[county].push_back(
158 1188 : model[county].parameters.template get<CriticalPerSevere<double>>()[(AgeGroup)group]);
159 : }
160 : }
161 :
162 99 : BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
163 : num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
164 : t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
165 : t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf));
166 :
167 198 : for (size_t county = 0; county < model.size(); county++) {
168 : // if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) > 0) {
169 99 : size_t num_groups = (size_t)model[county].parameters.get_num_groups();
170 693 : for (size_t i = 0; i < num_groups; i++) {
171 594 : model[county].populations[{AgeGroup(i), InfectionState::ExposedNaive}] = num_Exposed[county][i];
172 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsNaive}] =
173 : num_InfectedNoSymptoms[county][i];
174 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsNaiveConfirmed}] = 0;
175 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaive}] =
176 : num_InfectedSymptoms[county][i];
177 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaiveConfirmed}] = 0;
178 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereNaive}] =
179 : num_InfectedSevere[county][i];
180 : // Only set the number of ICU patients here, if the date is not available in the data.
181 594 : if (!is_divi_data_available(date)) {
182 216 : model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalNaive}] = num_icu[county][i];
183 : }
184 594 : model[county].populations[{AgeGroup(i), InfectionState::SusceptibleImprovedImmunity}] = num_rec[county][i];
185 594 : if (set_death) {
186 : // in set_confirmed_cases_data initilization, deaths are now set to 0. In order to visualize
187 : // the extrapolated real number of deaths, they have to be set here. In the comparison of data
188 : // it has to be paid attention to the fact, the the simulation starts with deaths=0
189 : // while this method starts with deaths=number of reported deaths so far...
190 : // Additionally, we set the number of reported deaths to DeadNaive since no information on that is
191 : // available here.
192 : // Do only add deaths after substraction.
193 216 : model[county].populations[{AgeGroup(i), InfectionState::DeadNaive}] = num_death[county][i];
194 : }
195 : }
196 :
197 : // }
198 99 : if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) {
199 27 : log_warning(
200 : "No infections for unvaccinated reported on date {} for region {}. Population data has not been set.",
201 : date, region[county]);
202 : }
203 : }
204 :
205 : /*----------- PARTIALLY VACCINATED -----------*/
206 198 : for (size_t county = 0; county < model.size(); county++) {
207 99 : t_InfectedNoSymptoms[county].clear();
208 99 : t_Exposed[county].clear();
209 99 : t_InfectedSymptoms[county].clear();
210 99 : t_InfectedSevere[county].clear();
211 99 : t_InfectedCritical[county].clear();
212 :
213 99 : mu_C_R[county].clear();
214 99 : mu_I_H[county].clear();
215 99 : mu_H_U[county].clear();
216 :
217 99 : num_InfectedSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
218 99 : num_death[county] = std::vector<double>(num_age_groups, 0.0);
219 99 : num_rec[county] = std::vector<double>(num_age_groups, 0.0);
220 99 : num_Exposed[county] = std::vector<double>(num_age_groups, 0.0);
221 99 : num_InfectedNoSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
222 99 : num_InfectedSevere[county] = std::vector<double>(num_age_groups, 0.0);
223 99 : num_icu[county] = std::vector<double>(num_age_groups, 0.0);
224 693 : for (size_t group = 0; group < num_age_groups; group++) {
225 594 : double reduc_t = model[0].parameters.template get<ReducTimeInfectedMild<double>>()[(AgeGroup)group];
226 1188 : t_Exposed[county].push_back(static_cast<int>(
227 1188 : std::round(model[county].parameters.template get<TimeExposed<double>>()[(AgeGroup)group])));
228 1188 : t_InfectedNoSymptoms[county].push_back(static_cast<int>(std::round(
229 1188 : model[county].parameters.template get<TimeInfectedNoSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
230 1188 : t_InfectedSymptoms[county].push_back(static_cast<int>(std::round(
231 1188 : model[county].parameters.template get<TimeInfectedSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
232 1188 : t_InfectedSevere[county].push_back(static_cast<int>(
233 1188 : std::round(model[county].parameters.template get<TimeInfectedSevere<double>>()[(AgeGroup)group])));
234 1188 : t_InfectedCritical[county].push_back(static_cast<int>(
235 1188 : std::round(model[county].parameters.template get<TimeInfectedCritical<double>>()[(AgeGroup)group])));
236 :
237 594 : double exp_fac_part_immune =
238 1188 : model[county].parameters.template get<ReducExposedPartialImmunity<double>>()[(AgeGroup)group];
239 594 : double inf_fac_part_immune =
240 1188 : model[county].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[(AgeGroup)group];
241 594 : double hosp_fac_part_immune =
242 594 : model[county]
243 1188 : .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[(AgeGroup)group];
244 594 : double icu_fac_part_immune =
245 594 : model[county]
246 1188 : .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[(AgeGroup)group];
247 1188 : mu_C_R[county].push_back((
248 1188 : 1 - inf_fac_part_immune / exp_fac_part_immune *
249 1188 : (1 - model[county]
250 1188 : .parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[(AgeGroup)group])));
251 1188 : mu_I_H[county].push_back(
252 1782 : hosp_fac_part_immune / inf_fac_part_immune *
253 1188 : model[county].parameters.template get<SeverePerInfectedSymptoms<double>>()[(AgeGroup)group]);
254 : // transfer from H to U, D unchanged.
255 1188 : mu_H_U[county].push_back(
256 1782 : icu_fac_part_immune / hosp_fac_part_immune *
257 1188 : model[county].parameters.template get<CriticalPerSevere<double>>()[(AgeGroup)group]);
258 : }
259 : }
260 :
261 99 : BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
262 : num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
263 : t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
264 : t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf));
265 :
266 198 : for (size_t county = 0; county < model.size(); county++) {
267 : // if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) > 0) {
268 99 : size_t num_groups = (size_t)model[county].parameters.get_num_groups();
269 693 : for (size_t i = 0; i < num_groups; i++) {
270 594 : model[county].populations[{AgeGroup(i), InfectionState::ExposedPartialImmunity}] = num_Exposed[county][i];
271 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsPartialImmunity}] =
272 : num_InfectedNoSymptoms[county][i];
273 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] = 0;
274 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunity}] =
275 : num_InfectedSymptoms[county][i];
276 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunityConfirmed}] = 0;
277 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSeverePartialImmunity}] =
278 : num_InfectedSevere[county][i];
279 : // Only set the number of ICU patients here, if the date is not available in the data.
280 594 : if (!is_divi_data_available(date)) {
281 216 : model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalPartialImmunity}] =
282 : num_icu[county][i];
283 : }
284 : }
285 : // }
286 99 : if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) {
287 27 : log_warning("No infections for partially vaccinated reported on date {} for region {}. "
288 : "Population data has not been set.",
289 : date, region[county]);
290 : }
291 : }
292 :
293 : /*----------- FULLY VACCINATED -----------*/
294 198 : for (size_t county = 0; county < model.size(); county++) {
295 99 : t_InfectedNoSymptoms[county].clear();
296 99 : t_Exposed[county].clear();
297 99 : t_InfectedSymptoms[county].clear();
298 99 : t_InfectedSevere[county].clear();
299 99 : t_InfectedCritical[county].clear();
300 :
301 99 : mu_C_R[county].clear();
302 99 : mu_I_H[county].clear();
303 99 : mu_H_U[county].clear();
304 :
305 99 : num_InfectedSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
306 99 : num_death[county] = std::vector<double>(num_age_groups, 0.0);
307 99 : num_rec[county] = std::vector<double>(num_age_groups, 0.0);
308 99 : num_Exposed[county] = std::vector<double>(num_age_groups, 0.0);
309 99 : num_InfectedNoSymptoms[county] = std::vector<double>(num_age_groups, 0.0);
310 99 : num_InfectedSevere[county] = std::vector<double>(num_age_groups, 0.0);
311 99 : num_icu[county] = std::vector<double>(num_age_groups, 0.0);
312 693 : for (size_t group = 0; group < num_age_groups; group++) {
313 594 : double reduc_t = model[0].parameters.template get<ReducTimeInfectedMild<double>>()[(AgeGroup)group];
314 1188 : t_Exposed[county].push_back(static_cast<int>(
315 1188 : std::round(model[county].parameters.template get<TimeExposed<double>>()[(AgeGroup)group])));
316 1188 : t_InfectedNoSymptoms[county].push_back(static_cast<int>(std::round(
317 1188 : model[county].parameters.template get<TimeInfectedNoSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
318 1188 : t_InfectedSymptoms[county].push_back(static_cast<int>(std::round(
319 1188 : model[county].parameters.template get<TimeInfectedSymptoms<double>>()[(AgeGroup)group] * reduc_t)));
320 1188 : t_InfectedSevere[county].push_back(static_cast<int>(
321 1188 : std::round(model[county].parameters.template get<TimeInfectedSevere<double>>()[(AgeGroup)group])));
322 1188 : t_InfectedCritical[county].push_back(static_cast<int>(
323 1188 : std::round(model[county].parameters.template get<TimeInfectedCritical<double>>()[(AgeGroup)group])));
324 :
325 594 : double reduc_immune_exp =
326 1188 : model[county].parameters.template get<ReducExposedImprovedImmunity<double>>()[(AgeGroup)group];
327 594 : double reduc_immune_inf =
328 1188 : model[county].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[(AgeGroup)group];
329 594 : double reduc_immune_hosp =
330 1188 : model[county].parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[(
331 : AgeGroup)group];
332 594 : double reduc_immune_icu =
333 1188 : model[county].parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[(
334 : AgeGroup)group];
335 1188 : mu_C_R[county].push_back((
336 1188 : 1 - reduc_immune_inf / reduc_immune_exp *
337 1188 : (1 - model[county]
338 1188 : .parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[(AgeGroup)group])));
339 1188 : mu_I_H[county].push_back(
340 1782 : reduc_immune_hosp / reduc_immune_inf *
341 1188 : model[county].parameters.template get<SeverePerInfectedSymptoms<double>>()[(AgeGroup)group]);
342 : // transfer from H to U, D unchanged.
343 1188 : mu_H_U[county].push_back(
344 1782 : reduc_immune_icu / reduc_immune_hosp *
345 1188 : model[county].parameters.template get<CriticalPerSevere<double>>()[(AgeGroup)group]);
346 : }
347 : }
348 :
349 99 : BOOST_OUTCOME_TRY(read_confirmed_cases_data(case_data, region, date, num_Exposed, num_InfectedNoSymptoms,
350 : num_InfectedSymptoms, num_InfectedSevere, num_icu, num_death, num_rec,
351 : t_Exposed, t_InfectedNoSymptoms, t_InfectedSymptoms, t_InfectedSevere,
352 : t_InfectedCritical, mu_C_R, mu_I_H, mu_H_U, scaling_factor_inf));
353 :
354 198 : for (size_t county = 0; county < model.size(); county++) {
355 99 : size_t num_groups = (size_t)model[county].parameters.get_num_groups();
356 693 : for (size_t i = 0; i < num_groups; i++) {
357 594 : model[county].populations[{AgeGroup(i), InfectionState::ExposedImprovedImmunity}] = num_Exposed[county][i];
358 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsImprovedImmunity}] =
359 : num_InfectedNoSymptoms[county][i];
360 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] = 0;
361 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunity}] =
362 : num_InfectedSymptoms[county][i];
363 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] = 0;
364 594 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereImprovedImmunity}] =
365 : num_InfectedSevere[county][i];
366 : // Only set the number of ICU patients here, if the date is not available in the data.
367 594 : if (!is_divi_data_available(date)) {
368 216 : model[county].populations[{AgeGroup(i), InfectionState::InfectedCriticalImprovedImmunity}] =
369 : num_icu[county][i];
370 : }
371 : }
372 99 : if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) {
373 27 : log_warning("No infections for vaccinated reported on date {} for region {}. "
374 : "Population data has not been set.",
375 : date, region[county]);
376 : }
377 : }
378 :
379 198 : return success();
380 99 : }
381 :
382 : /**
383 : * @brief sets populations data from a transformed RKI cases file into a Model.
384 : * @param[in, out] model vector of objects in which the data is set
385 : * @param[in] path Path to transformed RKI cases file
386 : * @param[in] region vector of keys of the region of interest
387 : * @param[in] date Date for which the arrays are initialized
388 : * @param[in] scaling_factor_inf factors by which to scale the confirmed cases of
389 : * rki data
390 : * @param set_death[in] If true, set the number of deaths.
391 : */
392 : template <class Model>
393 54 : IOResult<void> set_confirmed_cases_data(std::vector<Model>& model, const std::string& path,
394 : std::vector<int> const& region, Date date,
395 : const std::vector<double>& scaling_factor_inf, bool set_death = false)
396 : {
397 54 : BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path));
398 54 : BOOST_OUTCOME_TRY(set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, set_death));
399 108 : return success();
400 54 : }
401 :
402 : /**
403 : * @brief reads number of ICU patients from DIVI register into Parameters
404 : * @param[in] path Path to transformed DIVI file
405 : * @param[in] vregion Keys of the region of interest
406 : * @param[in] date Date for which the arrays are initialized
407 : * @param[in, out] vnum_icu number of ICU patients
408 : * @see mio::read_divi_data
409 : * @{
410 : */
411 : IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& vregion, Date date,
412 : std::vector<double>& vnum_icu);
413 : IOResult<void> read_divi_data(const std::vector<DiviEntry>& divi_data, const std::vector<int>& vregion, Date date,
414 : std::vector<double>& vnum_icu);
415 : /**@}*/
416 :
417 : /**
418 : * @brief sets populations data from DIVI register into Model
419 : * @param[in, out] model vector of objects in which the data is set
420 : * @param[in] path Path to transformed DIVI file
421 : * @param[in] vregion vector of keys of the regions of interest
422 : * @param[in] date Date for which the arrays are initialized
423 : * @param[in] scaling_factor_icu factor by which to scale the icu cases of divi data
424 : */
425 : template <class Model>
426 108 : IOResult<void> set_divi_data(std::vector<Model>& model, const std::string& path, const std::vector<int>& vregion,
427 : Date date, double scaling_factor_icu)
428 : {
429 : // DIVI dataset will no longer be updated from CW29 2024 on.
430 108 : if (!is_divi_data_available(date)) {
431 45 : log_warning("No DIVI data available for date: {}. "
432 : "ICU compartment will be set based on Case data.",
433 : date);
434 90 : return success();
435 : }
436 126 : std::vector<double> sum_mu_I_U(vregion.size(), 0);
437 126 : std::vector<std::vector<double>> mu_I_U{model.size()};
438 126 : for (size_t region = 0; region < vregion.size(); region++) {
439 63 : auto num_groups = model[region].parameters.get_num_groups();
440 441 : for (auto i = AgeGroup(0); i < num_groups; i++) {
441 756 : sum_mu_I_U[region] += model[region].parameters.template get<CriticalPerSevere<double>>()[i] *
442 378 : model[region].parameters.template get<SeverePerInfectedSymptoms<double>>()[i];
443 756 : mu_I_U[region].push_back(model[region].parameters.template get<CriticalPerSevere<double>>()[i] *
444 378 : model[region].parameters.template get<SeverePerInfectedSymptoms<double>>()[i]);
445 : }
446 : }
447 126 : std::vector<double> num_icu(model.size(), 0.0);
448 63 : BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu));
449 :
450 126 : for (size_t region = 0; region < vregion.size(); region++) {
451 63 : auto num_groups = model[region].parameters.get_num_groups();
452 441 : for (auto i = AgeGroup(0); i < num_groups; i++) {
453 378 : model[region].populations[{i, InfectionState::InfectedCriticalNaive}] =
454 378 : scaling_factor_icu * num_icu[region] * mu_I_U[region][(size_t)i] / sum_mu_I_U[region];
455 : }
456 : }
457 :
458 126 : return success();
459 63 : }
460 :
461 : /**
462 : * @brief reads population data from census data.
463 : * @param[in] path Path to population data file.
464 : * @param[in] vregion vector of keys of the regions of interest
465 : * @see mio::read_population_data
466 : * @{
467 : */
468 : IOResult<std::vector<std::vector<double>>> read_population_data(const std::string& path,
469 : const std::vector<int>& vregion);
470 : IOResult<std::vector<std::vector<double>>> read_population_data(const std::vector<PopulationDataEntry>& population_data,
471 : const std::vector<int>& vregion);
472 : /**@}*/
473 :
474 : /**
475 : * @brief sets population data from census data which has been read into num_population
476 : * @param[in, out] model vector of objects in which the data is set
477 : * @param[in] num_population vector of population data
478 : * @param[in] case_data vector of case data. Each inner vector represents a different region
479 : * @param[in] vregion vector of keys of the regions of interest
480 : * @param[in] date Date for which the arrays are initialized
481 : */
482 : template <class Model>
483 108 : IOResult<void> set_population_data(std::vector<Model>& model, const std::vector<std::vector<double>>& num_population,
484 : const std::vector<ConfirmedCasesDataEntry>& case_data,
485 : const std::vector<int>& vregion, Date date)
486 : {
487 108 : auto num_age_groups = ConfirmedCasesDataEntry::age_group_names.size();
488 432 : std::vector<std::vector<double>> vnum_rec(model.size(), std::vector<double>(num_age_groups, 0.0));
489 :
490 108 : BOOST_OUTCOME_TRY(read_confirmed_cases_data_fix_recovered(case_data, vregion, date, vnum_rec, 14.));
491 :
492 216 : for (size_t region = 0; region < vregion.size(); region++) {
493 108 : if (std::accumulate(num_population[region].begin(), num_population[region].end(), 0.0) > 0) {
494 99 : auto num_groups = model[region].parameters.get_num_groups();
495 693 : for (auto i = AgeGroup(0); i < num_groups; i++) {
496 :
497 1188 : double S_v = std::min(
498 1188 : model[region].parameters.template get<DailyFullVaccinations<double>>()[{i, SimulationDay(0)}] +
499 : vnum_rec[region][size_t(i)],
500 : num_population[region][size_t(i)]);
501 594 : double S_pv = std::max(
502 2376 : model[region].parameters.template get<DailyPartialVaccinations<double>>()[{i, SimulationDay(0)}] -
503 1188 : model[region].parameters.template get<DailyFullVaccinations<double>>()[{i, SimulationDay(0)}],
504 1188 : 0.0); // use std::max with 0
505 : double S;
506 594 : if (num_population[region][size_t(i)] - S_pv - S_v < 0.0) {
507 9 : log_warning("Number of vaccinated persons greater than population in county {}, age group {}.",
508 18 : region, size_t(i));
509 9 : S = 0.0;
510 9 : S_v = num_population[region][size_t(i)] - S_pv;
511 : }
512 : else {
513 585 : S = num_population[region][size_t(i)] - S_pv - S_v;
514 : }
515 :
516 594 : double denom_E =
517 594 : 1 / (S + S_pv * model[region].parameters.template get<ReducExposedPartialImmunity<double>>()[i] +
518 594 : S_v * model[region].parameters.template get<ReducExposedImprovedImmunity<double>>()[i]);
519 594 : double denom_C =
520 594 : 1 / (S + S_pv * model[region].parameters.template get<ReducExposedPartialImmunity<double>>()[i] +
521 594 : S_v * model[region].parameters.template get<ReducExposedImprovedImmunity<double>>()[i]);
522 594 : double denom_I =
523 : 1 /
524 594 : (S +
525 594 : S_pv * model[region].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i] +
526 594 : S_v * model[region].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i]);
527 594 : double denom_HU =
528 : 1 /
529 594 : (S +
530 1188 : S_pv * model[region]
531 594 : .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] +
532 1188 : S_v * model[region]
533 594 : .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i]);
534 :
535 594 : model[region].populations[{i, InfectionState::ExposedNaive}] =
536 1188 : S * model[region].populations[{i, InfectionState::ExposedNaive}] * denom_E;
537 594 : model[region].populations[{i, InfectionState::ExposedPartialImmunity}] =
538 1188 : S_pv * model[region].parameters.template get<ReducExposedPartialImmunity<double>>()[i] *
539 1188 : model[region].populations[{i, InfectionState::ExposedPartialImmunity}] * denom_E;
540 594 : model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] =
541 1188 : S_v * model[region].parameters.template get<ReducExposedImprovedImmunity<double>>()[i] *
542 1188 : model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] * denom_E;
543 :
544 594 : model[region].populations[{i, InfectionState::InfectedNoSymptomsNaive}] =
545 1188 : S * model[region].populations[{i, InfectionState::InfectedNoSymptomsNaive}] * denom_C;
546 594 : model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] =
547 1188 : S_pv * model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] * denom_C;
548 594 : model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] =
549 1188 : S_v * model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] * denom_C;
550 :
551 594 : model[region].populations[{i, InfectionState::InfectedNoSymptomsNaiveConfirmed}] =
552 1188 : S * model[region].populations[{i, InfectionState::InfectedNoSymptomsNaiveConfirmed}] * denom_C;
553 594 : model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] =
554 1188 : S_pv * model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] *
555 : denom_C;
556 594 : model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] =
557 1188 : S_v * model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] *
558 : denom_C;
559 :
560 594 : model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] =
561 1188 : S * model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] * denom_I;
562 594 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] =
563 1188 : S_pv * model[region].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i] *
564 1188 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] * denom_I;
565 594 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] =
566 1188 : S_v * model[region].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i] *
567 1188 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] * denom_I;
568 :
569 594 : model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] =
570 1188 : S * model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] * denom_I;
571 594 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] =
572 1188 : S_pv * model[region].parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i] *
573 1188 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] * denom_I;
574 594 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] =
575 1188 : S_v * model[region].parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i] *
576 1188 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] * denom_I;
577 :
578 594 : model[region].populations[{i, InfectionState::InfectedSevereNaive}] =
579 1188 : S * model[region].populations[{i, InfectionState::InfectedSevereNaive}] * denom_HU;
580 594 : model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] =
581 594 : S_pv *
582 1188 : model[region].parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] *
583 1188 : model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] * denom_HU;
584 594 : model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] =
585 594 : S_v *
586 594 : model[region]
587 1188 : .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i] *
588 1188 : model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] * denom_HU;
589 :
590 594 : model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] =
591 594 : S_pv *
592 1188 : model[region].parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] *
593 1188 : model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU;
594 594 : model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] =
595 594 : S_v *
596 594 : model[region]
597 1188 : .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i] *
598 1188 : model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU;
599 594 : model[region].populations[{i, InfectionState::InfectedCriticalNaive}] =
600 1188 : S * model[region].populations[{i, InfectionState::InfectedCriticalNaive}] * denom_HU;
601 :
602 594 : model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] =
603 1782 : model[region].parameters.template get<DailyFullVaccinations<double>>()[{i, SimulationDay(0)}] +
604 1188 : model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] -
605 1782 : (model[region].populations[{i, InfectionState::InfectedSymptomsNaive}] +
606 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] +
607 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] +
608 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsNaiveConfirmed}] +
609 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] +
610 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] +
611 1782 : model[region].populations[{i, InfectionState::InfectedSevereNaive}] +
612 1782 : model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] +
613 1782 : model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] +
614 1782 : model[region].populations[{i, InfectionState::InfectedCriticalNaive}] +
615 1782 : model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] +
616 1782 : model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] +
617 1782 : model[region].populations[{i, InfectionState::DeadNaive}] +
618 1782 : model[region].populations[{i, InfectionState::DeadPartialImmunity}] +
619 1188 : model[region].populations[{i, InfectionState::DeadImprovedImmunity}]);
620 :
621 1188 : model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = std::min(
622 1782 : S_v - model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] -
623 1782 : model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] -
624 1782 : model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] -
625 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] -
626 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] -
627 1782 : model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] -
628 2970 : model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] -
629 1188 : model[region].populations[{i, InfectionState::DeadImprovedImmunity}],
630 1188 : std::max(0.0, double(model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}])));
631 :
632 594 : model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = std::max(
633 1188 : 0.0,
634 1782 : S_pv - model[region].populations[{i, InfectionState::ExposedPartialImmunity}] -
635 1782 : model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] -
636 1782 : model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] -
637 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] -
638 1782 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] -
639 1782 : model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] -
640 2970 : model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] -
641 1188 : model[region].populations[{i, InfectionState::DeadPartialImmunity}]);
642 :
643 594 : model[region].populations.template set_difference_from_group_total<AgeGroup>(
644 : {i, InfectionState::SusceptibleNaive}, num_population[region][size_t(i)]);
645 : }
646 :
647 693 : for (auto i = AgeGroup(0); i < AgeGroup(6); i++) {
648 16632 : for (auto j = Index<InfectionState>(0); j < InfectionState::Count; ++j) {
649 16038 : if (model[region].populations[{i, j}] < 0) {
650 0 : log_warning("Compartment at age group {}, infection state {}, is negative: {}", size_t(i),
651 0 : size_t(j), model[region].populations[{i, j}] / num_population[region][size_t(i)]);
652 : }
653 : }
654 : }
655 : }
656 : else {
657 9 : log_warning("No population data available for region " + std::to_string(region) +
658 : ". Population data has not been set.");
659 : }
660 : }
661 :
662 216 : return success();
663 108 : }
664 :
665 : /**
666 : * @brief Sets population data from census data which has been read into num_population.
667 : * @param[in, out] model Vector of objects in which the data is set.
668 : * @param[in] path Path to population data file.
669 : * @param[in] path_rki Path to RKI cases data file.
670 : * @param[in] vregion Vector of keys of the regions of interest.
671 : * @param[in] date Date for which the arrays are initialized.
672 : */
673 : template <class Model>
674 72 : IOResult<void> set_population_data(std::vector<Model>& model, const std::string& path, const std::string& path_rki,
675 : const std::vector<int>& vregion, Date date)
676 : {
677 72 : BOOST_OUTCOME_TRY(auto&& num_population, details::read_population_data(path, vregion));
678 72 : BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path_rki));
679 :
680 72 : BOOST_OUTCOME_TRY(set_population_data(model, num_population, rki_data, vregion, date));
681 144 : return success();
682 72 : }
683 :
684 : /**
685 : * @brief Sets vaccination data for models stored in a vector.
686 : *
687 : * @tparam FP Floating point type used in the Model objects.
688 : * @param[in, out] model Vector of Model objects in which the vaccination data is set.
689 : * @param[in] vacc_data Vector of VaccinationDataEntry objects containing the vaccination data.
690 : * @param[in] date Start date for the simulation.
691 : * @param[in] vregion Vector of region identifiers.
692 : * @param[in] num_days Number of days for which the simulation is run.
693 : */
694 : template <typename FP = double>
695 90 : IOResult<void> set_vaccination_data(std::vector<Model<FP>>& model, const std::vector<VaccinationDataEntry>& vacc_data,
696 : Date date, const std::vector<int>& vregion, int num_days)
697 : {
698 90 : auto num_groups = model[0].parameters.get_num_groups();
699 :
700 : // type conversion from UncertainValue -> FP -> int
701 90 : auto days_until_effective1 = static_cast<int>(
702 90 : static_cast<FP>(model[0].parameters.template get<DaysUntilEffectivePartialImmunity<FP>>()[AgeGroup(0)]));
703 90 : auto days_until_effective2 = static_cast<int>(
704 90 : static_cast<FP>(model[0].parameters.template get<DaysUntilEffectiveImprovedImmunity<FP>>()[AgeGroup(0)]));
705 90 : auto vaccination_distance =
706 90 : static_cast<int>(static_cast<FP>(model[0].parameters.template get<VaccinationGap<FP>>()[AgeGroup(0)]));
707 :
708 : // iterate over regions (e.g., counties)
709 180 : for (size_t i = 0; i < model.size(); ++i) {
710 : // iterate over age groups in region
711 630 : for (auto g = AgeGroup(0); g < num_groups; ++g) {
712 :
713 540 : model[i].parameters.template get<DailyPartialVaccinations<FP>>().resize(SimulationDay(num_days + 1));
714 540 : model[i].parameters.template get<DailyFullVaccinations<FP>>().resize(SimulationDay(num_days + 1));
715 3132 : for (auto d = SimulationDay(0); d < SimulationDay(num_days + 1); ++d) {
716 2592 : model[i].parameters.template get<DailyPartialVaccinations<FP>>()[{g, d}] = 0.0;
717 2592 : model[i].parameters.template get<DailyFullVaccinations<FP>>()[{g, d}] = 0.0;
718 : }
719 : }
720 : }
721 :
722 90 : auto max_date_entry = std::max_element(vacc_data.begin(), vacc_data.end(), [](auto&& a, auto&& b) {
723 46728 : return a.date < b.date;
724 : });
725 90 : if (max_date_entry == vacc_data.end()) {
726 0 : return failure(StatusCode::InvalidFileFormat, "Vaccination data file is empty.");
727 : }
728 90 : auto max_date = max_date_entry->date;
729 :
730 46908 : for (auto&& vacc_data_entry : vacc_data) {
731 118674 : auto it = std::find_if(vregion.begin(), vregion.end(), [&vacc_data_entry](auto&& r) {
732 114183 : return r == 0 || (vacc_data_entry.county_id && vacc_data_entry.county_id == regions::CountyId(r)) ||
733 120582 : (vacc_data_entry.state_id && vacc_data_entry.state_id == regions::StateId(r)) ||
734 116091 : (vacc_data_entry.district_id && vacc_data_entry.district_id == regions::DistrictId(r));
735 : });
736 46818 : auto date_df = vacc_data_entry.date;
737 46818 : if (it != vregion.end()) {
738 46818 : auto region_idx = size_t(it - vregion.begin());
739 46818 : auto age = vacc_data_entry.age_group;
740 :
741 267156 : for (size_t d = 0; d < (size_t)num_days + 1; ++d) {
742 : int days_plus;
743 : // In the following, second dose means previous 'full immunization', now 'Grundimmunisierung'.
744 : // ---
745 : // date: start_date of the simulation (Input from IO call read_input_data_county_vaccmodel())
746 : // d: day of simulation, counted from 0 to num_days (for which we need (approximated) vaccination numbers)
747 : // root[i]["Vacc_completed"]: accumulated number of total second doses up to day date_df;
748 : // taken from input dataframe, single value, per county and age group
749 : // ----
750 : // An averaged distance between first and second doses (vaccination_distance) is assumed in the following
751 : // and the first doses are computed based on the second doses given 'vaccination_distance' days later.
752 : // ----
753 : // a person whose second dose is reported at start_date + simulation_day - days_until_effective1 + vaccination_distance
754 : // had the first dose on start_date + simulation_day - days_until_effective1. Furthermore, he/she has the full protection
755 : // of the first dose at day X = start_date + simulation_day
756 : // Storing its value in get<DailyPartialVaccinations>() will eventually (in the simulation)
757 : // transfer the difference (between get<DailyPartialVaccinations>() at d and d-1) of
758 : // N susceptible individuals to 'Susceptible Partially Vaccinated' state at day d; see secir_vaccinated.h
759 220338 : auto offset_first_date =
760 220338 : offset_date_by_days(date, (int)d - days_until_effective1 + vaccination_distance);
761 220338 : if (max_date >= offset_first_date) {
762 : // Option 1: considered offset_first_date is available in input data frame
763 13950 : if (date_df == offset_first_date) {
764 0 : model[region_idx]
765 0 : .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] =
766 0 : vacc_data_entry.num_vaccinations_completed;
767 : }
768 : }
769 : else { // offset_first_date > max_date
770 : // Option 2: considered offset_first_date is NOT available in input data frame
771 : // Here, a constant number of first and second doses is assumed, i.e.,
772 : // the the number of vaccinationes at day d (N days after max_date) will be:
773 : // total number of vaccinations up to day max_date + N * number of vaccinations ON max_date
774 : // (where the latter is computed as the difference between the total number at max_date and max_date-1)
775 206388 : days_plus = get_offset_in_days(offset_first_date, max_date);
776 206388 : if (date_df == offset_date_by_days(max_date, -1)) {
777 2430 : model[region_idx]
778 4860 : .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] -=
779 2430 : days_plus * vacc_data_entry.num_vaccinations_completed;
780 : }
781 203958 : else if (date_df == max_date) {
782 2430 : model[region_idx]
783 4860 : .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] +=
784 2430 : (days_plus + 1) * vacc_data_entry.num_vaccinations_completed;
785 : }
786 : }
787 :
788 : // a person whose second dose is reported at start_date + simulation_day - days_until_effective2
789 : // has the full protection of the second dose at day X = start_date + simulation_day
790 : // Storing its value in get<DailyFullVaccinations>() will eventually (in the simulation)
791 : // transfer the difference (between get<DailyFullVaccinations>() at d and d-1) of
792 : // N susceptible, partially vaccinated individuals to 'SusceptibleImprovedImmunity' state at day d; see secir_vaccinated.h
793 220338 : auto offset_full_date = offset_date_by_days(date, (int)d - days_until_effective2);
794 220338 : if (max_date >= offset_full_date) {
795 : // Option 1: considered offset_full_date is available in input data frame
796 220338 : if (date_df == offset_full_date) {
797 2430 : model[region_idx]
798 2430 : .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] =
799 2430 : vacc_data_entry.num_vaccinations_completed;
800 : }
801 : }
802 : else { // offset_full_date > max_full_date
803 : // Option 2: considered offset_full_date is NOT available in input data frame
804 0 : days_plus = get_offset_in_days(offset_full_date, max_date);
805 0 : if (date_df == offset_date_by_days(max_date, -1)) {
806 0 : model[region_idx]
807 0 : .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] -=
808 0 : days_plus * vacc_data_entry.num_vaccinations_completed;
809 : }
810 0 : else if (date_df == max_date) {
811 0 : model[region_idx]
812 0 : .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] +=
813 0 : (days_plus + 1) * vacc_data_entry.num_vaccinations_completed;
814 : }
815 : }
816 : }
817 : }
818 : }
819 180 : return success();
820 : }
821 :
822 : /**
823 : * @brief Reads vaccination data from a file and sets it for each model.
824 : *
825 : * @tparam FP Floating point type used in the Model objects.
826 : * @param[in, out] model Vector of Model objects in which the vaccination data is set.
827 : * @param[in] path Path to vaccination data file.
828 : * @param[in] date Start date for the simulation.
829 : * @param[in] vregion Vector of region identifiers.
830 : * @param[in] num_days Number of days for which the simulation is run.
831 : */
832 : template <typename FP = double>
833 54 : IOResult<void> set_vaccination_data(std::vector<Model<FP>>& model, const std::string& path, Date date,
834 : const std::vector<int>& vregion, int num_days)
835 : {
836 54 : BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path));
837 54 : BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days));
838 108 : return success();
839 54 : }
840 :
841 : } // namespace details
842 :
843 : #ifdef MEMILIO_HAS_HDF5
844 :
845 : /**
846 : * @brief Uses the initialisation method, which uses the reported data to set the initial conditions for the model for a given day.
847 : * The initialisation is applied for a predefined number of days and finally saved in a timeseries for each region. In the end,
848 : * we save the files "Results_rki.h5" and "Results_rki_sum.h5" in the results_dir.
849 : * Results_rki.h5 contains a time series for each region and Results_rki_sum.h5 contains the sum of all regions.
850 : * @param[in] models Vector of models in which the data is set. Copy is made to avoid changing the original model.
851 : * @param[in] results_dir Path to result files.
852 : * @param[in] counties Vector of keys of the counties of interest.
853 : * @param[in] date Date for which the data should be read.
854 : * @param[in] scaling_factor_inf Factors by which to scale the confirmed cases of rki data.
855 : * @param[in] scaling_factor_icu Factor by which to scale the icu cases of divi data.
856 : * @param[in] num_days Number of days to be simulated/initialized.
857 : * @param[in] divi_data_path Path to DIVI file.
858 : * @param[in] confirmed_cases_path Path to confirmed cases file.
859 : * @param[in] population_data_path Path to population data file.
860 : * @param[in] vaccination_data_path Path to vaccination data file.
861 : */
862 : template <class Model>
863 18 : IOResult<void> export_input_data_county_timeseries(
864 : std::vector<Model> models, const std::string& results_dir, const std::vector<int>& counties, Date date,
865 : const std::vector<double>& scaling_factor_inf, const double scaling_factor_icu, const int num_days,
866 : const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path,
867 : const std::string& vaccination_data_path = "")
868 : {
869 18 : const auto num_groups = (size_t)models[0].parameters.get_num_groups();
870 18 : assert(scaling_factor_inf.size() == num_groups);
871 18 : assert(num_groups == ConfirmedCasesDataEntry::age_group_names.size());
872 18 : assert(models.size() == counties.size());
873 72 : std::vector<TimeSeries<double>> extrapolated_data(
874 18 : models.size(), TimeSeries<double>::zero(num_days + 1, (size_t)InfectionState::Count * num_groups));
875 :
876 18 : BOOST_OUTCOME_TRY(auto&& case_data, read_confirmed_cases_data(confirmed_cases_path));
877 18 : BOOST_OUTCOME_TRY(auto&& population_data, details::read_population_data(population_data_path, counties));
878 :
879 : // empty vector if set_vaccination_data is not set
880 18 : std::vector<VaccinationDataEntry> vacc_data;
881 18 : if (!vaccination_data_path.empty()) {
882 18 : BOOST_OUTCOME_TRY(vacc_data, read_vaccination_data(vaccination_data_path));
883 18 : }
884 :
885 90 : for (int t = 0; t <= num_days; ++t) {
886 36 : auto offset_day = offset_date_by_days(date, t);
887 :
888 36 : if (!vaccination_data_path.empty()) {
889 36 : BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days));
890 36 : }
891 :
892 : // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
893 : // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
894 36 : BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu));
895 :
896 36 : BOOST_OUTCOME_TRY(
897 : details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf, true));
898 36 : BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, case_data, counties, offset_day));
899 :
900 72 : for (size_t r = 0; r < counties.size(); r++) {
901 36 : extrapolated_data[r][t] = models[r].get_initial_values();
902 : // in set_population_data the number of death individuals is subtracted from the SusceptibleImprovedImmunity compartment.
903 : // Since we should be independent whether we consider them or not, we add them back here before we save the data.
904 252 : for (size_t age = 0; age < num_groups; age++) {
905 216 : extrapolated_data[r][t][(size_t)InfectionState::SusceptibleImprovedImmunity +
906 432 : age * (size_t)InfectionState::Count] +=
907 216 : extrapolated_data[r][t][(size_t)InfectionState::DeadNaive + age * (size_t)InfectionState::Count];
908 : }
909 : }
910 : }
911 36 : BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast<int>(num_groups),
912 : path_join(results_dir, "Results_rki.h5")));
913 :
914 90 : auto extrapolated_rki_data_sum = sum_nodes(std::vector<std::vector<TimeSeries<double>>>{extrapolated_data});
915 144 : BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast<int>(num_groups),
916 : path_join(results_dir, "Results_rki_sum.h5")));
917 :
918 36 : return success();
919 18 : }
920 : #else
921 : template <class Model>
922 : IOResult<void> export_input_data_county_timeseries(std::vector<Model>, const std::string&, const std::vector<int>&,
923 : Date, const std::vector<double>&, const double, const int,
924 : const std::string&, const std::string&, const std::string&,
925 : const std::string&)
926 : {
927 : mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
928 : return success();
929 : }
930 :
931 : #endif //MEMILIO_HAS_HDF5
932 :
933 : /**
934 : * Reads compartments for German counties at a specified date from data files.
935 : * Estimates all compartments from available data using the model parameters, so the
936 : * model parameters must be set before calling this function.
937 : * Uses data files that contain centered 7-day moving average.
938 : * @param[in, out] model Vector of SECIRVVS models, one per county.
939 : * @param[in] date Date for which the data should be read.
940 : * @param[in] county Ids of the counties.
941 : * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county.
942 : * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting.
943 : * @param[in] dir Directory that contains the data files.
944 : * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation.
945 : * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files.
946 : */
947 : template <class Model>
948 27 : IOResult<void> read_input_data_county(std::vector<Model>& model, Date date, const std::vector<int>& county,
949 : const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
950 : const std::string& dir, int num_days, bool export_time_series = false)
951 : {
952 54 : BOOST_OUTCOME_TRY(details::set_vaccination_data(
953 : model, path_join(dir, "pydata/Germany", "all_county_ageinf_vacc_ma7.json"), date, county, num_days));
954 :
955 : // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
956 : // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
957 54 : BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"), county,
958 : date, scaling_factor_icu));
959 :
960 54 : BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(
961 : model, path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"), county, date, scaling_factor_inf));
962 81 : BOOST_OUTCOME_TRY(
963 : details::set_population_data(model, path_join(dir, "pydata/Germany", "county_current_population.json"),
964 : path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"), county, date));
965 :
966 27 : if (export_time_series) {
967 : // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
968 : // Run time equals run time of the previous functions times the num_days !
969 : // (This only represents the vectorization of the previous function over all simulation days...)
970 0 : log_warning("Exporting time series of extrapolated real data. This may take some minutes. "
971 : "For simulation runs over the same time period, deactivate it.");
972 0 : BOOST_OUTCOME_TRY(
973 : export_input_data_county_timeseries(model, dir, county, date, scaling_factor_inf, scaling_factor_icu,
974 : num_days, path_join(dir, "pydata/Germany", "county_divi_ma7.json"),
975 : path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"),
976 : path_join(dir, "pydata/Germany", "county_current_population.json"),
977 : path_join(dir, "pydata/Germany", "all_county_ageinf_vacc_ma7.json")));
978 0 : }
979 :
980 54 : return success();
981 27 : }
982 :
983 : /**
984 : * Reads compartments for German counties at a specified date from data files.
985 : * Estimates all compartments from available data using the model parameters, so the
986 : * model parameters must be set before calling this function.
987 : * Uses data files that contain centered 7-day moving average.
988 : * @param[in, out] model Vector of SECIRVVS models, one per county.
989 : * @param[in] date Date for which the data should be read.
990 : * @param[in] node_ids Ids of the nodes.
991 : * @param[in] scaling_factor_inf Factor of confirmed cases to account for undetected cases in each county.
992 : * @param[in] scaling_factor_icu Factor of ICU cases to account for underreporting.
993 : * @param[in] dir Directory that contains the data files.
994 : * @param[in] num_days Number of days to be simulated; required to load data for vaccinations during the simulation.
995 : * @param[in] export_time_series If true, reads data for each day of simulation and writes it in the same directory as the input files.
996 : */
997 : template <class Model>
998 27 : IOResult<void> read_input_data(std::vector<Model>& model, Date date, const std::vector<int>& node_ids,
999 : const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
1000 : const std::string& data_dir, int num_days, bool export_time_series = false)
1001 : {
1002 :
1003 54 : BOOST_OUTCOME_TRY(
1004 : details::set_vaccination_data(model, path_join(data_dir, "vaccination_data.json"), date, node_ids, num_days));
1005 :
1006 : // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
1007 : // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
1008 54 : BOOST_OUTCOME_TRY(
1009 : details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date, scaling_factor_icu));
1010 :
1011 54 : BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(data_dir, "confirmed_cases.json"), node_ids,
1012 : date, scaling_factor_inf));
1013 81 : BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(data_dir, "population_data.json"),
1014 : path_join(data_dir, "confirmed_cases.json"), node_ids, date));
1015 :
1016 27 : if (export_time_series) {
1017 : // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
1018 : // Run time equals run time of the previous functions times the num_days !
1019 : // (This only represents the vectorization of the previous function over all simulation days...)
1020 0 : log_warning("Exporting time series of extrapolated real data. This may take some minutes. "
1021 : "For simulation runs over the same time period, deactivate it.");
1022 0 : BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
1023 : model, data_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days,
1024 : path_join(data_dir, "divi_data.json"), path_join(data_dir, "confirmed_cases.json"),
1025 : path_join(data_dir, "population_data.json"), path_join(data_dir, "vaccination_data.json")));
1026 0 : }
1027 :
1028 54 : return success();
1029 27 : }
1030 :
1031 : } // namespace osecirvvs
1032 : } // namespace mio
1033 :
1034 : #endif // MEMILIO_HAS_JSONCPP
1035 :
1036 : #endif // MIO_ODE_SECIRVVS_PARAMETERS_IO_H
|
