Skip to content

Commit 35e1c2a

Browse files
HenrZureneSchm
HenrZu
and
reneSchm
authored
1102 divNj creates NaN values if subpopulation is 0 (#1104)
-floating type specific zero tolerance check to prevent NaN values if a subpopulation is zero Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>
1 parent 32cd836 commit 35e1c2a

File tree

10 files changed

+168
-44
lines changed

10 files changed

+168
-44
lines changed

‎cpp/benchmarks/flow_simulation_ode_secirvvs.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -203,7 +203,7 @@ class FlowlessModel : public CompartmentalModel<ScalarType, osecirvvs::Infection
203203
pop[ICrPIj] + pop[INSPICj] + pop[ISyPICj] + pop[SIIj] + pop[EIIj] + pop[INSIIj] +
204204
pop[ISyIIj] + pop[ISevIIj] + pop[ICrIIj] + pop[INSIICj] + pop[ISyIICj];
205205

206-
double divNj = 1.0 / Nj; // precompute 1.0/Nj
206+
const double divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
207207

208208
double ext_inf_force_dummy =
209209
cont_freq_eff * divNj *

‎cpp/memilio/config.h‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,48 @@
2626
#define MIO_CONFIG_H
2727

2828
#include "memilio/config_internal.h"
29+
#include <type_traits>
2930

3031
using ScalarType = double;
3132

33+
namespace mio
34+
{
35+
/**
36+
* @brief Type specific limits for floating-points.
37+
*
38+
* Specializations are provided for `float` and `double` types. In order to use
39+
* other floating point types, a new specialization for that type must be added.
40+
*
41+
* Member functions of the unspecialized class are disabled. Using them results
42+
* in compile time errors.
43+
* @{
44+
*/
45+
/**
46+
* @tparam FP A floating point type.
47+
*/
48+
template <typename FP>
49+
struct Limits {
50+
static constexpr FP zero_tolerance() = delete;
51+
};
52+
53+
template <>
54+
struct Limits<float> {
55+
/// @brief Returns the limit under which a float may be rounded down to zero.
56+
static constexpr float zero_tolerance()
57+
{
58+
return 1e-6f;
59+
}
60+
};
61+
62+
template <>
63+
struct Limits<double> {
64+
/// @brief Returns the limit under which a double may be rounded down to zero.
65+
static constexpr double zero_tolerance()
66+
{
67+
return 1e-12;
68+
}
69+
};
70+
/** @} */
71+
72+
} // namespace mio
3273
#endif

‎cpp/models/ode_secir/model.h‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -122,31 +122,31 @@ class Model : public FlowModel<FP, InfectionState, Populations<FP, AgeGroup, Inf
122122
params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[j]);
123123

124124
// effective contact rate by contact rate between groups i and j and damping j
125-
double season_val =
125+
ScalarType season_val =
126126
(1 + params.template get<Seasonality<FP>>() *
127127
sin(3.141592653589793 * ((params.template get<StartDay>() + t) / 182.5 + 0.5)));
128-
double cont_freq_eff =
128+
ScalarType cont_freq_eff =
129129
season_val * contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
130130
static_cast<Eigen::Index>((size_t)j));
131-
double Nj =
131+
ScalarType Nj =
132132
pop[Sj] + pop[Ej] + pop[INSj] + pop[ISyj] + pop[ISevj] + pop[ICrj] + pop[Rj]; // without died people
133-
double divNj = 1.0 / Nj; // precompute 1.0/Nj
134-
double dummy_S = y[Si] * cont_freq_eff * divNj *
135-
params.template get<TransmissionProbabilityOnContact<FP>>()[i] *
136-
(params.template get<RelativeTransmissionNoSymptoms<FP>>()[j] * pop[INSj] +
137-
riskFromInfectedSymptomatic * pop[ISyj]);
133+
const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
134+
ScalarType dummy_S = y[Si] * cont_freq_eff * divNj *
135+
params.template get<TransmissionProbabilityOnContact<FP>>()[i] *
136+
(params.template get<RelativeTransmissionNoSymptoms<FP>>()[j] * pop[INSj] +
137+
riskFromInfectedSymptomatic * pop[ISyj]);
138138

139139
// Susceptible -> Exposed
140140
flows[this->template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Exposed>({i})] +=
141141
dummy_S;
142142
}
143143

