SciCompMod · annawendler · Nov 21, 2022 · Nov 28, 2022 · Nov 28, 2022 · Dec 2, 2022
diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt
@@ -112,6 +112,9 @@ add_executable(history_example history.cpp)
 target_link_libraries(history_example PRIVATE memilio)
 target_compile_options(history_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ide_ode_secir_example ide_ode_secir.cpp)
+target_link_libraries(ide_ode_secir_example PRIVATE memilio ide_secir)
+
 if(MEMILIO_HAS_JSONCPP)
     add_executable(ode_secir_read_graph_example ode_secir_read_graph.cpp)
     target_link_libraries(ode_secir_read_graph_example PRIVATE memilio ode_secir)

diff --git a/cpp/examples/ide_initialization.cpp b/cpp/examples/ide_initialization.cpp
@@ -18,7 +18,7 @@
 * limitations under the License.
 */
 
-#include "ide_secir/model.h"
+#include "ide_secir/model_ide.h"
 #include "ide_secir/infection_state.h"
 #include "ide_secir/simulation.h"
 #include "ide_secir/parameters_io.h"

diff --git a/cpp/examples/ide_ode_secir.cpp b/cpp/examples/ide_ode_secir.cpp
diff --git a/cpp/examples/ide_secir.cpp b/cpp/examples/ide_secir.cpp
@@ -18,7 +18,7 @@
 * limitations under the License.
 */
 
-#include "ide_secir/model.h"
+#include "ide_secir/model_ide.h"
 #include "ide_secir/infection_state.h"
 #include "ide_secir/simulation.h"
 #include "memilio/config.h"

diff --git a/cpp/models/ide_secir/CMakeLists.txt b/cpp/models/ide_secir/CMakeLists.txt
@@ -1,10 +1,12 @@
 add_library(ide_secir
     infection_state.h
-    model.h
-    model.cpp
+    model_ide.h
+    model_ide.cpp
     simulation.h
     simulation.cpp
     parameters.h
+    initialize_from_ode.h
+    initialize_from_ode.cpp
     parameters_io.h
     parameters_io.cpp
 )

