Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Lena Ploetzke, Anna Wendler
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 "ide_secir/parameters_io.h"
22 : #include "memilio/config.h"
23 : #include "memilio/epidemiology/age_group.h"
24 : #include <cstddef>
25 : #include <iostream>
26 :
27 : #ifdef MEMILIO_HAS_JSONCPP
28 :
29 : #include "ide_secir/model.h"
30 : #include "ide_secir/infection_state.h"
31 : #include "memilio/math/eigen.h"
32 : #include "memilio/io/epi_data.h"
33 : #include "memilio/io/io.h"
34 : #include "memilio/utils/date.h"
35 :
36 : #include <string>
37 : #include <cmath>
38 :
39 : namespace mio
40 : {
41 : namespace isecir
42 : {
43 :
44 5 : IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date,
45 : ScalarType scale_confirmed_cases)
46 : {
47 : //--- Preparations ---
48 : // Try to get RKI data from path.
49 5 : BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_data(path));
50 5 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
51 2204 : return a.date < b.date;
52 : });
53 5 : if (max_date_entry == rki_data.end()) {
54 1 : log_error("RKI data file is empty.");
55 1 : return failure(StatusCode::InvalidFileFormat, path + ", file is empty.");
56 : }
57 4 : auto max_date = max_date_entry->date;
58 4 : if (max_date < date) {
59 1 : log_error("Specified date does not exist in RKI data.");
60 1 : return failure(StatusCode::OutOfRange, path + ", specified date does not exist in RKI data.");
61 : }
62 3 : auto min_date_entry = std::min_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
63 1653 : return a.date < b.date;
64 : });
65 3 : auto min_date = min_date_entry->date;
66 :
67 : // Get (global) support_max to determine how many flows in the past we have to compute.
68 3 : ScalarType global_support_max = model.get_global_support_max(dt);
69 3 : Eigen::Index global_support_max_index = Eigen::Index(std::ceil(global_support_max / dt));
70 :
71 : // Get the number of AgeGroups.
72 3 : const size_t num_age_groups = ConfirmedCasesDataEntry::age_group_names.size();
73 :
74 : // m_transitions should be empty at the beginning.
75 :
76 3 : if (model.m_transitions.get_num_time_points() > 0) {
77 2 : model.m_transitions = TimeSeries<ScalarType>(Eigen::Index(InfectionTransition::Count) * num_age_groups);
78 : }
79 3 : if (model.m_populations.get_time(0) != 0) {
80 1 : model.m_populations.remove_last_time_point();
81 2 : model.m_populations.add_time_point<Eigen::VectorXd>(
82 2 : 0, TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count * num_age_groups, 0));
83 : }
84 :
85 : // The first time we need is -4 * global_support_max.
86 3 : Eigen::Index start_shift = 4 * global_support_max_index;
87 : // The last time needed is dependent on the mean stay time in the Exposed compartment and
88 : // the mean stay time of asymptomatic individuals in InfectedNoSymptoms.
89 : // The mean stay time in a compartment may be dependent on the AgeGroup.
90 3 : CustomIndexArray<ScalarType, AgeGroup> mean_ExposedToInfectedNoSymptoms =
91 3 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
92 3 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedNoSymptomsToInfectedSymptoms =
93 3 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
94 3 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedSymptomsToInfectedSevere =
95 3 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
96 3 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedSevereToInfectedCritical =
97 3 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
98 3 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedCriticalToDead =
99 3 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
100 3 : Eigen::Index last_time_index_needed = 0;
101 :
102 21 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) {
103 : // Set the Dead compartment to 0 so that RKI data can be added correctly.
104 18 : int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
105 18 : model.m_populations[0][Di] = 0;
106 :
107 18 : mean_ExposedToInfectedNoSymptoms[group] =
108 : model.parameters
109 18 : .get<TransitionDistributions>()[group][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)]
110 18 : .get_mean(dt);
111 18 : mean_InfectedNoSymptomsToInfectedSymptoms[group] =
112 : model.parameters
113 18 : .get<TransitionDistributions>()[group]
114 18 : [Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]
115 18 : .get_mean(dt);
116 18 : mean_InfectedSymptomsToInfectedSevere[group] =
117 : model.parameters
118 18 : .get<TransitionDistributions>()[group]
119 18 : [Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere)]
120 18 : .get_mean(dt);
121 18 : mean_InfectedSevereToInfectedCritical[group] =
122 : model.parameters
123 18 : .get<TransitionDistributions>()[group]
124 18 : [Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical)]
125 18 : .get_mean(dt);
126 18 : mean_InfectedCriticalToDead[group] =
127 : model.parameters
128 18 : .get<TransitionDistributions>()[group][Eigen::Index(InfectionTransition::InfectedCriticalToDead)]
129 18 : .get_mean(dt);
130 18 : if (last_time_index_needed <
131 18 : Eigen::Index(std::ceil(
132 18 : (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt))) {
133 3 : last_time_index_needed = Eigen::Index(std::ceil(
134 3 : (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt));
135 : }
136 : }
137 :
138 : // Create TimeSeries with zeros. The index of time zero is start_shift.
139 66 : for (Eigen::Index i = -start_shift; i <= last_time_index_needed; i++) {
140 : // Add time point.
141 126 : model.m_transitions.add_time_point(
142 126 : i * dt, TimeSeries<ScalarType>::Vector::Constant((size_t)InfectionTransition::Count * num_age_groups, 0.));
143 : }
144 :
145 : // Setting the entries in m_total_confirmed_cases to zero before overwriting it with the RKI data.
146 3 : model.m_total_confirmed_cases = CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
147 : //--- Calculate the flow InfectedNoSymptomsToInfectedSymptoms using the RKI data and store in the m_transitions object.---
148 3 : ScalarType min_offset_needed = std::ceil(
149 3 : model.m_transitions.get_time(0) -
150 : 1); // Need -1 if first time point is integer and just the floor value if not, therefore use ceil and -1
151 3 : ScalarType max_offset_needed = std::ceil(model.m_transitions.get_last_time());
152 3 : bool min_offset_needed_avail = false;
153 3 : bool max_offset_needed_avail = false;
154 :
155 : // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled.
156 : // Define dummy variables to store the first and the last index of the TimeSeries where the considered entry of
157 : // rki_data is potentially needed.
158 3 : Eigen::Index idx_needed_first = 0;
159 3 : Eigen::Index idx_needed_last = 0;
160 3 : ScalarType time_idx = 0;
161 1659 : for (auto&& entry : rki_data) {
162 1656 : int offset = get_offset_in_days(entry.date, date);
163 1656 : AgeGroup group = entry.age_group;
164 1656 : if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) {
165 150 : if (offset == min_offset_needed) {
166 12 : min_offset_needed_avail = true;
167 : }
168 150 : if (offset == max_offset_needed) {
169 12 : max_offset_needed_avail = true;
170 : }
171 : // Smallest index for which the entry is needed.
172 150 : idx_needed_first =
173 150 : Eigen::Index(std::max(std::floor((offset - model.m_transitions.get_time(0) - 1) / dt), 0.));
174 : // Biggest index for which the entry is needed.
175 150 : idx_needed_last = Eigen::Index(std::min(std::ceil((offset - model.m_transitions.get_time(0) + 1) / dt),
176 300 : double(model.m_transitions.get_num_time_points() - 1)));
177 :
178 300 : int INStISyi = model.get_transition_flat_index(
179 150 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group);
180 :
181 804 : for (Eigen::Index i = idx_needed_first; i <= idx_needed_last; i++) {
182 :
183 654 : time_idx = model.m_transitions.get_time(i);
184 654 : if (offset == int(std::floor(time_idx))) {
185 528 : model.m_transitions[i][INStISyi] +=
186 264 : (1 - (time_idx - std::floor(time_idx))) * scale_confirmed_cases * entry.num_confirmed;
187 : }
188 654 : if (offset == int(std::ceil(time_idx))) {
189 528 : model.m_transitions[i][INStISyi] +=
190 264 : (time_idx - std::floor(time_idx)) * scale_confirmed_cases * entry.num_confirmed;
191 : }
192 654 : if (offset == int(std::floor(time_idx - dt))) {
193 528 : model.m_transitions[i][INStISyi] -=
194 264 : (1 - (time_idx - dt - std::floor(time_idx - dt))) * scale_confirmed_cases * entry.num_confirmed;
195 : }
196 654 : if (offset == int(std::ceil(time_idx - dt))) {
197 528 : model.m_transitions[i][INStISyi] -=
198 264 : (time_idx - dt - std::floor(time_idx - dt)) * scale_confirmed_cases * entry.num_confirmed;
199 : }
200 : }
201 :
202 : // Compute Dead by shifting RKI data according to mean stay times.
203 : // This is done because the RKI reports death with the date of positive test instead of the date of deaths.
204 150 : int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
205 300 : if (offset ==
206 300 : std::floor(-mean_InfectedSymptomsToInfectedSevere[group] -
207 150 : mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) {
208 12 : model.m_populations[0][Di] +=
209 12 : (1 -
210 12 : (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] -
211 12 : mean_InfectedCriticalToDead[group] -
212 12 : std::floor(-mean_InfectedSymptomsToInfectedSevere[group] -
213 12 : mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group]))) *
214 12 : entry.num_deaths;
215 : }
216 300 : if (offset ==
217 300 : std::ceil(-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] -
218 150 : mean_InfectedCriticalToDead[group])) {
219 12 : model.m_populations[0][Di] +=
220 12 : (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] -
221 12 : mean_InfectedCriticalToDead[group] -
222 12 : std::floor(-mean_InfectedSymptomsToInfectedSevere[group] -
223 12 : mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) *
224 12 : entry.num_deaths;
225 : }
226 :
227 150 : if (offset == 0) {
228 18 : model.m_total_confirmed_cases[group] = scale_confirmed_cases * entry.num_confirmed;
229 : }
230 : }
231 : }
232 :
233 3 : if (!max_offset_needed_avail) {
234 1 : log_error("Necessary range of dates needed to compute initial values does not exist in RKI data.");
235 1 : return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data.");
236 : }
237 2 : if (!min_offset_needed_avail) {
238 1 : std::string min_date_string =
239 7 : std::to_string(min_date.day) + "." + std::to_string(min_date.month) + "." + std::to_string(min_date.year);
240 : // Get first date that is needed.
241 1 : mio::Date min_offset_date = offset_date_by_days(date, int(min_offset_needed));
242 6 : std::string min_offset_date_string = std::to_string(min_offset_date.day) + "." +
243 5 : std::to_string(min_offset_date.month) + "." +
244 3 : std::to_string(min_offset_date.year);
245 4 : log_warning("RKI data is needed from " + min_offset_date_string +
246 3 : " to compute initial values. RKI data is only available from " + min_date_string +
247 : ". Missing dates were set to 0.");
248 1 : }
249 :
250 : //--- Calculate the flows "after" InfectedNoSymptomsToInfectedSymptoms. ---
251 : // Set m_support_max_vector and m_derivative_vector in the model which is needed for the following computations.
252 2 : model.set_transitiondistributions_support_max(dt);
253 2 : model.set_transitiondistributions_derivative(dt);
254 :
255 14 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) {
256 : //--- Calculate the flows "after" InfectedNoSymptomsToInfectedSymptoms. ---
257 : // Compute flow InfectedSymptomsToInfectedSevere for -3 * global_support_max, ..., 0.
258 168 : for (Eigen::Index i = -3 * global_support_max_index; i <= 0; i++) {
259 156 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
260 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt,
261 : i + start_shift, group);
262 : }
263 : // Compute flow InfectedSevereToInfectedCritical for -2 * global_support_max, ..., 0.
264 120 : for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
265 108 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
266 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift,
267 : group);
268 : }
269 : // Compute flows from InfectedSymptoms, InfectedSevere and InfectedCritical to Recovered and
270 : // flow InfectedCriticalToDead for -global_support_max, ..., 0.
271 72 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
272 : // Compute flow InfectedSymptomsToRecovered.
273 60 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered),
274 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt,
275 : i + start_shift, group);
276 : // Compute flow InfectedSevereToRecovered.
277 60 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered),
278 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift,
279 : group);
280 : // Compute flow InfectedCriticalToRecovered.
281 60 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered),
282 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift,
283 : group);
284 : // Compute flow InfectedCriticalToDead.
285 60 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead),
286 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift,
287 : group);
288 : }
289 :
290 : //--- Calculate the remaining flows. ---
291 : // Compute flow ExposedToInfectedNoSymptoms for -2 * global_support_max, ..., 0.
292 : // Use mean value of the TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
293 :
294 12 : Eigen::Index index_shift_mean = Eigen::Index(std::round(mean_InfectedNoSymptomsToInfectedSymptoms[group] / dt));
295 : int EtINSi =
296 12 : model.get_transition_flat_index(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), group);
297 24 : int INStISyi = model.get_transition_flat_index(
298 12 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group);
299 120 : for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
300 108 : model.m_transitions[i + start_shift][EtINSi] =
301 216 : (1 / model.parameters.get<TransitionProbabilities>()[group][Eigen::Index(
302 216 : InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
303 216 : model.m_transitions[i + start_shift + index_shift_mean][INStISyi];
304 : }
305 :
306 : // Compute flow SusceptibleToExposed for -global_support_max, ..., 0.
307 : // Use mean values of the TransitionDistribution ExposedToInfectedNoSymptoms and of the
308 : // TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
309 12 : index_shift_mean = Eigen::Index(std::round(
310 12 : (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt));
311 12 : int StEi = model.get_transition_flat_index(Eigen::Index(InfectionTransition::SusceptibleToExposed), group);
312 :
313 72 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
314 60 : model.m_transitions[i + start_shift][StEi] =
315 120 : (1 / model.parameters.get<TransitionProbabilities>()[group][Eigen::Index(
316 120 : InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
317 120 : model.m_transitions[i + start_shift + index_shift_mean][INStISyi];
318 : }
319 :
320 : // InfectedNoSymptomsToRecovered for -global_support_max, ..., 0.
321 : // If we previously calculated the transition ExposedToInfectedNoSymptoms, we can calculate this transition
322 : // using the standard formula.
323 72 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
324 60 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered),
325 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, i + start_shift,
326 : group);
327 : }
328 : }
329 :
330 : // At the end of the calculation, delete all time points that are not required for the simulation.
331 2 : auto transition_copy(model.m_transitions);
332 2 : model.m_transitions = TimeSeries<ScalarType>(Eigen::Index(InfectionTransition::Count) * num_age_groups);
333 12 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
334 10 : model.m_transitions.add_time_point(i * dt, transition_copy.get_value(i + start_shift));
335 : }
336 :
337 4 : return mio::success();
338 5 : }
339 :
340 : } // namespace isecir
341 : } // namespace mio
342 :
343 : #endif // MEMILIO_HAS_JSONCPP
|