144144
// ICU capacity shortage is close
145-
double criticalPerSevereAdjusted = smoother_cosine(
145+
ScalarType criticalPerSevereAdjusted = smoother_cosine(
146146
icu_occupancy, 0.90 * params.template get<ICUCapacity<FP>>(), params.template get<ICUCapacity<FP>>(),
147147
params.template get<CriticalPerSevere<FP>>()[i], 0);
148148

149-
double deathsPerSevereAdjusted =
149+
ScalarType deathsPerSevereAdjusted =
150150
params.template get<CriticalPerSevere<FP>>()[i] - criticalPerSevereAdjusted;
151151

152152
// Exposed -> InfectedNoSymptoms

‎cpp/models/ode_secirvvs/model.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -242,8 +242,7 @@ class Model
242242
pop[ISyNCj] + pop[SPIj] + pop[EPIj] + pop[INSPIj] + pop[ISyPIj] + pop[ISevPIj] + pop[ICrPIj] +
243243
pop[INSPICj] + pop[ISyPICj] + pop[SIIj] + pop[EIIj] + pop[INSIIj] + pop[ISyIIj] + pop[ISevIIj] +
244244
pop[ICrIIj] + pop[INSIICj] + pop[ISyIICj];
245-
246-
FP divNj = 1.0 / Nj; // precompute 1.0/Nj
245+
const FP divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
247246

248247
FP ext_inf_force_dummy = cont_freq_eff * divNj *
249248
params.template get<TransmissionProbabilityOnContact<FP>>()[(AgeGroup)i] *

‎cpp/models/ode_seir/model.h‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -85,10 +85,11 @@ class Model
8585
const size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected});
8686
const size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
8787

88-
const double Nj_inv = 1.0 / (pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj]);
89-
const double coeffStoE =
88+
const ScalarType Nj = pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj];
89+
const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
90+
const ScalarType coeffStoE =
9091
params.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_matrix_at(t)(i.get(), j.get()) *
91-
params.template get<TransmissionProbabilityOnContact<FP>>()[i] * Nj_inv;
92+
params.template get<TransmissionProbabilityOnContact<FP>>()[i] * divNj;
9293

9394
flows[Base::template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Exposed>(i)] +=
9495
coeffStoE * y[Si] * pop[Ij];
@@ -127,16 +128,16 @@ class Model
127128
size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible});
128129
for (auto j = AgeGroup(0); j < AgeGroup(num_groups); j++) {
129130

130-
double Nj = this->populations.get_group_total(j);
131-
double divNj = 1.0 / Nj;
131+
const ScalarType Nj = this->populations.get_group_total(j);
132+
const ScalarType divNj = (Nj < 1e-12) ? 0.0 : 1.0 / Nj;
132133

133-
double coeffStoE = contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) *
134-
params.template get<TransmissionProbabilityOnContact<ScalarType>>()[i] * divNj;
134+
ScalarType coeffStoE = contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) *
135+
params.template get<TransmissionProbabilityOnContact<ScalarType>>()[i] * divNj;
135136
F((size_t)i, (size_t)j + num_groups) = coeffStoE * y.get_value(t_idx)[Si];
136137
}
137138

138-
double T_Ei = params.template get<mio::oseir::TimeExposed<ScalarType>>()[i];
139-
double T_Ii = params.template get<mio::oseir::TimeInfected<ScalarType>>()[i];
139+
ScalarType T_Ei = params.template get<mio::oseir::TimeExposed<ScalarType>>()[i];
140+
ScalarType T_Ii = params.template get<mio::oseir::TimeInfected<ScalarType>>()[i];
140141
V((size_t)i, (size_t)i) = 1.0 / T_Ei;
141142
V((size_t)i + num_groups, (size_t)i) = -1.0 / T_Ei;
142143
V((size_t)i + num_groups, (size_t)i + num_groups) = 1.0 / T_Ii;

‎cpp/models/ode_sir/model.h‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -77,11 +77,12 @@ class Model
7777
size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected});
7878
size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
7979

80-
double Nj = pop[Sj] + pop[Ij] + pop[Rj];
80+
const ScalarType Nj = pop[Sj] + pop[Ij] + pop[Rj];
81+
const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
8182

