1102 divNj creates NaN values if subpopulation is 0

cpp/benchmarks/flow_simulation_ode_secirvvs.h

@@ -203,7 +203,7 @@ class FlowlessModel : public CompartmentalModel<ScalarType, osecirvvs::Infection
                             pop[ICrPIj] + pop[INSPICj] + pop[ISyPICj] + pop[SIIj] + pop[EIIj] + pop[INSIIj] +
                             pop[ISyIIj] + pop[ISevIIj] + pop[ICrIIj] + pop[INSIICj] + pop[ISyIICj];
 
-                double divNj = 1.0 / Nj; // precompute 1.0/Nj
+                const double divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
 
                 double ext_inf_force_dummy =
                     cont_freq_eff * divNj *

cpp/memilio/config.h

@@ -26,7 +26,48 @@
 #define MIO_CONFIG_H
 
 #include "memilio/config_internal.h"
+#include <type_traits>
 
 using ScalarType = double;
 
+namespace mio
+{
+/**
+ * @brief Type specific limits for floating-points.
+ *
+ * Specializations are provided for `float` and `double` types. In order to use
+ * other floating point types, a new specialization for that type must be added.
+ * 
+ * Member functions of the unspecialized class are disabled. Using them results
+ * in compile time errors.
+ * @{
+ */
+/**
+ * @tparam FP A floating point type.
+ */
+template <typename FP>
+struct Limits {
+    static constexpr FP zero_tolerance() = delete;
+};
+
+template <>
+struct Limits<float> {
+    /// @brief Returns the limit under which a float may be rounded down to zero.
+    static constexpr float zero_tolerance()
+    {
+        return 1e-6f;
+    }
+};
+
+template <>
+struct Limits<double> {
+    /// @brief Returns the limit under which a double may be rounded down to zero.
+    static constexpr double zero_tolerance()
+    {
+        return 1e-12;
+    }
+};
+/** @} */
+
+} // namespace mio
 #endif

cpp/models/ode_secir/model.h

@@ -122,31 +122,31 @@ class Model : public FlowModel<FP, InfectionState, Populations<FP, AgeGroup, Inf
                                     params.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[j]);
 
                 // effective contact rate by contact rate between groups i and j and damping j
-                double season_val =
+                ScalarType season_val =
                     (1 + params.template get<Seasonality<FP>>() *
                              sin(3.141592653589793 * ((params.template get<StartDay>() + t) / 182.5 + 0.5)));
-                double cont_freq_eff =
+                ScalarType cont_freq_eff =
                     season_val * contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
                                                                  static_cast<Eigen::Index>((size_t)j));
-                double Nj =
+                ScalarType Nj =
                     pop[Sj] + pop[Ej] + pop[INSj] + pop[ISyj] + pop[ISevj] + pop[ICrj] + pop[Rj]; // without died people
-                double divNj   = 1.0 / Nj; // precompute 1.0/Nj
-                double dummy_S = y[Si] * cont_freq_eff * divNj *
-                                 params.template get<TransmissionProbabilityOnContact<FP>>()[i] *
-                                 (params.template get<RelativeTransmissionNoSymptoms<FP>>()[j] * pop[INSj] +
-                                  riskFromInfectedSymptomatic * pop[ISyj]);
+                const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
+                ScalarType dummy_S     = y[Si] * cont_freq_eff * divNj *
+                                     params.template get<TransmissionProbabilityOnContact<FP>>()[i] *
+                                     (params.template get<RelativeTransmissionNoSymptoms<FP>>()[j] * pop[INSj] +
+                                      riskFromInfectedSymptomatic * pop[ISyj]);
 
                 // Susceptible -> Exposed
                 flows[this->template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Exposed>({i})] +=
                     dummy_S;
             }
 
             // ICU capacity shortage is close
