1451 interpolate_simulation_result does not yield consistent timepoints

cpp/memilio/data/analyze_result.h

@@ -20,12 +20,12 @@
 #ifndef MEMILIO_DATA_ANALYZE_RESULT_H
 #define MEMILIO_DATA_ANALYZE_RESULT_H
 
+#include "memilio/config.h"
 #include "memilio/utils/time_series.h"
 #include "memilio/mobility/metapopulation_mobility_instant.h"
 #include "memilio/math/interpolation.h"
 #include "memilio/io/io.h"
 
-#include <functional>
 #include <vector>
 
 namespace mio
@@ -42,7 +42,8 @@ namespace mio
  * @return interpolated time series
  */
 template <typename FP>
-TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result, const FP abs_tol = 1e-14);
+TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result,
+                                             const FP abs_tol = FP{100.} * Limits<FP>::zero_tolerance());
 
 /**
  * @brief interpolate time series with freely chosen time points that lie in between the time points of the given time series up to a given tolerance.
@@ -185,9 +186,9 @@ TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_re
     const auto t0    = simulation_result.get_time(0);
     const auto t_max = simulation_result.get_last_time();
     // add another day if the first time point is equal to day_0 up to absolute tolerance tol
-    const auto day0 = (t0 - abs_tol < ceil(t0) - 1) ? floor(t0) : ceil(t0);
+    const auto day0 = (t0 - floor(t0) < abs_tol) ? floor(t0) : ceil(t0);
     // add another day if the last time point is equal to day_max up to absolute tolerance tol
-    const auto day_max = (t_max + abs_tol > floor(t_max) + 1) ? ceil(t_max) : floor(t_max);
+    const auto day_max = (ceil(t_max) - t_max < abs_tol) ? ceil(t_max) : floor(t_max);
 
     // create interpolation_times vector with all days between day0 and day_max
     std::vector<FP> tps(static_cast<int>(day_max) - static_cast<int>(day0) + 1);

cpp/memilio/math/integrator.h

@@ -188,11 +188,12 @@ class SystemIntegrator
         using std::fabs;
         using std::max;
         using std::min;
-        const FP t0 = results.get_last_time();
+        const FP t0  = results.get_last_time();
+        const FP eps = Limits<FP>::zero_tolerance();
         assert(tmax > t0);
         assert(dt > 0);
 
-        const size_t num_steps =
+        const auto num_steps =
             static_cast<size_t>(ceil((tmax - t0) / dt)); // estimated number of time steps (if equidistant)
 
         results.reserve(results.get_num_time_points() + num_steps);
@@ -204,7 +205,7 @@ class SystemIntegrator
         FP dt_min_restore = m_core->get_dt_min(); // used to restore dt_min, if it was decreased to reach tmax
         FP t              = t0;
 
-        for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > 1e-10; ++i) {
+        for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > eps; ++i) {
             // We don't make time steps too small as the error estimator of an adaptive integrator
             //may not be able to handle it. this is very conservative and maybe unnecessary,
             //but also unlikely to happen. may need to be reevaluated.
@@ -225,7 +226,7 @@ class SystemIntegrator
             results.get_last_time() = t;
 
             // if dt has been changed by step, register the current m_core as adaptive.
-            m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, Limits<FP>::zero_tolerance());
+            m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, eps);
         }
         m_core->get_dt_min() = dt_min_restore; // restore dt_min
         // if dt was decreased to reach tmax in the last time iteration,
@@ -236,7 +237,7 @@ class SystemIntegrator
             if (!step_okay) {
                 log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
             }
-            else if (fabs((tmax - t) / (tmax - t0)) > 1e-14) {
+            else if (fabs((tmax - t) / (tmax - t0)) > eps) {
                 log_warning("Last time step too small. Could not reach tmax exactly.");
             }
             else {

pycode/memilio-simulation/memilio/simulation/bindings/data/analyze_result.h

@@ -0,0 +1,60 @@
+/* 
+* Copyright (C) 2020-2026 MEmilio
+*
+* Authors: Rene Schmieding
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* 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 PYMIO_DATA_ANALYZE_RESULT_H
+#define PYMIO_DATA_ANALYZE_RESULT_H
+
+#include "memilio/data/analyze_result.h"
+#include "pybind_util.h"
+
+#include "pybind11/pybind11.h"
+
+namespace py = pybind11;
+
+namespace pymio
+{
+
+/**
+ * @brief Bind functions interpolate_simulation_result and interpolate_ensemble_results.
+ */
+void bind_interpolate_result_methods(py::module_& m)
+{
+
+    m.def(
+        "interpolate_simulation_result",
+        [](const mio::TimeSeries<double>& ts) {
+            return mio::interpolate_simulation_result(ts);
+        },
+        py::arg("ts"));
+    m.def("interpolate_simulation_result",
+          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
+              &mio::interpolate_simulation_result),
+          py::arg("ts"), py::arg("abs_tol"));
+
+    m.def("interpolate_simulation_result",
+          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
+              &mio::interpolate_simulation_result),
+          py::arg("ts"), py::arg("interpolation_times"));
+
+    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+}
+
+} // namespace pymio
+
+#endif // PYMIO_DATA_ANALYZE_RESULT_H

pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp

@@ -25,6 +25,7 @@
 #include "epidemiology/populations.h"
 #include "utils/parameter_set.h"
 #include "utils/index.h"
+#include "data/analyze_result.h"
 
 // Includes from MEmilio
 #include "ode_mseirs4/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::omseirs4::InfectionState>()
 PYBIND11_MODULE(_simulation_omseirs4, m)
 {
     // interpolation helpers
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     // InfectionState enum
     pymio::iterable_enum<mio::omseirs4::InfectionState>(m, "InfectionState")

pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp

@@ -33,6 +33,7 @@
 #include "mobility/metapopulation_mobility_instant.h"
 #include "io/mobility_io.h"
 #include "io/result_io.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "ode_secir/model.h"
@@ -87,17 +88,7 @@ PYBIND11_MAKE_OPAQUE(std::vector<MobilityGraph>);
 PYBIND11_MODULE(_simulation_osecir, m)
 {
     // https://github.com/pybind/pybind11/issues/1153
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     m.def("ensemble_mean", &mio::ensemble_mean<double>);
     m.def("ensemble_percentile", &mio::ensemble_percentile<double>);

pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp

@@ -31,6 +31,7 @@
 #include "epidemiology/populations.h"
 #include "io/mobility_io.h"
 #include "io/result_io.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "ode_secirvvs/model.h"
@@ -79,15 +80,7 @@ PYBIND11_MAKE_OPAQUE(std::vector<MobilityGraph>);
 
 PYBIND11_MODULE(_simulation_osecirvvs, m)
 {
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
+    pymio::bind_interpolate_result_methods(m);
 
     pymio::iterable_enum<mio::osecirvvs::InfectionState>(m, "InfectionState")
         .value("SusceptibleNaive", mio::osecirvvs::InfectionState::SusceptibleNaive)

pycode/memilio-simulation/memilio/simulation/bindings/models/oseir.cpp

@@ -28,6 +28,7 @@
 #include "compartments/compartmental_model.h"
 #include "epidemiology/age_group.h"
 #include "epidemiology/populations.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "ode_seir/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::oseir::InfectionState>()
 
 PYBIND11_MODULE(_simulation_oseir, m)
 {
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     pymio::iterable_enum<mio::oseir::InfectionState>(m, "InfectionState")
         .value("Susceptible", mio::oseir::InfectionState::Susceptible)

pycode/memilio-simulation/memilio/simulation/bindings/models/oseirdb.cpp

@@ -28,6 +28,7 @@
 #include "compartments/compartmental_model.h"
 #include "epidemiology/age_group.h"
 #include "epidemiology/populations.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "ode_seirdb/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::oseirdb::InfectionState>()
 
 PYBIND11_MODULE(_simulation_oseirdb, m)
 {
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     pymio::iterable_enum<mio::oseirdb::InfectionState>(m, "InfectionState")
         .value("Susceptible", mio::oseirdb::InfectionState::Susceptible)

pycode/memilio-simulation/memilio/simulation/bindings/models/osir.cpp

@@ -28,6 +28,7 @@
 #include "compartments/compartmental_model.h"
 #include "epidemiology/age_group.h"
 #include "epidemiology/populations.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "ode_sir/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::osir::InfectionState>()
 
 PYBIND11_MODULE(_simulation_osir, m)
 {
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     pymio::iterable_enum<mio::osir::InfectionState>(m, "InfectionState")
         .value("Susceptible", mio::osir::InfectionState::Susceptible)

pycode/memilio-simulation/memilio/simulation/bindings/models/ssir.cpp

@@ -25,6 +25,7 @@
 #include "utils/custom_index_array.h"
 #include "utils/parameter_set.h"
 #include "epidemiology/populations.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "sde_sir/model.h"
@@ -52,17 +53,7 @@ inline std::string pretty_name<mio::ssir::InfectionState>()
 
 PYBIND11_MODULE(_simulation_ssir, m)
 {
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     pymio::iterable_enum<mio::ssir::InfectionState>(m, "InfectionState")
         .value("Susceptible", mio::ssir::InfectionState::Susceptible)

pycode/memilio-simulation/memilio/simulation/bindings/models/ssirs.cpp

@@ -25,6 +25,7 @@
 #include "utils/custom_index_array.h"
 #include "utils/parameter_set.h"
 #include "epidemiology/populations.h"
+#include "data/analyze_result.h"
 
 //Includes from MEmilio
 #include "sde_sirs/model.h"
@@ -52,17 +53,7 @@ inline std::string pretty_name<mio::ssirs::InfectionState>()
 
 PYBIND11_MODULE(_simulation_ssirs, m)
 {
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("abs_tol") = 1e-14);
-
-    m.def("interpolate_simulation_result",
-          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
-              &mio::interpolate_simulation_result),
-          py::arg("ts"), py::arg("interpolation_times"));
-
-    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
+    pymio::bind_interpolate_result_methods(m);
 
     pymio::iterable_enum<mio::ssirs::InfectionState>(m, "InfectionState")
         .value("Susceptible", mio::ssirs::InfectionState::Susceptible)