diff --git a/cpp/models/ide_secir/initialize_from_ode.cpp b/cpp/models/ide_secir/initialize_from_ode.cpp
@@ -0,0 +1,160 @@
+/* 
+* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
+*
+* Authors: Anna Wendler
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+#include "ide_secir/model_ide.h"
+#include "ide_secir/parameters.h"
+#include "ode_secir/model.h"
+#include "infection_state.h"
+#include "memilio/config.h"
+#include "memilio/utils/logging.h"
+#include "memilio/math/eigen.h"
+
+namespace mio
+{
+namespace isecir
+{
+void get_flows_from_ode_compartments(mio::osecir::Model<ScalarType>& model_ode,
+                                     mio::TimeSeries<ScalarType> compartments, mio::TimeSeries<ScalarType>& flows,
+                                     ScalarType t_max, ScalarType t_window, ScalarType dt_reference,
+                                     ScalarType dt_comparison)
+{ //TODO: Check that flows is an empty TimeSeries at the beginning.
+    int num_transitions = (int)mio::isecir::InfectionTransition::Count;
+
+    // Use scale_timesteps to get from index wrt ODE timestep to index wrt IDE timestep.
+    // Here we assume that we solve the ODE model on a finer scale (or equal scale) than the IDE model.
+
+    if (dt_comparison < 1e-10) {
+        dt_comparison = dt_reference;
+    }
+
+    ScalarType scale_timesteps = dt_comparison / dt_reference;
+
+    // Compute index variables with respect to dt_comparison.
+    Eigen::Index t_window_index = Eigen::Index(std::ceil(t_window / dt_comparison));
+    Eigen::Index t_max_index    = Eigen::Index(std::ceil(t_max / dt_comparison));
+
+    Eigen::Index flows_start_index = t_max_index - t_window_index + 1;
+
+    // Add time points to TimeSeries containing flows.
+    for (Eigen::Index i = flows_start_index; i <= t_max_index; i++) {
+        flows.add_time_point(i * dt_comparison, mio::TimeSeries<ScalarType>::Vector::Constant(num_transitions, 0));
+        flows.get_last_value()[Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)] +=
+            compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::Susceptible)] -
+            compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::Susceptible)];
+    }
+
+    // compute resulting flows as combination of change in compartments and previously computed flows
+
+    // flow from E to C
+    for (int i = flows_start_index; i <= t_max_index; i++) {
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)] =
+            compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::Exposed)] -
+            compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::Exposed)] +
+            flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::SusceptibleToExposed)];
+    }
+
+    // flow from C to I and from C to R
+    for (int i = flows_start_index; i <= t_max_index; i++) {
+        flows[i -
+              flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] =
+            (1 -
+             model_ode.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<ScalarType>>()[(mio::AgeGroup)0]) *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedNoSymptoms)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedNoSymptoms)] +
+             flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)]);
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToRecovered)] =
+            (model_ode.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<ScalarType>>()[(mio::AgeGroup)0]) *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedNoSymptoms)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedNoSymptoms)] +
+             flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::ExposedToInfectedNoSymptoms)]);
+    }
+
+    // flow from I to H and from I to R
+    for (int i = flows_start_index; i <= t_max_index; i++) {
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere)] =
+            model_ode.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[(mio::AgeGroup)0] *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedSymptoms)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedSymptoms)] +
+             flows[i - flows_start_index]
+                  [Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]);
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToRecovered)] =
+            (1 - model_ode.parameters.get<mio::osecir::SeverePerInfectedSymptoms<ScalarType>>()[(mio::AgeGroup)0]) *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedSymptoms)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedSymptoms)] +
+             flows[i - flows_start_index]
+                  [Eigen::Index(mio::isecir::InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]);
+    }
+
+    // flow from H to U and from H to R
+    for (int i = flows_start_index; i <= t_max_index; i++) {
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical)] =
+            model_ode.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[(mio::AgeGroup)0] *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedSevere)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedSevere)] +
+             flows[i - flows_start_index]
+                  [Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere)]);
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToRecovered)] =
+            (1 - model_ode.parameters.get<mio::osecir::CriticalPerSevere<ScalarType>>()[(mio::AgeGroup)0]) *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedSevere)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedSevere)] +
+             flows[i - flows_start_index]
+                  [Eigen::Index(mio::isecir::InfectionTransition::InfectedSymptomsToInfectedSevere)]);
+    }
+
+    // flow from U to D and from U to R
+    for (int i = flows_start_index; i <= t_max_index; i++) {
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToDead)] =
+            model_ode.parameters.get<mio::osecir::DeathsPerCritical<ScalarType>>()[(mio::AgeGroup)0] *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedCritical)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedCritical)] +
+             flows[i - flows_start_index]
+                  [Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical)]);
+        flows[i - flows_start_index][Eigen::Index(mio::isecir::InfectionTransition::InfectedCriticalToRecovered)] =
+            (1 - model_ode.parameters.get<mio::osecir::DeathsPerCritical<ScalarType>>()[(mio::AgeGroup)0]) *
+            (compartments[scale_timesteps * (i - 1)][Eigen::Index(mio::osecir::InfectionState::InfectedCritical)] -
+             compartments[scale_timesteps * i][Eigen::Index(mio::osecir::InfectionState::InfectedCritical)] +
+             flows[i - flows_start_index]
+                  [Eigen::Index(mio::isecir::InfectionTransition::InfectedSevereToInfectedCritical)]);
+    }
+}
+
+void compute_initial_flows_for_ide_from_ode(mio::osecir::Model<ScalarType>& model_ode, mio::isecir::Model& model_ide,
+                                            mio::TimeSeries<ScalarType> secihurd_ode, ScalarType t0_ide,
+                                            ScalarType dt_ode, ScalarType dt_ide)
+{
+    std::cout << "Computing initial flows. \n";
+
+    // Use t_window=t0_ide to get flows from t0 onwards.
+    get_flows_from_ode_compartments(model_ode, secihurd_ode, model_ide.m_transitions, t0_ide, t0_ide, dt_ode, dt_ide);
+
+    // Correct values in populations.
+    if (model_ide.m_populations.get_last_time() != model_ide.m_transitions.get_last_time()) {
+        model_ide.m_populations.remove_last_time_point();
+        model_ide.m_populations.add_time_point<Eigen::VectorXd>(
+            model_ide.m_transitions.get_last_time(),
+            TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count, 0));
+        model_ide.m_populations[0][Eigen::Index(InfectionState::Dead)] =
+            secihurd_ode[(Eigen::Index)secihurd_ode.get_num_time_points() -
+                         (secihurd_ode.get_last_time() - t0_ide) / dt_ode - 1]
+                        [(Eigen::Index)mio::osecir::InfectionState::Dead];
+    }
+}
+
+} // namespace isecir
+} // namespace mio
diff --git a/cpp/models/ide_secir/initialize_from_ode.h b/cpp/models/ide_secir/initialize_from_ode.h
@@ -0,0 +1,81 @@
+/* 
+* Copyright (C) 2020-2023 German Aerospace Center (DLR-SC)
+*
+* Authors: Anna Wendler
+*
+* Contact: Martin J. Kuehn <[email protected]>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+#ifndef IDE_SECIR_INITIALIZE_FROM_ODE_H
+#define IDE_SECIR_INITIALIZE_FROM_ODE_H
+
+#include "ide_secir/model_ide.h"
+#include "ide_secir/parameters.h"
+#include "ode_secir/model.h"
+#include "infection_state.h"
+#include "memilio/config.h"
+#include "memilio/utils/logging.h"
+#include "memilio/math/eigen.h"
+
+namespace mio
+{
+namespace isecir
+{
+/**
+* @brief Takes the resulting compartments of an ODE simulation and computes the respective flows and stores them in a TimeSeries. 
+*
+* With t_max and t_window it can be determined for which time window the flows will be computed. 
+* dt_reference is the time step size of the ODE simulation. By default, we compute the flows with the same time step size.
+* It is also possible to compute the corresponding flows for a bigger time step which can be given by dt_comparison. Here we assume, that
+* dt_comparison is a multiple of dt_reference. 
+*
+* @param[in] model_ode ODE model that is used.
+* @param[in] compartments TimeSeries containing compartments from ODE simulation
+* @param[out] flows TimeSeries where the computed flows will be stored. 
+* @param[in] t_max Maximal time for which the flows are computed.
+* @param[in] t_window Time window before t_max for which flows will be computed.
+* @param[in] dt_reference Reference time step size (coming from ODE simulation).
+* @param[in] dt_comparison Time step size. Default is dt_reference. 
+*/
+void get_flows_from_ode_compartments(mio::osecir::Model<ScalarType>& model_ode,
+                                     mio::TimeSeries<ScalarType> compartments, mio::TimeSeries<ScalarType>& flows,
+                                     ScalarType t_max, ScalarType t_window, ScalarType dt_reference,
+                                     ScalarType dt_comparison = 0.);
+
+/**
+* @brief Computes the inital flows that are needed for an IDE simulation given we have the compartments from an ODE 
+* simulation for an adequate time window before t0_ide. 
+*
+* Using model_ode we get simulated results for the compartments stored in secihurd_ode. Then we compute the needed iniital flows 
+* for the IDE model model_ide until time t0_ide. Here we assumen, that model_ode and model_ide are matching, i.e. that the parameters
+* of model_ideare chosen such that model_ide shopuld reduce to model_ode by choosing exponentially distributed transitions with the
+* the corresponding mean stay times etc. The time step size of ODE and IDE simulation can be chosen independently. However, we assume
+* that dt_is a multiple of dt_ide. 
+*
+* @param[in] model_ode ODE model that is used.
+* @param[in] model_ide IDE model that is used.
+* @param[in] secihurd_ode TimeSeries containing compartments from ODE simulation. 
+* @param[in] t0_ide Start time of IDE simulation until which we need to compute flows. 
+* @param[in] dt_ode Time step size of ODE simulation.
+* @param[in] dt_ide Time step size of IDE simulation. 
+*/
+void compute_initial_flows_for_ide_from_ode(mio::osecir::Model<ScalarType>& model_ode, mio::isecir::Model& model_ide,
+                                            mio::TimeSeries<ScalarType> secihurd_ode, ScalarType t0_ide,
+                                            ScalarType dt_ode, ScalarType dt_ide);
+
+} // namespace isecir
+} // namespace mio
+
+#endif //IDE_SECIR_INITIALIZE_FROM_ODE_H
diff --git a/cpp/models/ide_secir/model.cpp → cpp/models/ide_secir/model_ide.cpp b/cpp/models/ide_secir/model.cpp → cpp/models/ide_secir/model_ide.cpp
@@ -17,8 +17,9 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#include "ide_secir/model.h"
+#include "ide_secir/model_ide.h"
 #include "ide_secir/parameters.h"