-            double criticalPerSevereAdjusted = smoother_cosine(
+            ScalarType criticalPerSevereAdjusted = smoother_cosine(
                 icu_occupancy, 0.90 * params.template get<ICUCapacity<FP>>(), params.template get<ICUCapacity<FP>>(),
                 params.template get<CriticalPerSevere<FP>>()[i], 0);
 
-            double deathsPerSevereAdjusted =
+            ScalarType deathsPerSevereAdjusted =
                 params.template get<CriticalPerSevere<FP>>()[i] - criticalPerSevereAdjusted;
 
             // Exposed -> InfectedNoSymptoms

cpp/models/ode_secirvvs/model.h

@@ -242,8 +242,7 @@ class Model
                         pop[ISyNCj] + pop[SPIj] + pop[EPIj] + pop[INSPIj] + pop[ISyPIj] + pop[ISevPIj] + pop[ICrPIj] +
                         pop[INSPICj] + pop[ISyPICj] + pop[SIIj] + pop[EIIj] + pop[INSIIj] + pop[ISyIIj] + pop[ISevIIj] +
                         pop[ICrIIj] + pop[INSIICj] + pop[ISyIICj];
-
-                FP divNj = 1.0 / Nj; // precompute 1.0/Nj
+                const FP divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
 
                 FP ext_inf_force_dummy = cont_freq_eff * divNj *
                                          params.template get<TransmissionProbabilityOnContact<FP>>()[(AgeGroup)i] *

cpp/models/ode_seir/model.h

@@ -85,10 +85,11 @@ class Model
                 const size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected});
                 const size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
 
-                const double Nj_inv = 1.0 / (pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj]);
-                const double coeffStoE =
+                const ScalarType Nj    = pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj];
+                const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
+                const ScalarType coeffStoE =
                     params.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_matrix_at(t)(i.get(), j.get()) *
-                    params.template get<TransmissionProbabilityOnContact<FP>>()[i] * Nj_inv;
+                    params.template get<TransmissionProbabilityOnContact<FP>>()[i] * divNj;
 
                 flows[Base::template get_flat_flow_index<InfectionState::Susceptible, InfectionState::Exposed>(i)] +=
                     coeffStoE * y[Si] * pop[Ij];
@@ -127,16 +128,16 @@ class Model
             size_t Si = this->populations.get_flat_index({i, InfectionState::Susceptible});
             for (auto j = AgeGroup(0); j < AgeGroup(num_groups); j++) {
 
-                double Nj    = this->populations.get_group_total(j);
-                double divNj = 1.0 / Nj;
+                const ScalarType Nj    = this->populations.get_group_total(j);
+                const ScalarType divNj = (Nj < 1e-12) ? 0.0 : 1.0 / Nj;
 
-                double coeffStoE = contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) *
-                                   params.template get<TransmissionProbabilityOnContact<ScalarType>>()[i] * divNj;
+                ScalarType coeffStoE = contact_matrix.get_matrix_at(y.get_time(t_idx))(i.get(), j.get()) *
+                                       params.template get<TransmissionProbabilityOnContact<ScalarType>>()[i] * divNj;
                 F((size_t)i, (size_t)j + num_groups) = coeffStoE * y.get_value(t_idx)[Si];
             }
 
-            double T_Ei                          = params.template get<mio::oseir::TimeExposed<ScalarType>>()[i];
-            double T_Ii                          = params.template get<mio::oseir::TimeInfected<ScalarType>>()[i];
+            ScalarType T_Ei                      = params.template get<mio::oseir::TimeExposed<ScalarType>>()[i];
+            ScalarType T_Ii                      = params.template get<mio::oseir::TimeInfected<ScalarType>>()[i];
             V((size_t)i, (size_t)i)              = 1.0 / T_Ei;
             V((size_t)i + num_groups, (size_t)i) = -1.0 / T_Ei;
             V((size_t)i + num_groups, (size_t)i + num_groups) = 1.0 / T_Ii;

cpp/models/ode_sir/model.h

@@ -77,11 +77,12 @@ class Model
                 size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected});
                 size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
 
-                double Nj = pop[Sj] + pop[Ij] + pop[Rj];
+                const ScalarType Nj    = pop[Sj] + pop[Ij] + pop[Rj];
+                const ScalarType divNj = (Nj < Limits<ScalarType>::zero_tolerance()) ? 0.0 : 1.0 / Nj;
 
