SciCompMod · HenrZu · Jul 18, 2025 · Jun 10, 2025 · Jun 10, 2025 · Jun 10, 2025
diff --git a/cpp/benchmarks/flow_simulation_ode_secirvvs.cpp b/cpp/benchmarks/flow_simulation_ode_secirvvs.cpp
@@ -38,7 +38,7 @@ void flowless_sim(::benchmark::State& state)
     Model model(cfg.num_agegroups);
     mio::benchmark::setup_model(model);
     // create simulation
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
         std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
             cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
     // run benchmark
@@ -63,7 +63,7 @@ void flow_sim_comp_only(::benchmark::State& state)
     Model model(cfg.num_agegroups);
     mio::benchmark::setup_model(model);
     // create simulation
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
         std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
             cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
     // run benchmark
@@ -88,7 +88,7 @@ void flow_sim(::benchmark::State& state)
     Model model(cfg.num_agegroups);
     mio::benchmark::setup_model(model);
     // create simulation
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
         std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
             cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
     // run benchmark

diff --git a/cpp/benchmarks/flow_simulation_ode_seir.cpp b/cpp/benchmarks/flow_simulation_ode_seir.cpp
@@ -125,7 +125,7 @@ void flowless_sim(::benchmark::State& state)
     Model model(cfg.num_agegroups);
     mio::benchmark::setup_model(model);
     // create simulation
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
         std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
             cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
     mio::TimeSeries<ScalarType> results(static_cast<size_t>(Model::Compartments::Count));
@@ -148,7 +148,7 @@ void flow_sim_comp_only(::benchmark::State& state)
     Model model(cfg.num_agegroups);
     mio::benchmark::setup_model(model);
     // create simulation
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
         std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
             cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
     mio::TimeSeries<ScalarType> results(static_cast<size_t>(Model::Compartments::Count));
