added check after initialization

Author: lenaploetzke

cpp/README.md

@@ -4,7 +4,7 @@ The MEmilio C++ library contains the implementation of the epidemiological model
 
 Directory structure:
 - memilio: framework for developing epidemiological models with, e.g., interregional mobility implementations, nonpharmaceutical interventions (NPIs), and  mathematical, programming, and IO utilities.
-- models: implementation of concrete models (ODE and ABM)
+- models: implementation of concrete models (ODE, IDE, LCT and ABM)
 - simulations: simulation applications that were used to generate the scenarios and data for publications
 - examples: small applications that help with using the framework and models
 - tests: unit tests for framework and models.

cpp/models/ide_secir/model.cpp

@@ -130,6 +130,19 @@ void Model::initialize(ScalarType dt)
                       "larger 0.");
         }
     }
+
+    // Verify that the initialization produces an appropriate result.
+    // Another check would be if the sum of the compartments is equal to N, but in all initialization methods one of the compartments is initialized via N - the others.
+    // This also means that if a compartment is greater than N, we will always have one or more compartments less than zero.
+    // Check if all compartments are non negative.
+    for (int i = 0; i < (int)InfectionState::Count; i++) {
+        if (m_populations[0][i] < 0) {
+            m_initialization_method = -2;
+            log_error("Initialization failed. One or more initial values for populations are less than zero. It may "
+                      "help to use a different method for initialization.");
+        }
+    }
+
     // Compute m_forceofinfection at time 0 needed for further simulation.
     update_forceofinfection(dt);
 }

cpp/models/ide_secir/model.h