-                double coeffStoI = contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
-                                                                   static_cast<Eigen::Index>((size_t)j)) *
-                                   params.template get<TransmissionProbabilityOnContact<FP>>()[i] / Nj;
+                ScalarType coeffStoI = contact_matrix.get_matrix_at(t)(static_cast<Eigen::Index>((size_t)i),
+                                                                       static_cast<Eigen::Index>((size_t)j)) *
+                                       params.template get<TransmissionProbabilityOnContact<FP>>()[i] * divNj;
 
                 dydt[Si] += -coeffStoI * y[Si] * pop[Ij];
                 dydt[Ii] += coeffStoI * y[Si] * pop[Ij];

cpp/tests/test_odesecir.cpp

@@ -572,6 +572,25 @@ TEST(TestOdeSecir, testSettersAndGetters)
     EXPECT_EQ(vec[21], model.parameters.get<mio::osecir::Seasonality<double>>());
 }
 
+// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
+TEST(TestOdeSecir, population_zero_no_nan)
+{
+    // initialize simple model with total population 0
+    mio::osecir::Model<double> model(1);
+    model.populations.set_total(0.0);
+
+    // call the get_flows function
+    auto dydt_default = Eigen::VectorXd(15);
+    dydt_default.setZero();
+    auto y0 = model.get_initial_values();
+    model.get_flows(y0, y0, 0, dydt_default);
+
+    // check that there are now NaN values in dydt_default
+    for (int i = 0; i < dydt_default.size(); i++) {
+        EXPECT_FALSE(std::isnan(dydt_default[i]));
+    }
+}
+
 TEST(TestOdeSecir, testDamping)
 {
     // Test functionality of dampings

cpp/tests/test_odesecirvvs.cpp

@@ -901,18 +901,20 @@ TEST(TestOdeSECIRVVS, test_commuters)
     model.parameters.get_commuter_nondetection()    = non_detection_factor;
     auto sim                                        = mio::osecirvvs::Simulation<>(model);
     auto before_testing                             = sim.get_result().get_last_value().eval();
-    auto mobile_population                                      = (sim.get_result().get_last_value() * mobility_factor).eval();
-    auto mobile_population_tested                               = mobile_population.eval();
+    auto mobile_population                          = (sim.get_result().get_last_value() * mobility_factor).eval();
+    auto mobile_population_tested                   = mobile_population.eval();
 
     mio::osecirvvs::test_commuters<double>(sim, mobile_population_tested, 0.0);
 
     ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)],
-                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] * non_detection_factor,
+                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] *
+                    non_detection_factor,
                 1e-5);
     ASSERT_NEAR(
         sim.get_result().get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaiveConfirmed)],
         before_testing[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaiveConfirmed)] +
-            mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] * (1 - non_detection_factor),
+            mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsNaive)] *
+                (1 - non_detection_factor),
         1e-5);
     ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity)],
                 mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity)] *
@@ -925,10 +927,11 @@ TEST(TestOdeSECIRVVS, test_commuters)
             mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsPartialImmunity)] *
                 (1 - non_detection_factor),
         1e-5);
-    ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)],
-                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)] *
-                    non_detection_factor,
-                1e-5);
+    ASSERT_NEAR(
+        mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)],
+        mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunity)] *
+            non_detection_factor,
+        1e-5);
     ASSERT_NEAR(
         sim.get_result()
             .get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedSymptomsImprovedImmunityConfirmed)],
@@ -938,29 +941,32 @@ TEST(TestOdeSECIRVVS, test_commuters)
         1e-5);
 
     ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)],
-                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)] * non_detection_factor,
+                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)] *
+                    non_detection_factor,
                 1e-5);
     ASSERT_NEAR(sim.get_result()
                     .get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaiveConfirmed)],
                 before_testing[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaiveConfirmed)] +
                     mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsNaive)] *
                         (1 - non_detection_factor),
                 1e-5);
