Line data Source code
1 : /*
2 : * Copyright (C) 2020-2025 MEmilio
3 : *
4 : * Authors: Nils Wassmuth, Rene Schmieding, Martin J. Kuehn
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 MIO_SDE_SEIRVV_MODEL_H
22 : #define MIO_SDE_SEIRVV_MODEL_H
23 :
24 : #include "memilio/compartments/flow_model.h"
25 : #include "memilio/epidemiology/populations.h"
26 : #include "memilio/utils/random_number_generator.h"
27 : #include "sde_seirvv/infection_state.h"
28 : #include "sde_seirvv/parameters.h"
29 :
30 : namespace mio
31 : {
32 : namespace sseirvv
33 : {
34 :
35 : /********************
36 : * define the model *
37 : ********************/
38 :
39 : using Flows = TypeList<Flow<InfectionState::Susceptible, InfectionState::ExposedV1>,
40 : Flow<InfectionState::Susceptible, InfectionState::ExposedV2>,
41 : Flow<InfectionState::ExposedV1, InfectionState::InfectedV1>,
42 : Flow<InfectionState::ExposedV2, InfectionState::InfectedV2>,
43 : Flow<InfectionState::InfectedV1, InfectionState::RecoveredV1>,
44 : Flow<InfectionState::InfectedV2, InfectionState::RecoveredV2>,
45 : Flow<InfectionState::RecoveredV1, InfectionState::ExposedV1V2>,
46 : Flow<InfectionState::ExposedV1V2, InfectionState::InfectedV1V2>,
47 : Flow<InfectionState::InfectedV1V2, InfectionState::RecoveredV1V2>>;
48 :
49 : class Model : public FlowModel<ScalarType, InfectionState, Populations<ScalarType, InfectionState>, Parameters, Flows>
50 : {
51 : using Base = FlowModel<ScalarType, InfectionState, mio::Populations<ScalarType, InfectionState>, Parameters, Flows>;
52 :
53 : public:
54 9 : Model()
55 9 : : Base(Populations({InfectionState::Count}, 0.), ParameterSet())
56 : {
57 9 : }
58 :
59 63 : void get_flows(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop, Eigen::Ref<const Eigen::VectorX<ScalarType>> y,
60 : ScalarType t, Eigen::Ref<Eigen::VectorX<ScalarType>> flows) const
61 : {
62 63 : auto& params = this->parameters;
63 63 : params.get<ContactPatterns>().get_matrix_at(t)(0, 0);
64 126 : ScalarType coeffStoIV1 = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
65 189 : params.get<TransmissionProbabilityOnContactV1>() / populations.get_total();
66 126 : ScalarType coeffStoIV2 = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
67 189 : params.get<TransmissionProbabilityOnContactV2>() / populations.get_total();
68 :
69 : // Normal distributed values for the stochastic part of the flows, variables are encoded
70 : // in the following way: x_y is the stochastic part for the flow from x to y. Variant
71 : // specific compartments also get an addendum v1 or v2 denoting the relevant variant.
72 :
73 : ScalarType s_ev1 =
74 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
75 : ScalarType s_ev2 =
76 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
77 : ScalarType ev1_iv1 =
78 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
79 : ScalarType ev2_iv2 =
80 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
81 : ScalarType iv1_rv1 =
82 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
83 : ScalarType iv2_rv2 =
84 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
85 : ScalarType rv1_ev1v2 =
86 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
87 : ScalarType ev1v2_iv1v2 =
88 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
89 : ScalarType iv1v2_rv1v2 =
90 63 : mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
91 :
92 : // Assuming that no person can change its InfectionState twice in a single time step,
93 : // take the minimum of the calculated flow and the source compartment, to ensure that
94 : // no compartment attains negative values.
95 :
96 : // Calculate inv_step_size and inv_sqrt_step_size for optimization.
97 63 : ScalarType inv_step_size = 1.0 / step_size;
98 63 : ScalarType inv_sqrt_step_size = 1.0 / sqrt(step_size);
99 :
100 : // Two outgoing flows from S so will clamp their sum to S * inv_step_size to ensure non-negative S.
101 63 : const ScalarType outflow1 = std::clamp(
102 63 : coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1] +
103 63 : sqrt(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1]) *
104 63 : inv_sqrt_step_size * s_ev1,
105 126 : 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size);
106 :
107 : const ScalarType outflow2 =
108 126 : std::clamp(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] *
109 126 : (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) +
110 63 : sqrt(coeffStoIV2 * y[(size_t)InfectionState::Susceptible] *
111 63 : (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2])) *
112 63 : inv_sqrt_step_size * s_ev2,
113 126 : 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size);
114 :
115 63 : const ScalarType outflow_sum = outflow1 + outflow2;
116 63 : if (outflow_sum > 0) {
117 : const ScalarType scale =
118 45 : std::clamp(outflow_sum, 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size) / outflow_sum;
119 45 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()] = outflow1 * scale;
120 45 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()] = outflow2 * scale;
121 : }
122 : else {
123 18 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()] = 0;
124 18 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()] = 0;
125 : }
126 :
127 63 : flows[get_flat_flow_index<InfectionState::ExposedV1, InfectionState::InfectedV1>()] =
128 189 : std::clamp((1.0 / params.get<TimeExposedV1>()) * y[(size_t)InfectionState::ExposedV1] +
129 63 : sqrt((1.0 / params.get<TimeExposedV1>()) * y[(size_t)InfectionState::ExposedV1]) *
130 63 : inv_sqrt_step_size * ev1_iv1,
131 126 : 0.0, y[(size_t)InfectionState::ExposedV1] * inv_step_size);
132 :
133 63 : flows[get_flat_flow_index<InfectionState::ExposedV2, InfectionState::InfectedV2>()] =
134 189 : std::clamp((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV2] +
135 63 : sqrt((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV2]) *
136 63 : inv_sqrt_step_size * ev2_iv2,
137 126 : 0.0, y[(size_t)InfectionState::ExposedV2] * inv_step_size);
138 :
139 63 : flows[get_flat_flow_index<InfectionState::InfectedV1, InfectionState::RecoveredV1>()] =
140 189 : std::clamp((1.0 / params.get<TimeInfectedV1>()) * y[(size_t)InfectionState::InfectedV1] +
141 63 : sqrt((1.0 / params.get<TimeInfectedV1>()) * y[(size_t)InfectionState::InfectedV1]) *
142 63 : inv_sqrt_step_size * iv1_rv1,
143 126 : 0.0, y[(size_t)InfectionState::InfectedV1] * inv_step_size);
144 :
145 63 : flows[get_flat_flow_index<InfectionState::InfectedV2, InfectionState::RecoveredV2>()] =
146 189 : std::clamp((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV2] +
147 63 : sqrt((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV2]) *
148 63 : inv_sqrt_step_size * iv2_rv2,
149 126 : 0.0, y[(size_t)InfectionState::InfectedV2] * inv_step_size);
150 :
151 63 : flows[get_flat_flow_index<InfectionState::RecoveredV1, InfectionState::ExposedV1V2>()] =
152 126 : std::clamp(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] *
153 189 : (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) +
154 63 : sqrt(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] *
155 63 : (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2])) *
156 63 : inv_sqrt_step_size * rv1_ev1v2,
157 126 : 0.0, y[(size_t)InfectionState::RecoveredV1] * inv_step_size);
158 :
159 63 : flows[get_flat_flow_index<InfectionState::ExposedV1V2, InfectionState::InfectedV1V2>()] =
160 189 : std::clamp((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV1V2] +
161 63 : sqrt((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV1V2]) /
162 63 : sqrt(step_size) * ev1v2_iv1v2,
163 126 : 0.0, y[(size_t)InfectionState::ExposedV1V2] * inv_step_size);
164 :
165 63 : flows[get_flat_flow_index<InfectionState::InfectedV1V2, InfectionState::RecoveredV1V2>()] =
166 189 : std::clamp((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV1V2] +
167 63 : sqrt((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV1V2]) /
168 63 : sqrt(step_size) * iv1v2_rv1v2,
169 126 : 0.0, y[(size_t)InfectionState::InfectedV1V2] * inv_step_size);
170 63 : }
171 :
172 : ScalarType step_size; ///< A step size of the model with which the stochastic process is realized.
173 : mutable RandomNumberGenerator rng;
174 :
175 : private:
176 : };
177 :
178 : } // namespace sseirvv
179 : } // namespace mio
180 :
181 : #endif // MIO_SDE_SEIRVV_MODEL_H
|
