Skip to content

1451 interpolate_simulation_result does not yield consistent timepoints #1482

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
reneSchm merged 5 commits into from
Jan 27, 2026
Merged

1451 interpolate_simulation_result does not yield consistent timepoints #1482

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
4ee4d84
update tolerances, removing magic numbers
reneSchm Jan 20, 2026
949af55
make conditionals on adding days easier to read
reneSchm Jan 20, 2026
c0d7d3b
update interpolation bindings, using a central method
reneSchm Jan 20, 2026
a46386e
fix new file attribution
reneSchm Jan 20, 2026
e34590e
fix documentation according to review.
reneSchm Jan 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 25 additions & 19 deletions cpp/memilio/data/analyze_result.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

/**
Expand Down Expand Up @@ -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);
Expand Down
11 changes: 6 additions & 5 deletions cpp/memilio/math/integrator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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.
Expand All @@ -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,
Expand All @@ -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 {
Expand Down
62 changes: 62 additions & 0 deletions pycode/memilio-simulation/memilio/simulation/bindings/data/analyze_result.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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")
Expand Down
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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>);
Expand Down
11 changes: 2 additions & 9 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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)
Expand Down
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/oseir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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)
Expand Down
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/oseirdb.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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)
Expand Down
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/osir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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)
Expand Down
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/ssir.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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)
Expand Down
13 changes: 2 additions & 11 deletions pycode/memilio-simulation/memilio/simulation/bindings/models/ssirs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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)
Expand Down