Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: 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 : #ifndef LCTSECIR_PARAMETERS_IO_H
22 : #define LCTSECIR_PARAMETERS_IO_H
23 :
24 : #include "memilio/config.h"
25 :
26 : #ifdef MEMILIO_HAS_JSONCPP
27 :
28 : #include "lct_secir/model.h"
29 : #include "lct_secir/infection_state.h"
30 : #include "lct_secir/parameters.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 : #include "memilio/utils/metaprogramming.h"
36 : #include "memilio/math/eigen.h"
37 : #include "memilio/math/floating_point.h"
38 : #include "memilio/epidemiology/lct_populations.h"
39 :
40 : #include <string>
41 : #include <vector>
42 : #include <type_traits>
43 : namespace mio
44 : {
45 : namespace lsecir
46 : {
47 :
48 : namespace details
49 : { // Use namespace to hide functions that are not intended to be used outside this file.
50 :
51 : /**
52 : * @brief Processes one entry of an RKI data set for the definition of an initial value vector for an LCT population.
53 : *
54 : * Takes one entry of an RKI data vector and changes the value in populations accordingly.
55 : * This function provides sub-functionality of the set_initial_values_from_confirmed_cases() function.
56 : *
57 : * @param[out] populations The populations for which the inital data should be computed and set.
58 : * @param[in] entry The entry of the RKI data set.
59 : * @param[in] offset The offset between the date of the entry and the date for which the
60 : * initial value vector is calculated.
61 : * @param[in] staytimes A vector of the average time spent in each compartment defined in InfectionState
62 : * (for the group under consideration) for which the initial value vector is calculated.
63 : * @param[in] inv_prob_SymptomsPerNoSymptoms The inverse of the probability InfectedSymptomsPerInfectedNoSymptoms
64 : * (for the group under consideration) for which the initial value vector is calculated.
65 : * @param[in] severePerInfectedSymptoms Probability (for the group under consideration)
66 : * for which the initial value vector is calculated.
67 : * @param[in] criticalPerSevere Probability (for the group under consideration)
68 : * for which the initial value vector is calculated.
69 : * @param[in] scale_confirmed_cases Factor by which to scale the confirmed cases of RKI data to consider unreported cases.
70 : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations.
71 : * This defines the number of age groups and the numbers of subcompartments.
72 : * @tparam EntryType The type of the data entry of the RKI data.
73 : * @tparam Group The age group of the entry the should be processed.
74 : */
75 : template <class Populations, class EntryType, size_t Group>
76 394 : void process_entry(Populations& populations, const EntryType& entry, int offset,
77 : const std::vector<ScalarType> staytimes, ScalarType inv_prob_SymptomsPerNoSymptoms,
78 : ScalarType severePerInfectedSymptoms, ScalarType criticalPerSevere, ScalarType scale_confirmed_cases)
79 : {
80 : using LctStateGroup = type_at_index_t<Group, typename Populations::LctStatesGroups>;
81 394 : size_t first_index_group = populations.template get_first_index_of_group<Group>();
82 :
83 : // Add confirmed cases at date to compartment Recovered.
84 394 : if (offset == 0) {
85 72 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] +=
86 36 : scale_confirmed_cases * entry.num_confirmed;
87 : }
88 :
89 : // Compute initial values for compartment InfectedNoSymptoms.
90 394 : if (offset >= 0 && offset <= std::ceil(staytimes[(size_t)InfectionState::InfectedNoSymptoms])) {
91 : // Mean stay time in each subcompartment.
92 108 : ScalarType time_INS_per_subcomp =
93 108 : staytimes[(size_t)InfectionState::InfectedNoSymptoms] /
94 : (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
95 : // Index of the last subcompartment of InfectedNoSymptoms.
96 108 : size_t last_index_INS =
97 108 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedNoSymptoms>() +
98 108 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() - 1;
99 360 : for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
100 : i++) {
101 252 : if (offset == std::floor(i * time_INS_per_subcomp)) {
102 168 : populations[last_index_INS - i] -=
103 84 : (1 - (i * time_INS_per_subcomp - std::floor(i * time_INS_per_subcomp))) *
104 84 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
105 : }
106 252 : if (offset == std::ceil(i * time_INS_per_subcomp)) {
107 168 : populations[last_index_INS - i] -= (i * time_INS_per_subcomp - std::floor(i * time_INS_per_subcomp)) *
108 168 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases *
109 84 : entry.num_confirmed;
110 : }
111 252 : if (offset == std::floor((i + 1) * time_INS_per_subcomp)) {
112 168 : populations[last_index_INS - i] +=
113 84 : (1 - ((i + 1) * time_INS_per_subcomp - std::floor((i + 1) * time_INS_per_subcomp))) *
114 84 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
115 : }
116 252 : if (offset == std::ceil((i + 1) * time_INS_per_subcomp)) {
117 168 : populations[last_index_INS - i] +=
118 84 : ((i + 1) * time_INS_per_subcomp - std::floor((i + 1) * time_INS_per_subcomp)) *
119 84 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
120 : }
121 : }
122 : }
123 :
124 : // Compute initial values for compartment Exposed.
125 536 : if (offset >= std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms]) &&
126 142 : offset <= std::ceil(staytimes[(size_t)InfectionState::Exposed] +
127 : staytimes[(size_t)InfectionState::InfectedNoSymptoms])) {
128 : // Mean stay time in each subcompartment.
129 142 : ScalarType time_Exposed_per_subcomp =
130 142 : staytimes[(size_t)InfectionState::Exposed] /
131 : (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>();
132 : // Index of the last subcompartment of Exposed.
133 284 : size_t last_index_Exposed = first_index_group +
134 142 : LctStateGroup::template get_first_index<InfectionState::Exposed>() +
135 142 : LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>() - 1;
136 378 : for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>(); i++) {
137 472 : if (offset ==
138 236 : std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp)) {
139 60 : populations[last_index_Exposed - i] -=
140 60 : (1 - (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp -
141 60 : std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
142 60 : i * time_Exposed_per_subcomp))) *
143 60 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
144 : }
145 472 : if (offset ==
146 236 : std::ceil(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp)) {
147 60 : populations[last_index_Exposed - i] -=
148 60 : (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp -
149 60 : std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp)) *
150 60 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
151 : }
152 236 : if (offset == std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
153 236 : (i + 1) * time_Exposed_per_subcomp)) {
154 60 : populations[last_index_Exposed - i] +=
155 60 : (1 - (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + (i + 1) * time_Exposed_per_subcomp -
156 60 : std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
157 60 : (i + 1) * time_Exposed_per_subcomp))) *
158 60 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
159 : }
160 472 : if (offset ==
161 236 : std::ceil(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + (i + 1) * time_Exposed_per_subcomp)) {
162 58 : populations[last_index_Exposed - i] +=
163 58 : (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + (i + 1) * time_Exposed_per_subcomp -
164 58 : std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
165 58 : (i + 1) * time_Exposed_per_subcomp)) *
166 58 : inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
167 : }
168 : }
169 : }
170 :
171 : // Compute initial values for compartment InfectedSymptoms.
172 394 : if (offset >= std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms]) && offset <= 0) {
173 : // Mean stay time in each subcompartment.
174 144 : ScalarType time_IS_per_subcomp =
175 144 : staytimes[(size_t)InfectionState::InfectedSymptoms] /
176 : (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>();
177 : // Index of the first subcompartment of InfectedSymptoms.
178 144 : size_t first_index_IS =
179 144 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>();
180 384 : for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>();
181 : i++) {
182 240 : if (offset == std::floor(-time_IS_per_subcomp * (i + 1))) {
183 120 : populations[first_index_IS + i] -=
184 60 : (1 - (-(i + 1.) * time_IS_per_subcomp - std::floor(-(i + 1) * time_IS_per_subcomp))) *
185 60 : scale_confirmed_cases * entry.num_confirmed;
186 : }
187 240 : if (offset == std::ceil(-time_IS_per_subcomp * (i + 1))) {
188 120 : populations[first_index_IS + i] -=
189 60 : (-(i + 1) * time_IS_per_subcomp - std::floor(-(i + 1) * time_IS_per_subcomp)) *
190 60 : scale_confirmed_cases * entry.num_confirmed;
191 : }
192 240 : if (offset == std::floor(-time_IS_per_subcomp * i)) {
193 120 : populations[first_index_IS + i] +=
194 120 : (1 - (-i * time_IS_per_subcomp - std::floor(-i * time_IS_per_subcomp))) * scale_confirmed_cases *
195 60 : entry.num_confirmed;
196 : }
197 240 : if (offset == std::ceil(-time_IS_per_subcomp * i)) {
198 120 : populations[first_index_IS + i] += (-i * time_IS_per_subcomp - std::floor(-i * time_IS_per_subcomp)) *
199 60 : scale_confirmed_cases * entry.num_confirmed;
200 : }
201 : }
202 : }
203 :
204 : // Compute initial values for compartment InfectedSevere.
205 394 : if (offset >= std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
206 752 : staytimes[(size_t)InfectionState::InfectedSevere]) &&
207 358 : offset <= std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms])) {
208 : // Mean stay time in each subcompartment.
209 144 : ScalarType time_ISev_per_subcomp =
210 144 : staytimes[(size_t)InfectionState::InfectedSevere] /
211 : (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>();
212 : // Index of the first subcompartment of InfectedSevere.
213 144 : size_t first_index_ISev =
214 144 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSevere>();
215 384 : for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>();
216 : i++) {
217 480 : if (offset ==
218 240 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * (i + 1))) {
219 60 : populations[first_index_ISev + i] -=
220 60 : severePerInfectedSymptoms *
221 60 : (1 - (-staytimes[(size_t)InfectionState::InfectedSymptoms] - (i + 1) * time_ISev_per_subcomp -
222 60 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
223 60 : (i + 1) * time_ISev_per_subcomp))) *
224 60 : scale_confirmed_cases * entry.num_confirmed;
225 : }
226 480 : if (offset ==
227 240 : std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * (i + 1))) {
228 60 : populations[first_index_ISev + i] -=
229 60 : severePerInfectedSymptoms *
230 60 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] - (i + 1) * time_ISev_per_subcomp -
231 60 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
232 60 : (i + 1) * time_ISev_per_subcomp)) *
233 60 : scale_confirmed_cases * entry.num_confirmed;
234 : }
235 480 : if (offset ==
236 240 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * i)) {
237 60 : populations[first_index_ISev + i] +=
238 60 : severePerInfectedSymptoms *
239 60 : (1 -
240 60 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp -
241 60 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp))) *
242 60 : scale_confirmed_cases * entry.num_confirmed;
243 : }
244 240 : if (offset == std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * i)) {
245 60 : populations[first_index_ISev + i] +=
246 60 : severePerInfectedSymptoms *
247 60 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp -
248 60 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp)) *
249 60 : scale_confirmed_cases * entry.num_confirmed;
250 : }
251 : }
252 : }
253 :
254 : // Compute initial values for compartment InfectedCritical.
255 788 : if (offset >= std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
256 394 : staytimes[(size_t)InfectionState::InfectedSevere] -
257 788 : staytimes[(size_t)InfectionState::InfectedCritical]) &&
258 394 : offset <= std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
259 : staytimes[(size_t)InfectionState::InfectedSevere])) {
260 : // Mean stay time in each subcompartment.
261 108 : ScalarType time_ICri_per_subcomp =
262 108 : staytimes[(size_t)InfectionState::InfectedCritical] /
263 : (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
264 : // Index of the first subcompartment of InfectedCritical.
265 108 : size_t first_index_ICri =
266 108 : first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedCritical>();
267 276 : for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
268 : i++) {
269 336 : if (offset ==
270 168 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
271 168 : staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * (i + 1))) {
272 56 : populations[first_index_ICri + i] -=
273 112 : severePerInfectedSymptoms * criticalPerSevere *
274 56 : (1 - (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
275 112 : staytimes[(size_t)InfectionState::InfectedSevere] - (i + 1) * time_ICri_per_subcomp -
276 56 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
277 : staytimes[(size_t)InfectionState::InfectedSevere] -
278 56 : (i + 1) * time_ICri_per_subcomp))) *
279 56 : scale_confirmed_cases * entry.num_confirmed;
280 : }
281 336 : if (offset ==
282 168 : std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
283 168 : staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * (i + 1))) {
284 56 : populations[first_index_ICri + i] -=
285 112 : severePerInfectedSymptoms * criticalPerSevere *
286 56 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
287 112 : staytimes[(size_t)InfectionState::InfectedSevere] - (i + 1) * time_ICri_per_subcomp -
288 56 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
289 56 : staytimes[(size_t)InfectionState::InfectedSevere] - (i + 1) * time_ICri_per_subcomp)) *
290 56 : scale_confirmed_cases * entry.num_confirmed;
291 : }
292 168 : if (offset == std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
293 168 : staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * i)) {
294 56 : populations[first_index_ICri + i] +=
295 112 : severePerInfectedSymptoms * criticalPerSevere *
296 56 : (1 - (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
297 112 : staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp -
298 56 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
299 56 : staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp))) *
300 56 : scale_confirmed_cases * entry.num_confirmed;
301 : }
302 168 : if (offset == std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
303 168 : staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * i)) {
304 56 : populations[first_index_ICri + i] +=
305 112 : severePerInfectedSymptoms * criticalPerSevere *
306 56 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
307 112 : staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp -
308 56 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
309 56 : staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp)) *
310 56 : scale_confirmed_cases * entry.num_confirmed;
311 : }
312 : }
313 : }
314 :
315 : // Compute Dead by shifting RKI data according to mean stay times.
316 : // This is because the RKI reports deaths with the date of the positive test, not the date of death.
317 788 : if (offset == std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
318 394 : staytimes[(size_t)InfectionState::InfectedSevere] -
319 : staytimes[(size_t)InfectionState::InfectedCritical])) {
320 36 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>()] +=
321 36 : (1 -
322 36 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] - staytimes[(size_t)InfectionState::InfectedSevere] -
323 36 : staytimes[(size_t)InfectionState::InfectedCritical] -
324 72 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
325 36 : staytimes[(size_t)InfectionState::InfectedSevere] -
326 36 : staytimes[(size_t)InfectionState::InfectedCritical]))) *
327 36 : entry.num_deaths;
328 : }
329 788 : if (offset == std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
330 394 : staytimes[(size_t)InfectionState::InfectedSevere] -
331 : staytimes[(size_t)InfectionState::InfectedCritical])) {
332 36 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>()] +=
333 36 : (-staytimes[(size_t)InfectionState::InfectedSymptoms] - staytimes[(size_t)InfectionState::InfectedSevere] -
334 36 : staytimes[(size_t)InfectionState::InfectedCritical] -
335 72 : std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
336 36 : staytimes[(size_t)InfectionState::InfectedSevere] -
337 36 : staytimes[(size_t)InfectionState::InfectedCritical])) *
338 36 : entry.num_deaths;
339 : }
340 394 : }
341 :
342 : /**
343 : * @brief Computes an initialization vector for an LCT population with case data from RKI recursively for each age group
344 : * (or for one age group in the case without age resolution).
345 : *
346 : * Please use the set_initial_values_from_reported_data() function, which calls this function automatically!
347 : * This function calculates a segment referring to the defined age group of the initial value vector with
348 : * subcompartments using the rki_data and the parameters.
349 : * The values for the whole initial value vector stored in populations are calculated recursively.
350 : *
351 : * @param[in] rki_data Vector with the RKI data.
352 : * @param[out] populations The populations for which the inital data should be computed and set.
353 : * @param[in] parameters The parameters that should be used to calculate the initial values.
354 : * Probabilities and mean stay times are used.
355 : * @param[in] date Date for which the initial values should be computed. date is the start date of the simulation.
356 : * @param[in] total_population Total size of the population of Germany or of every age group.
357 : * @param[in] scale_confirmed_cases Factor(s for each age group) by which to scale the confirmed cases of the rki data
358 : * to consider unreported cases.
359 : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations.
360 : * This defines the number of age groups and the numbers of subcompartments.
361 : * @tparam EntryType is expected to be ConfirmedCasesNoAgeEntry for data that is not age resolved and
362 : * ConfirmedCasesDataEntry for age resolved data. See also epi_data.h.
363 : * @tparam Group The age group for which the initial values should be calculated. The function is called recursively
364 : * such that the initial values are calculated for every age group if Group is zero at the beginning.
365 : * @returns Any io errors that happen during data processing.
366 : */
367 : template <class Populations, class EntryType, size_t Group = 0>
368 36 : IOResult<void> set_initial_values_from_confirmed_cases(Populations& populations, const std::vector<EntryType>& rki_data,
369 : const Parameters& parameters, const Date date,
370 : const std::vector<ScalarType>& total_population,
371 : const std::vector<ScalarType>& scale_confirmed_cases)
372 : {
373 : static_assert((Group < Populations::num_groups) && (Group >= 0), "The template parameter Group should be valid.");
374 : using LctStateGroup = type_at_index_t<Group, typename Populations::LctStatesGroups>;
375 36 : size_t first_index_group = populations.template get_first_index_of_group<Group>();
376 :
377 : // Define variables for parameters.
378 72 : std::vector<ScalarType> staytimes((size_t)InfectionState::Count, -1.);
379 36 : staytimes[(size_t)InfectionState::Exposed] = parameters.template get<TimeExposed>()[Group];
380 36 : staytimes[(size_t)InfectionState::InfectedNoSymptoms] = parameters.template get<TimeInfectedNoSymptoms>()[Group];
381 36 : staytimes[(size_t)InfectionState::InfectedSymptoms] = parameters.template get<TimeInfectedSymptoms>()[Group];
382 36 : staytimes[(size_t)InfectionState::InfectedSevere] = parameters.template get<TimeInfectedSevere>()[Group];
383 36 : staytimes[(size_t)InfectionState::InfectedCritical] = parameters.template get<TimeInfectedCritical>()[Group];
384 36 : ScalarType inv_prob_SymptomsPerNoSymptoms =
385 36 : 1 / (1 - parameters.template get<RecoveredPerInfectedNoSymptoms>()[Group]);
386 36 : ScalarType prob_SeverePerInfectedSymptoms = parameters.template get<SeverePerInfectedSymptoms>()[Group];
387 36 : ScalarType prob_CriticalPerSevere = parameters.template get<CriticalPerSevere>()[Group];
388 :
389 72 : ScalarType min_offset_needed = std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
390 36 : staytimes[(size_t)InfectionState::InfectedSevere] -
391 : staytimes[(size_t)InfectionState::InfectedCritical]);
392 : ScalarType max_offset_needed =
393 36 : std::ceil(staytimes[(size_t)InfectionState::Exposed] + staytimes[(size_t)InfectionState::InfectedNoSymptoms]);
394 :
395 36 : bool min_offset_needed_avail = false;
396 36 : bool max_offset_needed_avail = false;
397 :
398 : // Go through the entries of rki_data and check if the entry is age resolved and is referring to
399 : // the age_group Group in the case with age resolution. If the date is
400 : // needed for calculation, another function to handle the entry is called. Confirmed cases are scaled by scale_confirmed_cases.
401 2388 : for (auto&& entry : rki_data) {
402 : if constexpr (std::is_same_v<EntryType, ConfirmedCasesDataEntry>) {
403 2304 : if (!((size_t)entry.age_group == Group)) {
404 1920 : continue;
405 : }
406 : }
407 432 : int offset = get_offset_in_days(entry.date, date);
408 432 : if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) {
409 394 : if (offset == max_offset_needed) {
410 34 : max_offset_needed_avail = true;
411 : }
412 394 : if (offset == min_offset_needed) {
413 36 : min_offset_needed_avail = true;
414 : }
415 394 : process_entry<Populations, EntryType, Group>(populations, entry, offset, staytimes,
416 : inv_prob_SymptomsPerNoSymptoms, prob_SeverePerInfectedSymptoms,
417 : prob_CriticalPerSevere, scale_confirmed_cases[Group]);
418 : }
419 : }
420 :
421 : // Compute Recovered.
422 36 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] -=
423 : populations.get_compartments()
424 108 : .segment(first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>(),
425 36 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() +
426 36 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() +
427 36 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
428 72 : .sum();
429 36 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] -=
430 36 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>()];
431 :
432 : // Compute Susceptibles.
433 72 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Susceptible>()] =
434 36 : total_population[Group] -
435 72 : populations.get_compartments().segment(first_index_group + 1, LctStateGroup::Count - 1).sum();
436 :
437 : // Check whether all necessary dates are available.
438 36 : if (!max_offset_needed_avail || !min_offset_needed_avail) {
439 2 : log_error(
440 : "Necessary range of dates needed to compute initial values does not exist in RKI data for group {:d}.",
441 : Group);
442 : return failure(StatusCode::OutOfRange,
443 2 : "Necessary range of dates needed to compute initial values does not exist in RKI data.");
444 : }
445 :
446 : // Check if all values for populations are valid.
447 130 : for (size_t i = first_index_group; i < LctStateGroup::Count; i++) {
448 98 : if (floating_point_less<ScalarType>((ScalarType)populations[i], 0., Limits<ScalarType>::zero_tolerance())) {
449 2 : log_error("Something went wrong in the initialization of group {:d}. At least one entry is negative.",
450 : Group);
451 : return failure(StatusCode::InvalidValue,
452 2 : "Something went wrong in the initialization as at least one entry is negative.");
453 : }
454 : }
455 :
456 : if constexpr (Group + 1 < Populations::num_groups) {
457 : return set_initial_values_from_confirmed_cases<Populations, EntryType, Group + 1>(
458 25 : populations, rki_data, parameters, date, total_population, scale_confirmed_cases);
459 : }
460 : else {
461 14 : return success();
462 : }
463 36 : }
464 :
465 : /**
466 : * @brief Computes the total number of patients in Intensive Care Units (in all groups and subcompartments)
467 : * in the provided Population.
468 : *
469 : * This function calculates the total number of individuals within the compartment InfectedCritical
470 : * for the Populations-data provided, irrespective of their subcompartment or group.
471 : * This total number can be used to scale the entries so that the total number in InfectedCritical is equal to
472 : * the number of ICU patients reported in the DIVI data.
473 : * Please use the set_initial_values_from_reported_data() function, which calls this function automatically!
474 : *
475 : * @param[in] populations The populations for which the total number in InfectedCritical should be computed.
476 : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations.
477 : * This defines the number of age groups and the numbers of subcompartments.
478 : * @tparam Group The age group for which the total number should be calculated. The function is called recursively
479 : * such that the total number in InfectedCritical within all groups is calculated if Group is zero at the beginning.
480 : * @returns The total number of patients in Intensive Care Units (in all groups and subcompartments).
481 : */
482 : template <class Populations, size_t Group = 0>
483 24 : ScalarType get_total_InfectedCritical_from_populations(const Populations& populations)
484 : {
485 : using LctStateGroup = type_at_index_t<Group, typename Populations::LctStatesGroups>;
486 24 : size_t first_index_group = populations.template get_first_index_of_group<Group>();
487 :
488 24 : ScalarType infectedCritical_Group =
489 : populations.get_compartments()
490 48 : .segment(first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedCritical>(),
491 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
492 24 : .sum();
493 :
494 : if constexpr (Group + 1 < Populations::num_groups) {
495 : return infectedCritical_Group +
496 20 : get_total_InfectedCritical_from_populations<Populations, Group + 1>(populations);
497 : }
498 : else {
499 4 : return infectedCritical_Group;
500 : }
501 : }
502 :
503 : /**
504 : * @brief Extract the reported number of patients in ICU for a specific date from DIVI data.
505 : *
506 : * @param[in] divi_data Vector with reported DIVI data.
507 : * @param[in] date Date for which the reported number of patients in ICU should be extracted.
508 : * @returns The reported number of patients in ICU or any io errors that happen during data processing.
509 : */
510 4 : IOResult<ScalarType> get_icu_from_divi_data(const std::vector<DiviEntry>& divi_data, const Date date)
511 : {
512 28 : for (auto&& entry : divi_data) {
513 27 : int offset = get_offset_in_days(entry.date, date);
514 27 : if (offset == 0) {
515 3 : return success(entry.num_icu);
516 : }
517 : }
518 1 : log_error("Specified date does not exist in DIVI data.");
519 1 : return failure(StatusCode::OutOfRange, "Specified date does not exist in DIVI data.");
520 : }
521 :
522 : /**
523 : * @brief Rescales the entries for InfectedCritical in populations such that the total number
524 : * equals the reported number.
525 : *
526 : * This function rescales the entries for InfectedCritical in the given population for every group and subcompartment
527 : * such that the total number in all InfectedCritical compartments equals the reported number infectedCritical_reported.
528 : *
529 : * If the total number of individuals in InfectedCritical in populations is zero and the reported number is not,
530 : * the reported number is distributed uniformly across the groups.
531 : * Within the groups, the number is distributed uniformly to the subcompartments.
532 : * Note that especially the uniform distribution across groups is not necessarily realistic,
533 : * because the need for intensive care can differ by group.
534 : *
535 : * @param[in,out] populations The populations for which the entries of the InfectedCritical compartments are rescaled.
536 : * @param[in] infectedCritical_reported The reported number for patients in ICU. The total number of individuals in the
537 : * InfectedCritical compartment in populations will be equal to this value afterward.
538 : * You can calculate this value with the get_icu_from_divi_data() function.
539 : * @param[in] infectedCritical_populations The current total number of individuals in the InfectedCritical compartment
540 : * in populations. You can calculate this value with the get_total_InfectedCritical_from_populations() function.
541 : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations.
542 : * This defines the number of age groups and the numbers of subcompartments.
543 : * @tparam Group The age group for which the entries of InfectedCritical should be scaled.
544 : * The function is called recursively for the groups. The total number in the InfectedCritical compartments is only
545 : * equal to infectedCritical_reported after the function call if Group is set to zero in the beginning.
546 : * @returns Any io errors that happen during the scaling.
547 : */
548 : template <class Populations, size_t Group = 0>
549 14 : IOResult<void> rescale_to_divi_data(Populations& populations, const ScalarType infectedCritical_reported,
550 : const ScalarType infectedCritical_populations)
551 : {
552 14 : if (floating_point_less<ScalarType>(infectedCritical_reported, 0., Limits<ScalarType>::zero_tolerance())) {
553 1 : log_error("The provided reported number of InfectedCritical is negative. Please check the data.");
554 : return failure(StatusCode::InvalidValue,
555 1 : "The provided reported number of InfectedCritical is negative. Please check the data.");
556 : }
557 :
558 : using LctStateGroup = type_at_index_t<Group, typename Populations::LctStatesGroups>;
559 13 : size_t first_index_group = populations.template get_first_index_of_group<Group>();
560 :
561 13 : if (floating_point_equal<ScalarType>(infectedCritical_populations, 0., Limits<ScalarType>::zero_tolerance())) {
562 7 : if (!(floating_point_equal<ScalarType>(infectedCritical_reported, 0., Limits<ScalarType>::zero_tolerance()))) {
563 7 : log_info("The calculated number of patients in intensive care is zero, although the reported number is "
564 : "not. The reported number is uniformly distributed across groups and subcompartments. Note "
565 : "that this is not necessarily realistic.");
566 7 : size_t num_InfectedCritical =
567 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
568 7 : ScalarType num_age_groups = (ScalarType)Populations::num_groups;
569 : // Distribute reported number uniformly to age groups and subcompartments.
570 30 : for (size_t subcompartment = 0;
571 30 : subcompartment < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
572 : subcompartment++) {
573 23 : populations[first_index_group +
574 23 : LctStateGroup::template get_first_index<InfectionState::InfectedCritical>() +
575 23 : subcompartment] =
576 23 : infectedCritical_reported / ((ScalarType)num_InfectedCritical * num_age_groups);
577 : }
578 : // Adjust Recovered compartment.
579 7 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] -=
580 7 : infectedCritical_reported / num_age_groups;
581 : // Number of Susceptibles is not affected because Recovered is adjusted accordingly.
582 : }
583 : }
584 : else {
585 : // Adjust number of Recovered by adding the old number in InfectedCritical
586 : // and subtracting the new number (= scaling_factor * old number).
587 6 : ScalarType scaling_factor_infectedCritical = infectedCritical_reported / infectedCritical_populations;
588 6 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] +=
589 6 : (1 - scaling_factor_infectedCritical) *
590 : populations.get_compartments()
591 18 : .segment(first_index_group +
592 6 : LctStateGroup::template get_first_index<InfectionState::InfectedCritical>(),
593 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
594 6 : .sum();
595 : // Adjust InfectedCritical.
596 6 : for (size_t subcompartment = 0;
597 15 : subcompartment < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
598 : subcompartment++) {
599 18 : populations[first_index_group +
600 18 : LctStateGroup::template get_first_index<InfectionState::InfectedCritical>() + subcompartment] *=
601 : scaling_factor_infectedCritical;
602 : }
603 : // Number of Susceptibles is not affected because Recovered is adjusted accordingly.
604 : }
605 26 : if (floating_point_less<ScalarType>(
606 13 : (ScalarType)
607 13 : populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()],
608 : 0., Limits<ScalarType>::zero_tolerance())) {
609 1 : log_error(
610 : "Scaling with reported DIVI data led to a negative entry in the Recovered compartment for group {:d}.",
611 : Group);
612 : return failure(StatusCode::InvalidValue,
613 1 : "Scaling with reported DIVI data led to a negative entry in a Recovered compartment.");
614 : }
615 : if constexpr (Group + 1 < Populations::num_groups) {
616 : return rescale_to_divi_data<Populations, Group + 1>(populations, infectedCritical_reported,
617 10 : infectedCritical_populations);
618 : }
619 : else {
620 4 : return success();
621 : }
622 : }
623 : } // namespace details
624 :
625 : /**
626 : * @brief Computes an initialization vector for an LCT population with case data from RKI (and possibly DIVI data).
627 : *
628 : * Use just one group in the definition of the populations to not divide between age groups.
629 : * Otherwise, the number of groups has to match the number of RKI age groups.
630 : * The function calculates an initial value vector referring to an LCT population and updates the initial value vector
631 : * in the populations class.
632 : * For the computation expected stay times in the subcompartments defined in the parameters variable are used.
633 : * To calculate the initial values, we assume for simplicity that individuals stay in the subcompartment
634 : * for exactly the expected time.
635 : * The RKI data are linearly interpolated within one day to match the expected stay time in a subcompartment.
636 : * The RKI data should contain data for each needed day with or without division of age groups,
637 : * the completeness of the dates is not verified.
638 : * Data can be downloaded e.g. with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, which creates files
639 : * named e.g. cases_all_germany.json for no groups or cases_all_age.json with division in age groups or similar names.
640 : * One should set impute_dates=True so that missing dates are imputed.
641 : * To read the data into a vector, use the functionality from epi_data.h.
642 : * The data and the number of entries in the total_population and scale_confirmed_cases vectors have to match the
643 : * number of groups used in Populations.
644 : *
645 : * Additionally, one can scale the result from the calculation with the RKI data to match the reported number of
646 : * patients in ICUs. The patient numbers are provided by DIVI and can be downloaded e.g. using
647 : * pycode/memilio-epidata/memilio/epidata/getDIVIData.py (One should also set impute_dates=True so that missing dates
648 : * are imputed.). Again, to read the data into a vector, use the functionality from epi_data.h.
649 : *
650 : * @param[in] rki_data Vector with the RKI data.
651 : * @param[out] populations The populations for which the inital data should be computed and set.
652 : * @param[in] parameters The parameters that should be used to calculate the initial values.
653 : * Probabilities and mean stay times are used.
654 : * @param[in] date Date for which the initial values should be computed. date is the start date of the simulation.
655 : * @param[in] total_population Total size of the population of Germany or of every age group.
656 : * @param[in] scale_confirmed_cases Factor(s for each age group) by which to scale the confirmed cases of the rki data
657 : * to consider unreported cases.
658 : * @param[in] divi_data Vector with DIVI data used to scale the number of individuals in the InfectedCritical
659 : * compartments in populations so that the total number match the reported number.
660 : * For the default value (an empty vector), the calculated populations using the RKI data is not scaled.
661 : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations.
662 : * This defines the number of age groups and the numbers of subcompartments.
663 : * @tparam EntryType is expected to be ConfirmedCasesNoAgeEntry for data that is not age resolved and
664 : * ConfirmedCasesDataEntry for age resolved data. See also epi_data.h.
665 : * @returns Any io errors that happen during data processing.
666 : */
667 : template <class Populations, class EntryType>
668 15 : IOResult<void> set_initial_values_from_reported_data(const std::vector<EntryType>& rki_data, Populations& populations,
669 : const Parameters& parameters, const Date date,
670 : const std::vector<ScalarType>& total_population,
671 : const std::vector<ScalarType>& scale_confirmed_cases,
672 : const std::vector<DiviEntry>& divi_data = std::vector<DiviEntry>())
673 : { // Check if the inputs are matching.
674 15 : assert(total_population.size() == Populations::num_groups);
675 15 : assert(scale_confirmed_cases.size() == Populations::num_groups);
676 : if constexpr (Populations::num_groups > 1) {
677 : static_assert(std::is_same_v<EntryType, ConfirmedCasesDataEntry>);
678 9 : assert(ConfirmedCasesDataEntry::age_group_names.size() == Populations::num_groups);
679 : }
680 : else {
681 : static_assert(std::is_same_v<EntryType, ConfirmedCasesNoAgeEntry>);
682 : }
683 : // Check if RKI data vector is valid.
684 638 : auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
685 623 : return a.date < b.date;
686 : });
687 15 : if (max_date_entry == rki_data.end()) {
688 2 : log_error("RKI data file is empty.");
689 2 : return failure(StatusCode::InvalidFileFormat, "RKI data is empty.");
690 : }
691 13 : auto max_date = max_date_entry->date;
692 13 : if (max_date < date) {
693 2 : log_error("Specified date does not exist in RKI data.");
694 2 : return failure(StatusCode::OutOfRange, "Specified date does not exist in RKI data.");
695 : }
696 : // Initially set populations to zero.
697 579 : for (size_t i = 0; i < populations.get_num_compartments(); i++) {
698 568 : populations[i] = 0;
699 : }
700 : // Set populations using the RKI data.
701 11 : IOResult<void> ioresult_confirmedcases = details::set_initial_values_from_confirmed_cases<Populations, EntryType>(
702 : populations, rki_data, parameters, date, total_population, scale_confirmed_cases);
703 11 : if (!(ioresult_confirmedcases)) {
704 4 : return ioresult_confirmedcases;
705 : }
706 :
707 : // Check if DIVI data is provided and scale the result in populations accordingly.
708 7 : if (!divi_data.empty()) {
709 : ScalarType infectedCritical_populations =
710 2 : details::get_total_InfectedCritical_from_populations<Populations>(populations);
711 2 : auto infectedCritical_reported = details::get_icu_from_divi_data(divi_data, date);
712 2 : if (!(infectedCritical_reported)) {
713 1 : return infectedCritical_reported.error();
714 : }
715 1 : return details::rescale_to_divi_data<Populations>(populations, infectedCritical_reported.value(),
716 1 : infectedCritical_populations);
717 2 : }
718 :
719 10 : return success();
720 11 : }
721 :
722 : } // namespace lsecir
723 : } // namespace mio
724 : #endif // MEMILIO_HAS_JSONCPP
725 :
726 : #endif // LCTSECIR_PARAMETERS_IO_H
|