82-
double coeffStoI = contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
83-
static_cast<Eigen::Index>((size_t)j)) *
84-
params.template get<TransmissionProbabilityOnContact<FP>>()[i] / Nj;
83+
ScalarType coeffStoI = contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
84+
static_cast<Eigen::Index>((size_t)j)) *
85+
params.template get<TransmissionProbabilityOnContact<FP>>()[i] * divNj;
8586

8687
dydt[Si] += -coeffStoI * y[Si] * pop[Ij];
8788
dydt[Ii] += coeffStoI * y[Si] * pop[Ij];

‎cpp/tests/test_odesecir.cpp‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -573,6 +573,25 @@ TEST(TestOdeSecir, testSettersAndGetters)
573573
EXPECT_EQ(vec[21], model.parameters.get<mio::osecir::Seasonality<double>>());
574574
}
575575

576+
// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
577+
TEST(TestOdeSecir, population_zero_no_nan)
578+
{
579+
// initialize simple model with total population 0
580+
mio::osecir::Model<double> model(1);
581+
model.populations.set_total(0.0);
582+
583+
// call the get_flows function
584+
auto dydt_default = Eigen::VectorXd(15);
585+
dydt_default.setZero();
586+
auto y0 = model.get_initial_values();
587+
model.get_flows(y0, y0, 0, dydt_default);
588+
589+
// check that there are now NaN values in dydt_default
590+
for (int i = 0; i < dydt_default.size(); i++) {
591+
EXPECT_FALSE(std::isnan(dydt_default[i]));
592+
}
593+
}
594+
576595
TEST(TestOdeSecir, testDamping)
577596
{
578597
// Test functionality of dampings

‎cpp/tests/test_odesecirvvs.cpp‎

Lines changed: 42 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1091,18 +1091,20 @@ TEST(TestOdeSECIRVVS, test_commuters)
10911091
model.parameters.get_commuter_nondetection() = non_detection_factor;
10921092
auto sim = mio::osecirvvs::Simulation<>(model);
10931093
auto before_testing = sim.get_result().get_last_value().eval();
1094-
auto mobile_population = (sim.get_result().get_last_value() * mobility_factor).eval();
1095-
auto mobile_population_tested = mobile_population.eval();
1094+
auto mobile_population = (sim.get_result().get_last_value() * mobility_factor).eval();
1095+
auto mobile_population_tested = mobile_population.eval();
10961096

10971097
mio::osecirvvs::test_commuters<double>(sim, mobile_population_tested, 0.0);
10981098

10991099
ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)],
1100-
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] * non_detection_factor,
1100+
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] *
1101+
non_detection_factor,
11011102
1e-5);
11021103
ASSERT_NEAR(
11031104
sim.get_result().get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaiveConfirmed)],
11041105
before_testing[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaiveConfirmed)] +
1105-
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] * (1 - non_detection_factor),
1106+
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] *
1107+
(1 - non_detection_factor),
11061108
1e-5);
11071109
ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity)],
11081110
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity)] *
@@ -1115,10 +1117,11 @@ TEST(TestOdeSECIRVVS, test_commuters)
11151117
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity)] *
11161118
(1 - non_detection_factor),
11171119
1e-5);
1118-
ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)],
1119-
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)] *
1120-
non_detection_factor,
1121-
1e-5);
1120+
ASSERT_NEAR(
1121+
mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)],
1122+
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)] *
1123+
non_detection_factor,
1124+
1e-5);
11221125
ASSERT_NEAR(
11231126
sim.get_result()
11241127
.get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunityConfirmed)],
@@ -1128,29 +1131,32 @@ TEST(TestOdeSECIRVVS, test_commuters)
11281131
1e-5);
11291132

