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