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_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(Eigen::VectorX<ScalarType> const& total_population,
89 : Eigen::VectorX<ScalarType> const& deaths,
90 : Eigen::VectorX<ScalarType> const& total_confirmed_cases)
91 : {
92 :
93 13 : Eigen::VectorX<ScalarType> init(m_model.populations.get_compartments());
94 :
95 13 : bool error = compute_initialization_vector_impl(init, total_population, deaths, total_confirmed_cases);
96 13 : if (error) {
97 6 : return true;
98 : }
99 :
100 : // Update initial value vector of the model.
101 217 : for (size_t i = 0; i < m_model.populations.get_num_compartments(); i++) {
102 210 : m_model.populations[i] = init[i];
103 : }
104 :
105 : // Check if constraints are fulfilled.
106 7 : return check_constraints();
107 13 : }
108 :
109 : /**
110 : * @brief Setter for the tolerance used to calculate the maximum support of ErlangDensity%s.
111 : *
112 : * @param[in] new_tol New tolerance.
113 : */
114 3 : void set_tol_for_support_max(ScalarType new_tol)
115 : {
116 3 : m_tol = new_tol;
117 3 : }
118 :
119 : private:
120 : /**
121 : * @brief Implementation of the calculation of the initial value vector slice that corresponds to a specified group.
122 : *
123 : * Computes a vector that can be used for the initialization of an LCT model stratified by groups with the number
124 : * of persons for each subcompartment. The groups are calculated recursively.
125 : *
126 : * @tparam Group The group for which the corresponding slice of the initial value vector is calculated.
127 : * @param[out] init The initial value vector under consideration.
128 : * @param[in] total_population The total size of the considered population.
129 : * @param[in] deaths Number of deceased people from the disease at time 0.
130 : * @param[in] total_confirmed_cases Total number of confirmed cases at time 0.
131 : * @return Returns true if one (or more) constraint(s) of the computed initial value vector are not satisfied,
132 : * otherwise false.
133 : */
134 : template <size_t Group = 0>
135 22 : bool compute_initialization_vector_impl(Eigen::VectorX<ScalarType>& init,
136 : Eigen::VectorX<ScalarType> const& total_population,
137 : Eigen::VectorX<ScalarType> const& deaths,
138 : Eigen::VectorX<ScalarType> const& total_confirmed_cases)
139 : {
140 : static_assert((Group < Model::num_groups) && (Group >= 0), "The template parameter Group should be valid.");
141 : using LctStateGroup = type_at_index_t<Group, LctStatesGroups>;
142 22 : Eigen::Index first_index = m_model.populations.template get_first_index_of_group<Group>();
143 22 : Eigen::Index first_index_flows = Group * (Eigen::Index)InfectionTransition::Count;
144 :
145 22 : bool error =
146 : // Exposed.
147 44 : (compute_compartment<InfectionState::Exposed, Group>(
148 : init, first_index_flows + (Eigen::Index)InfectionTransition::SusceptibleToExposed,
149 22 : 1 / m_model.parameters.template get<TimeExposed>()[Group])) ||
150 : // InfectedNoSymptoms.
151 42 : (compute_compartment<InfectionState::InfectedNoSymptoms, Group>(
152 : init, first_index_flows + (Eigen::Index)InfectionTransition::ExposedToInfectedNoSymptoms,
153 21 : 1 / m_model.parameters.template get<TimeInfectedNoSymptoms>()[Group])) ||
154 : // InfectedSymptoms.
155 40 : (compute_compartment<InfectionState::InfectedSymptoms, Group>(
156 : init, first_index_flows + (Eigen::Index)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms,
157 20 : 1 / m_model.parameters.template get<TimeInfectedSymptoms>()[Group])) ||
158 : // InfectedSevere.
159 38 : (compute_compartment<InfectionState::InfectedSevere, Group>(
160 : init, first_index_flows + (Eigen::Index)InfectionTransition::InfectedSymptomsToInfectedSevere,
161 62 : 1 / m_model.parameters.template get<TimeInfectedSevere>()[Group])) ||
162 : // InfectedCritical.
163 36 : (compute_compartment<InfectionState::InfectedCritical, Group>(
164 : init, first_index_flows + (Eigen::Index)InfectionTransition::InfectedSevereToInfectedCritical,
165 18 : 1 / m_model.parameters.template get<TimeInfectedCritical>()[Group]));
166 22 : if (error) {
167 5 : return true;
168 : }
169 :
170 : // Recovered.
171 : // Number of recovered is equal to the cumulative number of confirmed cases minus the number of people
172 : // in the group who are infected at the moment.
173 17 : init[first_index + LctStateGroup::template get_first_index<InfectionState::Recovered>()] =
174 17 : total_confirmed_cases[Group] -
175 51 : init.segment(first_index + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>(),
176 17 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() +
177 17 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() +
178 17 : LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
179 17 : .sum() -
180 : deaths[Group];
181 :
182 : // Susceptibles.
183 17 : init[first_index + LctStateGroup::template get_first_index<InfectionState::Susceptible>()] =
184 17 : total_population[Group] -
185 34 : init.segment(first_index + LctStateGroup::template get_first_index<InfectionState::Exposed>(),
186 : LctStateGroup::Count - 2)
187 17 : .sum() -
188 : deaths[Group];
189 :
190 : // Dead.
191 17 : init[first_index + LctStateGroup::template get_first_index<InfectionState::Dead>()] = deaths[Group];
192 :
193 67 : for (size_t state_idx : {LctStateGroup::template get_first_index<InfectionState::Susceptible>(),
194 : LctStateGroup::template get_first_index<InfectionState::Recovered>(),
195 : LctStateGroup::template get_first_index<InfectionState::Dead>()}) {
196 51 : if (floating_point_less(init[first_index + state_idx], 0.0, 1e-8)) {
197 1 : log_error("Initialization failed. Values of total_confirmed_cases, deaths and total_population do not "
198 : "match the specified values for the transitions, leading to at least one negative result.");
199 1 : return true;
200 : }
201 : }
202 : // Function call for next group if applicable.
203 : if constexpr (Group + 1 < Model::num_groups) {
204 9 : return compute_initialization_vector_impl<Group + 1>(init, total_population, deaths, total_confirmed_cases);
205 : }
206 : else {
207 7 : return false;
208 : }
209 : }
210 :
211 : /**
212 : * @brief Checks constraints of the Initializer including checks for the model.
213 : * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false.
214 : */
215 7 : bool check_constraints() const
216 : {
217 7 : if (!((Eigen::Index)InfectionTransition::Count * Model::num_groups == m_flows.get_num_elements())) {
218 1 : log_error("Initial condition size does not match subcompartments.");
219 1 : return true;
220 : }
221 :
222 6 : if (!(floating_point_equal(0., m_flows.get_last_time(), 1e-8, 1e-14))) {
223 1 : log_error("Last time point in flows has to be 0.");
224 1 : return true;
225 : }
226 :
227 725 : for (int i = 1; i < m_flows.get_num_time_points(); i++) {
228 721 : if (!(floating_point_equal(m_dt, m_flows.get_time(i) - m_flows.get_time(i - 1), 1e-8, 1e-14))) {
229 1 : log_error("Time points in flows have to be equidistant.");
230 1 : return true;
231 : }
232 : }
233 :
234 4 : if (!(m_dt < 1)) {
235 1 : log_warning("Step size was set very large. The result could be distorted.");
236 1 : return true;
237 : }
238 :
239 : // Check if model is valid and the calculated initial value vector is valid.
240 3 : return m_model.check_constraints();
241 : }
242 :
243 : /**
244 : * @brief Computes a slice of the initial value vector for each subcompartment of one InfectionState
245 : * for a specified group.
246 : *
247 : * @tparam Group The group for which the corresponding slice of the initial value vector is calculated.
248 : * @tparam State The InfectionState for which the corresponding slice of the initial value vector is calculated
249 : * @param[out] init The initial value vector under consideration.
250 : * @param[in] idx_incoming_flow Index of the flow which is relevant for the calculation, so the flow
251 : * to the InfectionState.
252 : * @param[in] transition_rate Specifies the transition rate of the InfectionState.
253 : * 'This is equal to 1 / (expected time in the InfectionState).
254 :
255 : * @return Returns true if one (or more) constraint(s) of the computed slice of the initial value vector
256 : * are not satisfied, otherwise false.
257 : */
258 : template <InfectionState State, size_t Group>
259 100 : bool compute_compartment(Eigen::VectorX<ScalarType>& init, Eigen::Index idx_incoming_flow,
260 : ScalarType transition_rate) const
261 : {
262 100 : size_t first_index = m_model.populations.template get_first_index_of_group<Group>() +
263 100 : type_at_index_t<Group, LctStatesGroups>::template get_first_index<State>();
264 100 : size_t num_subcompartments = type_at_index_t<Group, LctStatesGroups>::template get_num_subcompartments<State>();
265 :
266 : // Initialize relevant density for the InfectionState.
267 : // For the first subcompartment a shape parameter of one is needed.
268 100 : ErlangDensity erlang(1, 1. / (num_subcompartments * transition_rate));
269 :
270 : // Initialize other relevant parameters.
271 100 : ScalarType calc_time{0};
272 100 : Eigen::Index calc_time_index{0};
273 100 : ScalarType state_age{0};
274 100 : ScalarType sum{0};
275 100 : Eigen::Index num_time_points = m_flows.get_num_time_points();
276 :
277 : // Calculate number of people in each subcompartment.
278 328 : for (size_t j = 0; j < num_subcompartments; j++) {
279 : // For subcompartment number j+1, shape parameter j+1 is needed.
280 233 : erlang.set_distribution_parameter((ScalarType)j + 1);
281 :
282 : // Determine relevant calculation area and corresponding index.
283 233 : calc_time = erlang.get_support_max(m_dt, m_tol);
284 233 : calc_time_index = (Eigen::Index)std::ceil(calc_time / m_dt) - 1;
285 :
286 233 : if (num_time_points < calc_time_index) {
287 1 : log_error("Initialization failed. Not enough time points for the transitions are given. More than {} "
288 : "are needed but just {} are given.",
289 : calc_time_index, num_time_points);
290 1 : return true;
291 : }
292 : else {
293 : // Approximate integral with non-standard scheme.
294 8042 : for (Eigen::Index i = num_time_points - calc_time_index; i < num_time_points; i++) {
295 7810 : state_age = (num_time_points - i) * m_dt;
296 7810 : sum += erlang.eval(state_age) * m_flows[i][idx_incoming_flow];
297 : }
298 232 : init[first_index + j] = 1. / (num_subcompartments * transition_rate) * sum;
299 232 : if (init[first_index + j] < 0) {
300 4 : log_error(
301 : "Initialization failed. Result for at least one subcompartment is less than zero. Please check "
302 : "the input values.");
303 4 : return true;
304 : }
305 : }
306 228 : sum = 0;
307 : }
308 95 : return false;
309 100 : }
310 :
311 : TimeSeries<ScalarType> m_flows; ///< TimeSeries with the flows which are used to calculate the initial vector.
312 : Model& m_model; ///< The LCT-SECIR model for which the initialization should be performed.
313 : ScalarType m_dt{}; ///< Step size of the times in m_flows and time step for the approximation of the integral.
314 : ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the ErlangDensity%s.
315 : };
316 :
317 : } // namespace lsecir
318 : } // namespace mio
319 :
320 : #endif
|