+#include "ode_secir/model.h"
 #include "infection_state.h"
 #include "memilio/config.h"
 #include "memilio/utils/logging.h"

diff --git a/cpp/models/ide_secir/model.h → cpp/models/ide_secir/model_ide.h b/cpp/models/ide_secir/model.h → cpp/models/ide_secir/model_ide.h
@@ -24,6 +24,7 @@
 #include "abm/virus_variant.h"
 #include "ide_secir/parameters.h"
 #include "ide_secir/infection_state.h"
+#include "ode_secir/model.h"
 #include "memilio/config.h"
 #include "memilio/utils/time_series.h"
 
@@ -53,6 +54,7 @@ class Model
     * @param[in] total_confirmed_cases Total confirmed cases at time t0 can be set if it should be used for initialization.
     * @param[in, out] Parameterset_init Used Parameters for simulation. 
     */
+
     Model(TimeSeries<ScalarType>&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases = 0,
           const ParameterSet& Parameterset_init = ParameterSet());
 
@@ -205,8 +207,9 @@ class Model
 
         // It may be possible to run the simulation with fewer time points, but this number ensures that it is possible.
         if (m_transitions.get_num_time_points() < (Eigen::Index)std::ceil(get_global_support_max(dt) / dt)) {
-            log_error(
-                "Initialization failed. Not enough time points for transitions given before start of simulation.");
+            log_error("Initialization failed. Not enough time points for transitions given before start of simulation. "
+                      "{} time points are required.",
+                      (Eigen::Index)std::ceil(get_global_support_max(dt) / dt));
             return true;
         }
 
