Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Wadim Koslow, Daniel Abele, Martin J. Kühn
5 : *
6 : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
7 : *
8 : * Licensed under the Apache License, Version 2.0 (the "License");
9 : * you may not use this file except in compliance with the License.
10 : * You may obtain a copy of the License at
11 : *
12 : * http://www.apache.org/licenses/LICENSE-2.0
13 : *
14 : * Unless required by applicable law or agreed to in writing, software
15 : * distributed under the License is distributed on an "AS IS" BASIS,
16 : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 : * See the License for the specific language governing permissions and
18 : * limitations under the License.
19 : */
20 : #include "memilio/config.h"
21 :
22 : #ifdef MEMILIO_HAS_JSONCPP
23 :
24 : #include "ode_secirvvs/parameters_io.h"
25 :
26 : namespace mio
27 : {
28 : namespace osecirvvs
29 : {
30 : namespace details
31 : {
32 : //gets the county or state id of the entry if available, 0 (for whole country) otherwise
33 : //used for comparisons of entry to integer region id
34 : template <class EpiDataEntry>
35 86940 : int get_region_id(const EpiDataEntry& rki_entry)
36 : {
37 194580 : return rki_entry.county_id ? rki_entry.county_id->get()
38 62100 : : (rki_entry.state_id ? rki_entry.state_id->get()
39 128340 : : (rki_entry.district_id ? rki_entry.district_id->get() : 0));
40 : }
41 :
42 9 : IOResult<void> read_confirmed_cases_data(
43 : std::string const& path, std::vector<int> const& vregion, Date date, std::vector<std::vector<double>>& vnum_Exposed,
44 : std::vector<std::vector<double>>& vnum_InfectedNoSymptoms, std::vector<std::vector<double>>& vnum_InfectedSymptoms,
45 : std::vector<std::vector<double>>& vnum_InfectedSevere, std::vector<std::vector<double>>& vnum_icu,
46 : std::vector<std::vector<double>>& vnum_death, std::vector<std::vector<double>>& vnum_rec,
47 : const std::vector<std::vector<int>>& vt_Exposed, const std::vector<std::vector<int>>& vt_InfectedNoSymptoms,
48 : const std::vector<std::vector<int>>& vt_InfectedSymptoms, const std::vector<std::vector<int>>& vt_InfectedSevere,
49 : const std::vector<std::vector<int>>& vt_InfectedCritical, const std::vector<std::vector<double>>& vmu_C_R,
50 : const std::vector<std::vector<double>>& vmu_I_H, const std::vector<std::vector<double>>& vmu_H_U,
51 : const std::vector<double>& scaling_factor_inf)
52 : {
53 9 : BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path));
54 : return read_confirmed_cases_data(rki_data, vregion, date, vnum_Exposed, vnum_InfectedNoSymptoms,
55 : vnum_InfectedSymptoms, vnum_InfectedSevere, vnum_icu, vnum_death, vnum_rec,
56 : vt_Exposed, vt_InfectedNoSymptoms, vt_InfectedSymptoms, vt_InfectedSevere,
57 9 : vt_InfectedCritical, vmu_C_R, vmu_I_H, vmu_H_U, scaling_factor_inf);
58 9 : }
59 :
60 306 : IOResult<void> read_confirmed_cases_data(
61 : const std::vector<ConfirmedCasesDataEntry>& rki_data, std::vector<int> const& vregion, Date date,
62 : std::vector<std::vector<double>>& vnum_Exposed, std::vector<std::vector<double>>& vnum_InfectedNoSymptoms,
63 : std::vector<std::vector<double>>& vnum_InfectedSymptoms, std::vector<std::vector<double>>& vnum_InfectedSevere,
64 : std::vector<std::vector<double>>& vnum_icu, std::vector<std::vector<double>>& vnum_death,
65 : std::vector<std::vector<double>>& vnum_rec, const std::vector<std::vector<int>>& vt_Exposed,
66 : const std::vector<std::vector<int>>& vt_InfectedNoSymptoms,
67 : const std::vector<std::vector<int>>& vt_InfectedSymptoms, const std::vector<std::vector<int>>& vt_InfectedSevere,
68 : const std::vector<std::vector<int>>& vt_InfectedCritical, const std::vector<std::vector<double>>& vmu_C_R,
69 : const std::vector<std::vector<double>>& vmu_I_H, const std::vector<std::vector<double>>& vmu_H_U,
70 : const std::vector<double>& scaling_factor_inf)
71 : {
72 306 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
73 168606 : return a.date < b.date;
74 : });
75 306 : if (max_date_entry == rki_data.end()) {
76 0 : log_error("RKI data file is empty.");
77 0 : return failure(StatusCode::InvalidValue, "RKI data is empty.");
78 : }
79 306 : auto max_date = max_date_entry->date;
80 306 : if (max_date < date) {
81 0 : log_error("Specified date does not exist in RKI data");
82 0 : return failure(StatusCode::OutOfRange, "RKI data does not contain specified date.");
83 : }
84 :
85 : // shifts the initilization to the recent past if simulation starts
86 : // around current day and data of the future would be required.
87 : // Only needed for preinfection compartments, exposed and InfectedNoSymptoms.
88 306 : auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E + T_C
89 306 : if (days_surplus > 0) {
90 306 : days_surplus = 0;
91 : }
92 :
93 169218 : for (auto&& entry : rki_data) {
94 233496 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
95 233496 : return r == 0 || get_region_id(entry) == r;
96 : });
97 168912 : if (it != vregion.end()) {
98 168912 : auto region_idx = size_t(it - vregion.begin());
99 :
100 168912 : auto& t_Exposed = vt_Exposed[region_idx];
101 168912 : auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx];
102 168912 : auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx];
103 168912 : auto& t_InfectedSevere = vt_InfectedSevere[region_idx];
104 168912 : auto& t_InfectedCritical = vt_InfectedCritical[region_idx];
105 :
106 168912 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
107 168912 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
108 168912 : auto& num_rec = vnum_rec[region_idx];
109 168912 : auto& num_Exposed = vnum_Exposed[region_idx];
110 168912 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
111 168912 : auto& num_death = vnum_death[region_idx];
112 168912 : auto& num_icu = vnum_icu[region_idx];
113 :
114 168912 : auto& mu_C_R = vmu_C_R[region_idx];
115 168912 : auto& mu_I_H = vmu_I_H[region_idx];
116 168912 : auto& mu_H_U = vmu_H_U[region_idx];
117 :
118 168912 : auto age = (size_t)entry.age_group;
119 168912 : if (entry.date == offset_date_by_days(date, 0)) {
120 1215 : num_InfectedSymptoms[age] += scaling_factor_inf[age] * entry.num_confirmed;
121 1215 : num_rec[age] += entry.num_confirmed;
122 : }
123 168912 : if (entry.date == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) {
124 1215 : num_InfectedNoSymptoms[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed;
125 1215 : num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed;
126 : }
127 168912 : if (entry.date == offset_date_by_days(date, days_surplus)) {
128 1215 : num_InfectedNoSymptoms[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed;
129 : }
130 168912 : if (entry.date == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) {
131 1215 : num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * entry.num_confirmed;
132 : }
133 168912 : if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age])) {
134 1215 : num_InfectedSymptoms[age] -= scaling_factor_inf[age] * entry.num_confirmed;
135 1215 : num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed;
136 : }
137 168912 : if (entry.date == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) {
138 1215 : num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * entry.num_confirmed;
139 1215 : num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed;
140 : }
141 337824 : if (entry.date ==
142 337824 : offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) {
143 1215 : num_death[age] += entry.num_deaths;
144 1215 : num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * entry.num_confirmed;
145 : }
146 : }
147 : }
148 :
149 612 : for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) {
150 306 : auto region = vregion[region_idx];
151 :
152 306 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
153 306 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
154 306 : auto& num_rec = vnum_rec[region_idx];
155 306 : auto& num_Exposed = vnum_Exposed[region_idx];
156 306 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
157 306 : auto& num_death = vnum_death[region_idx];
158 306 : auto& num_icu = vnum_icu[region_idx];
159 :
160 306 : size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size();
161 2142 : for (size_t i = 0; i < num_groups; i++) {
162 : // subtract infected confirmed cases which are not yet recovered
163 : // and remove dark figure scaling factor
164 1836 : num_rec[i] -= num_InfectedSymptoms[i] / scaling_factor_inf[i];
165 1836 : num_rec[i] -= num_InfectedSevere[i] / scaling_factor_inf[i];
166 1836 : num_rec[i] -=
167 3672 : num_icu[i] /
168 1836 : scaling_factor_inf[i]; // TODO: this has to be adapted for scaling_factor_inf != 1 or != ***_icu
169 1836 : num_rec[i] -= num_death[i] / scaling_factor_inf[i];
170 1863 : auto try_fix_constraints = [region, i](double& value, double error, auto str) {
171 12852 : if (value < error) {
172 : // this should probably return a failure
173 : // but the algorithm is not robust enough to avoid large negative
174 : // values and there are tests that rely on it
175 27 : log_error("{:s} for age group {:s} is {:.4f} for region {:d}, "
176 : "exceeds expected negative value.",
177 27 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
178 27 : value = 0.0;
179 : }
180 12825 : else if (value < 0) {
181 0 : log_info("{:s} for age group {:s} is {:.4f} for region {:d}, "
182 : "automatically corrected",
183 0 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
184 0 : value = 0.0;
185 : }
186 12852 : };
187 :
188 1836 : try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms");
189 1836 : try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms");
190 1836 : try_fix_constraints(num_Exposed[i], -5, "Exposed");
191 1836 : try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere");
192 1836 : try_fix_constraints(num_death[i], -5, "Dead");
193 1836 : try_fix_constraints(num_icu[i], -5, "InfectedCritical");
194 1836 : try_fix_constraints(num_rec[i], -20, "Recovered or vaccinated");
195 : }
196 : }
197 :
198 612 : return success();
199 : }
200 :
201 0 : IOResult<void> read_confirmed_cases_data_fix_recovered(std::string const& path, std::vector<int> const& vregion,
202 : Date date, std::vector<std::vector<double>>& vnum_rec,
203 : double delay)
204 : {
205 0 : BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path));
206 0 : return read_confirmed_cases_data_fix_recovered(rki_data, vregion, date, vnum_rec, delay);
207 0 : }
208 :
209 108 : IOResult<void> read_confirmed_cases_data_fix_recovered(const std::vector<ConfirmedCasesDataEntry>& rki_data,
210 : std::vector<int> const& vregion, Date date,
211 : std::vector<std::vector<double>>& vnum_rec, double delay)
212 : {
213 108 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
214 59508 : return a.date < b.date;
215 : });
216 108 : if (max_date_entry == rki_data.end()) {
217 0 : log_error("RKI data is empty.");
218 0 : return failure(StatusCode::InvalidValue, "RKI data is empty.");
219 : }
220 108 : auto max_date = max_date_entry->date;
221 108 : if (max_date < date) {
222 0 : log_error("Specified date does not exist in RKI data");
223 0 : return failure(StatusCode::OutOfRange, "RKI data does not contain specified date.");
224 : }
225 :
226 : // shifts the initilization to the recent past if simulation starts
227 : // around current day and data of the future would be required.
228 : // Only needed for preinfection compartments, exposed and InfectedNoSymptoms.
229 108 : auto days_surplus = get_offset_in_days(max_date, date) - 6; // 6 > T_E^C + T_C^I
230 108 : if (days_surplus > 0) {
231 108 : days_surplus = 0;
232 : }
233 :
234 59724 : for (auto&& rki_entry : rki_data) {
235 79488 : auto it = std::find_if(vregion.begin(), vregion.end(), [&rki_entry](auto r) {
236 79488 : return r == 0 || get_region_id(rki_entry) == r;
237 : });
238 59616 : if (it != vregion.end()) {
239 54648 : auto region_idx = size_t(it - vregion.begin());
240 54648 : if (rki_entry.date == offset_date_by_days(date, int(-delay))) {
241 378 : vnum_rec[region_idx][size_t(rki_entry.age_group)] = rki_entry.num_confirmed;
242 : }
243 : }
244 : }
245 :
246 216 : for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) {
247 108 : auto region = vregion[region_idx];
248 108 : auto& num_rec = vnum_rec[region_idx];
249 :
250 108 : size_t num_groups = ConfirmedCasesDataEntry::age_group_names.size();
251 756 : for (size_t i = 0; i < num_groups; i++) {
252 648 : auto try_fix_constraints = [region, i](double& value, double error, auto str) {
253 648 : if (value < error) {
254 : // this should probably return a failure
255 : // but the algorithm is not robust enough to avoid large negative
256 : // values and there are tests that rely on it
257 0 : log_error("{:s} for age group {:s} is {:.4f} for region {:d}, "
258 : "exceeds expected negative value.",
259 0 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
260 0 : value = 0.0;
261 : }
262 648 : else if (value < 0) {
263 0 : log_info("{:s} for age group {:s} is {:.4f} for region {:d}, "
264 : "automatically corrected",
265 0 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
266 0 : value = 0.0;
267 : }
268 648 : };
269 648 : try_fix_constraints(num_rec[i], 0, "Recovered");
270 : }
271 : }
272 :
273 216 : return success();
274 : }
275 :
276 63 : IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& vregion, Date date,
277 : std::vector<double>& vnum_icu)
278 : {
279 63 : BOOST_OUTCOME_TRY(auto&& divi_data, mio::read_divi_data(path));
280 63 : return read_divi_data(divi_data, vregion, date, vnum_icu);
281 63 : }
282 :
283 63 : IOResult<void> read_divi_data(const std::vector<DiviEntry>& divi_data, const std::vector<int>& vregion, Date date,
284 : std::vector<double>& vnum_icu)
285 : {
286 63 : auto max_date_entry = std::max_element(divi_data.begin(), divi_data.end(), [](auto&& a, auto&& b) {
287 5733 : return a.date < b.date;
288 : });
289 63 : if (max_date_entry == divi_data.end()) {
290 0 : log_error("DIVI data is empty.");
291 0 : return failure(StatusCode::InvalidValue, "DIVI data is empty.");
292 : }
293 63 : auto max_date = max_date_entry->date;
294 63 : if (max_date < date) {
295 0 : log_error("DIVI data does not contain the specified date.");
296 0 : return failure(StatusCode::OutOfRange, "DIVI data does not contain the specified date.");
297 : }
298 :
299 5859 : for (auto&& entry : divi_data) {
300 8280 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
301 8280 : return r == 0 || r == get_region_id(entry);
302 : });
303 5796 : auto date_df = entry.date;
304 5796 : if (it != vregion.end() && date_df == date) {
305 63 : auto region_idx = size_t(it - vregion.begin());
306 63 : vnum_icu[region_idx] = entry.num_icu;
307 : }
308 : }
309 :
310 126 : return success();
311 : }
312 :
313 135 : IOResult<std::vector<std::vector<double>>> read_population_data(const std::string& path,
314 : const std::vector<int>& vregion)
315 : {
316 135 : BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path));
317 135 : return read_population_data(population_data, vregion);
318 135 : }
319 :
320 135 : IOResult<std::vector<std::vector<double>>> read_population_data(const std::vector<PopulationDataEntry>& population_data,
321 : const std::vector<int>& vregion)
322 : {
323 135 : std::vector<std::vector<double>> vnum_population(
324 540 : vregion.size(), std::vector<double>(ConfirmedCasesDataEntry::age_group_names.size(), 0.0));
325 :
326 54270 : for (auto&& county_entry : population_data) {
327 : //accumulate population of states or country from population of counties
328 54135 : if (!county_entry.county_id && !county_entry.district_id) {
329 0 : return failure(StatusCode::InvalidFileFormat, "File with county population expected.");
330 : }
331 : //find region that this county belongs to
332 : //all counties belong to the country (id = 0)
333 191241 : auto it = std::find_if(vregion.begin(), vregion.end(), [&county_entry](auto r) {
334 14436 : return r == 0 ||
335 54135 : (county_entry.county_id &&
336 97443 : regions::StateId(r) == regions::get_state_id(int(*county_entry.county_id))) ||
337 187650 : (county_entry.county_id && regions::CountyId(r) == *county_entry.county_id) ||
338 158742 : (county_entry.district_id && regions::DistrictId(r) == *county_entry.district_id);
339 : });
340 54135 : if (it != vregion.end()) {
341 39726 : auto region_idx = size_t(it - vregion.begin());
342 39726 : auto& num_population = vnum_population[region_idx];
343 278082 : for (size_t age = 0; age < num_population.size(); age++) {
344 238356 : num_population[age] += county_entry.population[AgeGroup(age)];
345 : }
346 : }
347 : }
348 :
349 135 : return success(vnum_population);
350 135 : }
351 :
352 : } // namespace details
353 : } // namespace osecirvvs
354 : } // namespace mio
355 :
356 : #endif // MEMILIO_HAS_JSONCPP
|
