Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Henrik Zunker, 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_SECIRTS_PARAMETERS_IO_H
21 : #define MIO_ODE_SECIRTS_PARAMETERS_IO_H
22 :
23 : #include "memilio/config.h"
24 :
25 : #ifdef MEMILIO_HAS_JSONCPP
26 :
27 : #include "ode_secirts/model.h"
28 : #include "ode_secirts/analyze_result.h"
29 : #include "memilio/math/eigen_util.h"
30 : #include "memilio/mobility/graph.h"
31 : #include "memilio/mobility/metapopulation_mobility_instant.h"
32 : #include "memilio/io/epi_data.h"
33 : #include "memilio/io/io.h"
34 : #include "memilio/io/json_serializer.h"
35 : #include "memilio/io/result_io.h"
36 : #include "memilio/utils/logging.h"
37 : #include "memilio/utils/date.h"
38 :
39 : namespace mio
40 : {
41 : namespace osecirts
42 : {
43 :
44 : namespace details
45 : {
46 :
47 : /**
48 : * @brief Gets the region ID (county, state, or district) of an EpiDataEntry.
49 : *
50 : * If none are available, it defaults to 0 which is representing the whole country.
51 : *
52 : * @tparam EpiDataEntry The type of the data entry.
53 : * @param data_entry The (RKI) data entry to extract the region ID from.
54 : * @return The region ID as integer, or 0 if no specific region information is available.
55 : */
56 : template <class EpiDataEntry>
57 52164 : int get_region_id(const EpiDataEntry& data_entry)
58 : {
59 104328 : if (data_entry.county_id) {
60 36432 : return data_entry.county_id->get();
61 : }
62 31464 : if (data_entry.state_id) {
63 0 : return data_entry.state_id->get();
64 : }
65 31464 : if (data_entry.district_id) {
66 15732 : return data_entry.district_id->get();
67 : }
68 0 : return 0;
69 : }
70 :
71 : /**
72 : * @brief Computes the distribution of confirmed cases across infection states based on Case (RKI) data.
73 : *
74 : * This function processes case data for given regions and distributes the cases across different
75 : * infection states, considering the corresponding transition times and probabilities defined in the model.
76 : *
77 : * @tparam Model The type of the model used.
78 : * @tparam FP Floating point type (default: double).
79 : *
80 : * @param[in] case_data Vector of confirmed case data entries (defined in epi_data.h).
81 : * @param[out] vnum_Exposed Output vector for the number of exposed individuals per age group and region.
82 : * @param[out] vnum_InfectedNoSymptoms Output vector for the number of infected individuals without symptoms.
83 : * @param[out] vnum_InfectedSymptoms Output vector for the number of infected individuals with symptoms.
84 : * @param[out] vnum_InfectedSevere Output vector for the number of severely infected individuals.
85 : * @param[out] vnum_icu Output vector for the number of individuals in critical condition (ICU).
86 : * @param[out] vnum_death Output vector for the number of deaths.
87 : * @param[out] vnum_timm_i Output vector for the number of individuals in temporary immunity state.
88 : * @param[in] vregion Vector of region IDs representing the regions in the model vector.
89 : * @param[in] date Date for which the simulation starts.
90 : * @param[in] model Vector of models, each representing a region and containing the parameters.
91 : * @param[in] scaling_factor_inf Vector of scaling factors for confirmed cases for
92 : * @param[in] layer Specifies the immunity layer: 0 (Naive), 1 (Partial Immunity), 2 (Improved Immunity).
93 : *
94 : * @return An IOResult showing success or failure.
95 : */
96 : template <class Model, typename FP = double>
97 216 : IOResult<void> compute_confirmed_cases_data(
98 : const std::vector<ConfirmedCasesDataEntry>& case_data, std::vector<std::vector<FP>>& vnum_Exposed,
99 : std::vector<std::vector<FP>>& vnum_InfectedNoSymptoms, std::vector<std::vector<FP>>& vnum_InfectedSymptoms,
100 : std::vector<std::vector<FP>>& vnum_InfectedSevere, std::vector<std::vector<FP>>& vnum_icu,
101 : std::vector<std::vector<FP>>& vnum_death, std::vector<std::vector<FP>>& vnum_timm_i,
102 : std::vector<int> const& vregion, Date date, const std::vector<Model>& model,
103 : const std::vector<FP>& scaling_factor_inf, const size_t layer, bool read_icu = false)
104 : {
105 216 : auto max_date_entry = std::max_element(case_data.begin(), case_data.end(), [](auto&& a, auto&& b) {
106 114057 : return a.date < b.date;
107 : });
108 216 : if (max_date_entry == case_data.end()) {
109 9 : log_error("Case data file is empty.");
110 9 : return failure(StatusCode::InvalidValue, "Case data is empty.");
111 : }
112 207 : auto max_date = max_date_entry->date;
113 207 : if (max_date < date) {
114 9 : log_error("Specified date does not exist in case data");
115 9 : return failure(StatusCode::OutOfRange, "Case data does not contain specified date.");
116 : }
117 :
118 : // shifts the initilization to the recent past if simulation starts
119 : // around current day and data of the future would be required.
120 : // Only needed for preinfection compartments, exposed and InfectedNoSymptoms.
121 198 : auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E + T_C
122 198 : if (days_surplus > 0) {
123 198 : days_surplus = 0;
124 : }
125 :
126 109494 : for (auto&& entry : case_data) {
127 158976 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
128 158976 : return r == 0 || get_region_id(entry) == r;
129 : });
130 109296 : if (it != vregion.end()) {
131 109296 : auto region_idx = size_t(it - vregion.begin());
132 :
133 109296 : auto params_region = model[region_idx].parameters;
134 109296 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
135 109296 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
136 109296 : auto& num_Exposed = vnum_Exposed[region_idx];
137 109296 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
138 109296 : auto& num_death = vnum_death[region_idx];
139 109296 : auto& num_icu = vnum_icu[region_idx];
140 109296 : auto& num_imm = vnum_timm_i[region_idx];
141 :
142 109296 : auto age = (size_t)entry.age_group;
143 : // (rounded) transition times
144 109296 : const int t_exposed =
145 109296 : static_cast<int>(std::round(params_region.template get<TimeExposed<FP>>()[entry.age_group]));
146 109296 : int t_InfectedNoSymptoms =
147 109296 : static_cast<int>(std::round(params_region.template get<TimeInfectedNoSymptoms<FP>>()[entry.age_group]));
148 109296 : int t_InfectedSymptoms =
149 109296 : static_cast<int>(std::round(params_region.template get<TimeInfectedSymptoms<FP>>()[entry.age_group]));
150 109296 : const int t_InfectedSevere =
151 109296 : static_cast<int>(std::round(params_region.template get<TimeInfectedSevere<FP>>()[entry.age_group]));
152 109296 : const int t_InfectedCritical =
153 109296 : static_cast<int>(std::round(params_region.template get<TimeInfectedCritical<FP>>()[entry.age_group]));
154 109296 : const int t_imm_interval_i = static_cast<int>(
155 109296 : std::round(params_region.template get<TimeTemporaryImmunityPI<FP>>()[entry.age_group]));
156 :
157 : // transition probabilities
158 109296 : FP recoveredPerInfectedNoSymptoms =
159 109296 : params_region.template get<RecoveredPerInfectedNoSymptoms<FP>>()[entry.age_group];
160 109296 : FP severePerInfectedSymptoms = params_region.template get<SeverePerInfectedSymptoms<FP>>()[entry.age_group];
161 109296 : FP criticalPerSevere = params_region.template get<CriticalPerSevere<FP>>()[entry.age_group];
162 :
163 : // if we select a layer with better immunity (layer > 0), we need to adjust the times and transition rates
164 109296 : if (layer > 0) {
165 69552 : t_InfectedNoSymptoms = static_cast<int>(std::round(
166 69552 : t_InfectedNoSymptoms * params_region.template get<ReducTimeInfectedMild<FP>>()[entry.age_group]));
167 69552 : t_InfectedSymptoms = static_cast<int>(std::round(
168 69552 : t_InfectedSymptoms * params_region.template get<ReducTimeInfectedMild<FP>>()[entry.age_group]));
169 :
170 69552 : const FP red_fact_exp =
171 69552 : layer == 1 ? params_region.template get<ReducExposedPartialImmunity<FP>>()[entry.age_group]
172 34776 : : params_region.template get<ReducExposedImprovedImmunity<FP>>()[entry.age_group];
173 :
174 69552 : const FP red_fact_inf =
175 : layer == 1
176 69552 : ? params_region.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[entry.age_group]
177 34776 : : params_region.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[entry.age_group];
178 :
179 69552 : const FP red_fact_sev =
180 : layer == 1
181 69552 : ? params_region
182 34776 : .template get<ReducInfectedSevereCriticalDeadPartialImmunity<FP>>()[entry.age_group]
183 : : params_region
184 34776 : .template get<ReducInfectedSevereCriticalDeadImprovedImmunity<FP>>()[entry.age_group];
185 :
186 69552 : recoveredPerInfectedNoSymptoms = 1 - red_fact_inf / red_fact_exp * (1 - recoveredPerInfectedNoSymptoms);
187 69552 : severePerInfectedSymptoms = red_fact_sev / red_fact_inf * severePerInfectedSymptoms;
188 : }
189 :
190 109296 : if (entry.date == offset_date_by_days(date, 0)) {
191 1188 : num_InfectedSymptoms[age] += scaling_factor_inf[age] * entry.num_confirmed;
192 1188 : num_imm[age] += entry.num_confirmed;
193 : }
194 109296 : if (entry.date == offset_date_by_days(date, t_InfectedNoSymptoms + days_surplus)) {
195 1188 : num_InfectedNoSymptoms[age] +=
196 1188 : 1 / (1 - recoveredPerInfectedNoSymptoms) * scaling_factor_inf[age] * entry.num_confirmed;
197 1188 : num_Exposed[age] -=
198 1188 : 1 / (1 - recoveredPerInfectedNoSymptoms) * scaling_factor_inf[age] * entry.num_confirmed;
199 : }
200 109296 : if (entry.date == offset_date_by_days(date, days_surplus)) {
201 1188 : num_InfectedNoSymptoms[age] -=
202 1188 : 1 / (1 - recoveredPerInfectedNoSymptoms) * scaling_factor_inf[age] * entry.num_confirmed;
203 : }
204 109296 : if (entry.date == offset_date_by_days(date, t_exposed + t_InfectedNoSymptoms + days_surplus)) {
205 1188 : num_Exposed[age] +=
206 1188 : 1 / (1 - recoveredPerInfectedNoSymptoms) * scaling_factor_inf[age] * entry.num_confirmed;
207 : }
208 109296 : if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms)) {
209 1188 : num_InfectedSymptoms[age] -= scaling_factor_inf[age] * entry.num_confirmed;
210 1188 : num_InfectedSevere[age] += severePerInfectedSymptoms * scaling_factor_inf[age] * entry.num_confirmed;
211 : }
212 109296 : if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms - t_InfectedSevere)) {
213 1188 : num_InfectedSevere[age] -= severePerInfectedSymptoms * scaling_factor_inf[age] * entry.num_confirmed;
214 1188 : if (read_icu) {
215 54 : num_icu[age] +=
216 54 : severePerInfectedSymptoms * criticalPerSevere * scaling_factor_inf[age] * entry.num_confirmed;
217 : }
218 : }
219 109296 : if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms - t_InfectedSevere - t_InfectedCritical)) {
220 1188 : num_death[age] += entry.num_deaths;
221 1188 : if (read_icu) {
222 54 : num_icu[age] -=
223 54 : severePerInfectedSymptoms * criticalPerSevere * scaling_factor_inf[age] * entry.num_confirmed;
224 : }
225 : }
226 109296 : if (entry.date == offset_date_by_days(date, 0 - t_imm_interval_i)) {
227 1188 : num_imm[age] -= entry.num_confirmed;
228 : }
229 109296 : }
230 : }
231 :
232 396 : for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) {
233 198 : auto region = vregion[region_idx];
234 :
235 198 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
236 198 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
237 198 : auto& num_Exposed = vnum_Exposed[region_idx];
238 198 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
239 198 : auto& num_death = vnum_death[region_idx];
240 198 : auto& num_icu = vnum_icu[region_idx];
241 198 : auto& num_timm_i = vnum_timm_i[region_idx];
242 :
243 198 : size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size();
244 1386 : for (size_t i = 0; i < num_groups; i++) {
245 9504 : auto try_fix_constraints = [region, i](FP& value, FP error, auto str) {
246 8316 : if (value < error) {
247 : // this should probably return a failure
248 : // but the algorithm is not robust enough to avoid large negative
249 : // values and there are tests that rely on it
250 0 : log_error("{:s} for age group {:s} is {:.4f} for region {:d}, "
251 : "exceeds expected negative value.",
252 0 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
253 0 : value = 0.0;
254 : }
255 8316 : else if (value < 0) {
256 0 : log_info("{:s} for age group {:s} is {:.4f} for region {:d}, "
257 : "automatically corrected",
258 0 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
259 0 : value = 0.0;
260 : }
261 : };
262 :
263 1188 : const FP tol_error = -1e-8;
264 1188 : try_fix_constraints(num_InfectedSymptoms[i], tol_error, "InfectedSymptoms");
265 1188 : try_fix_constraints(num_InfectedNoSymptoms[i], tol_error, "InfectedNoSymptoms");
266 1188 : try_fix_constraints(num_Exposed[i], tol_error, "Exposed");
267 1188 : try_fix_constraints(num_InfectedSevere[i], tol_error, "InfectedSevere");
268 1188 : try_fix_constraints(num_death[i], tol_error, "Dead");
269 1188 : try_fix_constraints(num_icu[i], tol_error, "InfectedCritical");
270 1188 : try_fix_constraints(num_timm_i[i], tol_error, "Recently Recovered or Vaccinated");
271 : }
272 : }
273 :
274 396 : return success();
275 : }
276 :
277 : /**
278 : * @brief Reads confirmed case data from a file and computes the distribution of cases across infection states.
279 : *
280 : * This function reads transformed RKI data from a specified file and processes the confirmed cases
281 : * to distribute them across different infection states and age groups.
282 : *
283 : * @tparam Model The type of the model used.
284 : * @tparam FP Floating point type (default: double).
285 : *
286 : * @param[in] path Path to the file containing transformed case (rki) data.
287 : * @param[in] vregion Vector of region IDs to process.
288 : * @param[in] date Date for which the simulation starts.
289 : * @param[out] vnum_Exposed Output vector for the number of exposed individuals per age group and region.
290 : * @param[out] vnum_InfectedNoSymptoms Output vector for the number of infected individuals without symptoms.
291 : * @param[out] vnum_InfectedSymptoms Output vector for the number of infected individuals with symptoms.
292 : * @param[out] vnum_InfectedSevere Output vector for the number of severely infected individuals (Hospitalized).
293 : * @param[out] vnum_icu Output vector for the number of individuals in critical condition (ICU).
294 : * @param[out] vnum_death Output vector for the number of deaths.
295 : * @param[out] vnum_timm_i Output vector for the number of individuals in a temporary immunity state.
296 : * @param[in] model Vector of models, each representing a region and containing the parameters.
297 : * @param[in] scaling_factor_inf Vector of scaling factors for confirmed cases.
298 : * @param[in] layer Specifies the immunity layer: 0 (Naive), 1 (Partial Immunity), 2 (Improved Immunity).
299 : *
300 : * @return An IOResult indicating success or failure.
301 : */
302 : template <class Model, typename FP = double>
303 18 : IOResult<void> read_confirmed_cases_data(
304 : std::string const& path, std::vector<int> const& vregion, Date date, std::vector<std::vector<FP>>& vnum_Exposed,
305 : std::vector<std::vector<FP>>& vnum_InfectedNoSymptoms, std::vector<std::vector<FP>>& vnum_InfectedSymptoms,
306 : std::vector<std::vector<FP>>& vnum_InfectedSevere, std::vector<std::vector<FP>>& vnum_icu,
307 : std::vector<std::vector<FP>>& vnum_death, std::vector<std::vector<FP>>& vnum_timm_i,
308 : const std::vector<Model>& model, const std::vector<FP>& scaling_factor_inf, const size_t layer,
309 : bool read_icu = false)
310 : {
311 18 : BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path));
312 : return compute_confirmed_cases_data(case_data, vnum_Exposed, vnum_InfectedNoSymptoms, vnum_InfectedSymptoms,
313 : vnum_InfectedSevere, vnum_icu, vnum_death, vnum_timm_i, vregion, date, model,
314 18 : scaling_factor_inf, layer, read_icu);
315 18 : }
316 :
317 : /**
318 : * @brief Sets the confirmed cases data in the model considering different immunity layers.
319 : *
320 : * This function distributes confirmed case data across infection states for regions and age groups
321 : * in the model. It considers different levels of immunity (naive, partial, and improved).
322 : *
323 : * @tparam Model The type of the model used.
324 : * @tparam FP Floating point type (default: double).
325 : *
326 : * @param[in,out] model Vector of models, each representing a region, where the compartments are updated.
327 : * @param[in] case_data Vector of confirmed case data entries.
328 : * @param[in] region Vector of region IDs for which the data is processed.
329 : * @param[in] date Date for which the confirmed cases are set in the model.
330 : * @param[in] scaling_factor_inf Vector of scaling factors for confirmed cases.
331 : * @param[in] immunity_population Vector containing the immunity distribution for naive, partial, and improved immunity layers.
332 : *
333 : * @return An IOResult indicating success or failure.
334 : */
335 : template <class Model, typename FP = double>
336 : IOResult<void>
337 63 : set_confirmed_cases_data(std::vector<Model>& model, const std::vector<ConfirmedCasesDataEntry>& case_data,
338 : std::vector<int> const& region, Date date, const std::vector<FP>& scaling_factor_inf,
339 : const std::vector<std::vector<FP>> immunity_population)
340 : {
341 63 : auto num_age_groups = (size_t)model[0].parameters.get_num_groups();
342 63 : assert(scaling_factor_inf.size() == num_age_groups); //TODO: allow vector or scalar valued scaling factors
343 63 : assert(ConfirmedCasesDataEntry::age_group_names.size() == num_age_groups);
344 :
345 126 : std::vector<std::vector<FP>> num_InfectedSymptoms(model.size());
346 126 : std::vector<std::vector<FP>> num_death(model.size());
347 126 : std::vector<std::vector<FP>> num_Exposed(model.size());
348 126 : std::vector<std::vector<FP>> num_InfectedNoSymptoms(model.size());
349 126 : std::vector<std::vector<FP>> num_InfectedSevere(model.size());
350 126 : std::vector<std::vector<FP>> num_icu(model.size());
351 126 : std::vector<std::vector<FP>> num_timm1(model.size());
352 126 : std::vector<std::vector<FP>> num_timm2(model.size());
353 :
354 126 : std::vector<FP> denom_E(num_age_groups, 0.0);
355 126 : std::vector<FP> denom_I_NS(num_age_groups, 0.0);
356 126 : std::vector<FP> denom_I_Sy(num_age_groups, 0.0);
357 126 : std::vector<FP> denom_I_Sev(num_age_groups, 0.0);
358 :
359 : /*----------- Naive immunity -----------*/
360 126 : for (size_t county = 0; county < model.size(); county++) {
361 63 : num_InfectedSymptoms[county] = std::vector<FP>(num_age_groups, 0.0);
362 63 : num_death[county] = std::vector<FP>(num_age_groups, 0.0);
363 63 : num_Exposed[county] = std::vector<FP>(num_age_groups, 0.0);
364 63 : num_InfectedNoSymptoms[county] = std::vector<FP>(num_age_groups, 0.0);
365 63 : num_InfectedSevere[county] = std::vector<FP>(num_age_groups, 0.0);
366 63 : num_icu[county] = std::vector<FP>(num_age_groups, 0.0);
367 63 : num_timm1[county] = std::vector<FP>(num_age_groups, 0.0);
368 63 : num_timm2[county] = std::vector<FP>(num_age_groups, 0.0);
369 441 : for (size_t group = 0; group < num_age_groups; group++) {
370 : // calculate the denominators to split the reported case numbers to the different immunity layers.
371 378 : denom_E[group] =
372 378 : 1 / (immunity_population[0][group] +
373 756 : immunity_population[1][group] *
374 756 : model[county].parameters.template get<ReducExposedPartialImmunity<FP>>()[(AgeGroup)group] +
375 756 : immunity_population[2][group] *
376 756 : model[county].parameters.template get<ReducExposedImprovedImmunity<FP>>()[(AgeGroup)group]);
377 :
378 378 : denom_I_NS[group] =
379 378 : 1 / (immunity_population[0][group] +
380 756 : immunity_population[1][group] *
381 756 : model[county].parameters.template get<ReducExposedPartialImmunity<FP>>()[(AgeGroup)group] +
382 756 : immunity_population[2][group] *
383 756 : model[county].parameters.template get<ReducExposedImprovedImmunity<FP>>()[(AgeGroup)group]);
384 :
385 378 : denom_I_Sy[group] =
386 378 : 1 / (immunity_population[0][group] +
387 756 : immunity_population[1][group] *
388 378 : model[county]
389 756 : .parameters.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[(AgeGroup)group] +
390 756 : immunity_population[2][group] *
391 378 : model[county]
392 756 : .parameters.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[(AgeGroup)group]);
393 :
394 378 : denom_I_Sev[group] =
395 378 : 1 / (immunity_population[0][group] +
396 756 : immunity_population[1][group] *
397 756 : model[county].parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<FP>>()[(
398 378 : AgeGroup)group] +
399 756 : immunity_population[2][group] *
400 756 : model[county].parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<FP>>()[(
401 : AgeGroup)group]);
402 : }
403 : }
404 :
405 63 : BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms,
406 : num_InfectedSevere, num_icu, num_death, num_timm1, region, date,
407 : model, scaling_factor_inf, 0));
408 :
409 126 : for (size_t county = 0; county < model.size(); county++) {
410 63 : size_t num_groups = (size_t)model[county].parameters.get_num_groups();
411 441 : for (size_t i = 0; i < num_groups; i++) {
412 378 : model[county].populations[{AgeGroup(i), InfectionState::ExposedNaive}] =
413 378 : immunity_population[0][i] * denom_E[i] * num_Exposed[county][i];
414 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsNaive}] =
415 378 : immunity_population[0][i] * denom_I_NS[i] * num_InfectedNoSymptoms[county][i];
416 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsNaiveConfirmed}] = 0;
417 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaive}] =
418 378 : immunity_population[0][i] * denom_I_Sy[i] * num_InfectedSymptoms[county][i];
419 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsNaiveConfirmed}] = 0;
420 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereNaive}] =
421 378 : immunity_population[0][i] * denom_I_Sev[i] * num_InfectedSevere[county][i];
422 : }
423 63 : if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) {
424 0 : log_warning("No infections for unvaccinated reported on date " + std::to_string(date.year) + "-" +
425 : std::to_string(date.month) + "-" + std::to_string(date.day) + " for region " +
426 0 : std::to_string(region[county]) + ". Population data has not been set.");
427 : }
428 : }
429 :
430 : /*----------- PARTIAL Immunity -----------*/
431 126 : for (size_t county = 0; county < model.size(); county++) {
432 63 : num_InfectedSymptoms[county] = std::vector<FP>(num_age_groups, 0.0);
433 63 : num_death[county] = std::vector<FP>(num_age_groups, 0.0);
434 63 : num_Exposed[county] = std::vector<FP>(num_age_groups, 0.0);
435 63 : num_InfectedNoSymptoms[county] = std::vector<FP>(num_age_groups, 0.0);
436 63 : num_InfectedSevere[county] = std::vector<FP>(num_age_groups, 0.0);
437 63 : num_icu[county] = std::vector<FP>(num_age_groups, 0.0);
438 : }
439 :
440 63 : BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms,
441 : num_InfectedSevere, num_icu, num_death, num_timm1, region, date,
442 : model, scaling_factor_inf, 1));
443 :
444 126 : for (size_t county = 0; county < model.size(); county++) {
445 63 : size_t num_groups = (size_t)model[county].parameters.get_num_groups();
446 441 : for (size_t i = 0; i < num_groups; i++) {
447 378 : model[county].populations[{AgeGroup(i), InfectionState::ExposedPartialImmunity}] =
448 756 : immunity_population[1][i] *
449 1134 : model[county].parameters.template get<ReducExposedPartialImmunity<FP>>()[(AgeGroup)i] * denom_E[i] *
450 378 : num_Exposed[county][i];
451 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsPartialImmunity}] =
452 756 : immunity_population[1][i] *
453 1134 : model[county].parameters.template get<ReducExposedPartialImmunity<FP>>()[(AgeGroup)i] * denom_I_NS[i] *
454 378 : num_InfectedNoSymptoms[county][i];
455 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] = 0;
456 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunity}] =
457 756 : immunity_population[1][i] *
458 1134 : model[county].parameters.template get<ReducInfectedSymptomsPartialImmunity<FP>>()[(AgeGroup)i] *
459 756 : denom_I_Sy[i] * num_InfectedSymptoms[county][i];
460 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsPartialImmunityConfirmed}] = 0;
461 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSeverePartialImmunity}] =
462 756 : immunity_population[1][i] *
463 378 : model[county]
464 1134 : .parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<FP>>()[(AgeGroup)i] *
465 756 : denom_I_Sev[i] * num_InfectedSevere[county][i];
466 : // the += is necessary because we already set the previous vaccinated individuals
467 378 : model[county].populations[{AgeGroup(i), InfectionState::TemporaryImmunePartialImmunity}] +=
468 756 : immunity_population[1][i] *
469 1512 : model[county].parameters.template get<ReducExposedPartialImmunity<FP>>()[(AgeGroup)i] * denom_E[i] *
470 378 : num_timm1[county][i];
471 : }
472 63 : if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) {
473 0 : log_warning("No infections for partially vaccinated reported on date " + std::to_string(date.year) + "-" +
474 : std::to_string(date.month) + "-" + std::to_string(date.day) + " for region " +
475 0 : std::to_string(region[county]) + ". Population data has not been set.");
476 : }
477 : }
478 :
479 : /*----------- Improved Immunity -----------*/
480 126 : for (size_t county = 0; county < model.size(); county++) {
481 63 : num_InfectedSymptoms[county] = std::vector<FP>(num_age_groups, 0.0);
482 63 : num_death[county] = std::vector<FP>(num_age_groups, 0.0);
483 63 : num_Exposed[county] = std::vector<FP>(num_age_groups, 0.0);
484 63 : num_InfectedNoSymptoms[county] = std::vector<FP>(num_age_groups, 0.0);
485 63 : num_InfectedSevere[county] = std::vector<FP>(num_age_groups, 0.0);
486 63 : num_icu[county] = std::vector<FP>(num_age_groups, 0.0);
487 : }
488 :
489 63 : BOOST_OUTCOME_TRY(compute_confirmed_cases_data(case_data, num_Exposed, num_InfectedNoSymptoms, num_InfectedSymptoms,
490 : num_InfectedSevere, num_icu, num_death, num_timm2, region, date,
491 : model, scaling_factor_inf, 2));
492 :
493 126 : for (size_t county = 0; county < model.size(); county++) {
494 63 : size_t num_groups = (size_t)model[county].parameters.get_num_groups();
495 441 : for (size_t i = 0; i < num_groups; i++) {
496 378 : model[county].populations[{AgeGroup(i), InfectionState::ExposedImprovedImmunity}] =
497 756 : immunity_population[2][i] *
498 1134 : model[county].parameters.template get<ReducExposedImprovedImmunity<FP>>()[(AgeGroup)i] * denom_E[i] *
499 378 : num_Exposed[county][i];
500 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsImprovedImmunity}] =
501 756 : immunity_population[2][i] *
502 1134 : model[county].parameters.template get<ReducExposedImprovedImmunity<FP>>()[(AgeGroup)i] * denom_I_NS[i] *
503 378 : num_InfectedNoSymptoms[county][i];
504 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] = 0;
505 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunity}] =
506 756 : immunity_population[2][i] *
507 1134 : model[county].parameters.template get<ReducInfectedSymptomsImprovedImmunity<FP>>()[(AgeGroup)i] *
508 756 : denom_I_Sy[i] * num_InfectedSymptoms[county][i];
509 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] = 0;
510 378 : model[county].populations[{AgeGroup(i), InfectionState::InfectedSevereImprovedImmunity}] =
511 756 : immunity_population[2][i] *
512 378 : model[county]
513 1134 : .parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<FP>>()[(AgeGroup)i] *
514 756 : denom_I_Sev[i] * num_InfectedSevere[county][i];
515 : // the += is necessary because we already set the previous vaccinated individuals
516 378 : model[county].populations[{AgeGroup(i), InfectionState::TemporaryImmuneImprovedImmunity}] +=
517 756 : immunity_population[2][i] *
518 1512 : model[county].parameters.template get<ReducExposedImprovedImmunity<FP>>()[(AgeGroup)i] * denom_E[i] *
519 378 : num_timm2[county][i];
520 : }
521 63 : if (std::accumulate(num_InfectedSymptoms[county].begin(), num_InfectedSymptoms[county].end(), 0.0) == 0) {
522 0 : log_warning("No infections for vaccinated reported on date " + std::to_string(date.year) + "-" +
523 : std::to_string(date.month) + "-" + std::to_string(date.day) + " for region " +
524 0 : std::to_string(region[county]) + ". Population data has not been set.");
525 : }
526 : }
527 126 : return success();
528 63 : }
529 :
530 : /**
531 : * @brief Reads confirmed case data from a file and sets it in the model.
532 : *
533 : * This function reads transformed RKI data from the specified file and distributes the confirmed case data
534 : * across different infection states for regions and age groups in the model. It considers naive, partial,
535 : * and improved immunity layers.
536 : *
537 : * @tparam Model The type of the model used.
538 : * @tparam FP Floating point type (default: double).
539 : *
540 : * @param[in,out] model Vector of models, each representing a region, where the compartments are updated.
541 : * @param[in] path Path to the file containing case (RKI) data.
542 : * @param[in] region Vector of region IDs for which the data is processed.
543 : * @param[in] date Date for which the confirmed cases are set in the model.
544 : * @param[in] scaling_factor_inf Vector of scaling factors for confirmed cases.
545 : * @param[in] immunity_population Vector containing the immunity distribution for naive, partial, and improved immunity layers.
546 : *
547 : * @return An IOResult indicating success or failure.
548 : */
549 : template <class Model, typename FP = double>
550 36 : IOResult<void> set_confirmed_cases_data(std::vector<Model>& model, const std::string& path,
551 : std::vector<int> const& region, Date date,
552 : const std::vector<FP>& scaling_factor_inf,
553 : const std::vector<std::vector<FP>> immunity_population)
554 : {
555 36 : BOOST_OUTCOME_TRY(auto&& case_data, mio::read_confirmed_cases_data(path));
556 72 : BOOST_OUTCOME_TRY(
557 : set_confirmed_cases_data(model, case_data, region, date, scaling_factor_inf, immunity_population));
558 72 : return success();
559 36 : }
560 :
561 : /**
562 : * @brief Extracts the number of individuals in critical condition (ICU) for each region
563 : * on a specified date from the provided DIVI data.
564 : *
565 : * @tparam FP Floating point type (default: double).
566 : *
567 : * @param[in] divi_data Vector of DIVI data entries containing date, region, and ICU information.
568 : * @param[in] vregion Vector of region IDs for which the data is computed.
569 : * @param[in] date Date for which the ICU data is computed.
570 : * @param[out] vnum_icu Output vector containing the number of ICU cases for each region.
571 : *
572 : * @return An IOResult indicating success or failure.
573 : */
574 : template <typename FP = double>
575 63 : IOResult<void> compute_divi_data(const std::vector<DiviEntry>& divi_data, const std::vector<int>& vregion, Date date,
576 : std::vector<FP>& vnum_icu)
577 : {
578 63 : auto max_date_entry = std::max_element(divi_data.begin(), divi_data.end(), [](auto&& a, auto&& b) {
579 5733 : return a.date < b.date;
580 : });
581 63 : if (max_date_entry == divi_data.end()) {
582 0 : log_error("DIVI data is empty.");
583 0 : return failure(StatusCode::InvalidValue, "DIVI data is empty.");
584 : }
585 63 : auto max_date = max_date_entry->date;
586 63 : if (max_date < date) {
587 0 : log_error("DIVI data does not contain the specified date.");
588 0 : return failure(StatusCode::OutOfRange, "DIVI data does not contain the specified date.");
589 : }
590 :
591 5859 : for (auto&& entry : divi_data) {
592 8280 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
593 8280 : return r == 0 || r == get_region_id(entry);
594 : });
595 5796 : auto date_df = entry.date;
596 5796 : if (it != vregion.end() && date_df == date) {
597 63 : auto region_idx = size_t(it - vregion.begin());
598 63 : vnum_icu[region_idx] = entry.num_icu;
599 : }
600 : }
601 :
602 126 : return success();
603 : }
604 :
605 : /**
606 : * @brief Reads DIVI data from a file and computes the ICU data for specified regions and date.
607 : *
608 : * @tparam FP Floating point type (default: double).
609 : *
610 : * @param[in] path Path to the file containing DIVI data.
611 : * @param[in] vregion Vector of region IDs for which the data is computed.
612 : * @param[in] date Date for which the ICU data is computed.
613 : * @param[out] vnum_icu Output vector containing the number of ICU cases for each region.
614 : *
615 : * @return An IOResult indicating success or failure.
616 : */
617 : template <typename FP = double>
618 63 : IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& vregion, Date date,
619 : std::vector<FP>& vnum_icu)
620 : {
621 63 : BOOST_OUTCOME_TRY(auto&& divi_data, mio::read_divi_data(path));
622 63 : return compute_divi_data(divi_data, vregion, date, vnum_icu);
623 63 : }
624 :
625 : /**
626 : * @brief Sets ICU data from DIVI data into the a vector of models, distributed across age groups.
627 : *
628 : * This function reads DIVI data from a file, computes the number of individuals in critical condition (ICU)
629 : * for each region, and sets these values in the model. The ICU cases are distributed across age groups
630 : * using the transition probabilities from severe to critical.
631 : *
632 : * @tparam Model The type of the model used.
633 : * @tparam FP Floating point type (default: double).
634 : *
635 : * @param[in,out] model Vector of models, each representing a region, where the ICU population is updated.
636 : * @param[in] path Path to the file containing DIVI data.
637 : * @param[in] vregion Vector of region IDs for which the data is computed.
638 : * @param[in] date Date for which the ICU data is computed.
639 : * @param[in] scaling_factor_icu Scaling factor for reported ICU cases.
640 : *
641 : * @return An IOResult indicating success or failure.
642 : */
643 : template <class Model, typename FP = double>
644 63 : IOResult<void> set_divi_data(std::vector<Model>& model, const std::string& path, const std::vector<int>& vregion,
645 : Date date, FP scaling_factor_icu)
646 : {
647 126 : std::vector<FP> sum_mu_I_U(vregion.size(), 0);
648 126 : std::vector<std::vector<FP>> mu_I_U{model.size()};
649 126 : for (size_t region = 0; region < vregion.size(); region++) {
650 63 : auto num_groups = model[region].parameters.get_num_groups();
651 441 : for (auto i = AgeGroup(0); i < num_groups; i++) {
652 756 : sum_mu_I_U[region] += model[region].parameters.template get<CriticalPerSevere<FP>>()[i] *
653 378 : model[region].parameters.template get<SeverePerInfectedSymptoms<FP>>()[i];
654 756 : mu_I_U[region].push_back(model[region].parameters.template get<CriticalPerSevere<FP>>()[i] *
655 378 : model[region].parameters.template get<SeverePerInfectedSymptoms<FP>>()[i]);
656 : }
657 : }
658 126 : std::vector<FP> num_icu(model.size(), 0.0);
659 63 : BOOST_OUTCOME_TRY(read_divi_data(path, vregion, date, num_icu));
660 :
661 126 : for (size_t region = 0; region < vregion.size(); region++) {
662 63 : auto num_groups = model[region].parameters.get_num_groups();
663 441 : for (auto i = AgeGroup(0); i < num_groups; i++) {
664 378 : model[region].populations[{i, InfectionState::InfectedCriticalNaive}] =
665 378 : scaling_factor_icu * num_icu[region] * mu_I_U[region][(size_t)i] / sum_mu_I_U[region];
666 : }
667 : }
668 :
669 126 : return success();
670 63 : }
671 :
672 : /**
673 : * @brief Reads population data from census data.
674 : *
675 : * @param[in] path Path to the population data file.
676 : * @param[in] vregion Vector of keys representing the regions of interest.
677 : * @return An IOResult containing a vector of vectors, where each inner vector represents the population
678 : * distribution across age groups for a specific region, or an error if the function fails.
679 : * @see mio::read_population_data
680 : */
681 : IOResult<std::vector<std::vector<double>>> read_population_data(const std::string& path,
682 : const std::vector<int>& vregion);
683 :
684 : /**
685 : * @brief Reads population data from a vector of population data entries.
686 : *
687 : * @param[in] population_data Vector of population data entries.
688 : * @param[in] vregion Vector of keys representing the regions of interest.
689 : * @return An IOResult containing a vector of vectors, where each inner vector represents the population
690 : * distribution across age groups for a specific region, or an error if the function fails.
691 : * @see mio::read_population_data
692 : */
693 : IOResult<std::vector<std::vector<double>>> read_population_data(const std::vector<PopulationDataEntry>& population_data,
694 : const std::vector<int>& vregion);
695 :
696 : /**
697 : * @brief Sets the population data for the given models based on the provided population distribution and immunity levels.
698 : *
699 : * @tparam Model The type of the model used.
700 : * @tparam FP Floating point type (default: double).
701 : *
702 : * @param[in,out] model A vector of models for which population data will be set.
703 : * @param[in] num_population A 2D vector where each row represents the age group population distribution for a specific region.
704 : * @param[in] vregion A vector of region identifiers corresponding to the population data.
705 : * @param[in] immunity_population A 2D vector where each row represents the immunity distribution for a specific region
706 : * across different levels of immunity (e.g., naive, partial, improved immunity).
707 : *
708 : * @return An IOResult indicating success or failure.
709 : */
710 : template <class Model, typename FP = double>
711 63 : IOResult<void> set_population_data(std::vector<Model>& model, const std::vector<std::vector<FP>>& num_population,
712 : const std::vector<int>& vregion,
713 : const std::vector<std::vector<FP>> immunity_population)
714 : {
715 126 : for (size_t region = 0; region < vregion.size(); region++) {
716 63 : if (std::accumulate(num_population[region].begin(), num_population[region].end(), 0.0) > 0) {
717 63 : auto num_groups = model[region].parameters.get_num_groups();
718 441 : for (auto i = AgeGroup(0); i < num_groups; i++) {
719 :
720 378 : FP SN = num_population[region][size_t(i)] * immunity_population[0][size_t(i)];
721 378 : FP SPI = num_population[region][size_t(i)] * immunity_population[1][size_t(i)];
722 378 : FP SII = num_population[region][size_t(i)] - SN - SPI;
723 :
724 378 : model[region].populations[{i, InfectionState::SusceptibleImprovedImmunity}] = std::max(
725 756 : 0.0,
726 756 : FP(SII -
727 1134 : (model[region].populations[{i, InfectionState::ExposedImprovedImmunity}] +
728 1134 : model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunity}] +
729 1134 : model[region].populations[{i, InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] +
730 1134 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunity}] +
731 1134 : model[region].populations[{i, InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] +
732 1134 : model[region].populations[{i, InfectionState::InfectedSevereImprovedImmunity}] +
733 1134 : model[region].populations[{i, InfectionState::InfectedCriticalImprovedImmunity}] +
734 1134 : model[region].populations[{i, InfectionState::DeadImprovedImmunity}] +
735 756 : model[region].populations[{i, InfectionState::TemporaryImmuneImprovedImmunity}])));
736 :
737 378 : model[region].populations[{i, InfectionState::SusceptiblePartialImmunity}] = std::max(
738 756 : 0.0,
739 1134 : SPI - model[region].populations[{i, InfectionState::ExposedPartialImmunity}] -
740 1134 : model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunity}] -
741 1134 : model[region].populations[{i, InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] -
742 1134 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunity}] -
743 1134 : model[region].populations[{i, InfectionState::InfectedSymptomsPartialImmunityConfirmed}] -
744 1134 : model[region].populations[{i, InfectionState::InfectedSeverePartialImmunity}] -
745 1134 : model[region].populations[{i, InfectionState::InfectedCriticalPartialImmunity}] -
746 1890 : model[region].populations[{i, InfectionState::DeadPartialImmunity}] -
747 756 : model[region].populations[{i, InfectionState::TemporaryImmunePartialImmunity}]);
748 :
749 756 : model[region].populations.template set_difference_from_group_total<AgeGroup>(
750 378 : {i, InfectionState::SusceptibleNaive}, num_population[region][size_t(i)]);
751 : }
752 :
753 441 : for (auto i = AgeGroup(0); i < AgeGroup(6); i++) {
754 11340 : for (auto j = Index<InfectionState>(0); j < InfectionState::Count; ++j) {
755 10962 : if (model[region].populations[{i, j}] < 0) {
756 0 : log_warning("Compartment at age group {}, infection state {}, is negative: {}", size_t(i),
757 0 : size_t(j), model[region].populations[{i, j}]);
758 : }
759 : }
760 : }
761 : }
762 : else {
763 0 : log_warning("No population data available for region " + std::to_string(region) +
764 : ". Population data has not been set.");
765 : }
766 : }
767 :
768 126 : return success();
769 : }
770 :
771 : /**
772 : * @brief Reads population data from a file and sets it for the each given model.
773 : *
774 : * @tparam Model The type of the model used.
775 : *
776 : * @param[in,out] model A vector of models for which population data will be set.
777 : * @param[in] path The file path to the population data.
778 : * @param[in] vregion A vector of region identifiers corresponding to the population data.
779 : * @param[in] immunity_population A 2D vector where each row represents the immunity distribution for a specific region
780 : * across different levels of immunity (e.g., naive, partial, improved).
781 : *
782 : * @return An IOResult indicating success or failure.
783 : */
784 : template <class Model>
785 36 : IOResult<void> set_population_data(std::vector<Model>& model, const std::string& path, const std::vector<int>& vregion,
786 : const std::vector<std::vector<double>> immunity_population)
787 : {
788 36 : BOOST_OUTCOME_TRY(auto&& num_population, details::read_population_data(path, vregion));
789 72 : BOOST_OUTCOME_TRY(set_population_data(model, num_population, vregion, immunity_population));
790 72 : return success();
791 36 : }
792 :
793 : /**
794 : * @brief Sets vaccination data for the given models using provided vaccination (partial, full, and booster) data.
795 : *
796 : *
797 : * @tparam FP Floating point type (default: double).
798 : *
799 : * @param[in,out] model A vector of models for which vaccination data will be set.
800 : * @param[in] vacc_data A vector of VaccinationDataEntry objects containing the vaccination data.
801 : * @param[in] date The starting date for the simulation.
802 : * @param[in] vregion A vector of region identifiers corresponding to the vaccination data.
803 : * @param[in] num_days The number of days for which the simulation runs.
804 : *
805 : * @return An IOResult indicating success or failure.
806 : */
807 : template <typename FP = double>
808 64 : IOResult<void> set_vaccination_data(std::vector<Model<FP>>& model, const std::vector<VaccinationDataEntry>& vacc_data,
809 : Date date, const std::vector<int>& vregion, int num_days)
810 : {
811 64 : auto num_groups = model[0].parameters.get_num_groups();
812 :
813 64 : auto days_until_effective_n =
814 64 : (int)(double)model[0].parameters.template get<DaysUntilEffectivePartialVaccination<FP>>()[AgeGroup(0)];
815 64 : auto days_until_effective_pi =
816 64 : (int)(double)model[0].parameters.template get<DaysUntilEffectiveImprovedVaccination<FP>>()[AgeGroup(0)];
817 64 : auto days_until_effective_ii =
818 64 : (int)(double)model[0].parameters.template get<DaysUntilEffectiveBoosterImmunity<FP>>()[AgeGroup(0)];
819 : // iterate over regions (e.g., counties)
820 128 : for (size_t i = 0; i < model.size(); ++i) {
821 : // iterate over age groups in region
822 443 : for (auto g = AgeGroup(0); g < num_groups; ++g) {
823 :
824 379 : model[i].parameters.template get<DailyPartialVaccinations<FP>>().resize(SimulationDay(num_days + 1));
825 379 : model[i].parameters.template get<DailyFullVaccinations<FP>>().resize(SimulationDay(num_days + 1));
826 379 : model[i].parameters.template get<DailyBoosterVaccinations<FP>>().resize(SimulationDay(num_days + 1));
827 2819 : for (auto d = SimulationDay(0); d < SimulationDay(num_days + 1); ++d) {
828 2440 : model[i].parameters.template get<DailyPartialVaccinations<FP>>()[{g, d}] = 0.0;
829 2440 : model[i].parameters.template get<DailyFullVaccinations<FP>>()[{g, d}] = 0.0;
830 2440 : model[i].parameters.template get<DailyBoosterVaccinations<FP>>()[{g, d}] = 0.0;
831 : }
832 : }
833 : }
834 :
835 64 : auto max_date_entry = std::max_element(vacc_data.begin(), vacc_data.end(), [](auto&& a, auto&& b) {
836 32815 : return a.date < b.date;
837 : });
838 64 : if (max_date_entry == vacc_data.end()) {
839 0 : return failure(StatusCode::InvalidFileFormat, "Vaccination data file is empty.");
840 : }
841 64 : auto max_date = max_date_entry->date;
842 64 : if (max_date < offset_date_by_days(date, num_days)) {
843 1 : log_error("Vaccination data does not contain all dates. After the last day the data, vaccinations per day are "
844 : "set to 0.");
845 : }
846 :
847 32943 : for (auto&& vacc_data_entry : vacc_data) {
848 104779 : auto it = std::find_if(vregion.begin(), vregion.end(), [&vacc_data_entry](auto&& r) {
849 100310 : return r == 0 || (vacc_data_entry.county_id && vacc_data_entry.county_id == regions::CountyId(r)) ||
850 92715 : (vacc_data_entry.state_id && vacc_data_entry.state_id == regions::StateId(r)) ||
851 88213 : (vacc_data_entry.district_id && vacc_data_entry.district_id == regions::DistrictId(r));
852 : });
853 32879 : auto date_df = vacc_data_entry.date;
854 32879 : if (it != vregion.end()) {
855 32879 : auto region_idx = size_t(it - vregion.begin());
856 32879 : AgeGroup age = vacc_data_entry.age_group;
857 :
858 : // get daily vaccinations for each layer
859 239377 : for (size_t d = 0; d < (size_t)num_days + 1; ++d) {
860 206498 : auto offset_first_date = offset_date_by_days(date, (int)d - days_until_effective_n);
861 206498 : if (max_date >= offset_first_date) {
862 206487 : if (date_df == offset_first_date) {
863 2250 : model[region_idx]
864 2250 : .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] =
865 2250 : vacc_data_entry.num_vaccinations_partial;
866 : }
867 : }
868 : else {
869 11 : if (date_df == offset_first_date) {
870 0 : model[region_idx]
871 0 : .parameters.template get<DailyPartialVaccinations<FP>>()[{age, SimulationDay(d)}] = 0;
872 : }
873 : }
874 :
875 206498 : auto offset_full_date = offset_date_by_days(date, (int)d - days_until_effective_pi);
876 206498 : if (max_date >= offset_full_date) {
877 206498 : if (date_df == offset_full_date) {
878 2440 : model[region_idx]
879 2440 : .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] =
880 2440 : vacc_data_entry.num_vaccinations_completed;
881 : }
882 : }
883 : else {
884 0 : if (date_df == offset_first_date) {
885 0 : model[region_idx]
886 0 : .parameters.template get<DailyFullVaccinations<FP>>()[{age, SimulationDay(d)}] = 0;
887 : }
888 : }
889 :
890 206498 : auto offset_booster_date = offset_date_by_days(date, (int)d - days_until_effective_ii);
891 206498 : if (max_date >= offset_booster_date) {
892 206487 : if (date_df == offset_booster_date) {
893 2439 : model[region_idx]
894 2439 : .parameters.template get<DailyBoosterVaccinations<FP>>()[{age, SimulationDay(d)}] =
895 2439 : vacc_data_entry.num_vaccinations_refreshed_first +
896 2439 : vacc_data_entry.num_vaccinations_refreshed_additional;
897 : }
898 : }
899 : else {
900 11 : if (date_df == offset_first_date) {
901 0 : model[region_idx]
902 0 : .parameters.template get<DailyBoosterVaccinations<FP>>()[{age, SimulationDay(d)}] = 0;
903 : }
904 : }
905 : }
906 : }
907 : }
908 128 : return success();
909 : }
910 :
911 : /**
912 : * @brief Sets vaccination data for the given models using vaccination data from a file.
913 : *
914 : * This function reads vaccination data from a specified file, and assigns daily vaccination numbers
915 : * (partial, full, and booster) to each region and age group in the models.
916 : *
917 : * @tparam FP Floating point type (default: double).
918 : *
919 : * @param[in,out] model A vector of models for which vaccination data will be set.
920 : * @param[in] path The file path to the vaccination data.
921 : * @param[in] date The starting date for the simulation.
922 : * @param[in] vregion A vector of region identifiers corresponding to the vaccination data.
923 : * @param[in] num_days The number of days for which the simulation runs.
924 : *
925 : * @return An IOResult indicating success or failure.
926 : */
927 : template <typename FP = double>
928 37 : IOResult<void> set_vaccination_data(std::vector<Model<FP>>& model, const std::string& path, Date date,
929 : const std::vector<int>& vregion, int num_days)
930 : {
931 37 : BOOST_OUTCOME_TRY(auto&& vacc_data, read_vaccination_data(path));
932 37 : BOOST_OUTCOME_TRY(set_vaccination_data(model, vacc_data, date, vregion, num_days));
933 74 : return success();
934 37 : }
935 :
936 : } // namespace details
937 :
938 : #ifdef MEMILIO_HAS_HDF5
939 :
940 : /**
941 : * @brief Exports extrapolated real-world data time series for specified regions.
942 : *
943 : * This function generates and exports time series data based on the extrapolation and approximation methods
944 : * used to initialize the model with real-world data. The resulting data represents the initialized states of
945 : * the model over the specified time range.
946 : *
947 : * @tparam Model Type of the model used.
948 : *
949 : * @param[in] models A vector of models for which the extrapolated data is set.
950 : * @param[in] results_dir Path to the directory where the extrapolated results will be saved in a h5 file.
951 : * @param[in] counties A vector of region identifiers for which the time series will be exported.
952 : * @param[in] date The starting date of the time series.
953 : * @param[in] scaling_factor_inf A vector of scaling factors applied to confirmed cases.
954 : * @param[in] scaling_factor_icu A scaling factor applied to ICU cases.
955 : * @param[in] num_days The number of days for which will be extrapolated.
956 : * @param[in] divi_data_path Path to the DIVI ICU data file.
957 : * @param[in] confirmed_cases_path Path to the confirmed cases data file.
958 : * @param[in] population_data_path Path to the population data file.
959 : * @param[in] immunity_population A vector of vectors specifying immunity for each age group and immunity layer.
960 : * @param[in] vaccination_data_path Path to the vaccination data file (optional).
961 : *
962 : * @return An IOResult indicating success or failure.
963 : */
964 : template <class Model>
965 9 : IOResult<void> export_input_data_county_timeseries(
966 : std::vector<Model> models, const std::string& results_dir, const std::vector<int>& counties, Date date,
967 : const std::vector<double>& scaling_factor_inf, const double scaling_factor_icu, const int num_days,
968 : const std::string& divi_data_path, const std::string& confirmed_cases_path, const std::string& population_data_path,
969 : const std::vector<std::vector<double>> immunity_population, const std::string& vaccination_data_path = "")
970 : {
971 9 : const auto num_groups = (size_t)models[0].parameters.get_num_groups();
972 9 : assert(scaling_factor_inf.size() == num_groups);
973 9 : assert(num_groups == ConfirmedCasesDataEntry::age_group_names.size());
974 9 : assert(models.size() == counties.size());
975 36 : std::vector<TimeSeries<double>> extrapolated_data(
976 9 : models.size(), TimeSeries<double>::zero(num_days + 1, (size_t)InfectionState::Count * num_groups));
977 :
978 9 : BOOST_OUTCOME_TRY(auto&& case_data, read_confirmed_cases_data(confirmed_cases_path));
979 9 : BOOST_OUTCOME_TRY(auto&& population_data, details::read_population_data(population_data_path, counties));
980 :
981 : // empty vector if set_vaccination_data is not set
982 9 : std::vector<VaccinationDataEntry> vacc_data;
983 9 : if (!vaccination_data_path.empty()) {
984 9 : BOOST_OUTCOME_TRY(vacc_data, read_vaccination_data(vaccination_data_path));
985 9 : }
986 :
987 63 : for (int t = 0; t <= num_days; ++t) {
988 27 : auto offset_day = offset_date_by_days(date, t);
989 :
990 27 : if (!vaccination_data_path.empty()) {
991 27 : BOOST_OUTCOME_TRY(details::set_vaccination_data(models, vacc_data, offset_day, counties, num_days));
992 27 : }
993 :
994 : // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
995 : // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
996 27 : if (offset_day > Date(2020, 4, 23)) {
997 27 : BOOST_OUTCOME_TRY(details::set_divi_data(models, divi_data_path, counties, offset_day, scaling_factor_icu));
998 27 : }
999 : else {
1000 0 : log_warning("No DIVI data available for for date: {}-{}-{}", offset_day.day, offset_day.month,
1001 : offset_day.year);
1002 : }
1003 :
1004 54 : BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(models, case_data, counties, offset_day, scaling_factor_inf,
1005 : immunity_population));
1006 :
1007 54 : BOOST_OUTCOME_TRY(details::set_population_data(models, population_data, counties, immunity_population));
1008 :
1009 54 : for (size_t r = 0; r < counties.size(); r++) {
1010 27 : extrapolated_data[r][t] = models[r].get_initial_values();
1011 : // in set_population_data the number of death individuals is subtracted from the SusceptibleImprovedImmunity compartment.
1012 : // Since we should be independent whether we consider them or not, we add them back here before we save the data.
1013 189 : for (size_t age = 0; age < num_groups; age++) {
1014 162 : extrapolated_data[r][t][(size_t)InfectionState::SusceptibleImprovedImmunity +
1015 324 : age * (size_t)InfectionState::Count] +=
1016 162 : extrapolated_data[r][t][(size_t)InfectionState::DeadNaive + age * (size_t)InfectionState::Count];
1017 : }
1018 : }
1019 : }
1020 18 : BOOST_OUTCOME_TRY(save_result(extrapolated_data, counties, static_cast<int>(num_groups),
1021 : path_join(results_dir, "Results_rki.h5")));
1022 :
1023 45 : auto extrapolated_rki_data_sum = sum_nodes(std::vector<std::vector<TimeSeries<double>>>{extrapolated_data});
1024 72 : BOOST_OUTCOME_TRY(save_result({extrapolated_rki_data_sum[0][0]}, {0}, static_cast<int>(num_groups),
1025 : path_join(results_dir, "Results_rki_sum.h5")));
1026 :
1027 18 : return success();
1028 9 : }
1029 :
1030 : #else
1031 : template <class Model>
1032 : IOResult<void> export_input_data_county_timeseries(std::vector<Model>, const std::string&, const std::vector<int>&,
1033 : Date, const std::vector<double>&, const double, const int,
1034 : const std::string&, const std::string&, const std::string&,
1035 : const std::vector<std::vector<double>>, const std::string&)
1036 : {
1037 : mio::log_warning("HDF5 not available. Cannot export time series of extrapolated real data.");
1038 : return success();
1039 : }
1040 :
1041 : #endif //MEMILIO_HAS_HDF5
1042 :
1043 : /**
1044 : * @brief Reads input data for specified counties and initializes the model accordingly.
1045 : *
1046 : * This function loads real-world data for specified counties and initializes
1047 : * the corresponding model compartments. Optionally, it can extrapolate real-world data
1048 : * using the export_input_data_county_timeseries function.
1049 : *
1050 : * @tparam Model The model type used.
1051 : *
1052 : * @param[in,out] model Vector of model objects that will be initialized with the input data.
1053 : * @param[in] date The starting date for the simulation.
1054 : * @param[in] county Vector of county IDs for which the data is read.
1055 : * @param[in] scaling_factor_inf Vector of scaling factors for confirmed cases.
1056 : * @param[in] scaling_factor_icu Scaling factor for ICU data.
1057 : * @param[in] dir Path to the directory containing input data files.
1058 : * @param[in] num_days The number of days for which the simulation runs.
1059 : * @param[in] immunity_population Vector of vectors representing immunity proportions for each age group and immunity layer.
1060 : * @param[in] export_time_series Boolean flag indicating whether to export time series of extrapolated data (default: false).
1061 : *
1062 : * @return An IOResult indicating success or failure.
1063 : */
1064 : template <class Model>
1065 18 : IOResult<void> read_input_data_county(std::vector<Model>& model, Date date, const std::vector<int>& county,
1066 : const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
1067 : const std::string& dir, int num_days,
1068 : const std::vector<std::vector<double>> immunity_population,
1069 : bool export_time_series = false)
1070 : {
1071 36 : BOOST_OUTCOME_TRY(details::set_vaccination_data(
1072 : model, path_join(dir, "pydata/Germany", "all_county_ageinf_vacc_ma7.json"), date, county, num_days));
1073 :
1074 : // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
1075 : // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
1076 18 : if (date > Date(2020, 4, 23)) {
1077 36 : BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(dir, "pydata/Germany", "county_divi_ma7.json"),
1078 : county, date, scaling_factor_icu));
1079 18 : }
1080 : else {
1081 0 : log_warning("No DIVI data available for this date");
1082 : }
1083 :
1084 54 : BOOST_OUTCOME_TRY(
1085 : details::set_confirmed_cases_data(model, path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"),
1086 : county, date, scaling_factor_inf, immunity_population));
1087 54 : BOOST_OUTCOME_TRY(details::set_population_data(
1088 : model, path_join(dir, "pydata/Germany", "county_current_population.json"), county, immunity_population));
1089 :
1090 18 : if (export_time_series) {
1091 : // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
1092 : // Run time equals run time of the previous functions times the num_days !
1093 : // (This only represents the vectorization of the previous function over all simulation days...)
1094 0 : log_info("Exporting time series of extrapolated real data. This may take some minutes. "
1095 : "For simulation runs over the same time period, deactivate it.");
1096 0 : BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
1097 : model, dir, county, date, scaling_factor_inf, scaling_factor_icu, num_days,
1098 : path_join(dir, "pydata/Germany", "county_divi_ma7.json"),
1099 : path_join(dir, "pydata/Germany", "cases_all_county_age_ma7.json"),
1100 : path_join(dir, "pydata/Germany", "county_current_population.json"), immunity_population,
1101 : path_join(dir, "pydata/Germany", "all_county_ageinf_vacc_ma7.json")));
1102 0 : }
1103 :
1104 36 : return success();
1105 18 : }
1106 :
1107 : /**
1108 : * @brief Reads compartments for geographic units at a specified date from data files.
1109 : *
1110 : * This function estimates all compartments from available data using the provided model parameters.
1111 : *
1112 : * @tparam Model The type of tmodel used.
1113 : *
1114 : * @param[in,out] model Vector of models, one per county, to be initialized with data.
1115 : * @param[in] date Date for which the data should be read.
1116 : * @param[in] node_ids Vector of IDs of the units for which data is read.
1117 : * @param[in] scaling_factor_inf Vector of scaling factors for confirmed cases.
1118 : * @param[in] scaling_factor_icu Scaling factor for ICU cases.
1119 : * @param[in] data_dir Directory containing the input data files.
1120 : * @param[in] num_days Number of days to simulate.
1121 : * @param[in] immunity_population Matrix containing immunity proportions for each age group and immunity layer.
1122 : * @param[in] export_time_series Boolean flag indicating whether to export time series of extrapolated data (default: false).
1123 : *
1124 : * @return An IOResult indicating success or failure.
1125 : */
1126 : template <class Model>
1127 18 : IOResult<void> read_input_data(std::vector<Model>& model, Date date, const std::vector<int>& node_ids,
1128 : const std::vector<double>& scaling_factor_inf, double scaling_factor_icu,
1129 : const std::string& data_dir, int num_days,
1130 : const std::vector<std::vector<double>> immunity_population,
1131 : bool export_time_series = false)
1132 : {
1133 :
1134 36 : BOOST_OUTCOME_TRY(
1135 : details::set_vaccination_data(model, path_join(data_dir, "vaccination_data.json"), date, node_ids, num_days));
1136 :
1137 : // TODO: Reuse more code, e.g., set_divi_data (in secir) and a set_divi_data (here) only need a different ModelType.
1138 : // TODO: add option to set ICU data from confirmed cases if DIVI or other data is not available.
1139 18 : if (date > Date(2020, 4, 23)) {
1140 36 : BOOST_OUTCOME_TRY(details::set_divi_data(model, path_join(data_dir, "critical_cases.json"), node_ids, date,
1141 : scaling_factor_icu));
1142 18 : }
1143 : else {
1144 0 : log_warning("No DIVI data available for this date");
1145 : }
1146 :
1147 54 : BOOST_OUTCOME_TRY(details::set_confirmed_cases_data(model, path_join(data_dir, "confirmed_cases.json"), node_ids,
1148 : date, scaling_factor_inf, immunity_population));
1149 54 : BOOST_OUTCOME_TRY(details::set_population_data(model, path_join(data_dir, "population_data.json"), node_ids,
1150 : immunity_population));
1151 :
1152 18 : if (export_time_series) {
1153 : // Use only if extrapolated real data is needed for comparison. EXPENSIVE !
1154 : // Run time equals run time of the previous functions times the num_days !
1155 : // (This only represents the vectorization of the previous function over all simulation days...)
1156 0 : log_info("Exporting time series of extrapolated real data. This may take some minutes. "
1157 : "For simulation runs over the same time period, deactivate it.");
1158 0 : BOOST_OUTCOME_TRY(export_input_data_county_timeseries(
1159 : model, data_dir, node_ids, date, scaling_factor_inf, scaling_factor_icu, num_days,
1160 : path_join(data_dir, "critical_cases.json"), path_join(data_dir, "confirmed_cases.json"),
1161 : path_join(data_dir, "population_data.json"), immunity_population,
1162 : path_join(data_dir, "vaccination_data.json")));
1163 0 : }
1164 :
1165 36 : return success();
1166 18 : }
1167 :
1168 : } // namespace osecirts
1169 : } // namespace mio
1170 :
1171 : #endif // MEMILIO_HAS_JSONCPP
1172 :
1173 : #endif // MIO_ODE_SECIRTS_PARAMETERS_IO_H
|
