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