-    ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)],
-                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)] *
-                    non_detection_factor,
-                1e-5);
+    ASSERT_NEAR(
+        mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)],
+        mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)] *
+            non_detection_factor,
+        1e-5);
     ASSERT_NEAR(
         sim.get_result()
             .get_last_value()[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed)],
         before_testing[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed)] +
             mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsPartialImmunity)] *
                 (1 - non_detection_factor),
         1e-5);
-    ASSERT_NEAR(mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)],
-                mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)] *
-                    non_detection_factor,
-                1e-5);
+    ASSERT_NEAR(
+        mobile_population_tested[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)],
+        mobile_population[Eigen::Index(mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunity)] *
+            non_detection_factor,
+        1e-5);
     ASSERT_NEAR(
         sim.get_result().get_last_value()[Eigen::Index(
             mio::osecirvvs::InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed)],
@@ -970,6 +976,25 @@ TEST(TestOdeSECIRVVS, test_commuters)
         1e-5);
 }
 
+// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
+TEST(TestOdeSECIRVVS, population_zero_no_nan)
+{
+    // initialize simple model with total population 0
+    mio::osecirvvs::Model<double> model(1);
+    model.populations.set_total(0.0);
+
+    // call the get_flows function
+    auto dydt_default = Eigen::VectorXd(45);
+    dydt_default.setZero();
+    auto y0 = model.get_initial_values();
+    model.get_flows(y0, y0, 0, dydt_default);
+
+    // check that there are now NaN values in dydt_default
+    for (int i = 0; i < dydt_default.size(); i++) {
+        EXPECT_FALSE(std::isnan(dydt_default[i]));
+    }
+}
+
 TEST(TestOdeSECIRVVS, check_constraints_parameters)
 {
     auto model = mio::osecirvvs::Model<double>(1);

cpp/tests/test_odeseir.cpp

@@ -311,6 +311,25 @@ TEST(TestOdeSeir, get_reproduction_number)
     EXPECT_NEAR(model.get_reproduction_number(0.9, result).value(), 1.858670429549998504, 1e-12);
 }
 
+// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
+TEST(TestSeir, population_zero_no_nan)
+{
+    // initialize simple model with total population 0
+    mio::oseir::Model<double> model(1);
+    model.populations.set_total(0.0);
+
+    // call the get_flows function
+    auto dydt_default = Eigen::VectorXd(3);
+    dydt_default.setZero();
+    auto y0 = model.get_initial_values();
+    model.get_flows(y0, y0, 0, dydt_default);
+
+    // check that there are now NaN values in dydt_default
+    for (int i = 0; i < dydt_default.size(); i++) {
+        EXPECT_FALSE(std::isnan(dydt_default[i]));
+    }
+}
+
 TEST(TestSeir, get_flows)
 {
     mio::oseir::Model<double> model(1);

cpp/tests/test_odesir.cpp

@@ -211,6 +211,25 @@ TEST(Testsir, get_derivatives)
     EXPECT_NEAR(dydt_default[2], 25, 1e-12);
 }
 
+// Test model initialization with total population of 0 and ensure get_flows returns no NaN values
+TEST(Testsir, population_zero_no_nan)
+{
+    // initialize simple model with total population 0
+    mio::osir::Model<double> model(1);
+    model.populations.set_total(0.0);
+
+    // call the get_derivatives function
+    auto dydt_default = Eigen::VectorXd(Eigen::Index(mio::osir::InfectionState::Count));
+    dydt_default.setZero();
+    auto y0 = model.get_initial_values();
+    model.get_derivatives(y0, y0, 0, dydt_default);
+
+    // check that there are now NaN values in dydt_default
+    for (int i = 0; i < dydt_default.size(); i++) {
+        EXPECT_FALSE(std::isnan(dydt_default[i]));
+    }
+}
+
 TEST(Testsir, Simulation)
 {
     double t0   = 0;
@@ -242,4 +261,4 @@ TEST(Testsir, Simulation)
     EXPECT_NEAR(results_t1[0], 150, 1e-12);
     EXPECT_NEAR(results_t1[1], 125, 1e-12);
     EXPECT_NEAR(results_t1[2], 125, 1e-12);
-}
+}