Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Anna Wendler, 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 : #ifndef IDE_SECIR_PARAMS_H
21 : #define IDE_SECIR_PARAMS_H
22 :
23 : #include "memilio/config.h"
24 : #include "memilio/math/floating_point.h"
25 : #include "memilio/utils/custom_index_array.h"
26 : #include "memilio/utils/parameter_set.h"
27 : #include "ide_secir/infection_state.h"
28 : #include "memilio/math/eigen.h"
29 : #include "memilio/math/smoother.h"
30 : #include "memilio/epidemiology/state_age_function.h"
31 : #include "memilio/epidemiology/uncertain_matrix.h"
32 : #include "memilio/epidemiology/age_group.h"
33 :
34 : #include <memory>
35 : #include <cstddef>
36 : #include <vector>
37 :
38 : namespace mio
39 : {
40 : namespace isecir
41 : {
42 :
43 : /**********************************************
44 : * Define Parameters of the IDE-SECIHURD model *
45 : **********************************************/
46 :
47 : /**
48 : * @brief Transition distribution for each transition in #InfectionTransition.
49 : *
50 : * As a default we use SmootherCosine functions for all transitions with m_parameter=2.
51 : */
52 : struct TransitionDistributions {
53 :
54 : using Type = CustomIndexArray<std::vector<StateAgeFunctionWrapper>, AgeGroup>;
55 22 : static Type get_default(AgeGroup size)
56 : {
57 22 : SmootherCosine smoothcos(2.0);
58 22 : StateAgeFunctionWrapper delaydistribution(smoothcos);
59 22 : std::vector<StateAgeFunctionWrapper> state_age_function_vector((int)InfectionTransition::Count,
60 44 : delaydistribution);
61 44 : return Type(size, state_age_function_vector);
62 22 : }
63 :
64 : static std::string name()
65 : {
66 : return "TransitionDistributions";
67 : }
68 : };
69 :
70 : /**
71 : * @brief Defines the probability for each possible transition to take this flow/transition.
72 : */
73 : struct TransitionProbabilities {
74 : /*For consistency, also define TransitionProbabilities for each transition in #InfectionTransition.
75 : Transition Probabilities should be set to 1 if there is no possible other flow from starting compartment.*/
76 : using Type = CustomIndexArray<std::vector<ScalarType>, AgeGroup>;
77 22 : static Type get_default(AgeGroup size)
78 : {
79 44 : std::vector<ScalarType> probs((int)InfectionTransition::Count, 0.5);
80 : // Set the following probablities to 1 as there is no other option to go anywhere else.
81 22 : probs[Eigen::Index(InfectionTransition::SusceptibleToExposed)] = 1;
82 22 : probs[Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] = 1;
83 44 : return Type(size, probs);
84 22 : }
85 :
86 : static std::string name()
87 : {
88 : return "TransitionProbabilities";
89 : }
90 : };
91 :
92 : /**
93 : * @brief The contact patterns within the society are modelled using an UncertainContactMatrix.
94 : */
95 : struct ContactPatterns {
96 : using Type = UncertainContactMatrix<double>;
97 :
98 22 : static Type get_default(AgeGroup size)
99 : {
100 22 : ContactMatrixGroup contact_matrix = ContactMatrixGroup(1, static_cast<Eigen::Index>((size_t)size));
101 66 : contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(static_cast<Eigen::Index>((size_t)size),
102 66 : static_cast<Eigen::Index>((size_t)size), 10.));
103 44 : return Type(contact_matrix);
104 22 : }
105 : static std::string name()
106 : {
107 : return "ContactPatterns";
108 : }
109 : };
110 :
111 : /**
112 : * @brief Probability of getting infected from a contact.
113 : */
114 : struct TransmissionProbabilityOnContact {
115 : using Type = CustomIndexArray<StateAgeFunctionWrapper, AgeGroup>;
116 22 : static Type get_default(AgeGroup size)
117 : {
118 22 : ConstantFunction constfunc(1.0);
119 44 : return Type(size, constfunc);
120 22 : }
121 : static std::string name()
122 : {
123 : return "TransmissionProbabilityOnContact";
124 : }
125 : };
126 :
127 : /**
128 : * @brief The relative InfectedNoSymptoms infectability.
129 : */
130 : struct RelativeTransmissionNoSymptoms {
131 :
132 : using Type = CustomIndexArray<StateAgeFunctionWrapper, AgeGroup>;
133 22 : static Type get_default(AgeGroup size)
134 : {
135 22 : ConstantFunction constfunc(1.0);
136 44 : return Type(size, constfunc);
137 22 : }
138 : static std::string name()
139 : {
140 : return "RelativeTransmissionNoSymptoms";
141 : }
142 : };
143 :
144 : /**
145 : * @brief The risk of infection from symptomatic cases in the SECIR model.
146 : */
147 : struct RiskOfInfectionFromSymptomatic {
148 : using Type = CustomIndexArray<StateAgeFunctionWrapper, AgeGroup>;
149 22 : static Type get_default(AgeGroup size)
150 : {
151 22 : ConstantFunction constfunc(1.0);
152 44 : return Type(size, constfunc);
153 22 : }
154 : static std::string name()
155 : {
156 : return "RiskOfInfectionFromSymptomatic";
157 : }
158 : };
159 :
160 : /**
161 : * @brief Sets the day in a year at which a simulation with an IDE-SECIR model is started.
162 : *
163 : * The value 0.0 corresponds to the 1st of January at 0:00 am, 31.25 is the 1st of February at 6:00 am, and so on.
164 : * The start day defines in which season the simulation is started.
165 : * If the start day is 180 and simulation takes place from t0=0 to
166 : * tmax=100 the days 180 to 280 of the year are simulated.
167 : * The parameter is used in the formula of the seasonality in the model.
168 : */
169 : struct StartDay {
170 : using Type = ScalarType;
171 22 : static Type get_default(AgeGroup)
172 : {
173 22 : return 0.;
174 : }
175 : static std::string name()
176 : {
177 : return "StartDay";
178 : }
179 : };
180 :
181 : /**
182 : * @brief The seasonality parameter k in the IDE-SECIR model.
183 : * The formula for the seasonality used in the model is given as (1+k*sin()) where the sine
184 : * curve is below one in summer and above one in winter.
185 : */
186 : struct Seasonality {
187 : using Type = ScalarType;
188 22 : static Type get_default(AgeGroup)
189 : {
190 22 : return Type(0.);
191 : }
192 : static std::string name()
193 : {
194 : return "Seasonality";
195 : }
196 : };
197 :
198 : // Define Parameterset for IDE-SECIR model.
199 : using ParametersBase =
200 : ParameterSet<TransitionDistributions, TransitionProbabilities, ContactPatterns, TransmissionProbabilityOnContact,
201 : RelativeTransmissionNoSymptoms, RiskOfInfectionFromSymptomatic, StartDay, Seasonality>;
202 :
203 : /**
204 : * @brief Parameters of an age-resolved SECIR/SECIHURD model.
205 : */
206 : class Parameters : public ParametersBase
207 : {
208 : public:
209 22 : Parameters(AgeGroup num_agegroups)
210 22 : : ParametersBase(num_agegroups)
211 22 : , m_num_groups{num_agegroups}
212 : {
213 22 : }
214 :
215 : /**
216 : * @brief Checks whether all Parameters satisfy their corresponding constraints and logs an error.
217 : * @param dt Time step size which is used to get model specific StateAgeFunction%s support.
218 : * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false.
219 : */
220 17 : bool check_constraints() const
221 : {
222 25 : for (AgeGroup group = AgeGroup(0); group < m_num_groups; ++group) {
223 : // For parameters potentially depending on the infectious age, values are checked
224 : // equidistantly on a realistic maximum window.
225 : // Please note that this is an incomplete check on correctness.
226 19 : size_t infectious_window_check = 50; // parameter defining minimal window on x-axis
227 919 : for (size_t i = 0; i < infectious_window_check; i++) {
228 1801 : if (this->get<TransmissionProbabilityOnContact>()[group].eval((ScalarType)i) < 0.0 ||
229 900 : this->get<TransmissionProbabilityOnContact>()[group].eval((ScalarType)i) > 1.0) {
230 1 : log_error("Constraint check: TransmissionProbabilityOnContact smaller {:d} or larger {:d} at some "
231 : "time {:d}",
232 2 : 0, 1, i);
233 1 : return true;
234 : }
235 : }
236 :
237 868 : for (size_t i = 0; i < infectious_window_check; i++) {
238 1701 : if (this->get<RelativeTransmissionNoSymptoms>()[group].eval((ScalarType)i) < 0.0 ||
239 850 : this->get<RelativeTransmissionNoSymptoms>()[group].eval((ScalarType)i) > 1.0) {
240 1 : log_error("Constraint check: RelativeTransmissionNoSymptoms smaller {:d} or larger {:d} at some "
241 : "time {:d}",
242 2 : 0, 1, i);
243 1 : return true;
244 : }
245 : }
246 :
247 817 : for (size_t i = 0; i < infectious_window_check; i++) {
248 1601 : if (this->get<RiskOfInfectionFromSymptomatic>()[group].eval((ScalarType)i) < 0.0 ||
249 800 : this->get<RiskOfInfectionFromSymptomatic>()[group].eval((ScalarType)i) > 1.0) {
250 1 : log_error("Constraint check: RiskOfInfectionFromSymptomatic smaller {:d} or larger {:d} at some "
251 : "time {:d}",
252 2 : 0, 1, i);
253 1 : return true;
254 : }
255 : }
256 :
257 174 : for (size_t i = 0; i < (int)InfectionTransition::Count; i++) {
258 317 : if (this->get<TransitionProbabilities>()[group][i] < 0.0 ||
259 158 : this->get<TransitionProbabilities>()[group][i] > 1.0) {
260 1 : log_error("Constraint check: One parameter in TransitionProbabilities smaller {:d} or larger {:d}",
261 2 : 0, 1);
262 1 : return true;
263 : }
264 : }
265 :
266 30 : if (!floating_point_equal(
267 15 : this->get<TransitionProbabilities>()[group][(int)InfectionTransition::SusceptibleToExposed], 1.0,
268 : 1e-14)) {
269 1 : log_error("Constraint check: Parameter transition probability for SusceptibleToExposed unequal to {:d}",
270 2 : 1);
271 1 : return true;
272 : }
273 :
274 28 : if (!floating_point_equal(
275 14 : this->get<TransitionProbabilities>()[group][(int)InfectionTransition::ExposedToInfectedNoSymptoms],
276 : 1.0, 1e-14)) {
277 1 : log_error("Constraint check: Parameter transition probability for ExposedToInfectedNoSymptoms unequal "
278 : "to {:d}",
279 2 : 1);
280 1 : return true;
281 : }
282 :
283 39 : if (!floating_point_equal(this->get<TransitionProbabilities>()[group][(
284 13 : int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] +
285 13 : this->get<TransitionProbabilities>()[group][(
286 13 : int)InfectionTransition::InfectedNoSymptomsToRecovered],
287 : 1.0, 1e-14)) {
288 1 : log_error(
289 : "Constraint check: Sum of transition probability for InfectedNoSymptomsToInfectedSymptoms and "
290 : "InfectedNoSymptomsToRecovered not equal to {:d}",
291 2 : 1);
292 1 : return true;
293 : }
294 :
295 36 : if (!floating_point_equal(this->get<TransitionProbabilities>()[group][(
296 12 : int)InfectionTransition::InfectedSymptomsToInfectedSevere] +
297 12 : this->get<TransitionProbabilities>()[group][(
298 12 : int)InfectionTransition::InfectedSymptomsToRecovered],
299 : 1.0, 1e-14)) {
300 1 : log_error("Constraint check: Sum of transition probability for InfectedSymptomsToInfectedSevere and "
301 : "InfectedSymptomsToRecovered not equal to {:d}",
302 2 : 1);
303 1 : return true;
304 : }
305 :
306 33 : if (!floating_point_equal(this->get<TransitionProbabilities>()[group][(
307 11 : int)InfectionTransition::InfectedSevereToInfectedCritical] +
308 11 : this->get<TransitionProbabilities>()[group][(
309 11 : int)InfectionTransition::InfectedSevereToRecovered],
310 : 1.0, 1e-14)) {
311 1 : log_error("Constraint check: Sum of transition probability for InfectedSevereToInfectedCritical and "
312 : "InfectedSevereToRecovered not equal to {:d}",
313 2 : 1);
314 1 : return true;
315 : }
316 :
317 30 : if (!floating_point_equal(
318 10 : this->get<TransitionProbabilities>()[group][(int)InfectionTransition::InfectedCriticalToDead] +
319 10 : this->get<TransitionProbabilities>()[group]
320 10 : [(int)InfectionTransition::InfectedCriticalToRecovered],
321 : 1.0, 1e-14)) {
322 1 : log_error("Constraint check: Sum of transition probability for InfectedCriticalToDead and "
323 : "InfectedCriticalToRecovered not equal to {:d}",
324 2 : 1);
325 1 : return true;
326 : }
327 :
328 : /* The first entry of TransitionDistributions is not checked because the distribution S->E is never used
329 : (and it makes no sense to use the distribution). The support does not need to be valid.*/
330 81 : for (size_t i = 1; i < (int)InfectionTransition::Count; i++) {
331 73 : if (floating_point_less(this->get<TransitionDistributions>()[group][i].get_support_max(10), 0.0,
332 : 1e-14)) {
333 1 : log_error("Constraint check: One function in TransitionDistributions has invalid support.");
334 1 : return true;
335 : }
336 : }
337 : }
338 :
339 6 : if (this->get<Seasonality>() < 0.0 || this->get<Seasonality>() > 0.5) {
340 1 : log_warning("Constraint check: Parameter Seasonality should lie between {:0.4f} and {:.4f}", 0.0, 0.5);
341 1 : return true;
342 : }
343 :
344 5 : return false;
345 : }
346 :
347 : /**
348 : * deserialize an object of this class.
349 : * @see mio::deserialize
350 : */
351 : template <class IOContext>
352 : static IOResult<Parameters> deserialize(IOContext& io)
353 : {
354 : BOOST_OUTCOME_TRY(auto&& base, ParametersBase::deserialize(io));
355 : return success(Parameters(std::move(base)));
356 : }
357 :
358 : private:
359 : Parameters(ParametersBase&& base)
360 : : ParametersBase(std::move(base))
361 : , m_num_groups(get<ContactPatterns>().get_cont_freq_mat().get_num_groups())
362 : {
363 : }
364 :
365 : private:
366 : AgeGroup m_num_groups;
367 : };
368 : } // namespace isecir
369 :
370 : } // namespace mio
371 :
372 : #endif // IDE_SECIR_PARAMS_H
|