@@ -219,8 +222,10 @@ class Model
             }
         }
         if (m_transitions.get_last_time() != m_populations.get_last_time()) {
-            log_error("Last time point of TimeSeries for transitions does not match last time point of TimeSeries for "
-                      "compartments. Both of these time points have to agree for a sensible simulation.");
+            log_error(
+                "Last time point of TimeSeries for transitions {} does not match last time point of TimeSeries for "
+                "compartments {}. For a meaningful simulation, these two points must match.",
+                m_transitions.get_last_time(), m_populations.get_last_time());
             return true;
         }
 
@@ -304,7 +309,7 @@ class Model
     ScalarType m_N{0}; ///< Total population size of the considered region.
     ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the TransitionDistributions.
     int m_initialization_method{0}; ///< Gives the index of the method used for the initialization of the model.
-        //See also get_initialization_method_compartments() for the number code.
+        // See also get_initialization_method_compartments() for the number code.
 };
 
 } // namespace isecir

diff --git a/cpp/models/ide_secir/parameters.h b/cpp/models/ide_secir/parameters.h
@@ -40,7 +40,6 @@ namespace isecir
 /**********************************************
 * Define Parameters of the IDE-SECIHURD model *
 **********************************************/
-
 /**
  * @brief Transition distribution for each transition in #InfectionTransition.
  *
@@ -325,6 +324,7 @@ class Parameters : public ParametersBase
         return false;
     }
 
+public:
     /**
      * deserialize an object of this class.
      * @see mio::deserialize

diff --git a/cpp/models/ide_secir/parameters_io.cpp b/cpp/models/ide_secir/parameters_io.cpp
@@ -23,7 +23,7 @@
 
 #ifdef MEMILIO_HAS_JSONCPP
 
-#include "ide_secir/model.h"
+#include "ide_secir/model_ide.h"
 #include "ide_secir/infection_state.h"
 #include "memilio/math/eigen.h"
 #include "memilio/io/epi_data.h"

diff --git a/cpp/models/ide_secir/parameters_io.h b/cpp/models/ide_secir/parameters_io.h
@@ -24,7 +24,7 @@
 
 #ifdef MEMILIO_HAS_JSONCPP
 
-#include "ide_secir/model.h"
+#include "ide_secir/model_ide.h"
 #include "memilio/io/io.h"
 #include "memilio/utils/date.h"