793 add seasonality to ide model (#972)

lenaploetzke · reneSchm · web-flow · commit a5cbb91c68bf · 2024-04-02T12:17:05.000+02:00
- Functionality to simulate the IDE model with a seasonality has been added.
- The test for the function check_constraints() of the Parameters class was rewritten. The use of the class Model made the test unnecessarily complicated.

Co-authored-by: reneSchm &lt;49305466+reneSchm@users.noreply.github.com&gt;
diff --git a/cpp/examples/ide_secir.cpp b/cpp/examples/ide_secir.cpp
@@ -91,11 +91,14 @@ int main()
     model.parameters.set<mio::isecir::TransmissionProbabilityOnContact>(prob);
     model.parameters.set<mio::isecir::RelativeTransmissionNoSymptoms>(prob);
     model.parameters.set<mio::isecir::RiskOfInfectionFromSymptomatic>(prob);
+    model.parameters.set<mio::isecir::Seasonality>(0.1);
+    // Start the simulation on the 40th day of a year (i.e. in February).
+    model.parameters.set<mio::isecir::StartDay>(40);
 
     model.check_constraints(dt);
 
     // Carry out simulation.
-    mio::isecir::Simulation sim(model, 0, dt);
+    mio::isecir::Simulation sim(model, dt);
     sim.advance(tmax);
 
     auto interpolated_results = mio::interpolate_simulation_result(sim.get_result(), dt / 2);
diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp
@@ -275,9 +275,12 @@ void Model::update_forceofinfection(ScalarType dt, bool initialization)
     for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) {
 
         ScalarType state_age = (num_time_points - 1 - i) * dt;
-
+        ScalarType season_val =
+            1 +
+            parameters.get<Seasonality>() *
+                sin(3.141592653589793 * (std::fmod((parameters.get<StartDay>() + current_time), 365.0) / 182.5 + 0.5));
         m_forceofinfection +=
-            parameters.get<TransmissionProbabilityOnContact>().eval(state_age) *
+            season_val * parameters.get<TransmissionProbabilityOnContact>().eval(state_age) *
             parameters.get<ContactPatterns>().get_cont_freq_mat().get_matrix_at(current_time)(0, 0) *
             ((parameters
                       .get<TransitionProbabilities>()[(int)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms] *
diff --git a/cpp/models/ide_secir/parameters.h b/cpp/models/ide_secir/parameters.h
@@ -150,10 +150,48 @@ struct RiskOfInfectionFromSymptomatic {
     }
 };
 
-// Define Parameterset for IDE SECIR model.
+/**
+ * @brief Sets the day in a year at which a simulation with an IDE-SECIR model is started.
+ *
+ * The value 0.0 corresponds to the 1st of January at 0:00 am, 31.25 is the 1st of February at 6:00 am, and so on.
+ * The start day defines in which season the simulation is started.
+ * If the start day is 180 and simulation takes place from t0=0 to
+ * tmax=100 the days 180 to 280 of the year are simulated.
+ * The parameter is used in the formula of the seasonality in the model.
+ */
+struct StartDay {
+    using Type = ScalarType;
+    static Type get_default()
+    {
+        return 0.;
+    }
+    static std::string name()
+    {
+        return "StartDay";
+    }
+};
+
+/**
+ * @brief The seasonality parameter k in the IDE-SECIR model.
+ * The formula for the seasonality used in the model is given as (1+k*sin()) where the sine
+ * curve is below one in summer and above one in winter.
+ */
+struct Seasonality {
+    using Type = ScalarType;
+    static Type get_default()
+    {
+        return Type(0.);
+    }
+    static std::string name()
+    {
+        return "Seasonality";
+    }
+};
+
+// Define Parameterset for IDE-SECIR model.
 using ParametersBase =
     ParameterSet<TransitionDistributions, TransitionProbabilities, ContactPatterns, TransmissionProbabilityOnContact,
-                 RelativeTransmissionNoSymptoms, RiskOfInfectionFromSymptomatic>;
+                 RelativeTransmissionNoSymptoms, RiskOfInfectionFromSymptomatic, StartDay, Seasonality>;
 
 /**
  * @brief Parameters of an age-resolved SECIR/SECIHURD model.
@@ -269,6 +307,7 @@ class Parameters : public ParametersBase
                       1);
             return true;
         }
+
         /* The first entry of TransitionDistributions is not checked because the distribution S->E is never used 
         (and it makes no sense to use the distribution). The support does not need to be valid.*/
         for (size_t i = 1; i < (int)InfectionTransition::Count; i++) {
@@ -278,6 +317,11 @@ class Parameters : public ParametersBase
             }
         }
 
+        if (this->get<Seasonality>() < 0.0 || this->get<Seasonality>() > 0.5) {
+            log_warning("Constraint check: Parameter Seasonality should lie between {:0.4f} and {:.4f}", 0.0, 0.5);
+            return true;
+        }
+
         return false;
     }
 
