Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 MEmilio
3 : *
4 : * Authors: Wadim Koslow, Daniel Abele, Martin J. Kuehn, Lena Ploetzke
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 :
21 : #include "memilio/utils/compiler_diagnostics.h"
22 :
23 : //see below for line that causes this warning
24 : GCC_CLANG_DIAGNOSTIC(push)
25 : GCC_CLANG_DIAGNOSTIC(ignored "-Wmaybe-uninitialized")
26 :
27 : #include "memilio/config.h"
28 :
29 : #ifdef MEMILIO_HAS_JSONCPP
30 :
31 : #include "ode_secir/parameters_io.h"
32 : #include "memilio/io/epi_data.h"
33 : #include "memilio/io/io.h"
34 : #include "memilio/utils/stl_util.h"
35 : #include "memilio/utils/date.h"
36 :
37 : namespace mio
38 : {
39 :
40 : namespace osecir
41 : {
42 :
43 : namespace details
44 : {
45 : //district, county or state id of a data entry if available, 0 (for whole country) otherwise
46 : //used to compare data entries to integer ids in STL algorithms
47 : template <class EpiDataEntry>
48 1023604 : int get_region_id(const EpiDataEntry& entry)
49 : {
50 1023604 : return entry.county_id
51 2195352 : ? entry.county_id->get()
52 1455683 : : (entry.state_id ? entry.state_id->get() : (entry.district_id ? entry.district_id->get() : 0));
53 : }
54 : //overload for integers, so the comparison of data entry to integers is symmetric (required by e.g. equal_range)
55 1672 : int get_region_id(int id)
56 : {
57 1672 : return id;
58 : }
59 :
60 88 : IOResult<void> read_confirmed_cases_data(
61 : 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 88 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
73 48488 : return a.date < b.date;
74 : });
75 88 : if (max_date_entry == rki_data.end()) {
76 0 : log_error("RKI data file is empty.");
77 0 : return failure(StatusCode::InvalidFileFormat, "RKI file is empty.");
78 : }
79 88 : auto max_date = max_date_entry->date;
80 88 : if (max_date < date) {
81 0 : log_error("Specified date does not exist in RKI data");
82 0 : return failure(StatusCode::OutOfRange, "Specified date does not exist in RKI data.");
83 : }
84 88 : auto days_surplus = std::min(get_offset_in_days(max_date, date) - 6, 0);
85 :
86 : //this statement causes maybe-uninitialized warning for some versions of gcc.
87 : //the error is reported in an included header, so the warning is disabled for the whole file
88 88 : std::sort(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
89 507792 : return std::make_tuple(get_region_id(a), a.date) < std::make_tuple(get_region_id(b), b.date);
90 : });
91 :
92 176 : for (auto region_idx = size_t(0); region_idx < vregion.size(); ++region_idx) {
93 88 : auto region_entry_range_it =
94 1760 : std::equal_range(rki_data.begin(), rki_data.end(), vregion[region_idx], [](auto&& a, auto&& b) {
95 1672 : return get_region_id(a) < get_region_id(b);
96 : });
97 88 : auto region_entry_range = make_range(region_entry_range_it);
98 88 : if (region_entry_range.begin() == region_entry_range.end()) {
99 0 : log_error("No entries found for region {}", vregion[region_idx]);
100 0 : return failure(StatusCode::InvalidFileFormat,
101 0 : "No entries found for region " + std::to_string(vregion[region_idx]));
102 : }
103 48664 : for (auto&& region_entry : region_entry_range) {
104 :
105 48576 : auto& t_Exposed = vt_Exposed[region_idx];
106 48576 : auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx];
107 48576 : auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx];
108 48576 : auto& t_InfectedSevere = vt_InfectedSevere[region_idx];
109 48576 : auto& t_InfectedCritical = vt_InfectedCritical[region_idx];
110 :
111 48576 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
112 48576 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
113 48576 : auto& num_rec = vnum_rec[region_idx];
114 48576 : auto& num_Exposed = vnum_Exposed[region_idx];
115 48576 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
116 48576 : auto& num_death = vnum_death[region_idx];
117 48576 : auto& num_icu = vnum_icu[region_idx];
118 :
119 48576 : auto& mu_C_R = vmu_C_R[region_idx];
120 48576 : auto& mu_I_H = vmu_I_H[region_idx];
121 48576 : auto& mu_H_U = vmu_H_U[region_idx];
122 :
123 48576 : auto date_df = region_entry.date;
124 48576 : auto age = size_t(region_entry.age_group);
125 :
126 48576 : bool read_icu = false; //params.populations.get({age, SecirCompartments::U}) == 0;
127 :
128 48576 : if (date_df == offset_date_by_days(date, 0)) {
129 420 : num_InfectedSymptoms[age] += scaling_factor_inf[age] * region_entry.num_confirmed;
130 420 : num_rec[age] += region_entry.num_confirmed;
131 : }
132 48576 : if (date_df == offset_date_by_days(date, days_surplus)) {
133 420 : num_InfectedNoSymptoms[age] -=
134 420 : 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
135 : }
136 48576 : if (date_df == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) {
137 420 : num_InfectedNoSymptoms[age] +=
138 420 : 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
139 420 : num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
140 : }
141 48576 : if (date_df == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) {
142 420 : num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
143 : }
144 48576 : if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age])) {
145 420 : num_InfectedSymptoms[age] -= scaling_factor_inf[age] * region_entry.num_confirmed;
146 420 : num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
147 : }
148 48576 : if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) {
149 420 : num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
150 420 : if (read_icu) {
151 0 : num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
152 : }
153 : }
154 97152 : if (date_df ==
155 97152 : offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) {
156 420 : num_death[age] += region_entry.num_deaths;
157 420 : if (read_icu) {
158 0 : num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
159 : }
160 : }
161 : }
162 : }
163 :
164 176 : for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) {
165 88 : auto region = vregion[region_idx];
166 :
167 88 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
168 88 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
169 88 : auto& num_rec = vnum_rec[region_idx];
170 88 : auto& num_Exposed = vnum_Exposed[region_idx];
171 88 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
172 88 : auto& num_death = vnum_death[region_idx];
173 88 : auto& num_icu = vnum_icu[region_idx];
174 :
175 616 : for (size_t i = 0; i < ConfirmedCasesDataEntry::age_group_names.size(); i++) {
176 528 : auto try_fix_constraints = [region, i](double& value, double error, auto str) {
177 3696 : if (value < error) {
178 : //this should probably return a failure
179 : //but the algorithm is not robust enough to avoid large negative values and there are tests that rely on it
180 0 : log_error("{:s} for age group {:s} is {:.4f} for region {:d}, exceeds expected negative value.",
181 0 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
182 0 : value = 0.0;
183 : }
184 3696 : else if (value < 0) {
185 0 : log_info("{:s} for age group {:s} is {:.4f} for region {:d}, automatically corrected", str,
186 0 : ConfirmedCasesDataEntry::age_group_names[i], value, region);
187 0 : value = 0.0;
188 : }
189 3696 : };
190 :
191 528 : try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms");
192 528 : try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms");
193 528 : try_fix_constraints(num_Exposed[i], -5, "Exposed");
194 528 : try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere");
195 528 : try_fix_constraints(num_death[i], -5, "Dead");
196 528 : try_fix_constraints(num_icu[i], -5, "InfectedCritical");
197 528 : try_fix_constraints(num_rec[i], -20, "Recovered");
198 : }
199 : }
200 :
201 176 : return success();
202 : }
203 :
204 70 : IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& vregion, Date date,
205 : std::vector<double>& vnum_icu)
206 : {
207 70 : BOOST_OUTCOME_TRY(auto&& divi_data, mio::read_divi_data(path));
208 :
209 70 : auto max_date_entry = std::max_element(divi_data.begin(), divi_data.end(), [](auto&& a, auto&& b) {
210 6370 : return a.date < b.date;
211 : });
212 70 : if (max_date_entry == divi_data.end()) {
213 0 : log_error("DIVI data file is empty.");
214 0 : return failure(StatusCode::InvalidFileFormat, path + ", file is empty.");
215 : }
216 70 : auto max_date = max_date_entry->date;
217 70 : if (max_date < date) {
218 0 : log_error("Specified date does not exist in DIVI data.");
219 0 : return failure(StatusCode::OutOfRange, path + ", specified date does not exist in DIVI data.");
220 : }
221 :
222 6510 : for (auto&& entry : divi_data) {
223 19136 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
224 19136 : return r == 0 || r == get_region_id(entry);
225 : });
226 6440 : auto date_df = entry.date;
227 6440 : if (it != vregion.end() && date_df == date) {
228 70 : auto region_idx = size_t(it - vregion.begin());
229 70 : vnum_icu[region_idx] = entry.num_icu;
230 : }
231 : }
232 :
233 140 : return success();
234 70 : }
235 :
236 : IOResult<std::vector<std::vector<double>>>
237 105 : read_population_data(const std::string& path, const std::vector<int>& vregion, bool accumulate_age_groups)
238 : {
239 105 : BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path, !accumulate_age_groups));
240 : //if we set up the model for one age group, the population data should be read in with the
241 : //age groups given in the population data json file and are accumulated later
242 : //otherwise the populations are directly saved for the correct model age groups
243 105 : size_t age_group_size = accumulate_age_groups ? PopulationDataEntry::age_group_names.size()
244 96 : : ConfirmedCasesDataEntry::age_group_names.size();
245 420 : std::vector<std::vector<double>> vnum_population(vregion.size(), std::vector<double>(age_group_size, 0.0));
246 :
247 42210 : for (auto&& entry : population_data) {
248 246477 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
249 200901 : return r == 0 || (entry.county_id && regions::StateId(r) == regions::get_state_id(int(*entry.county_id))) ||
250 246878 : (entry.county_id && regions::CountyId(r) == *entry.county_id) ||
251 171412 : (entry.district_id && regions::DistrictId(r) == *entry.district_id);
252 : });
253 42105 : if (it != vregion.end()) {
254 519 : auto region_idx = size_t(it - vregion.begin());
255 519 : auto& num_population = vnum_population[region_idx];
256 3678 : for (size_t age = 0; age < num_population.size(); age++) {
257 3159 : num_population[age] += entry.population[AgeGroup(age)];
258 : }
259 : }
260 : }
261 105 : if (accumulate_age_groups) {
262 36 : std::vector<std::vector<double>> vnum_pop_acc(vregion.size(), std::vector<double>(1, 0));
263 18 : for (size_t region = 0; region < vregion.size(); ++region) {
264 9 : vnum_pop_acc[region][0] =
265 9 : std::accumulate(vnum_population[region].begin(), vnum_population[region].end(), 0.0);
266 : }
267 9 : return success(vnum_pop_acc);
268 9 : }
269 : else {
270 96 : return success(vnum_population);
271 : }
272 105 : }
273 :
274 : } // namespace details
275 : } // namespace osecir
276 : } // namespace mio
277 :
278 : #endif // MEMILIO_HAS_JSONCPP
279 :
280 : GCC_CLANG_DIAGNOSTIC(pop)
|