11301133
ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)],
1131-
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)] * non_detection_factor,
1134+
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)] *
1135+
non_detection_factor,
11321136
1e-5);
11331137
ASSERT_NEAR(sim.get_result()
11341138
.get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaiveConfirmed)],
11351139
before_testing[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaiveConfirmed)] +
11361140
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)] *
11371141
(1 - non_detection_factor),
11381142
1e-5);
1139-
ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)],
1140-
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)] *
1141-
non_detection_factor,
1142-
1e-5);
1143+
ASSERT_NEAR(
1144+
mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)],
1145+
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)] *
1146+
non_detection_factor,
1147+
1e-5);
11431148
ASSERT_NEAR(
11441149
sim.get_result()
11451150
.get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed)],
11461151
before_testing[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed)] +
11471152
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)] *
11481153
(1 - non_detection_factor),
11491154
1e-5);
1150-
ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)],
1151-
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)] *
1152-
non_detection_factor,
1153-
1e-5);
1155+
ASSERT_NEAR(
1156+
mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)],
1157+
mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)] *
1158+
non_detection_factor,
1159+
1e-5);
11541160
ASSERT_NEAR(
11551161
sim.get_result().get_last_value()[Eigen::Index(
11561162
mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed)],
@@ -1160,6 +1166,25 @@ TEST(TestOdeSECIRVVS, test_commuters)
11601166
1e-5);
11611167
}
11621168

1169+
// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
1170+
TEST(TestOdeSECIRVVS, population_zero_no_nan)
1171+
{
1172+
// initialize simple model with total population 0
1173+
mio::osecirvvs::Model<double> model(1);
1174+
model.populations.set_total(0.0);
1175+
1176+
// call the get_flows function
1177+
auto dydt_default = Eigen::VectorXd(45);
1178+
dydt_default.setZero();
1179+
auto y0 = model.get_initial_values();
1180+
model.get_flows(y0, y0, 0, dydt_default);
1181+
1182+
// check that there are now NaN values in dydt_default
1183+
for (int i = 0; i < dydt_default.size(); i++) {
1184+
EXPECT_FALSE(std::isnan(dydt_default[i]));
1185+
}
1186+
}
1187+
11631188
TEST(TestOdeSECIRVVS, check_constraints_parameters)
11641189
{
11651190
auto model = mio::osecirvvs::Model<double>(1);

‎cpp/tests/test_odeseir.cpp‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -311,6 +311,25 @@ TEST(TestOdeSeir, get_reproduction_number)
311311
EXPECT_NEAR(model.get_reproduction_number(0.9, result).value(), 1.858670429549998504, 1e-12);
312312
}
313313

314+
// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
315+
TEST(TestSeir, population_zero_no_nan)
316+
{
317+
// initialize simple model with total population 0
318+
mio::oseir::Model<double> model(1);
319+
model.populations.set_total(0.0);
320+
321+
// call the get_flows function
322+
auto dydt_default = Eigen::VectorXd(3);
323+
dydt_default.setZero();
324+
auto y0 = model.get_initial_values();
325+
model.get_flows(y0, y0, 0, dydt_default);
326+
327+
// check that there are now NaN values in dydt_default
328+
for (int i = 0; i < dydt_default.size(); i++) {
329+
EXPECT_FALSE(std::isnan(dydt_default[i]));
330+
}
331+
}
332+
314333
TEST(TestSeir, get_flows)
315334
{
316335
mio::oseir::Model<double> model(1);

‎cpp/tests/test_odesir.cpp‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -211,6 +211,25 @@ TEST(Testsir, get_derivatives)
211211
EXPECT_NEAR(dydt_default[2], 25, 1e-12);
212212
}
213213

214+
// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
215+
TEST(Testsir, population_zero_no_nan)
216+
{
217+
// initialize simple model with total population 0
218+
mio::osir::Model<double> model(1);
219+
model.populations.set_total(0.0);
220+
221+
// call the get_derivatives function
222+
auto dydt_default = Eigen::VectorXd(Eigen::Index(mio::osir::InfectionState::Count));
223+
dydt_default.setZero();
224+
auto y0 = model.get_initial_values();
225+
model.get_derivatives(y0, y0, 0, dydt_default);
226+
227+
// check that there are now NaN values in dydt_default
228+
for (int i = 0; i < dydt_default.size(); i++) {
229+
EXPECT_FALSE(std::isnan(dydt_default[i]));
230+
}
231+
}
232+
214233
TEST(Testsir, Simulation)
215234
{
216235
double t0 = 0;
@@ -242,4 +261,4 @@ TEST(Testsir, Simulation)
242261
EXPECT_NEAR(results_t1[0], 150, 1e-12);
243262
EXPECT_NEAR(results_t1[1], 125, 1e-12);
244263
EXPECT_NEAR(results_t1[2], 125, 1e-12);
245-
}
264+
}

0 commit comments

Comments
 (0)