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 : #ifndef IDE_INITIALFLOWS_H
21 : #define IDE_INITIALFLOWS_H
22 :
23 : #include "memilio/config.h"
24 :
25 : #ifdef MEMILIO_HAS_JSONCPP
26 :
27 : #include "ide_secir/model.h"
28 : #include "ide_secir/infection_state.h"
29 :
30 : #include "memilio/epidemiology/age_group.h"
31 : #include "memilio/io/epi_data.h"
32 : #include "memilio/io/io.h"
33 : #include "memilio/utils/date.h"
34 : #include "memilio/utils/logging.h"
35 :
36 : #include <string>
37 :
38 : namespace mio
39 : {
40 : namespace isecir
41 : {
42 :
43 : /**
44 : * @brief Computes a TimeSeries of flows to provide initial data for an IDE-SECIR model with data from RKI.
45 : *
46 : * The flows InfectedNoSymptomsToInfectedSymptoms are calculated using the confirmed cases in the RKI data.
47 : * If necessary, the RKI data are linearly interpolated within one day.
48 : * The RKI data should contain data for each needed day without division of age groups, the completeness of the dates is not verified.
49 : * Data can be downloaded e.g. with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py,
50 : * which creates a file named cases_all_germany.json or a similar name. One should set impute_dates=True so that missing dates are imputed.
51 : *
52 : * The flows InfectedSymptomsToInfectedSevere, InfectedSymptomsToRecovered, InfectedSevereToInfectedCritical,
53 : * InfectedSevereToRecovered, InfectedCriticalToDead and InfectedCriticalToRecovered can then be calculated using
54 : * the InfectedNoSymptomsToInfectedSymptoms flow with the standard formula from the IDE model.
55 : * The ExposedToInfectedNoSymptoms and InfectedNoSymptomsToInfectedSymptoms flows are calculated
56 : * using the means of the respective TransitionDistribution.
57 : * The flow InfectedNoSymptomsToInfectedSymptoms is calculated with the standard formula from the IDE model
58 : * using the results for ExposedToInfectedNoSymptoms.
59 : *
60 : * The number of deaths used in the model is computed by shifting the reported RKI data according to the mean values of the transitions
61 : * InfectedSymptomsToInfectedSevere, InfectedSevereToInfectedCritical and InfectedCriticalToDead.
62 : * We also set the number of total confirmed cases in the model.
63 : * Therefore the initialization method using the total confirmed cases is used in the model. See also the documentation of the model class.
64 : *
65 : * The start date of the model simulation is set to t0=0.
66 : *
67 : * @param[in, out] model The model for which the initial flows should be computed.
68 : * @param[in] dt Time step size.
69 : * @param[in] rki_data Vector containing RKI data.
70 : * @param[in] date The start date of the simulation and the last time point of the TimeSeries used for initialization.
71 : * @param[in] scale_confirmed_cases Vector with factor(s for each age group) by which to scale the confirmed cases of
72 : * rki_data to consider unreported cases.
73 : * @tparam EntryType is expected to be ConfirmedCasesNoAgeEntry for data that is not age resolved and
74 : * ConfirmedCasesDataEntry for age resolved data. See also epi_data.h.
75 : * @returns Any io errors that happen during reading of the files.
76 : */
77 :
78 : template <typename EntryType>
79 12 : IOResult<void> set_initial_flows(Model& model, const ScalarType dt, const std::vector<EntryType> rki_data,
80 : const Date date, const CustomIndexArray<ScalarType, AgeGroup> scale_confirmed_cases)
81 : {
82 : // Check if scale_confirmed_cases has the right size (= number of age groups).
83 12 : assert(model.get_num_agegroups() == (size_t)scale_confirmed_cases.size());
84 : // Check if the correct EntryType was used.
85 : if constexpr (std::is_same_v<EntryType, ConfirmedCasesDataEntry>) {
86 6 : assert(model.get_num_agegroups() == (size_t)EntryType::age_group_names.size());
87 : }
88 : else {
89 6 : assert(model.get_num_agegroups() == 1);
90 : }
91 :
92 : //--- Preparations ---
93 2837 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
94 2825 : return a.date < b.date;
95 : });
96 12 : if (max_date_entry == rki_data.end()) {
97 2 : log_error("RKI data file is empty.");
98 2 : return failure(StatusCode::InvalidFileFormat, "RKI data file is empty.");
99 : }
100 10 : auto max_date = max_date_entry->date;
101 10 : if (max_date < date) {
102 2 : log_error("Specified date does not exist in RKI data.");
103 2 : return failure(StatusCode::OutOfRange, "Specified date does not exist in RKI data.");
104 : }
105 :
106 : // Get (global) support_max to determine how many flows in the past we have to compute.
107 8 : ScalarType global_support_max = model.get_global_support_max(dt);
108 8 : Eigen::Index global_support_max_index = Eigen::Index(std::ceil(global_support_max / dt));
109 :
110 : // Get the number of AgeGroups.
111 8 : const size_t num_age_groups = model.get_num_agegroups();
112 :
113 : // m_transitions should be empty at the beginning.
114 8 : if (model.transitions.get_num_time_points() > 0) {
115 6 : model.transitions = TimeSeries<ScalarType>(Eigen::Index(InfectionTransition::Count) * num_age_groups);
116 : }
117 8 : if (model.populations.get_time(0) != 0) {
118 2 : model.populations.remove_last_time_point();
119 4 : model.populations.add_time_point<Eigen::VectorXd>(
120 2 : 0, TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count * num_age_groups, 0));
121 : }
122 :
123 : // The first time we need is -4 * global_support_max because we need values for
124 : // InfectedNoSymptomsToInfectedSymptoms on this time window to compute all consecutive transitions on the time
125 : // window from -global_support_max to 0.
126 8 : Eigen::Index start_shift = 4 * global_support_max_index;
127 : // The last time needed is dependent on the mean stay time in the Exposed compartment and
128 : // the mean stay time of asymptomatic individuals in InfectedNoSymptoms.
129 : // The mean stay time in a compartment may be dependent on the AgeGroup.
130 24 : CustomIndexArray<ScalarType, AgeGroup> mean_ExposedToInfectedNoSymptoms =
131 24 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
132 24 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedNoSymptomsToInfectedSymptoms =
133 24 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
134 24 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedSymptomsToInfectedSevere =
135 24 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
136 24 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedSevereToInfectedCritical =
137 24 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
138 24 : CustomIndexArray<ScalarType, AgeGroup> mean_InfectedCriticalToDead =
139 24 : CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
140 8 : Eigen::Index last_time_index_needed = 0;
141 :
142 36 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) {
143 : // Set the Dead compartment to 0 so that RKI data can be added correctly.
144 28 : int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
145 28 : model.populations[0][Di] = 0;
146 :
147 28 : mean_ExposedToInfectedNoSymptoms[group] =
148 : model.parameters
149 28 : .get<TransitionDistributions>()[group][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)]
150 28 : .get_mean(dt);
151 28 : mean_InfectedNoSymptomsToInfectedSymptoms[group] =
152 : model.parameters
153 28 : .get<TransitionDistributions>()[group]
154 : [Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]
155 28 : .get_mean(dt);
156 28 : mean_InfectedSymptomsToInfectedSevere[group] =
157 : model.parameters
158 28 : .get<TransitionDistributions>()[group]
159 : [Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere)]
160 28 : .get_mean(dt);
161 28 : mean_InfectedSevereToInfectedCritical[group] =
162 : model.parameters
163 28 : .get<TransitionDistributions>()[group]
164 : [Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical)]
165 28 : .get_mean(dt);
166 28 : mean_InfectedCriticalToDead[group] =
167 : model.parameters
168 28 : .get<TransitionDistributions>()[group][Eigen::Index(InfectionTransition::InfectedCriticalToDead)]
169 28 : .get_mean(dt);
170 28 : if (last_time_index_needed <
171 28 : Eigen::Index(std::ceil(
172 28 : (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt))) {
173 8 : last_time_index_needed = Eigen::Index(std::ceil(
174 8 : (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt));
175 : }
176 : }
177 :
178 : // Create TimeSeries with zeros. The index of time zero is start_shift.
179 176 : for (Eigen::Index i = -start_shift; i <= last_time_index_needed; i++) {
180 : // Add time point.
181 336 : model.transitions.add_time_point(
182 336 : i * dt, TimeSeries<ScalarType>::Vector::Constant((size_t)InfectionTransition::Count * num_age_groups, 0.));
183 : }
184 :
185 : // Setting the entries in m_total_confirmed_cases to zero before overwriting it with the RKI data.
186 8 : model.total_confirmed_cases = CustomIndexArray<ScalarType, AgeGroup>(AgeGroup(num_age_groups), 0.);
187 : //--- Calculate the flow InfectedNoSymptomsToInfectedSymptoms using the RKI data and store in the m_transitions object.---
188 16 : ScalarType min_offset_needed = std::ceil(
189 8 : model.transitions.get_time(0) -
190 : 1); // Need -1 if first time point is integer and just the floor value if not, therefore use ceil and -1
191 8 : ScalarType max_offset_needed = std::ceil(model.transitions.get_last_time());
192 8 : bool min_offset_needed_avail = false;
193 8 : bool max_offset_needed_avail = false;
194 :
195 : // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled.
196 : // Define dummy variables to store the first and the last index of the TimeSeries where the considered entry of
197 : // rki_data is potentially needed.
198 8 : Eigen::Index idx_needed_first = 0;
199 8 : Eigen::Index idx_needed_last = 0;
200 8 : ScalarType time_idx = 0;
201 :
202 2276 : for (auto&& entry : rki_data) {
203 2268 : int offset = get_offset_in_days(entry.date, date);
204 :
205 : // Get the index regarding the age group.
206 : // If we don't have age resolution and use EntryType=ConfirmedCasesNoAge, the index is set to 1.
207 : // If we consider multiple age groups and use EntryType=ConfirmedCasesDataEntry, it is determined accordingly.
208 2268 : AgeGroup group = AgeGroup(0);
209 : if constexpr (std::is_same_v<EntryType, ConfirmedCasesDataEntry>) {
210 2208 : group = entry.age_group;
211 : }
212 :
213 2268 : if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) {
214 259 : if (offset == min_offset_needed) {
215 21 : min_offset_needed_avail = true;
216 : }
217 259 : if (offset == max_offset_needed) {
218 21 : max_offset_needed_avail = true;
219 : }
220 : // Smallest index for which the entry is needed.
221 259 : idx_needed_first =
222 259 : Eigen::Index(std::max(std::floor((offset - model.transitions.get_time(0) - 1) / dt), 0.));
223 : // Biggest index for which the entry is needed.
224 259 : idx_needed_last = Eigen::Index(std::min(std::ceil((offset - model.transitions.get_time(0) + 1) / dt),
225 518 : double(model.transitions.get_num_time_points() - 1)));
226 :
227 259 : int INStISyi = model.get_transition_flat_index(
228 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group);
229 :
230 1386 : for (Eigen::Index i = idx_needed_first; i <= idx_needed_last; i++) {
231 :
232 1127 : time_idx = model.transitions.get_time(i);
233 1127 : if (offset == int(std::floor(time_idx))) {
234 455 : model.transitions[i][INStISyi] +=
235 455 : (1 - (time_idx - std::floor(time_idx))) * scale_confirmed_cases[group] * entry.num_confirmed;
236 : }
237 1127 : if (offset == int(std::ceil(time_idx))) {
238 455 : model.transitions[i][INStISyi] +=
239 455 : (time_idx - std::floor(time_idx)) * scale_confirmed_cases[group] * entry.num_confirmed;
240 : }
241 1127 : if (offset == int(std::floor(time_idx - dt))) {
242 910 : model.transitions[i][INStISyi] -= (1 - (time_idx - dt - std::floor(time_idx - dt))) *
243 455 : scale_confirmed_cases[group] * entry.num_confirmed;
244 : }
245 1127 : if (offset == int(std::ceil(time_idx - dt))) {
246 910 : model.transitions[i][INStISyi] -= (time_idx - dt - std::floor(time_idx - dt)) *
247 455 : scale_confirmed_cases[group] * entry.num_confirmed;
248 : }
249 : }
250 :
251 : // Compute Dead by shifting RKI data according to mean stay times.
252 : // This is done because the RKI reports death with the date of positive test instead of the date of deaths.
253 259 : int Di = model.get_state_flat_index(Eigen::Index(InfectionState::Dead), group);
254 518 : if (offset ==
255 518 : std::floor(-mean_InfectedSymptomsToInfectedSevere[group] -
256 259 : mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) {
257 21 : model.populations[0][Di] +=
258 21 : (1 -
259 21 : (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] -
260 21 : mean_InfectedCriticalToDead[group] -
261 21 : std::floor(-mean_InfectedSymptomsToInfectedSevere[group] -
262 21 : mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group]))) *
263 21 : entry.num_deaths;
264 : }
265 518 : if (offset ==
266 259 : std::ceil(-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] -
267 : mean_InfectedCriticalToDead[group])) {
268 21 : model.populations[0][Di] +=
269 21 : (-mean_InfectedSymptomsToInfectedSevere[group] - mean_InfectedSevereToInfectedCritical[group] -
270 21 : mean_InfectedCriticalToDead[group] -
271 21 : std::floor(-mean_InfectedSymptomsToInfectedSevere[group] -
272 21 : mean_InfectedSevereToInfectedCritical[group] - mean_InfectedCriticalToDead[group])) *
273 21 : entry.num_deaths;
274 : }
275 :
276 259 : if (offset == 0) {
277 28 : model.total_confirmed_cases[group] = scale_confirmed_cases[group] * entry.num_confirmed;
278 : }
279 : }
280 : }
281 :
282 8 : if (!max_offset_needed_avail) {
283 2 : log_error("Necessary range of dates needed to compute initial values does not exist in RKI data.");
284 2 : return failure(StatusCode::OutOfRange, "Necessary range of dates does not exist in RKI data.");
285 : }
286 6 : if (!min_offset_needed_avail) {
287 567 : auto min_date_entry = std::min_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
288 565 : return a.date < b.date;
289 : });
290 2 : auto min_date = min_date_entry->date;
291 : // Get first date that is needed.
292 2 : mio::Date min_offset_date = offset_date_by_days(date, int(min_offset_needed));
293 2 : log_warning("RKI data is needed from {} to compute initial values. RKI data is only available from {}. Missing "
294 : "dates were set to 0.",
295 : min_offset_date, min_date);
296 : }
297 :
298 : //--- Calculate the flows "after" InfectedNoSymptomsToInfectedSymptoms (that were set above using rki_data). ---
299 : // Set m_support_max_vector and m_derivative_vector in the model which is needed for the following computations.
300 6 : model.set_transitiondistributions_support_max(dt);
301 6 : model.set_transitiondistributions_derivative(dt);
302 :
303 27 : for (AgeGroup group = AgeGroup(0); group < AgeGroup(num_age_groups); group++) {
304 : //--- Calculate the flows "after" InfectedNoSymptomsToInfectedSymptoms. ---
305 : // Compute flow InfectedSymptomsToInfectedSevere for -3 * global_support_max, ..., 0.
306 294 : for (Eigen::Index i = -3 * global_support_max_index; i <= 0; i++) {
307 273 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
308 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt,
309 : i + start_shift, group);
310 : }
311 : // Compute flow InfectedSevereToInfectedCritical for -2 * global_support_max, ..., 0.
312 210 : for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
313 189 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
314 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift,
315 : group);
316 : }
317 : // Compute flows from InfectedSymptoms, InfectedSevere and InfectedCritical to Recovered and
318 : // flow InfectedCriticalToDead for -global_support_max, ..., 0.
319 126 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
320 : // Compute flow InfectedSymptomsToRecovered.
321 105 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered),
322 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt,
323 : i + start_shift, group);
324 : // Compute flow InfectedSevereToRecovered.
325 105 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered),
326 : Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift,
327 : group);
328 : // Compute flow InfectedCriticalToRecovered.
329 105 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered),
330 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift,
331 : group);
332 : // Compute flow InfectedCriticalToDead.
333 105 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead),
334 : Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift,
335 : group);
336 : }
337 :
338 : //--- Calculate the remaining flows. ---
339 : // Compute flow ExposedToInfectedNoSymptoms for -2 * global_support_max, ..., 0.
340 : // Use mean value of the TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
341 :
342 21 : Eigen::Index index_shift_mean = Eigen::Index(std::round(mean_InfectedNoSymptomsToInfectedSymptoms[group] / dt));
343 21 : int EtINSi =
344 42 : model.get_transition_flat_index(Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), group);
345 21 : int INStISyi = model.get_transition_flat_index(
346 : Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), group);
347 210 : for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
348 189 : model.transitions[i + start_shift][EtINSi] =
349 189 : (1 / model.parameters.get<TransitionProbabilities>()[group][Eigen::Index(
350 567 : InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
351 189 : model.transitions[i + start_shift + index_shift_mean][INStISyi];
352 : }
353 :
354 : // Compute flow SusceptibleToExposed for -global_support_max, ..., 0.
355 : // Use mean values of the TransitionDistribution ExposedToInfectedNoSymptoms and of the
356 : // TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
357 21 : index_shift_mean = Eigen::Index(std::round(
358 21 : (mean_ExposedToInfectedNoSymptoms[group] + mean_InfectedNoSymptomsToInfectedSymptoms[group]) / dt));
359 21 : int StEi = model.get_transition_flat_index(Eigen::Index(InfectionTransition::SusceptibleToExposed), group);
360 :
361 126 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
362 105 : model.transitions[i + start_shift][StEi] =
363 105 : (1 / model.parameters.get<TransitionProbabilities>()[group][Eigen::Index(
364 315 : InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
365 105 : model.transitions[i + start_shift + index_shift_mean][INStISyi];
366 : }
367 :
368 : // InfectedNoSymptomsToRecovered for -global_support_max, ..., 0.
369 : // If we previously calculated the transition ExposedToInfectedNoSymptoms, we can calculate this transition
370 : // using the standard formula.
371 126 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
372 105 : model.compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered),
373 : Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, i + start_shift,
374 : group);
375 : }
376 : }
377 :
378 : // At the end of the calculation, delete all time points that are not required for the simulation.
379 6 : auto transition_copy(model.transitions);
380 6 : model.transitions = TimeSeries<ScalarType>(Eigen::Index(InfectionTransition::Count) * num_age_groups);
381 36 : for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
382 30 : model.transitions.add_time_point(i * dt, transition_copy.get_value(i + start_shift));
383 : }
384 :
385 12 : return mio::success();
386 8 : }
387 :
388 : } // namespace isecir
389 : } // namespace mio
390 :
391 : #endif // MEMILIO_HAS_JSONCPP
392 :
393 : #endif // IDE_INITIALFLOWS_H
|