@@ -171,7 +171,7 @@ void flow_sim(::benchmark::State& state)
     Model model(cfg.num_agegroups);
     mio::benchmark::setup_model(model);
     // create simulation
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
         std::make_shared<mio::ControlledStepperWrapper<ScalarType, boost::numeric::odeint::runge_kutta_cash_karp54>>(
             cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
     mio::TimeSeries<ScalarType> results(static_cast<size_t>(Model::Compartments::Count));

diff --git a/cpp/benchmarks/simulation.cpp b/cpp/benchmarks/simulation.cpp
@@ -35,7 +35,7 @@ void simulation(::benchmark::State& state)
 
     for (auto _ : state) {
         // This code gets timed
-        std::shared_ptr<mio::IntegratorCore<ScalarType>> I =
+        std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> I =
             std::make_shared<Integrator>(cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
         simulate(cfg.t0, cfg.t_max, cfg.dt, model, I);
     }

diff --git a/cpp/examples/adapt_rk_test.cpp b/cpp/examples/adapt_rk_test.cpp
@@ -17,13 +17,11 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#include "memilio/math/euler.h"
 #include "memilio/math/adapt_rk.h"
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <vector>
-#include <iostream>
 #include <cmath>
 
 void init_vectors(std::vector<Eigen::VectorXd>& y, std::vector<Eigen::VectorXd>& sol, size_t n)
@@ -40,7 +38,7 @@ void integration_test(std::vector<Eigen::VectorXd>& y, std::vector<Eigen::Vector
         dydt[0] = std::cos(t);
     };
 
-    mio::RKIntegratorCore<> rkf45;
+    mio::RKIntegratorCore<double> rkf45;
     rkf45.set_abs_tolerance(1e-7);
     rkf45.set_rel_tolerance(1e-7);
     rkf45.set_dt_min(1e-3);

diff --git a/cpp/examples/euler_test.cpp b/cpp/examples/euler_test.cpp
@@ -22,7 +22,6 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <vector>
-#include <iostream>
 #include <cmath>
 
 void init_vectors(std::vector<Eigen::VectorXd>& y, std::vector<Eigen::VectorXd>& sol, size_t n)
@@ -46,7 +45,7 @@ void integration_test(std::vector<Eigen::VectorXd>& y, std::vector<Eigen::Vector
     for (size_t i = 0; i < n - 1; i++) {
         sol[i + 1][0] = std::sin((i + 1) * dt);
 
-        mio::EulerIntegratorCore().step(f, y[i], t, dt, y[i + 1]);
+        mio::EulerIntegratorCore<double>().step(f, y[i], t, dt, y[i + 1]);
 
         printf("\n %.8f\t %.8f", y[i + 1][0], sol[i + 1][0]);
         // printf("\n approx: %.4e, sol: %.4e, error %.4e", y[i+1][0], sol[i+1][0], err);

diff --git a/cpp/examples/ode_secir_contact_changes.cpp b/cpp/examples/ode_secir_contact_changes.cpp
@@ -18,7 +18,7 @@
 * limitations under the License.
 */
 #include "ode_secir/model.h"
-#include "memilio/compartments/simulation.h"
+#include "memilio/compartments/flow_simulation.h"
 #include "memilio/utils/logging.h"
 #include "memilio/math/euler.h"
 

diff --git a/cpp/examples/ode_sir.cpp b/cpp/examples/ode_sir.cpp
@@ -56,7 +56,7 @@ int main()
     contact_matrix[0].add_damping(0.6, mio::SimulationTime(12.5));
     model.check_constraints();
 
-    std::shared_ptr<mio::IntegratorCore<ScalarType>> integrator =
+    std::shared_ptr<mio::OdeIntegratorCore<ScalarType>> integrator =
         std::make_shared<mio::EulerIntegratorCore<ScalarType>>();
     auto sir = simulate(t0, tmax, dt, model, integrator);
 

diff --git a/cpp/examples/sde_seirvv.cpp b/cpp/examples/sde_seirvv.cpp
@@ -17,10 +17,10 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
+#include "memilio/compartments/stochastic_simulation.h"
 #include "memilio/utils/logging.h"
 #include "memilio/utils/uncertain_value.h"
 #include "sde_seirvv/model.h"
-#include "sde_seirvv/simulation.h"
 
 #include <vector>
 
@@ -73,13 +73,13 @@ int main()
     model.check_constraints();
 
     // Simulate the model up until tmid, with only the first variant.
-    auto sseirv = mio::sseirvv::simulate(t0, tmid, dt, model);
+    auto sseirv = mio::simulate_stochastic(t0, tmid, dt, model);
     // Set the model population to the simulation result, so it is used as initial value for the second simulation.
     model.populations.array() = sseirv.get_last_value().cast<mio::UncertainValue<ScalarType>>();
     // The second variant enters with 100 individuals. This increases the model population to total_population + 100.
     model.populations[{mio::sseirvv::InfectionState::InfectedV2}] = 100;
     // Simulate the model from tmid to tmax, now with both variants.
-    auto sseirv2 = mio::sseirvv::simulate(tmid, tmax, dt, model);
+    auto sseirv2 = mio::simulate_stochastic(tmid, tmax, dt, model);
     sseirv.print_table({"Susceptible", "ExposedV1", "InfectedV1", "RecoveredV1", "ExposedV2", "InfectedV2",
                         "RecoveredV2", "ExposedV1V2", "InfectedV1V2", "RecoveredV1V2"});
     sseirv2.print_table({"Susceptible", "ExposedV1", "InfectedV1", "RecoveredV1", "ExposedV2", "InfectedV2",

diff --git a/cpp/examples/sde_sir.cpp b/cpp/examples/sde_sir.cpp
@@ -18,9 +18,8 @@
 * limitations under the License.
 */
 
-#include "memilio/utils/logging.h"
+#include "memilio/compartments/stochastic_simulation.h"
 #include "sde_sir/model.h"
-#include "sde_sir/simulation.h"
 
 int main()
 {
@@ -49,7 +48,7 @@ int main()
 
     model.check_constraints();
 
-    auto sir = mio::ssir::simulate(t0, tmax, dt, model);
+    auto sir = mio::simulate_stochastic(t0, tmax, dt, model);
 
     sir.print_table();
 }
diff --git a/cpp/examples/sde_sirs.cpp b/cpp/examples/sde_sirs.cpp
@@ -17,10 +17,10 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
+#include "memilio/compartments/stochastic_simulation.h"
 #include "memilio/data/analyze_result.h"
 #include "memilio/utils/logging.h"
 #include "sde_sirs/model.h"
-#include "sde_sirs/simulation.h"
 
 int main()
 {
@@ -50,7 +50,7 @@ int main()
 
     model.check_constraints();
 
-    auto ssirs = mio::ssirs::simulate(t0, tmax, dt, model);
+    auto ssirs = mio::simulate_stochastic(t0, tmax, dt, model);
 
     // interpolate results
     auto interpolated_results = mio::interpolate_simulation_result(ssirs);

diff --git a/cpp/memilio/CMakeLists.txt b/cpp/memilio/CMakeLists.txt
@@ -25,8 +25,12 @@ add_library(memilio
     geography/holiday_data.ipp
     compartments/compartmental_model.h
     compartments/flow_model.h
-    compartments/simulation.h
+    compartments/flow_simulation_base.h
     compartments/flow_simulation.h
+    compartments/simulation_base.h
+    compartments/simulation.h
+    compartments/stochastic_simulation.h
+    compartments/stochastic_model.h
     compartments/parameter_studies.h
     io/default_serialize.h
     io/default_serialize.cpp
@@ -64,6 +68,7 @@ add_library(memilio
     math/interpolation.cpp
     math/time_series_functor.h
     math/time_series_functor.cpp
+    math/math_utils.h
     mobility/metapopulation_mobility_instant.h
     mobility/metapopulation_mobility_instant.cpp
     mobility/metapopulation_mobility_stochastic.h

diff --git a/cpp/memilio/compartments/README.md b/cpp/memilio/compartments/README.md
@@ -1,11 +1,31 @@
 This directory contains the framework used to build SIR-type compartment models.
 
 Classes:
-- CompartmentalModel: Template base class for compartmental models. Specialize the class template using a parameter set (e.g. using the [ParameterSet class](../utils/parameter_set.h)) and populations (e.g. using the [Populations class](../epidemiology/populations.h)). The population is divided into compartments (and optionally other subcategories, e.g. age groups). Derive from the class to define the change in populations.
-- FlowModel: A CompartmentalModel defined by the flows between populations. Requires the additional template parameter Flows, which is a [TypeList](../utils/type_list.h) of [Flow](../utils/flow.h). Instead of defining `get_derivatives`, the function `get_flows` has to be implemented for a FlowModel. 
-- Simulation: Template class that runs the simulation using a specified compartment model. Can be derived from to implement behavior that cannot be modeled inside the usual compartment flow structure.
-- FlowSimulation: Template class derived from Simulation, that additionally computes the flows between populations. It requires a FlowModel. While the compartment and flow results can be changed outside of the simulation, it is required that both still have the same number of time points to further advance this simulation.
-- FeedbackSimulation: Template class derived from Simulation (or FlowSimulation) that extends the standard simulation functionality with a feedback mechanism. FeedbackSimulation holds its own parameters (using a specialized ParameterSet) and a perceived risk time series. It overrides the simulation advance routine to compute and apply contact adjustments at regular intervals. The contact adjustments are computed by first deriving the perceived risk from the historical ICU occupancy—normalized by the nominal ICU capacity—and then transforming that risk via a softplus function into an effective damping factor to reduce contacts as described in Dönges et al. ([doi.org/10.3389/fphy.2022.842180](https://doi.org/10.3389/fphy.2022.842180)).
+- CompartmentalModel: Template base class for compartmental models. Specialize the class template using a parameter set
+  (e.g. using the [ParameterSet class](../utils/parameter_set.h)) and populations (e.g. using the
+  [Populations class](../epidemiology/populations.h)). The population is divided into compartments (and optionally other
+  subcategories, e.g. age groups). Derive from the class to define the change in populations.
+- FlowModel: A CompartmentalModel defined by the flows between populations. Requires the additional template parameter
+  Flows, which is a [TypeList](../utils/type_list.h) of [Flow](../utils/flow.h). Instead of defining `get_derivatives`,
+  the function `get_flows` has to be implemented for a FlowModel. 
+- StochasticModel: A CompartmentalModel (or FlowModel) defined through a stochastic differential equation (SDE). In
+  addition, a function `get_noise` needs to implemented that returns random changes to each compartment.
+- Simulation: Template class that runs the simulation using a specified compartment model. Can be derived from to
+  implement behavior that cannot be modeled inside the usual compartment flow structure.
+- FlowSimulation: Template class derived from Simulation, that additionally computes the flows between populations. It
+  requires a FlowModel. While the compartment and flow results can be changed outside of the simulation, it is required
+  that both still have the same number of time points to further advance this simulation.
+- FeedbackSimulation: Template class derived from Simulation (or FlowSimulation) that extends the standard simulation
+  functionality with a feedback mechanism. FeedbackSimulation holds its own parameters (using a specialized
+  ParameterSet) and a perceived risk time series. It overrides the simulation advance routine to compute and apply
+  contact adjustments at regular intervals. The contact adjustments are computed by first deriving the perceived risk
+  from the historical ICU occupancy—normalized by the nominal ICU capacity—and then transforming that risk via a
+  softplus function into an effective damping factor to reduce contacts as described in Dönges et al.
+  ([doi.org/10.3389/fphy.2022.842180](https://doi.org/10.3389/fphy.2022.842180)).
+- StochasticSimulation: A simulation specifically for running a StochasticModel. Cannot simulate flows, because the
+  integrator uses mitigations for negative values caused by random noise that (in general) only work on compartments.
+- SimulationBase and FlowSimulationBase: Base classes for the simulations that define most members and setters/getters,
+  so that only a public advance method needs to be defined. 
 
 Note that currently, only the last value of the simulation result (i.e. the TimeSeries obtained from `get_result()`) is used to advance the simulation. Thus, any changes to the population of the simulation's model will have no effect on the simulation result.
 

diff --git a/cpp/memilio/compartments/compartmental_model.h b/cpp/memilio/compartments/compartmental_model.h
@@ -17,8 +17,8 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#ifndef MIO_COMPARTMENTALMODEL_H
-#define MIO_COMPARTMENTALMODEL_H
+#ifndef MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H
+#define MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H
 
 #include "memilio/config.h"
 #include "memilio/math/eigen.h"
@@ -95,8 +95,9 @@ struct CompartmentalModel {
     virtual ~CompartmentalModel()                            = default;
 
     // REMARK: Not pure virtual for easier java/python bindings.
-    virtual void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> /*pop*/, Eigen::Ref<const Eigen::VectorX<FP>> /*y*/,
-                                 FP /*t*/, Eigen::Ref<Eigen::VectorX<FP>> /*dydt*/) const {};
+    virtual void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> /*pop*/,
+                                 Eigen::Ref<const Eigen::VectorX<FP>> /*y*/, FP /*t*/,
+                                 Eigen::Ref<Eigen::VectorX<FP>> /*dydt*/) const {};
 
     /**
      * @brief This function evaluates the right-hand-side f of the ODE dydt = f(y, t).
@@ -214,4 +215,4 @@ using is_compartment_model =
 
 } // namespace mio
 
-#endif // MIO_COMPARTMENTALMODEL_H
+#endif // MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H
diff --git a/cpp/memilio/compartments/feedback_simulation.h b/cpp/memilio/compartments/feedback_simulation.h
@@ -17,8 +17,8 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#ifndef FEEDBACK_SIMULATION_H
-#define FEEDBACK_SIMULATION_H
+#ifndef MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
+#define MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
 
 #include "memilio/compartments/simulation.h"
 #include "memilio/utils/time_series.h"
@@ -375,4 +375,4 @@ class FeedbackSimulation
 
 } // namespace mio
 
-#endif // FEEDBACK_SIMULATION_H
+#endif // MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
diff --git a/cpp/memilio/compartments/flow_model.h b/cpp/memilio/compartments/flow_model.h
@@ -18,8 +18,8 @@
 * limitations under the License.
 */
 
-#ifndef MIO_FLOW_MODEL_H_
-#define MIO_FLOW_MODEL_H_
+#ifndef MIO_COMPARTMENTS_FLOW_MODEL_H
+#define MIO_COMPARTMENTS_FLOW_MODEL_H
 
 #include "memilio/compartments/compartmental_model.h"
 #include "memilio/utils/index_range.h"
@@ -273,4 +273,4 @@ using is_flow_model = std::integral_constant<bool, (is_expression_valid<get_deri
 
 } // namespace mio
 
-#endif // MIO_FLOW_MODEL_H_
+#endif // MIO_COMPARTMENTS_FLOW_MODEL_H