Skip to content

Commit db770f7

Browse files
reneSchmreneSchm
authored
Remove debug output from SDE seirvv model (#1236)
1 parent 01af90b commit db770f7

File tree

1 file changed

+14
-75
lines changed

1 file changed

+14
-75
lines changed

‎cpp/models/sde_seirvv/model.h‎

Lines changed: 14 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -59,175 +59,114 @@ class Model : public FlowModel<ScalarType, InfectionState, Populations<ScalarTyp
5959
void get_flows(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop, Eigen::Ref<const Eigen::VectorX<ScalarType>> y,
6060
ScalarType t, Eigen::Ref<Eigen::VectorX<ScalarType>> flows) const
6161
{
62-
std::cout << "calc flows" << std::endl;
6362
auto& params = this->parameters;
63+
params.get<ContactPatterns>().get_matrix_at(t)(0, 0);
64+
ScalarType coeffStoIV1 = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
65+
params.get<TransmissionProbabilityOnContactV1>() / populations.get_total();
66+
ScalarType coeffStoIV2 = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
67+
params.get<TransmissionProbabilityOnContactV2>() / populations.get_total();
6468

65-
// Hole den Kontaktmuster-Wert aus der Matrix bei Zeit t
66-
ScalarType cp_val = params.get<ContactPatterns>().get_matrix_at(t)(0, 0);
67-
std::cout << "ContactPatterns value: " << cp_val << std::endl;
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.
6872

69-
// Berechne Koeffizienten für die Transmission
70-
ScalarType coeffStoIV1 = cp_val * params.get<TransmissionProbabilityOnContactV1>() / populations.get_total();
71-
std::cout << "coeffStoIV1: " << coeffStoIV1 << std::endl;
72-
73-
ScalarType coeffStoIV2 = cp_val * params.get<TransmissionProbabilityOnContactV2>() / populations.get_total();
74-
std::cout << "coeffStoIV2: " << coeffStoIV2 << std::endl;
75-
76-
// Normalverteilte Zufallszahlen für den stochastischen Anteil der Flows
7773
ScalarType s_ev1 =
7874
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
79-
std::cout << "s_ev1: " << s_ev1 << std::endl;
80-
8175
ScalarType s_ev2 =
8276
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
83-
std::cout << "s_ev2: " << s_ev2 << std::endl;
84-
8577
ScalarType ev1_iv1 =
8678
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
87-
std::cout << "ev1_iv1: " << ev1_iv1 << std::endl;
88-
8979
ScalarType ev2_iv2 =
9080
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
91-
std::cout << "ev2_iv2: " << ev2_iv2 << std::endl;
92-
9381
ScalarType iv1_rv1 =
9482
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
95-
std::cout << "iv1_rv1: " << iv1_rv1 << std::endl;
96-
9783
ScalarType iv2_rv2 =
9884
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
99-
std::cout << "iv2_rv2: " << iv2_rv2 << std::endl;
100-
10185
ScalarType rv1_ev1v2 =
10286
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
103-
std::cout << "rv1_ev1v2: " << rv1_ev1v2 << std::endl;
104-
10587
ScalarType ev1v2_iv1v2 =
10688
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
107-
std::cout << "ev1v2_iv1v2: " << ev1v2_iv1v2 << std::endl;
108-
10989
ScalarType iv1v2_rv1v2 =
11090
mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
111-
std::cout << "iv1v2_rv1v2: " << iv1v2_rv1v2 << std::endl;
11291

113-
// Berechne Optimierungsgrößen
114-
ScalarType inv_step_size = 1.0 / step_size;
115-
std::cout << "inv_step_size: " << inv_step_size << std::endl;
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.
11695

96+
// Calculate inv_step_size and inv_sqrt_step_size for optimization.
97+
ScalarType inv_step_size = 1.0 / step_size;
11798
ScalarType inv_sqrt_step_size = 1.0 / sqrt(step_size);
118-
std::cout << "inv_sqrt_step_size: " << inv_sqrt_step_size << std::endl;
11999

120-
// Berechne den ersten Outflow aus S
100+
// Two outgoing flows from S so will clamp their sum to S * inv_step_size to ensure non-negative S.
121101
const ScalarType outflow1 = std::clamp(
122102
coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1] +
123103
sqrt(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::InfectedV1]) *
124104
inv_sqrt_step_size * s_ev1,
125105
0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size);
126-
std::cout << "outflow1: " << outflow1 << std::endl;
127106

128-
// Berechne den zweiten Outflow aus S
129107
const ScalarType outflow2 =
130108
std::clamp(coeffStoIV1 * y[(size_t)InfectionState::Susceptible] *
131109
(pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) +
132110
sqrt(coeffStoIV2 * y[(size_t)InfectionState::Susceptible] *
133111
(pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2])) *
134112
inv_sqrt_step_size * s_ev2,
135113
0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size);
136-
std::cout << "outflow2: " << outflow2 << std::endl;
137114

138-
// Summe der Outflows
139115
const ScalarType outflow_sum = outflow1 + outflow2;
140-
std::cout << "outflow_sum: " << outflow_sum << std::endl;
141-
142116
if (outflow_sum > 0) {
143117
const ScalarType scale =
144118
std::clamp(outflow_sum, 0.0, y[(size_t)InfectionState::Susceptible] * inv_step_size) / outflow_sum;
145-
std::cout << "scale: " << scale << std::endl;
146119
flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()] = outflow1 * scale;
147-
std::cout << "flow Sus -> ExpV1: "
148-
<< flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()]
149-
<< std::endl;
150120
flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()] = outflow2 * scale;
151-
std::cout << "flow Sus -> ExpV2: "
152-
<< flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()]
153-
<< std::endl;
154121
}
155122
else {
156123
flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV1>()] = 0;
157-
std::cout << "flow Sus -> ExpV1 auf 0 gesetzt" << std::endl;
158124
flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::ExposedV2>()] = 0;
159-
std::cout << "flow Sus -> ExpV2 auf 0 gesetzt" << std::endl;
160125
}
161126

