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