958 implement a check for the results of the initialization (#959)

- Optimization of IDE SECIR initialization to avoid overwriting of beforehand set values.
- Check and return of initialization could be conducted properly.
- Added unit tests to check_constraints() function to increase code coverage.

Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com>
Co-authored-by: Martin J. Kühn <62713180+mknaranja@users.noreply.github.com>

Author: 3 people

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

@@ -66,14 +66,51 @@ void Model::initialize(ScalarType dt)
             m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
             m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
     }
-    else {
+    else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) {
+        // Take initialized value for Susceptibles if value can't be calculated via the standard formula.
+        m_initialization_method = 2;
+        // Calculate other compartment sizes for t=0.
+        other_compartments_current_timestep(dt);
 
+        // R; need an initial value for R, therefore do not calculate via compute_recovered()
+        m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
+            m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
+    }
+    else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) {
+        // If value for Recovered is initialized and standard method is not applicable (i.e., no value for total infections
+        //  or suceptibles given directly), calculate Susceptibles via other compartments.
+        // The calculation of other compartments' values is not dependent on Susceptibles at time 0, i.e., S(0), but only on the transitions of the past.
+        m_initialization_method = 3;
+        other_compartments_current_timestep(dt);
+
+        m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
+            m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
+            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
+    }
+    else {
         // compute Susceptibles at time 0  and m_forceofinfection at time -dt as initial values for discretization scheme
         // use m_forceofinfection at -dt to be consistent with further calculations of S (see compute_susceptibles()),
         // where also the value of m_forceofinfection for the previous timestep is used
         update_forceofinfection(dt, true);
         if (m_forceofinfection > 1e-12) {
-            m_initialization_method = 2;
+            m_initialization_method = 4;
+
+            /* Attention: With an inappropriate combination of parameters and initial conditions, it can happen that S 
+            becomes greater than N when this formula is used. In this case, at least one compartment is initialized 
+            with a number less than zero, so that a log_error is thrown.
+            However, this initialization method is consistent with the numerical solver of the model equations,
+            so it may sometimes make sense to use this method. */
             m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
                 m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] /
                 (dt * m_forceofinfection);
@@ -91,46 +128,27 @@ void Model::initialize(ScalarType dt)
                 m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
                 m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
         }
-        else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) {
-            //take initialized value for Susceptibles if value can't be calculated via the standard formula
-            m_initialization_method = 3;
-            //calculate other compartment sizes for t=0
-            other_compartments_current_timestep(dt);
-
-            //R; need an initial value for R, therefore do not calculate via compute_recovered()
-            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] =
-                m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
-        }
-        else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) {
-            //if value for Recovered is initialized and standard method is not applicable, calculate Susceptibles via other compartments
-            //determining other compartment sizes is not dependent of Susceptibles(0), just of the transitions of the past.
-            //calculate other compartment sizes for t=0
-            m_initialization_method = 4;
-            other_compartments_current_timestep(dt);
-
-            m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] =
-                m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] -
-                m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)];
-        }
         else {
             m_initialization_method = -1;
             log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither "
                       "Susceptibles nor Recovered or total confirmed cases for time 0 were set. One of them should be "
                       "larger 0.");
         }
     }
-    // compute m_forceofinfection at time 0 needed for further simulation
+
+    // 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

@@ -53,23 +53,15 @@ class Model
           const ParameterSet& Parameterset_init = ParameterSet());
 
     /**
-    * @brief Checks constraints on model parameters.
+    * @brief Checks constraints on model parameters and initial data.
+    * @return Returns true if one (or more) constraint(s) are not satisfied, otherwise false.
     */
-    void check_constraints(ScalarType dt) const
+    bool check_constraints(ScalarType dt) const
     {
-        if (!(m_populations.get_num_time_points() > 0)) {
-            log_error("Model construction failed. No initial time point for populations.");
-        }
-
-        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.");
-            }
-        }
-
         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;
         }
 
         ScalarType support_max = std::max(
@@ -95,17 +87,25 @@ class Model
         if (m_transitions.get_num_time_points() < (Eigen::Index)std::ceil(support_max / dt)) {
             log_error(
                 "Initialization failed. Not enough time points for transitions given before start of simulation.");
+            return true;
         }
 
-        parameters.check_constraints();
+        return parameters.check_constraints();
     }
 
     /**
-     * @brief Calculate the number of individuals in each compartment for time 0.
-     * 
-     * Initial transitions are used to calculate the initial compartment sizes.
-     * @param[in] dt Time discretization step size.         
-     */
+    * @brief Initializes the model and sets compartment values at time zero.
+    *
+    * The initialization method is selected automatically based on the different values that need to be set beforehand.
+    * Infection compartments are always computed through historic flow.
+    * Initialization methods for susceptibles and recovered are tested in the following order:
+    * 1.) If a positive number for the total number of confirmed cases is set, recovered is set according to that value and #Susceptible%s are derived.
+    * 2.) If #Susceptible%s are set, recovered will be derived.
+    * 3.) If #Recovered are set directly, #Susceptible%s are derived.
+    * 4.) If none of the above is set with positive value, the force of infection is used as in Messina et al (2021) to set the #Susceptible%s.
+    *
+    * @param[in] dt Time discretization step size.         
+    */
     void initialize(ScalarType dt);
 
     /**
@@ -208,14 +208,15 @@ class Model
     }
 
     /**
-     * @brief Specifies a number associated with the method used for initialization.
+     * @brief Returns the number associated with the method selected automatically for initialization.
      *
-     * @returns 0 if the initialization method has not yet been selected,
+    * @returns 0 if the model has not yet been initialized and the 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 force of infection method is used,
-     *      3 if the initialization is calculated using a prior set value for S,
-     *      4 if the initialization is calculated using a prior set value for R and
-     *      -1 if the initialization was not possible using any of the methods.
+     *      2 if the initialization is calculated using a beforehand set value for S,
+     *      3 if the initialization is calculated using a beforehand 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()
     {