162-
// Fluss von ExposedV1 zu InfectedV1
163127
flows[get_flat_flow_index<InfectionState::ExposedV1, InfectionState::InfectedV1>()] =
164128
std::clamp((1.0 / params.get<TimeExposedV1>()) * y[(size_t)InfectionState::ExposedV1] +
165129
sqrt((1.0 / params.get<TimeExposedV1>()) * y[(size_t)InfectionState::ExposedV1]) *
166130
inv_sqrt_step_size * ev1_iv1,
167131
0.0, y[(size_t)InfectionState::ExposedV1] * inv_step_size);
168-
std::cout << "flow ExpV1 -> InfV1: "
169-
<< flows[get_flat_flow_index<InfectionState::ExposedV1, InfectionState::InfectedV1>()] << std::endl;
170132

171-
// Fluss von ExposedV2 zu InfectedV2
172133
flows[get_flat_flow_index<InfectionState::ExposedV2, InfectionState::InfectedV2>()] =
173134
std::clamp((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV2] +
174135
sqrt((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV2]) *
175136
inv_sqrt_step_size * ev2_iv2,
176137
0.0, y[(size_t)InfectionState::ExposedV2] * inv_step_size);
177-
std::cout << "flow ExpV2 -> InfV2: "
178-
<< flows[get_flat_flow_index<InfectionState::ExposedV2, InfectionState::InfectedV2>()] << std::endl;
179138

180-
// Fluss von InfectedV1 zu RecoveredV1
181139
flows[get_flat_flow_index<InfectionState::InfectedV1, InfectionState::RecoveredV1>()] =
182140
std::clamp((1.0 / params.get<TimeInfectedV1>()) * y[(size_t)InfectionState::InfectedV1] +
183141
sqrt((1.0 / params.get<TimeInfectedV1>()) * y[(size_t)InfectionState::InfectedV1]) *
184142
inv_sqrt_step_size * iv1_rv1,
185143
0.0, y[(size_t)InfectionState::InfectedV1] * inv_step_size);
186-
std::cout << "flow InfV1 -> RecV1: "
187-
<< flows[get_flat_flow_index<InfectionState::InfectedV1, InfectionState::RecoveredV1>()] << std::endl;
188144

189-
// Fluss von InfectedV2 zu RecoveredV2
190145
flows[get_flat_flow_index<InfectionState::InfectedV2, InfectionState::RecoveredV2>()] =
191146
std::clamp((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV2] +
192147
sqrt((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV2]) *
193148
inv_sqrt_step_size * iv2_rv2,
194149
0.0, y[(size_t)InfectionState::InfectedV2] * inv_step_size);
195-
std::cout << "flow InfV2 -> RecV2: "
196-
<< flows[get_flat_flow_index<InfectionState::InfectedV2, InfectionState::RecoveredV2>()] << std::endl;
197150

198-
// Fluss von RecoveredV1 zu ExposedV1V2
199151
flows[get_flat_flow_index<InfectionState::RecoveredV1, InfectionState::ExposedV1V2>()] =
200152
std::clamp(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] *
201153
(pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2]) +
202154
sqrt(coeffStoIV2 * y[(size_t)InfectionState::RecoveredV1] *
203155
(pop[(size_t)InfectionState::InfectedV1V2] + pop[(size_t)InfectionState::InfectedV2])) *
204156
inv_sqrt_step_size * rv1_ev1v2,
205157
0.0, y[(size_t)InfectionState::RecoveredV1] * inv_step_size);
206-
std::cout << "flow RecV1 -> ExpV1V2: "
207-
<< flows[get_flat_flow_index<InfectionState::RecoveredV1, InfectionState::ExposedV1V2>()]
208-
<< std::endl;
209158

210-
// Fluss von ExposedV1V2 zu InfectedV1V2
211159
flows[get_flat_flow_index<InfectionState::ExposedV1V2, InfectionState::InfectedV1V2>()] =
212160
std::clamp((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV1V2] +
213161
sqrt((1.0 / params.get<TimeExposedV2>()) * y[(size_t)InfectionState::ExposedV1V2]) /
214162
sqrt(step_size) * ev1v2_iv1v2,
215163
0.0, y[(size_t)InfectionState::ExposedV1V2] * inv_step_size);
216-
std::cout << "flow ExpV1V2 -> InfV1V2: "
217-
<< flows[get_flat_flow_index<InfectionState::ExposedV1V2, InfectionState::InfectedV1V2>()]
218-
<< std::endl;
219164

220-
// Fluss von InfectedV1V2 zu RecoveredV1V2
221165
flows[get_flat_flow_index<InfectionState::InfectedV1V2, InfectionState::RecoveredV1V2>()] =
222166
std::clamp((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV1V2] +
223167
sqrt((1.0 / params.get<TimeInfectedV2>()) * y[(size_t)InfectionState::InfectedV1V2]) /
224168
sqrt(step_size) * iv1v2_rv1v2,
225169
0.0, y[(size_t)InfectionState::InfectedV1V2] * inv_step_size);
226-
std::cout << "flow InfV1V2 -> RecV1V2: "
227-
<< flows[get_flat_flow_index<InfectionState::InfectedV1V2, InfectionState::RecoveredV1V2>()]
228-
<< std::endl;
229-
230-
std::cout << "calc flows done" << std::endl;
231170
}
232171

233172
ScalarType step_size; ///< A step size of the model with which the stochastic process is realized.

0 commit comments

Comments
 (0)