Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 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_INITIALIZER_H
22 : #define LCTSECIR_INITIALIZER_H
23 :
24 : #include "memilio/config.h"
25 : #include "lct_secir/model.h"
26 : #include "lct_secir/infection_state.h"
27 : #include "memilio/math/eigen.h"
28 : #include "memilio/utils/time_series.h"
29 : #include "memilio/math/floating_point.h"
30 : #include "memilio/utils/logging.h"
31 : #include "memilio/utils/metaprogramming.h"
32 : #include "memilio/epidemiology/state_age_function.h"
33 :
34 : namespace mio
35 : {
36 : namespace lsecir
37 : {
38 :
39 : /**
40 : * @brief Class that can be used to compute an initialization vector out of flows for an LCT Model
41 : * with division in groups.
42 : *
43 : * The initialization method is based on the fact that the LCT model is a special case of an IDE model with special
44 : * choices of stay time distributions (so-called Erlang distributions). Accordingly, the method for calculating
45 : * initial values for the compartments from given flows/transitions is taken from the IDE-SECIR model.
46 : * We cannot use the functionality of the IDE model directly, as we have to calculate a division into sub-compartments
47 : * and not only the sizes of the compartments.
48 : * See also the IDE-SECIR model for the general method and for a better understanding of flows/transitions.
49 : *
50 : * @tparam Model is expected to be an LCT-SECIR model defined in models/lct_secir/model.h.
51 : */
52 : template <typename Model>
53 : class Initializer
54 : {
55 : public:
56 : using LctStatesGroups = typename Model::LctStatesGroups;
57 :
58 : /**
59 : * @brief Constructs a new Initializer object.
60 : *
61 : * @param[in] flows Initializing TimeSeries with flows fitting to these defined in InfectionTransition.
62 : * For each group of m_model, InfectionTransition::Count entries are required.
63 : * Timesteps should be equidistant and the values should be non-negative.
64 : * The time history has to be long enough so that it is possible to calculate the initial vector.
65 : * The length of the required time history depends on the Erlang densities used to compute the initial vector.
66 : * @param[in, out] model The LCT-SECIR model for which the initialization should be performed.
67 : */
68 13 : Initializer(TimeSeries<ScalarType>&& flows, Model& model)
69 13 : : m_flows(std::move(flows))
70 13 : , m_model(model)
71 : {
72 13 : m_dt = m_flows.get_time(1) - m_flows.get_time(0);
73 13 : }
74 :
75 : /**
76 : * @brief Core function of Initializer.
77 : *
78 : * Computes a vector that can be used for the initialization of an LCT model stratified by groups with the number
79 : * of persons for each subcompartment for each group.
80 : * The initial value vector is updated in the model.
81 : *
82 : * @param[in] total_population The total size of the considered population.
83 : * @param[in] deaths Number of deceased people from the disease at time 0.
84 : * @param[in] total_confirmed_cases Total number of confirmed cases at time 0.
85 : * @return Returns true if one (or more) constraint(s) of the model, the initial flows
86 : * or the computed initial value vector are not satisfied, otherwise false.
87 : */
88 13 : bool compute_initialization_vector(Vector<ScalarType> const& total_population, Vector<ScalarType> const& deaths,
89 : Vector<ScalarType> const& total_confirmed_cases)
90 : {
91 :
92 13 : Vector<ScalarType> init(m_model.populations.get_compartments());
93 :
94 13 : bool error = compute_initialization_vector_impl(init, total_population, deaths, total_confirmed_cases);
95 13 : if (error) {
96 6 : return true;
97 : }
98 :
99 : // Update initial value vector of the model.
100 217 : for (size_t i = 0; i < m_model.populations.get_num_compartments(); i++) {
101 210 : m_model.populations[i] = init[i];
102 : }
103 :
104 : // Check if constraints are fulfilled.
105 7 : return check_constraints();
106 13 : }
107 :
108 : /**
109 : * @brief Setter for the tolerance used to calculate the maximum support of ErlangDensity%s.
110 : *
111 : * @param[in] new_tol New tolerance.
112 : */
113 3 : void set_tol_for_support_max(ScalarType new_tol)
114 : {
115 3 : m_tol = new_tol;
116 3 : }
117 :
118 : private:
119 : /**
120 : * @brief Implementation of the calculation of the initial value vector slice that corresponds to a specified group.
121 : *
122 : * Computes a vector that can be used for the initialization of an LCT model stratified by groups with the number
123 : * of persons for each subcompartment. The groups are calculated recursively.
124 : *
125 : * @tparam Group The group for which the corresponding slice of the initial value vector is calculated.
126 : * @param[out] init The initial value vector under consideration.
127 : * @param[in] total_population The total size of the considered population.
128 : * @param[in] deaths Number of deceased people from the disease at time 0.
129 : * @param[in] total_confirmed_cases Total number of confirmed cases at time 0.
130 : * @return Returns true if one (or more) constraint(s) of the computed initial value vector are not satisfied,
131 : * otherwise false.
132 : */
133 : template <size_t Group = 0>
134 22 : bool compute_initialization_vector_impl(Vector<ScalarType>& init, Vector<ScalarType> const& total_population,
135 : Vector<ScalarType> const& deaths,
136 : Vector<ScalarType> const& total_confirmed_cases)
137 : {
138 : static_assert((Group < Model::num_groups) && (Group >= 0), "The template parameter Group should be valid.");
139 : using LctStateGroup = type_at_index_t<Group, LctStatesGroups>;
140 22 : Eigen::Index first_index = m_model.populations.template get_first_index_of_group<Group>();
141 22 : Eigen::Index first_index_flows = Group * (Eigen::Index)InfectionTransition::Count;
142 :
143 22 : bool error =
144 : // Exposed.
145 44 : (compute_compartment<InfectionState::Exposed, Group>(
146 : init, first_index_flows + (Eigen::Index)InfectionTransition::SusceptibleToExposed,
147 22 : 1 / m_model.parameters.template get<TimeExposed>()[Group])) ||
148 : // InfectedNoSymptoms.
149 42 : (compute_compartment<InfectionState::InfectedNoSymptoms, Group>(
150 : init, first_index_flows + (Eigen::Index)InfectionTransition::ExposedToInfectedNoSymptoms,
151 21 : 1 / m_model.parameters.template get<TimeInfectedNoSymptoms>()[Group])) ||
152 : // InfectedSymptoms.
153 40 : (compute_compartment<InfectionState::InfectedSymptoms, Group>(
154 : init, first_index_flows + (Eigen::Index)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
155 20 : 1 / m_model.parameters.template get<TimeInfectedSymptoms>()[Group])) ||
156 : // InfectedSevere.
157 38 : (compute_compartment<InfectionState::InfectedSevere, Group>(
158 : init, first_index_flows + (Eigen::Index)InfectionTransition::InfectedSymptomsToInfectedSevere,
159 62 : 1 / m_model.parameters.template get<TimeInfectedSevere>()[Group])) ||
160 : // InfectedCritical.
161 36 : (compute_compartment<InfectionState::InfectedCritical, Group>(
162 : init, first_index_flows + (Eigen::Index)InfectionTransition::InfectedSevereToInfectedCritical,
163 18 : 1 / m_model.parameters.template get<TimeInfectedCritical>()[Group]));
164 22 : if (error) {
165 5 : return true;
166 : }
167 :
168 : // Recovered.
169 : // Number of recovered is equal to the cumulative number of confirmed cases minus the number of people
170 : // in the group who are infected at the moment.
171 17 : init[first_index + LctStateGroup::template get_first_index<InfectionState::Recovered>()] =
172 17 : total_confirmed_cases[Group] -
173 34 : init.segment(first_index + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>(),
174 17 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() +
175 17 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() +
176 17 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
177 34 : .sum() -
178 : deaths[Group];
179 :
180 : // Susceptibles.
181 17 : init[first_index + LctStateGroup::template get_first_index<InfectionState::Susceptible>()] =
182 17 : total_population[Group] -
183 17 : init.segment(first_index + LctStateGroup::template get_first_index<InfectionState::Exposed>(),
184 : LctStateGroup::Count - 2)
185 34 : .sum() -
186 : deaths[Group];
187 :
188 : // Dead.
189 17 : init[first_index + LctStateGroup::template get_first_index<InfectionState::Dead>()] = deaths[Group];
190 :
191 67 : for (size_t state_idx : {LctStateGroup::template get_first_index<InfectionState::Susceptible>(),
192 : LctStateGroup::template get_first_index<InfectionState::Recovered>(),
193 : LctStateGroup::template get_first_index<InfectionState::Dead>()}) {
194 51 : if (floating_point_less(init[first_index + state_idx], 0.0, 1e-8)) {
195 2 : log_error("Initialization failed. Values of total_confirmed_cases, deaths and total_population do not "
196 : "match the specified values for the transitions, leading to at least one negative result.");
197 1 : return true;
198 : }
199 : }
200 : // Function call for next group if applicable.
201 : if constexpr (Group + 1 < Model::num_groups) {
202 9 : return compute_initialization_vector_impl<Group + 1>(init, total_population, deaths, total_confirmed_cases);
203 : }
204 : else {
205 7 : return false;
206 : }
207 : }
208 :
209 : /**
210 : * @brief Checks constraints of the Initializer including checks for the model.
211 : * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false.
212 : */
213 7 : bool check_constraints() const
214 : {
215 7 : if (!((Eigen::Index)InfectionTransition::Count * Model::num_groups == m_flows.get_num_elements())) {
216 2 : log_error("Initial condition size does not match subcompartments.");
217 1 : return true;
218 : }
219 :
220 6 : if (!(floating_point_equal(0., m_flows.get_last_time(), 1e-8, 1e-14))) {
221 2 : log_error("Last time point in flows has to be 0.");
222 1 : return true;
223 : }
224 :
225 725 : for (int i = 1; i < m_flows.get_num_time_points(); i++) {
226 721 : if (!(floating_point_equal(m_dt, m_flows.get_time(i) - m_flows.get_time(i - 1), 1e-8, 1e-14))) {
227 2 : log_error("Time points in flows have to be equidistant.");
228 1 : return true;
229 : }
230 : }
231 :
232 4 : if (!(m_dt < 1)) {
233 2 : log_warning("Step size was set very large. The result could be distorted.");
234 1 : return true;
235 : }
236 :
237 : // Check if model is valid and the calculated initial value vector is valid.
238 3 : return m_model.check_constraints();
239 : }
240 :
241 : /**
242 : * @brief Computes a slice of the initial value vector for each subcompartment of one InfectionState
243 : * for a specified group.
244 : *
245 : * @tparam Group The group for which the corresponding slice of the initial value vector is calculated.
246 : * @tparam State The InfectionState for which the corresponding slice of the initial value vector is calculated
247 : * @param[out] init The initial value vector under consideration.
248 : * @param[in] idx_incoming_flow Index of the flow which is relevant for the calculation, so the flow
249 : * to the InfectionState.
250 : * @param[in] transition_rate Specifies the transition rate of the InfectionState.
251 : * 'This is equal to 1 / (expected time in the InfectionState).
252 :
253 : * @return Returns true if one (or more) constraint(s) of the computed slice of the initial value vector
254 : * are not satisfied, otherwise false.
255 : */
256 : template <InfectionState State, size_t Group>
257 100 : bool compute_compartment(Vector<ScalarType>& init, Eigen::Index idx_incoming_flow, ScalarType transition_rate) const
258 : {
259 100 : size_t first_index = m_model.populations.template get_first_index_of_group<Group>() +
260 100 : type_at_index_t<Group, LctStatesGroups>::template get_first_index<State>();
261 100 : size_t num_subcompartments = type_at_index_t<Group, LctStatesGroups>::template get_num_subcompartments<State>();
262 :
263 : // Initialize relevant density for the InfectionState.
264 : // For the first subcompartment a shape parameter of one is needed.
265 100 : ErlangDensity erlang(1, 1. / (num_subcompartments * transition_rate));
266 :
267 : // Initialize other relevant parameters.
268 100 : ScalarType calc_time{0};
269 100 : Eigen::Index calc_time_index{0};
270 100 : ScalarType state_age{0};
271 100 : ScalarType sum{0};
272 100 : Eigen::Index num_time_points = m_flows.get_num_time_points();
273 :
274 : // Calculate number of people in each subcompartment.
275 328 : for (size_t j = 0; j < num_subcompartments; j++) {
276 : // For subcompartment number j+1, shape parameter j+1 is needed.
277 233 : erlang.set_distribution_parameter((ScalarType)j + 1);
278 :
279 : // Determine relevant calculation area and corresponding index.
280 233 : calc_time = erlang.get_support_max(m_dt, m_tol);
281 233 : calc_time_index = (Eigen::Index)std::ceil(calc_time / m_dt) - 1;
282 :
283 233 : if (num_time_points < calc_time_index) {
284 2 : log_error("Initialization failed. Not enough time points for the transitions are given. More than {} "
285 : "are needed but just {} are given.",
286 : calc_time_index, num_time_points);
287 1 : return true;
288 : }
289 : else {
290 : // Approximate integral with non-standard scheme.
291 9024 : for (Eigen::Index i = num_time_points - calc_time_index; i < num_time_points; i++) {
292 8792 : state_age = (num_time_points - i) * m_dt;
293 8792 : sum += erlang.eval(state_age) * m_flows[i][idx_incoming_flow];
294 : }
295 232 : init[first_index + j] = 1. / (num_subcompartments * transition_rate) * sum;
296 232 : if (init[first_index + j] < 0) {
297 8 : log_error(
298 : "Initialization failed. Result for at least one subcompartment is less than zero. Please check "
299 : "the input values.");
300 4 : return true;
301 : }
302 : }
303 228 : sum = 0;
304 : }
305 95 : return false;
306 100 : }
307 :
308 : TimeSeries<ScalarType> m_flows; ///< TimeSeries with the flows which are used to calculate the initial vector.
309 : Model& m_model; ///< The LCT-SECIR model for which the initialization should be performed.
310 : ScalarType m_dt{}; ///< Step size of the times in m_flows and time step for the approximation of the integral.
311 : ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the ErlangDensity%s.
312 : };
313 :
314 : } // namespace lsecir
315 : } // namespace mio
316 :
317 : #endif
|
