Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 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 1111885 : int get_region_id(const EpiDataEntry& entry)
49 : {
50 1111885 : return entry.county_id
51 2371914 : ? entry.county_id->get()
52 1543964 : : (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 1843 : int get_region_id(int id)
56 : {
57 1843 : return id;
58 : }
59 :
60 97 : 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 97 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
73 53447 : return a.date < b.date;
74 : });
75 97 : 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 97 : auto max_date = max_date_entry->date;
80 97 : 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 97 : 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 97 : std::sort(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
89 551847 : return std::make_tuple(get_region_id(a), a.date) < std::make_tuple(get_region_id(b), b.date);
90 : });
91 :
92 194 : for (auto region_idx = size_t(0); region_idx < vregion.size(); ++region_idx) {
93 97 : auto region_entry_range_it =
94 1940 : std::equal_range(rki_data.begin(), rki_data.end(), vregion[region_idx], [](auto&& a, auto&& b) {
95 1843 : return get_region_id(a) < get_region_id(b);
96 : });
97 97 : auto region_entry_range = make_range(region_entry_range_it);
98 97 : 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 53641 : for (auto&& region_entry : region_entry_range) {
104 :
105 53544 : auto& t_Exposed = vt_Exposed[region_idx];
106 53544 : auto& t_InfectedNoSymptoms = vt_InfectedNoSymptoms[region_idx];
107 53544 : auto& t_InfectedSymptoms = vt_InfectedSymptoms[region_idx];
108 53544 : auto& t_InfectedSevere = vt_InfectedSevere[region_idx];
109 53544 : auto& t_InfectedCritical = vt_InfectedCritical[region_idx];
110 :
111 53544 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
112 53544 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
113 53544 : auto& num_rec = vnum_rec[region_idx];
114 53544 : auto& num_Exposed = vnum_Exposed[region_idx];
115 53544 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
116 53544 : auto& num_death = vnum_death[region_idx];
117 53544 : auto& num_icu = vnum_icu[region_idx];
118 :
119 53544 : auto& mu_C_R = vmu_C_R[region_idx];
120 53544 : auto& mu_I_H = vmu_I_H[region_idx];
121 53544 : auto& mu_H_U = vmu_H_U[region_idx];
122 :
123 53544 : auto date_df = region_entry.date;
124 53544 : auto age = size_t(region_entry.age_group);
125 :
126 53544 : if (date_df == offset_date_by_days(date, 0)) {
127 429 : num_InfectedSymptoms[age] += scaling_factor_inf[age] * region_entry.num_confirmed;
128 429 : num_rec[age] += region_entry.num_confirmed;
129 : }
130 53544 : if (date_df == offset_date_by_days(date, days_surplus)) {
131 429 : num_InfectedNoSymptoms[age] -=
132 429 : 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
133 : }
134 53544 : if (date_df == offset_date_by_days(date, t_InfectedNoSymptoms[age] + days_surplus)) {
135 429 : num_InfectedNoSymptoms[age] +=
136 429 : 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
137 429 : num_Exposed[age] -= 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
138 : }
139 53544 : if (date_df == offset_date_by_days(date, t_Exposed[age] + t_InfectedNoSymptoms[age] + days_surplus)) {
140 429 : num_Exposed[age] += 1 / (1 - mu_C_R[age]) * scaling_factor_inf[age] * region_entry.num_confirmed;
141 : }
142 53544 : if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age])) {
143 429 : num_InfectedSymptoms[age] -= scaling_factor_inf[age] * region_entry.num_confirmed;
144 429 : num_InfectedSevere[age] += mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
145 : }
146 53544 : if (date_df == offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age])) {
147 429 : num_InfectedSevere[age] -= mu_I_H[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
148 429 : num_icu[age] += mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
149 : }
150 107088 : if (date_df ==
151 107088 : offset_date_by_days(date, -t_InfectedSymptoms[age] - t_InfectedSevere[age] - t_InfectedCritical[age])) {
152 429 : num_death[age] += region_entry.num_deaths;
153 429 : num_icu[age] -= mu_I_H[age] * mu_H_U[age] * scaling_factor_inf[age] * region_entry.num_confirmed;
154 : }
155 : }
156 : }
157 :
158 194 : for (size_t region_idx = 0; region_idx < vregion.size(); ++region_idx) {
159 97 : auto region = vregion[region_idx];
160 :
161 97 : auto& num_InfectedNoSymptoms = vnum_InfectedNoSymptoms[region_idx];
162 97 : auto& num_InfectedSymptoms = vnum_InfectedSymptoms[region_idx];
163 97 : auto& num_rec = vnum_rec[region_idx];
164 97 : auto& num_Exposed = vnum_Exposed[region_idx];
165 97 : auto& num_InfectedSevere = vnum_InfectedSevere[region_idx];
166 97 : auto& num_death = vnum_death[region_idx];
167 97 : auto& num_icu = vnum_icu[region_idx];
168 :
169 679 : for (size_t i = 0; i < ConfirmedCasesDataEntry::age_group_names.size(); i++) {
170 591 : auto try_fix_constraints = [region, i](double& value, double error, auto str) {
171 4074 : if (value < error) {
172 : //this should probably return a failure
173 : //but the algorithm is not robust enough to avoid large negative values and there are tests that rely on it
174 9 : log_error("{:s} for age group {:s} is {:.4f} for region {:d}, exceeds expected negative value.",
175 9 : str, ConfirmedCasesDataEntry::age_group_names[i], value, region);
176 9 : value = 0.0;
177 : }
178 4065 : else if (value < 0) {
179 0 : log_info("{:s} for age group {:s} is {:.4f} for region {:d}, automatically corrected", str,
180 0 : ConfirmedCasesDataEntry::age_group_names[i], value, region);
181 0 : value = 0.0;
182 : }
183 4074 : };
184 :
185 582 : try_fix_constraints(num_InfectedSymptoms[i], -5, "InfectedSymptoms");
186 582 : try_fix_constraints(num_InfectedNoSymptoms[i], -5, "InfectedNoSymptoms");
187 582 : try_fix_constraints(num_Exposed[i], -5, "Exposed");
188 582 : try_fix_constraints(num_InfectedSevere[i], -5, "InfectedSevere");
189 582 : try_fix_constraints(num_death[i], -5, "Dead");
190 582 : try_fix_constraints(num_icu[i], -5, "InfectedCritical");
191 582 : try_fix_constraints(num_rec[i], -20, "Recovered");
192 : }
193 : }
194 :
195 194 : return success();
196 : }
197 :
198 70 : IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& vregion, Date date,
199 : std::vector<double>& vnum_icu)
200 : {
201 70 : BOOST_OUTCOME_TRY(auto&& divi_data, mio::read_divi_data(path));
202 :
203 70 : auto max_date_entry = std::max_element(divi_data.begin(), divi_data.end(), [](auto&& a, auto&& b) {
204 6370 : return a.date < b.date;
205 : });
206 70 : if (max_date_entry == divi_data.end()) {
207 0 : log_error("DIVI data file is empty.");
208 0 : return failure(StatusCode::InvalidFileFormat, path + ", file is empty.");
209 : }
210 70 : auto max_date = max_date_entry->date;
211 70 : if (max_date < date) {
212 0 : log_error("Specified date does not exist in DIVI data.");
213 0 : return failure(StatusCode::OutOfRange, path + ", specified date does not exist in DIVI data.");
214 : }
215 :
216 6510 : for (auto&& entry : divi_data) {
217 19136 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
218 19136 : return r == 0 || r == get_region_id(entry);
219 : });
220 6440 : auto date_df = entry.date;
221 6440 : if (it != vregion.end() && date_df == date) {
222 70 : auto region_idx = size_t(it - vregion.begin());
223 70 : vnum_icu[region_idx] = entry.num_icu;
224 : }
225 : }
226 :
227 140 : return success();
228 70 : }
229 :
230 : IOResult<std::vector<std::vector<double>>>
231 105 : read_population_data(const std::string& path, const std::vector<int>& vregion, bool accumulate_age_groups)
232 : {
233 105 : BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path, !accumulate_age_groups));
234 : //if we set up the model for one age group, the population data should be read in with the
235 : //age groups given in the population data json file and are accumulated later
236 : //otherwise the populations are directly saved for the correct model age groups
237 105 : size_t age_group_size = accumulate_age_groups ? PopulationDataEntry::age_group_names.size()
238 96 : : ConfirmedCasesDataEntry::age_group_names.size();
239 420 : std::vector<std::vector<double>> vnum_population(vregion.size(), std::vector<double>(age_group_size, 0.0));
240 :
241 42210 : for (auto&& entry : population_data) {
242 246477 : auto it = std::find_if(vregion.begin(), vregion.end(), [&entry](auto r) {
243 200901 : return r == 0 || (entry.county_id && regions::StateId(r) == regions::get_state_id(int(*entry.county_id))) ||
244 246878 : (entry.county_id && regions::CountyId(r) == *entry.county_id) ||
245 171412 : (entry.district_id && regions::DistrictId(r) == *entry.district_id);
246 : });
247 42105 : if (it != vregion.end()) {
248 519 : auto region_idx = size_t(it - vregion.begin());
249 519 : auto& num_population = vnum_population[region_idx];
250 3678 : for (size_t age = 0; age < num_population.size(); age++) {
251 3159 : num_population[age] += entry.population[AgeGroup(age)];
252 : }
253 : }
254 : }
255 105 : if (accumulate_age_groups) {
256 36 : std::vector<std::vector<double>> vnum_pop_acc(vregion.size(), std::vector<double>(1, 0));
257 18 : for (size_t region = 0; region < vregion.size(); ++region) {
258 9 : vnum_pop_acc[region][0] =
259 9 : std::accumulate(vnum_population[region].begin(), vnum_population[region].end(), 0.0);
260 : }
261 9 : return success(vnum_pop_acc);
262 9 : }
263 : else {
264 96 : return success(vnum_population);
265 : }
266 105 : }
267 :
268 : } // namespace details
269 : } // namespace osecir
270 : } // namespace mio
271 :
272 : #endif // MEMILIO_HAS_JSONCPP
273 :
274 : GCC_CLANG_DIAGNOSTIC(pop)
|
