Line data Source code
1 : /*
2 : * Copyright (C) 2020-2024 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 Vector<ScalarType>> pop, Eigen::Ref<const Vector<ScalarType>> y, ScalarType t,
60 : Eigen::Ref<Vector<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 63 : ScalarType s_ev1 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
74 63 : ScalarType s_ev2 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
75 63 : ScalarType ev1_iv1 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
76 63 : ScalarType ev2_iv2 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
77 63 : ScalarType iv1_rv1 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
78 63 : ScalarType iv2_rv2 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
79 63 : ScalarType rv1_ev1v2 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
80 63 : ScalarType ev1v2_iv1v2 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
81 63 : ScalarType iv1v2_rv1v2 = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
82 :
83 : // Assuming that no person can change its InfectionState twice in a single time step,
84 : // take the minimum of the calculated flow and the source compartment, to ensure that
85 : // no compartment attains negative values.
86 :
87 : // Calculate inv_step_size and inv_sqrt_step_size for optimization.
88 63 : ScalarType inv_step_size = 1.0 / step_size;
89 63 : ScalarType inv_sqrt_step_size = 1.0 / sqrt(step_size);
90 :
91 : // Two outgoing flows from S so will clamp their sum to S * inv_step_size to ensure non-negative S.
92 63 : const ScalarType outflow1 = std::clamp(
93 63 : coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1] +
94 63 : sqrt(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1])
95 63 : * inv_sqrt_step_size * s_ev1,
96 126 : 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size);
97 :
98 63 : const ScalarType outflow2 = std::clamp(
99 63 : coeffStoIV1 * y[(size_t)InfectionState::Susceptible] *
100 126 : (pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) +
101 63 : sqrt(coeffStoIV2 * y[(size_t)InfectionState::Susceptible] * (pop[(size_t)InfectionState::InfectedV1V2]
102 63 : + pop[(size_t)InfectionState::InfectedV2])) * inv_sqrt_step_size * s_ev2,
103 126 : 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size);
104 :
105 63 : const ScalarType outflow_sum = outflow1 + outflow2;
106 63 : if (outflow_sum > 0)
107 : {
108 45 : const ScalarType scale = std::clamp(outflow_sum, 0.0 , y[(size_t)InfectionState::Susceptible] * inv_step_size) / outflow_sum;
109 45 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()] = outflow1 * scale;
110 45 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()] = outflow2 * scale;
111 : } else
112 : {
113 18 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()] = 0;
114 18 : flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()] = 0;
115 : }
116 :
117 63 : flows[get_flat_flow_index<InfectionState::ExposedV1, InfectionState::InfectedV1>()] = std::clamp(
118 126 : (1.0 / params.get<TimeExposedV1>()) * y[(size_t)InfectionState::ExposedV1] +
119 63 : sqrt((1.0 / params.get<TimeExposedV1>()) * y[(size_t)InfectionState::ExposedV1]) * inv_sqrt_step_size * ev1_iv1,
120 126 : 0.0, y[(size_t)InfectionState::ExposedV1] * inv_step_size);
121 :
122 63 : flows[get_flat_flow_index<InfectionState::ExposedV2, InfectionState::InfectedV2>()] = std::clamp(
123 126 : (1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV2] +
124 63 : sqrt((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV2]) * inv_sqrt_step_size * ev2_iv2,
125 126 : 0.0, y[(size_t)InfectionState::ExposedV2] * inv_step_size);
126 :
127 63 : flows[get_flat_flow_index<InfectionState::InfectedV1, InfectionState::RecoveredV1>()] = std::clamp(
128 126 : (1.0 / params.get<TimeInfectedV1>()) * y[(size_t)InfectionState::InfectedV1] +
129 63 : sqrt((1.0 / params.get<TimeInfectedV1>()) * y[(size_t)InfectionState::InfectedV1]) * inv_sqrt_step_size * iv1_rv1,
130 126 : 0.0, y[(size_t)InfectionState::InfectedV1] * inv_step_size);
131 :
132 63 : flows[get_flat_flow_index<InfectionState::InfectedV2, InfectionState::RecoveredV2>()] = std::clamp(
133 126 : (1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV2] +
134 63 : sqrt((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV2]) * inv_sqrt_step_size * iv2_rv2,
135 126 : 0.0, y[(size_t)InfectionState::InfectedV2] * inv_step_size);
136 :
137 63 : flows[get_flat_flow_index<InfectionState::RecoveredV1, InfectionState::ExposedV1V2>()] = std::clamp(
138 63 : coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1]
139 189 : * (pop[(size_t)InfectionState::InfectedV1V2]+ pop[(size_t)InfectionState::InfectedV2]) +
140 63 : sqrt(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] * (pop[(size_t)InfectionState::InfectedV1V2]
141 63 : + pop[(size_t)InfectionState::InfectedV2])) * inv_sqrt_step_size * rv1_ev1v2,
142 126 : 0.0, y[(size_t)InfectionState::RecoveredV1] * inv_step_size);
143 :
144 63 : flows[get_flat_flow_index<InfectionState::ExposedV1V2, InfectionState::InfectedV1V2>()] = std::clamp(
145 126 : (1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV1V2] +
146 63 : sqrt((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV1V2]) /
147 63 : sqrt(step_size) * ev1v2_iv1v2,
148 126 : 0.0, y[(size_t)InfectionState::ExposedV1V2] * inv_step_size);
149 :
150 63 : flows[get_flat_flow_index<InfectionState::InfectedV1V2, InfectionState::RecoveredV1V2>()] = std::clamp(
151 126 : (1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV1V2] +
152 63 : sqrt((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV1V2]) /
153 63 : sqrt(step_size) * iv1v2_rv1v2,
154 126 : 0.0, y[(size_t)InfectionState::InfectedV1V2] * inv_step_size);
155 63 : }
156 :
157 : ScalarType step_size; ///< A step size of the model with which the stochastic process is realized.
158 : mutable RandomNumberGenerator rng;
159 :
160 : private:
161 : };
162 :
163 : } // namespace sseirvv
164 : } // namespace mio
165 :
166 : #endif // MIO_SDE_SEIRVV_MODEL_H
|
