1018 use floating point types consistently (#1307)

- Extended automatic differentiation (AD) support to the Secirvvs model.
- Templated floating-point types (FP) across key classes and functions for AD compatibility.
- Added math operations for UncertainValue to work correctly with AD types.
- Replaced most occurrences of double with FP or ScalarType in the core library.
- Added a function eval_intermediate to consolidate usage of Eigen and AD intermediates.

Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>
Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com>
Co-authored-by: jubicker <113909589+jubicker@users.noreply.github.com>
Co-authored-by: Henrik Zunker <69154294+HenrZu@users.noreply.github.com>
Co-authored-by: MaxBetz <104758467+MaxBetzDLR@users.noreply.github.com>

Author: 6 people

cpp/benchmarks/abm.cpp

@@ -71,7 +71,7 @@ mio::abm::Simulation<> make_simulation(size_t num_persons, std::initializer_list
         auto prng = mio::abm::PersonalRandomNumberGenerator(person);
         //some % of people are infected, large enough to have some infection activity without everyone dying
         auto pct_infected = 0.05;
-        if (mio::UniformDistribution<double>::get_instance()(prng, 0.0, 1.0) < pct_infected) {
+        if (mio::UniformDistribution<ScalarType>::get_instance()(prng, 0.0, 1.0) < pct_infected) {
             auto state = mio::abm::InfectionState(
                 mio::UniformIntDistribution<int>::get_instance()(prng, 1, int(mio::abm::InfectionState::Count) - 1));
             auto infection = mio::abm::Infection(prng, mio::abm::VirusVariant::Wildtype, person.get_age(),
@@ -91,7 +91,7 @@ mio::abm::Simulation<> make_simulation(size_t num_persons, std::initializer_list
         //skip homes so persons always have a place to go, simulation might break otherwise
         auto pct_require_mask = 0.2;
         if (loc.get_type() != mio::abm::LocationType::Home &&
-            mio::UniformDistribution<double>::get_instance()(model.get_rng()) < pct_require_mask) {
+            mio::UniformDistribution<ScalarType>::get_instance()(model.get_rng()) < pct_require_mask) {
             loc.set_required_mask(mio::abm::MaskType::Community);
         }
     }

cpp/benchmarks/flow_simulation_ode_secirvvs.h

@@ -1,4 +1,4 @@
-/* 
+/*
 * Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Rene Schmieding, Daniel Abele, Martin J. Kuehn
@@ -52,8 +52,9 @@ class FlowlessModel
     {
     }
 
-    void get_derivatives(Eigen::Ref<const Eigen::VectorXd> pop, Eigen::Ref<const Eigen::VectorXd> y, double t,
-                         Eigen::Ref<Eigen::VectorXd> dydt) const override
+    void get_derivatives(Eigen::Ref<const Eigen::VectorX<ScalarType>> pop,
+                         Eigen::Ref<const Eigen::VectorX<ScalarType>> y, ScalarType t,
+                         Eigen::Ref<Eigen::VectorX<ScalarType>> dydt) const override
     {
         auto& params                     = this->parameters;
         const Index<AgeGroup> age_groups = reduce_index<Index<AgeGroup>>(this->populations.size());
@@ -71,10 +72,10 @@ class FlowlessModel
                 size_t Ij = this->populations.get_flat_index({j, InfectionState::Infected});
                 size_t Rj = this->populations.get_flat_index({j, InfectionState::Recovered});
 
-                const double Nj_inv = 1.0 / (pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj]);
-                const double coeffStoE =
-                    params.template get<ContactPatterns<ScalarType>>().get_cont_freq_mat().get_matrix_at(t)(i.get(),
-                                                                                                            j.get()) *
+                const ScalarType Nj_inv = 1.0 / (pop[Sj] + pop[Ej] + pop[Ij] + pop[Rj]);
+                const ScalarType coeffStoE =
+                    params.template get<ContactPatterns<ScalarType>>().get_cont_freq_mat().get_matrix_at(
+                        SimulationTime<ScalarType>(t))(i.get(), j.get()) *
                     params.template get<TransmissionProbabilityOnContact<ScalarType>>()[i] * Nj_inv;
 
                 dydt[Si] -= y[Si] * pop[Ij] * coeffStoE;
@@ -91,8 +92,8 @@ class FlowlessModel
 template <class Model>
 void setup_model(Model& model)
 {
-    const double total_population = 10000.0;
-    const auto num_groups         = model.parameters.get_num_groups();
+    const ScalarType total_population = 10000.0;
+    const auto num_groups             = model.parameters.get_num_groups();
     for (AgeGroup i = 0; i < num_groups; i++) {
         model.populations[{i, oseir::InfectionState::Exposed}]   = 100.0 / static_cast<size_t>(num_groups);
         model.populations[{i, oseir::InfectionState::Infected}]  = 100.0 / static_cast<size_t>(num_groups);
@@ -106,8 +107,9 @@ void setup_model(Model& model)
     model.parameters.template set<mio::oseir::TimeExposed<ScalarType>>(5.2);
     model.parameters.template set<mio::oseir::TimeInfected<ScalarType>>(6);
     model.parameters.template set<mio::oseir::TransmissionProbabilityOnContact<ScalarType>>(0.04);
-    mio::ContactMatrixGroup& contact_matrix = model.parameters.template get<mio::oseir::ContactPatterns<ScalarType>>();
-    contact_matrix[0].get_baseline().setConstant(10.);
+    mio::ContactMatrixGroup<ScalarType>& contact_matrix =
+        model.parameters.template get<mio::oseir::ContactPatterns<ScalarType>>();
+    contact_matrix[0].get_baseline().setConstant(10.0);
 }
 
 } // namespace benchmark

cpp/benchmarks/flow_simulation_ode_seir.cpp

@@ -1,4 +1,4 @@
-/* 
+/*
 * Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Henrik Zunker
@@ -78,7 +78,7 @@ mio::osecirvvs::Model<ScalarType> create_model(size_t num_agegroups, const Scala
     auto& contact_matrix = contacts.get_cont_freq_mat();
     contact_matrix[0].get_baseline().setConstant(0.5);
     contact_matrix[0].get_baseline().diagonal().setConstant(5.0);
-    contact_matrix[0].add_damping(0.3, mio::SimulationTime(5.0));
+    contact_matrix[0].add_damping(0.3, mio::SimulationTime<ScalarType>(5.0));
 
     for (mio::AgeGroup i = 0; i < (mio::AgeGroup)num_agegroups; i++) {
         //times
@@ -115,9 +115,9 @@ auto create_simulation()
 {
     auto cfg = mio::benchmark::GraphConfig::initialize(config_path);
 
-    mio::osecirvvs::Model model = create_model(cfg.num_agegroups, cfg.t_max);
+    mio::osecirvvs::Model<ScalarType> model = create_model(cfg.num_agegroups, cfg.t_max);
 
-    mio::Graph<mio::SimulationNode<mio::Simulation<ScalarType, mio::osecirvvs::Model<ScalarType>>>,
+    mio::Graph<mio::SimulationNode<ScalarType, mio::Simulation<ScalarType, mio::osecirvvs::Model<ScalarType>>>,
                mio::MobilityEdge<ScalarType>>
         g;
     for (size_t county_id = 0; county_id < cfg.num_regions; county_id++) {
@@ -130,9 +130,9 @@ auto create_simulation()
         for (size_t county_idx_j = 0; county_idx_j < g.nodes().size(); ++county_idx_j) {
             if (county_idx_i == county_idx_j)
                 continue;
-            g.add_edge(
-                county_idx_i, county_idx_j,
-                Eigen::VectorXd::Constant((size_t)mio::osecirvvs::InfectionState::Count * cfg.num_agegroups, 0.01));
+            g.add_edge(county_idx_i, county_idx_j,
+                       Eigen::VectorX<ScalarType>::Constant(
+                           (size_t)mio::osecirvvs::InfectionState::Count * cfg.num_agegroups, 0.01));
         }
     }
 

cpp/benchmarks/graph_simulation.cpp

@@ -1,4 +1,4 @@
-/* 
+/*
 * Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Henrik Zunker
@@ -32,7 +32,7 @@ namespace benchmark
 /// @brief parameters for simulation benchmark
 struct GraphConfig {
     int num_agegroups, num_regions;
-    double t0, t_max, dt;
+    ScalarType t0, t_max, dt;
     /**
      * @brief creates configuration with default parameters for a secir model
      * @param num_agegroups number of agegroups
@@ -65,9 +65,9 @@ struct GraphConfig {
         auto obj              = io.expect_object("bench_graph_simulation");
         auto num_agegroups_io = obj.expect_element("num_agegroups", mio::Tag<int>{});
         auto num_regions_io   = obj.expect_element("num_regions", mio::Tag<int>{});
-        auto t_io             = obj.expect_element("t0", mio::Tag<double>{});
-        auto t_max_io         = obj.expect_element("t_max", mio::Tag<double>{});
-        auto dt_io            = obj.expect_element("dt", mio::Tag<double>{});
+        auto t_io             = obj.expect_element("t0", mio::Tag<ScalarType>{});
+        auto t_max_io         = obj.expect_element("t_max", mio::Tag<ScalarType>{});
+        auto dt_io            = obj.expect_element("dt", mio::Tag<ScalarType>{});
         return mio::apply(
             io,
             [](auto&& num_agegroups, auto&& num_regions, auto&& t0, auto&& t_max, auto&& dt) {

cpp/benchmarks/graph_simulation.h

@@ -1,4 +1,4 @@
-/* 
+/*
 * Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Rene Schmieding
@@ -36,13 +36,13 @@ void integrator_step(::benchmark::State& state)
     //auto cfg = mio::benchmark::IntegratorStepConfig::initialize();
     auto model = mio::benchmark::model::SecirAgeres(cfg.num_agegroups);
     // set deriv function and integrator
-    mio::DerivFunction<ScalarType> f = [model](Eigen::Ref<const Eigen::VectorXd> x, double s,
-                                               Eigen::Ref<Eigen::VectorXd> dxds) {
+    mio::DerivFunction<ScalarType> f = [model](Eigen::Ref<const Eigen::VectorX<ScalarType>> x, ScalarType s,
+                                               Eigen::Ref<Eigen::VectorX<ScalarType>> dxds) {
         model.eval_right_hand_side(x, x, s, dxds);
     };
     auto I = Integrator(cfg.abs_tol, cfg.rel_tol, cfg.dt_min, cfg.dt_max);
 
-    double t, dt;
+    ScalarType t, dt;
     for (auto _ : state) {
         // This code gets timed
         t  = cfg.t_init;

cpp/benchmarks/integrator_step.cpp

@@ -1,4 +1,4 @@
-/* 
+/*
 * Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Rene Schmieding
@@ -32,31 +32,31 @@ namespace benchmark
 /// @brief parameters for integrator-step benchmark
 struct IntegratorStepConfig {
     int num_agegroups;
-    double t_init, dt_init, abs_tol, rel_tol, dt_min, dt_max;
-    Eigen::VectorXd yt, ytp1;
+    ScalarType t_init, dt_init, abs_tol, rel_tol, dt_min, dt_max;
+    Eigen::VectorX<ScalarType> yt, ytp1;
     /**
-         * @brief creates configuration with default parameters for a secir model
-         * @return configuration for integrator-step benchmark
-         */
+     * @brief creates configuration with default parameters for a secir model
+     * @return configuration for integrator-step benchmark
+     */
     static IntegratorStepConfig initialize()
     {
-        const double vals[8] = {6377.873644, 35.249156, 30.029611,   182.145865,
-                                66.153059,   79.530621, 3069.383604, 159.634440};
+        const ScalarType vals[8] = {6377.873644, 35.249156, 30.029611,   182.145865,
+                                    66.153059,   79.530621, 3069.383604, 159.634440};
         return IntegratorStepConfig{1,
                                     50,
                                     1,
                                     1e-10,
                                     1e-5,
-                                    std::numeric_limits<double>::min(),
-                                    std::numeric_limits<double>::max(),
-                                    Eigen::Matrix<double, 8, 1>(vals),
-                                    Eigen::VectorXd::Zero(8)};
+                                    std::numeric_limits<ScalarType>::min(),
+                                    std::numeric_limits<ScalarType>::max(),
+                                    Eigen::Matrix<ScalarType, 8, 1>(vals),
+                                    Eigen::VectorX<ScalarType>::Zero(8)};
     }
     /**
-         * @brief reads configuration from json file
-         * @param path the path of the configfile
-         * @return configuration for integrator-step benchmark
-         */
+     * @brief reads configuration from json file
+     * @param path the path of the configfile
+     * @return configuration for integrator-step benchmark
+     */
     static IntegratorStepConfig initialize(std::string path)
     {
         auto result = mio::read_json(path, mio::Tag<IntegratorStepConfig>{});
@@ -72,13 +72,13 @@ struct IntegratorStepConfig {
     {
         auto obj              = io.expect_object("integrator_step");
         auto num_agegroups_io = obj.expect_element("num_agegroups", mio::Tag<int>{});
-        auto t_init_io        = obj.expect_element("t_init", mio::Tag<double>{});
-        auto dt_init_io       = obj.expect_element("dt_init", mio::Tag<double>{});
-        auto abs_tol_io       = obj.expect_element("abs_tol", mio::Tag<double>{});
-        auto rel_tol_io       = obj.expect_element("rel_tol", mio::Tag<double>{});
-        auto dt_min_io        = obj.expect_element("dt_min", mio::Tag<double>{});
-        auto dt_max_io        = obj.expect_element("dt_max", mio::Tag<double>{});
-        auto yt_io            = obj.expect_list("yt", mio::Tag<double>{});
+        auto t_init_io        = obj.expect_element("t_init", mio::Tag<ScalarType>{});
+        auto dt_init_io       = obj.expect_element("dt_init", mio::Tag<ScalarType>{});
+        auto abs_tol_io       = obj.expect_element("abs_tol", mio::Tag<ScalarType>{});
+        auto rel_tol_io       = obj.expect_element("rel_tol", mio::Tag<ScalarType>{});
+        auto dt_min_io        = obj.expect_element("dt_min", mio::Tag<ScalarType>{});
+        auto dt_max_io        = obj.expect_element("dt_max", mio::Tag<ScalarType>{});
+        auto yt_io            = obj.expect_list("yt", mio::Tag<ScalarType>{});
         return mio::apply(
             io,
             [](auto&& num_agegroups, auto&& t_init, auto&& dt_init, auto&& abs_tol, auto&& rel_tol, auto&& dt_min,
@@ -90,8 +90,8 @@ struct IntegratorStepConfig {
                                          rel_tol,
                                          dt_min,
                                          dt_max,
-                                         Eigen::VectorXd::Zero(yt.size()),
-                                         Eigen::VectorXd::Zero(yt.size())};
+                                         Eigen::VectorX<ScalarType>::Zero(yt.size()),
+                                         Eigen::VectorX<ScalarType>::Zero(yt.size())};
                 for (size_t i = 0; i < yt.size(); i++) {
                     cfg.yt[i] = yt[i];
                 }