@@ -59,18 +59,11 @@ class Model
     bool check_constraints(ScalarType dt) const
     {
         if (!((int)m_transitions.get_num_elements() == (int)InfectionTransition::Count)) {
-            log_error(
-                "Initialization failed. Number of elements in transition vector does not match the required number.");
+            log_error("A variable given for model construction is not valid. Number of elements in transition vector "
+                      "does not match the required number.");
             return true;
         }
 
-        for (int i = 0; i < (int)InfectionState::Count; i++) {
-            if (m_populations[0][i] < 0) {
-                log_error("Initialization failed. Initial values for populations are less than zero.");
-                return true;
-            }
-        }
-
         ScalarType support_max = std::max(
             {parameters.get<TransitionDistributions>()[(int)InfectionTransition::ExposedToInfectedNoSymptoms]
                  .get_support_max(dt, m_tol),
@@ -213,9 +206,10 @@ class Model
      * @returns 0 if the initialization method has not yet been selected,
      *      1 if the method using the total number of confirmed cases at time 0 is used,
      *      2 if the initialization is calculated using a prior set value for S,
-     *      3 if the initialization is calculated using a prior set value for R 
-     *      4 if the force of infection method is used,and
-     *      -1 if the initialization was not possible using any of the methods.
+     *      3 if the initialization is calculated using a prior set value for R,
+     *      4 if the force of infection method is used, 
+     *      -1 if the initialization was not possible using any of the methods and
+     *      -2 if the initialization was possible using one of the provided methods but the result is not appropriate.
      */
     int get_initialization_method()
     {

cpp/tests/test_ide_secir.cpp

@@ -289,50 +289,50 @@ TEST(IdeSecir, checkInitializations)
     contact_matrix[0]                              = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 0));
     parameters.get<mio::isecir::ContactPatterns>() = mio::UncertainContactMatrix(contact_matrix);
 
-    mio::TimeSeries<ScalarType> init_copy3(init);
-    mio::isecir::Model model3(std::move(init_copy3), N, deaths, 0, std::move(parameters));
+    mio::TimeSeries<ScalarType> init_copy2(init);
+    mio::isecir::Model model2(std::move(init_copy2), N, deaths, 0, std::move(parameters));
 
-    model3.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 5000;
-    model3.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 0;
+    model2.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 5000;
+    model2.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 0;
 
     // Carry out simulation.
-    mio::isecir::Simulation sim3(model3, 0, dt);
-    sim3.advance(tmax);
+    mio::isecir::Simulation sim2(model2, 0, dt);
+    sim2.advance(tmax);
 
     // Verify that the expected initialization method was used.
-    EXPECT_EQ(2, sim3.get_model().get_initialization_method());
+    EXPECT_EQ(2, sim2.get_model().get_initialization_method());
 
     // --- Case with R.
-    mio::TimeSeries<ScalarType> init_copy4(init);
-    mio::isecir::Model model4(std::move(init_copy4), N, deaths, 0, std::move(parameters));
+    mio::TimeSeries<ScalarType> init_copy3(init);
+    mio::isecir::Model model3(std::move(init_copy3), N, deaths, 0, std::move(parameters));
 
-    model4.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 0;
-    model4.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 1000;
+    model3.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 0;
+    model3.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 1000;
 
     // Carry out simulation.
-    mio::isecir::Simulation sim4(model4, 0, dt);
-    sim4.advance(tmax);
+    mio::isecir::Simulation sim3(model3, 0, dt);
+    sim3.advance(tmax);
 
     // Verify that the expected initialization method was used.
-    EXPECT_EQ(3, sim4.get_model().get_initialization_method());
+    EXPECT_EQ(3, sim3.get_model().get_initialization_method());
 
     // --- Case with forceofinfection.
-    mio::TimeSeries<ScalarType> init_copy2(init);
-    mio::isecir::Model model2(std::move(init_copy2), N, deaths, 0);
+    mio::TimeSeries<ScalarType> init_copy4(init);
+    mio::isecir::Model model4(std::move(init_copy4), N, deaths, 0);
 
     // Carry out simulation.
-    mio::isecir::Simulation sim2(model2, 0, dt);
-    sim2.advance(tmax);
+    mio::isecir::Simulation sim4(model4, 0, dt);
+    sim4.advance(tmax);
 
     // Verify that the expected initialization method was used.
-    EXPECT_EQ(4, sim2.get_model().get_initialization_method());
+    EXPECT_EQ(4, sim4.get_model().get_initialization_method());
 
     // --- Case without fitting initialization method.
-    // Deactivate temporarily log output for next test. Error is expected here.
+    // Deactivate temporarily log output for next tests. Errors are expected here.
     mio::set_log_level(mio::LogLevel::off);
 
-    // Here we do not need a copy of init as this is the last use of the vector. We can apply move directly.
-    mio::isecir::Model model5(std::move(init), N, deaths, 0, std::move(parameters));
+    mio::TimeSeries<ScalarType> init_copy5(init);
+    mio::isecir::Model model5(std::move(init_copy5), N, deaths, 0, std::move(parameters));
 
     model5.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 0;
     model5.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 0;
@@ -344,6 +344,22 @@ TEST(IdeSecir, checkInitializations)
     // Verify that initialization was not possible with one of the models methods.
     EXPECT_EQ(-1, sim5.get_model().get_initialization_method());
 
+    // --- Test with negative number of deaths.
+    deaths = -10;
+
+    // Here we do not need a copy of init as this is the last use of the vector. We can apply move directly.
+    mio::isecir::Model model6(std::move(init), N, deaths, 0, std::move(parameters));
+
+    model6.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 0;
+    model6.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered]   = 0;
+
+    // Carry out simulation.
+    mio::isecir::Simulation sim6(model6, 0, dt);
+    sim6.advance(tmax);
+
+    // Verify that initialization was possible but the result is not appropriate.
+    EXPECT_EQ(-2, sim6.get_model().get_initialization_method());
+
     // Reactive log output.
     mio::set_log_level(mio::LogLevel::warn);
 }
@@ -383,8 +399,7 @@ TEST(IdeSecir, testModelConstraints)
     auto constraint_check = model_wrong_size.check_constraints(dt);
     EXPECT_TRUE(constraint_check);
 
-    // --- Test with negative number of deaths.
-    deaths = -10;
+    // --- Test with too few time points.
     // Create TimeSeries with num_transitions elements.
     mio::TimeSeries<ScalarType> init(num_transitions);
     // Add time points for initialization of transitions.
@@ -395,16 +410,6 @@ TEST(IdeSecir, testModelConstraints)
         init.add_time_point(init.get_last_time() + dt, vec_init);
     }
 
-    // Initialize a model.
-    mio::TimeSeries<ScalarType> init_copy(init);
-    mio::isecir::Model model_negative_deaths(std::move(init_copy), N, deaths);
-
-    // Return true for negative entry in m_populations.
-    constraint_check = model_negative_deaths.check_constraints(dt);
-    EXPECT_TRUE(constraint_check);
-
-    // --- Test with too few time points.
-    deaths = 10;
     // Initialize a model.
     mio::isecir::Model model_few_timepoints(std::move(init), N, deaths);