1451 interpolate_simulation_result does not yield consistent timepoints

cpp/memilio/data/analyze_result.h

@@ -20,40 +20,45 @@
 #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 <cassert>
 #include <vector>
 
 namespace mio
 {
 
 /**
- * @brief interpolate time series with evenly spaced, integer time points that represent whole days.
- * time points [t0, t1, t2, ..., tmax] interpolated as days [ceil(t0), floor(t0) + 1,...,floor(tmax)].
- * tolerances in the first and last time point (t0 and t_max) are accounted for.
- * values at new time points are linearly interpolated from their immediate neighbors from the old time points.
+ * @brief Interpolate a given time series with evenly spaced, integer time points that represent whole days.
+ * We choose the time points for the interpolated time series using the first and last time points, t0 and tmax, as
+ * [ceil(t0 - abs_tol), floor(t0) + 1, ..., floor(tmax + abs_tol)].
+ * The tolerances in the first and last time point account for inaccuracies from the integration scheme or
+ * floating point arithmetic. Avoid using large(r) tolerances, as this function can not extrapolate results.
+ * The values at new time points are linearly interpolated from their immediate neighbors within the given time series.
  * @see interpolate_simulation_result
- * @param simulation_result time series to interpolate
- * @param abs_tol  absolute tolerance given for doubles t0 and tmax to account for small deviations from whole days.
- * @return interpolated time series
+ * @param simulation_result A time series to interpolate.
+ * @param abs_tol Optional parameter to set the absolute tolerance used to account for small deviations from whole days.
+ *   Must be less then 1. The default tolerance is chosen to minimize the chances of "loosing" days due to rounding.
+ * @return An interpolated time series with integer valued times.
  */
 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.
- * values at new time points are linearly interpolated from their immediate neighbors from the old time points.
- * @param simulation_result time series to interpolate
- * @param interpolations_times std::vector of time points at which simulation results are interpolated.
- * @return interpolated time series at given interpolation points
+ * @brief Interpolate a time series at the given time points.
+ * New time points must be monotonic increasing, and lie in between time points of the given time series.
+ * The values at new time points are linearly interpolated from their immediate neighbors within the given time series.
+ * @param simulation_result The time series to interpolate.
+ * @param interpolation_times An std::vector of time points at which simulation results are interpolated.
+ * @return The interpolated time series at given interpolation points.
  */
 template <typename FP>
 TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result,
-
                                              const std::vector<FP>& interpolation_times);
 
 /**
@@ -180,14 +185,15 @@ FP result_distance_2norm(const std::vector<mio::TimeSeries<FP>>& result1,
 template <typename FP>
 TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result, const FP abs_tol)
 {
+    assert(abs_tol < 1);
     using std::ceil;
     using std::floor;
     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);
-    // 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);
+    // add an additional day, if the first time point is within tolerance of floor(t0)
+    const auto day0 = ceil(t0 - abs_tol);
+    // add an additional day, if the last time point is within tolerance of ceil(tmax)
+    const auto day_max = floor(t_max + abs_tol);
 
     // 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,62 @@
+/* 
+* 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);
+        },
+        "Interpolate a given time series with evenly spaced, integer time points that represent whole days.",
+        py::arg("ts"));
+    m.def("interpolate_simulation_result",
+          static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
+              &mio::interpolate_simulation_result),
+          "Interpolate a given time series with evenly spaced, integer time points that represent whole days.",
+          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),
+          "Interpolate a time series at the given time points.", py::arg("ts"), py::arg("interpolation_times"));
+
+    m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>,
+          "Interpolates results of all runs with evenly spaced, integer time points that represent whole days.");
+}
+
+} // 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)