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