diff --git a/cpp/models/ide_secir/simulation.cpp b/cpp/models/ide_secir/simulation.cpp
@@ -59,10 +59,10 @@ void Simulation::advance(ScalarType tmax)
     }
 }
 
-TimeSeries<ScalarType> simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& m_model)
+TimeSeries<ScalarType> simulate(ScalarType tmax, ScalarType dt, Model const& m_model)
 {
     m_model.check_constraints(dt);
-    Simulation sim(m_model, t0, dt);
+    Simulation sim(m_model, dt);
     sim.advance(tmax);
     return sim.get_result();
 }
diff --git a/cpp/models/ide_secir/simulation.h b/cpp/models/ide_secir/simulation.h
@@ -44,12 +44,10 @@ class Simulation
     /**
      * @brief setup the Simulation for an IDE model.
      * @param[in] model An instance of the IDE model.
-     * @param[in] t0 Start time.
      * @param[in] dt Step size of numerical solver.
      */
-    Simulation(Model const& model, ScalarType t0 = 0., ScalarType dt = 0.1)
+    Simulation(Model const& model, ScalarType dt = 0.1)
         : m_model(std::make_unique<Model>(model))
-        , m_t0(t0)
         , m_dt(dt)
     {
     }
@@ -106,15 +104,6 @@ class Simulation
         return *m_model;
     }
 
-    /**
-     * @brief get the starting time of the simulation.
-     * 
-     */
-    ScalarType get_t0()
-    {
-        return m_t0;
-    }
-
     /**
      * @brief get the time step of the simulation.
      * 
@@ -126,19 +115,18 @@ class Simulation
 
 private:
     std::unique_ptr<Model> m_model; ///< Unique pointer to the Model simulated.
-    ScalarType m_t0; ///< Start time used for simulation.
     ScalarType m_dt; ///< Time step used for numerical computations in simulation.
 };
 
 /**
- * @brief simulates a compartmental model
- * @param[in] t0 start time
- * @param[in] tmax end time
- * @param[in] dt initial step size of integration
- * @param[in] model an instance of a compartmental model
- * @return a TimeSeries to represent the final simulation result
+ * @brief Run a Simulation of an IDE-SECIR model.
+ *
+ * @param[in] tmax End time.
+ * @param[in] dt Initial step size of integration.
+ * @param[in] model An instance of an IDE-SECIR model.
+ * @return A TimeSeries to represent the final simulation result.
  */
-TimeSeries<ScalarType> simulate(double t0, double tmax, double dt, Model const& model);
+TimeSeries<ScalarType> simulate(double tmax, double dt, Model const& model);
 
 } // namespace isecir
 } // namespace mio
diff --git a/cpp/tests/test_ide_secir.cpp b/cpp/tests/test_ide_secir.cpp

Original file line number	Diff line number	Diff line change
`@@ -59,10 +59,10 @@ void Simulation::advance(ScalarType tmax)`
`59`	`59`	`}`
`60`	`60`	`}`
`61`	`61`
`62`		`-TimeSeries<ScalarType> simulate(ScalarType t0, ScalarType tmax, ScalarType dt, Model const& m_model)`
	`62`	`+TimeSeries<ScalarType> simulate(ScalarType tmax, ScalarType dt, Model const& m_model)`
`63`	`63`	`{`
`64`	`64`	`m_model.check_constraints(dt);`
`65`		`- Simulation sim(m_model, t0, dt);`
	`65`	`+ Simulation sim(m_model, dt);`
`66`	`66`	`sim.advance(tmax);`
`67`	`67`	`return sim.get_result();`
`68`	`68`	`}`