Skip to content

Commit 564fbf1

Browse files
reneSchmHenrZu
reneSchm
and
HenrZu
authored
1451 interpolate_simulation_result does not yield consistent timepoints (#1482)
- Replace magic numbers used as tolerances in SystemIntegrator and interpolate_simulation_result by Limits::zero_tolerance. - Relax default tolerance in interpolate_simulation_result. - Fix bindings of interpolate_simulation_result, replacing duplicated code by a new bind_interpolate_result_methods function. Co-authored-by: Henrik Zunker <[email protected]>
1 parent 54aa7db commit 564fbf1

11 files changed

Lines changed: 109 additions & 110 deletions

File tree

‎cpp/memilio/data/analyze_result.h‎

Lines changed: 25 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -20,40 +20,45 @@
2020
#ifndef MEMILIO_DATA_ANALYZE_RESULT_H
2121
#define MEMILIO_DATA_ANALYZE_RESULT_H
2222

23+
#include "memilio/config.h"
2324
#include "memilio/utils/time_series.h"
2425
#include "memilio/mobility/metapopulation_mobility_instant.h"
2526
#include "memilio/math/interpolation.h"
2627
#include "memilio/io/io.h"
2728

28-
#include <functional>
29+
#include <cassert>
2930
#include <vector>
3031

3132
namespace mio
3233
{
3334

3435
/**
35-
* @brief interpolate time series with evenly spaced, integer time points that represent whole days.
36-
* time points [t0, t1, t2, ..., tmax] interpolated as days [ceil(t0), floor(t0) + 1,...,floor(tmax)].
37-
* tolerances in the first and last time point (t0 and t_max) are accounted for.
38-
* values at new time points are linearly interpolated from their immediate neighbors from the old time points.
36+
* @brief Interpolate a given time series with evenly spaced, integer time points that represent whole days.
37+
* We choose the time points for the interpolated time series using the first and last time points, t0 and tmax, as
38+
* [ceil(t0 - abs_tol), floor(t0) + 1, ..., floor(tmax + abs_tol)].
39+
* The tolerances in the first and last time point account for inaccuracies from the integration scheme or
40+
* floating point arithmetic. Avoid using large(r) tolerances, as this function can not extrapolate results.
41+
* The values at new time points are linearly interpolated from their immediate neighbors within the given time series.
3942
* @see interpolate_simulation_result
40-
* @param simulation_result time series to interpolate
41-
* @param abs_tol absolute tolerance given for doubles t0 and tmax to account for small deviations from whole days.
42-
* @return interpolated time series
43+
* @param simulation_result A time series to interpolate.
44+
* @param abs_tol Optional parameter to set the absolute tolerance used to account for small deviations from whole days.
45+
* Must be less then 1. The default tolerance is chosen to minimize the chances of "loosing" days due to rounding.
46+
* @return An interpolated time series with integer valued times.
4347
*/
4448
template <typename FP>
45-
TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result, const FP abs_tol = 1e-14);
49+
TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result,
50+
const FP abs_tol = FP{100.} * Limits<FP>::zero_tolerance());
4651

4752
/**
48-
* @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.
49-
* values at new time points are linearly interpolated from their immediate neighbors from the old time points.
50-
* @param simulation_result time series to interpolate
51-
* @param interpolations_times std::vector of time points at which simulation results are interpolated.
52-
* @return interpolated time series at given interpolation points
53+
* @brief Interpolate a time series at the given time points.
54+
* New time points must be monotonic increasing, and lie in between time points of the given time series.
55+
* The values at new time points are linearly interpolated from their immediate neighbors within the given time series.
56+
* @param simulation_result The time series to interpolate.
57+
* @param interpolation_times An std::vector of time points at which simulation results are interpolated.
58+
* @return The interpolated time series at given interpolation points.
5359
*/
5460
template <typename FP>
5561
TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result,
56-
5762
const std::vector<FP>& interpolation_times);
5863

5964
/**
@@ -180,14 +185,15 @@ FP result_distance_2norm(const std::vector<mio::TimeSeries<FP>>& result1,
180185
template <typename FP>
181186
TimeSeries<FP> interpolate_simulation_result(const TimeSeries<FP>& simulation_result, const FP abs_tol)
182187
{
188+
assert(abs_tol < 1);
183189
using std::ceil;
184190
using std::floor;
185191
const auto t0 = simulation_result.get_time(0);
186192
const auto t_max = simulation_result.get_last_time();
187-
// add another day if the first time point is equal to day_0 up to absolute tolerance tol
188-
const auto day0 = (t0 - abs_tol < ceil(t0) - 1) ? floor(t0) : ceil(t0);
189-
// add another day if the last time point is equal to day_max up to absolute tolerance tol
190-
const auto day_max = (t_max + abs_tol > floor(t_max) + 1) ? ceil(t_max) : floor(t_max);
193+
// add an additional day, if the first time point is within tolerance of floor(t0)
194+
const auto day0 = ceil(t0 - abs_tol);
195+
// add an additional day, if the last time point is within tolerance of ceil(tmax)
196+
const auto day_max = floor(t_max + abs_tol);
191197

192198
// create interpolation_times vector with all days between day0 and day_max
193199
std::vector<FP> tps(static_cast<int>(day_max) - static_cast<int>(day0) + 1);

‎cpp/memilio/math/integrator.h‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -188,11 +188,12 @@ class SystemIntegrator
188188
using std::fabs;
189189
using std::max;
190190
using std::min;
191-
const FP t0 = results.get_last_time();
191+
const FP t0 = results.get_last_time();
192+
const FP eps = Limits<FP>::zero_tolerance();
192193
assert(tmax > t0);
193194
assert(dt > 0);
194195

195-
const size_t num_steps =
196+
const auto num_steps =
196197
static_cast<size_t>(ceil((tmax - t0) / dt)); // estimated number of time steps (if equidistant)
197198

198199
results.reserve(results.get_num_time_points() + num_steps);
@@ -204,7 +205,7 @@ class SystemIntegrator
204205
FP dt_min_restore = m_core->get_dt_min(); // used to restore dt_min, if it was decreased to reach tmax
205206
FP t = t0;
206207

207-
for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > 1e-10; ++i) {
208+
for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > eps; ++i) {
208209
// We don't make time steps too small as the error estimator of an adaptive integrator
209210
//may not be able to handle it. this is very conservative and maybe unnecessary,
210211
//but also unlikely to happen. may need to be reevaluated.
@@ -225,7 +226,7 @@ class SystemIntegrator
225226
results.get_last_time() = t;
226227

227228
// if dt has been changed by step, register the current m_core as adaptive.
228-
m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, Limits<FP>::zero_tolerance());
229+
m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, eps);
229230
}
230231
m_core->get_dt_min() = dt_min_restore; // restore dt_min
231232
// if dt was decreased to reach tmax in the last time iteration,
@@ -236,7 +237,7 @@ class SystemIntegrator
236237
if (!step_okay) {
237238
log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
238239
}
239-
else if (fabs((tmax - t) / (tmax - t0)) > 1e-14) {
240+
else if (fabs((tmax - t) / (tmax - t0)) > eps) {
240241
log_warning("Last time step too small. Could not reach tmax exactly.");
241242
}
242243
else {

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

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
/*
2+
* Copyright (C) 2020-2026 MEmilio
3+
*
4+
* Authors: Rene Schmieding
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
7+
*
8+
* Licensed under the Apache License, Version 2.0 (the "License");
9+
* you may not use this file except in compliance with the License.
10+
* You may obtain a copy of the License at
11+
*
12+
* http://www.apache.org/licenses/LICENSE-2.0
13+
*
14+
* Unless required by applicable law or agreed to in writing, software
15+
* distributed under the License is distributed on an "AS IS" BASIS,
16+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
* See the License for the specific language governing permissions and
18+
* limitations under the License.
19+
*/
20+
#ifndef PYMIO_DATA_ANALYZE_RESULT_H
21+
#define PYMIO_DATA_ANALYZE_RESULT_H
22+
23+
#include "memilio/data/analyze_result.h"
24+
#include "pybind_util.h"
25+
26+
#include "pybind11/pybind11.h"
27+
28+
namespace py = pybind11;
29+
30+
namespace pymio
31+
{
32+
33+
/**
34+
* @brief Bind functions interpolate_simulation_result and interpolate_ensemble_results.
35+
*/
36+
void bind_interpolate_result_methods(py::module_& m)
37+
{
38+
m.def(
39+
"interpolate_simulation_result",
40+
[](const mio::TimeSeries<double>& ts) {
41+
return mio::interpolate_simulation_result(ts);
42+
},
43+
"Interpolate a given time series with evenly spaced, integer time points that represent whole days.",
44+
py::arg("ts"));
45+
m.def("interpolate_simulation_result",
46+
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
47+
&mio::interpolate_simulation_result),
48+
"Interpolate a given time series with evenly spaced, integer time points that represent whole days.",
49+
py::arg("ts"), py::arg("abs_tol"));
50+
51+
m.def("interpolate_simulation_result",
52+
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
53+
&mio::interpolate_simulation_result),
54+
"Interpolate a time series at the given time points.", py::arg("ts"), py::arg("interpolation_times"));
55+
56+
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>,
57+
"Interpolates results of all runs with evenly spaced, integer time points that represent whole days.");
58+
}
59+
60+
} // namespace pymio
61+
62+
#endif // PYMIO_DATA_ANALYZE_RESULT_H

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

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include "epidemiology/populations.h"
2626
#include "utils/parameter_set.h"
2727
#include "utils/index.h"
28+
#include "data/analyze_result.h"
2829

2930
// Includes from MEmilio
3031
#include "ode_mseirs4/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::omseirs4::InfectionState>()
5152
PYBIND11_MODULE(_simulation_omseirs4, m)
5253
{
5354
// interpolation helpers
54-
m.def("interpolate_simulation_result",
55-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
56-
&mio::interpolate_simulation_result),
57-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
58-
59-
m.def("interpolate_simulation_result",
60-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
61-
&mio::interpolate_simulation_result),
62-
py::arg("ts"), py::arg("interpolation_times"));
63-
64-
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
55+
pymio::bind_interpolate_result_methods(m);
6556

6657
// InfectionState enum
6758
pymio::iterable_enum<mio::omseirs4::InfectionState>(m, "InfectionState")

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

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@
3333
#include "mobility/metapopulation_mobility_instant.h"
3434
#include "io/mobility_io.h"
3535
#include "io/result_io.h"
36+
#include "data/analyze_result.h"
3637

3738
//Includes from MEmilio
3839
#include "ode_secir/model.h"
@@ -87,17 +88,7 @@ PYBIND11_MAKE_OPAQUE(std::vector<MobilityGraph>);
8788
PYBIND11_MODULE(_simulation_osecir, m)
8889
{
8990
// https://github.com/pybind/pybind11/issues/1153
90-
m.def("interpolate_simulation_result",
91-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
92-
&mio::interpolate_simulation_result),
93-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
94-
95-
m.def("interpolate_simulation_result",
96-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
97-
&mio::interpolate_simulation_result),
98-
py::arg("ts"), py::arg("interpolation_times"));
99-
100-
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
91+
pymio::bind_interpolate_result_methods(m);
10192

10293
m.def("ensemble_mean", &mio::ensemble_mean<double>);
10394
m.def("ensemble_percentile", &mio::ensemble_percentile<double>);

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

Lines changed: 2 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@
3131
#include "epidemiology/populations.h"
3232
#include "io/mobility_io.h"
3333
#include "io/result_io.h"
34+
#include "data/analyze_result.h"
3435

3536
//Includes from MEmilio
3637
#include "ode_secirvvs/model.h"
@@ -79,15 +80,7 @@ PYBIND11_MAKE_OPAQUE(std::vector<MobilityGraph>);
7980

8081
PYBIND11_MODULE(_simulation_osecirvvs, m)
8182
{
82-
m.def("interpolate_simulation_result",
83-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
84-
&mio::interpolate_simulation_result),
85-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
86-
87-
m.def("interpolate_simulation_result",
88-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
89-
&mio::interpolate_simulation_result),
90-
py::arg("ts"), py::arg("interpolation_times"));
83+
pymio::bind_interpolate_result_methods(m);
9184

9285
pymio::iterable_enum<mio::osecirvvs::InfectionState>(m, "InfectionState")
9386
.value("SusceptibleNaive", mio::osecirvvs::InfectionState::SusceptibleNaive)

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

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
#include "compartments/compartmental_model.h"
2929
#include "epidemiology/age_group.h"
3030
#include "epidemiology/populations.h"
31+
#include "data/analyze_result.h"
3132

3233
//Includes from MEmilio
3334
#include "ode_seir/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::oseir::InfectionState>()
5152

5253
PYBIND11_MODULE(_simulation_oseir, m)
5354
{
54-
m.def("interpolate_simulation_result",
55-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
56-
&mio::interpolate_simulation_result),
57-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
58-
59-
m.def("interpolate_simulation_result",
60-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
61-
&mio::interpolate_simulation_result),
62-
py::arg("ts"), py::arg("interpolation_times"));
63-
64-
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
55+
pymio::bind_interpolate_result_methods(m);
6556

6657
pymio::iterable_enum<mio::oseir::InfectionState>(m, "InfectionState")
6758
.value("Susceptible", mio::oseir::InfectionState::Susceptible)

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

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
#include "compartments/compartmental_model.h"
2929
#include "epidemiology/age_group.h"
3030
#include "epidemiology/populations.h"
31+
#include "data/analyze_result.h"
3132

3233
//Includes from MEmilio
3334
#include "ode_seirdb/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::oseirdb::InfectionState>()
5152

5253
PYBIND11_MODULE(_simulation_oseirdb, m)
5354
{
54-
m.def("interpolate_simulation_result",
55-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
56-
&mio::interpolate_simulation_result),
57-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
58-
59-
m.def("interpolate_simulation_result",
60-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
61-
&mio::interpolate_simulation_result),
62-
py::arg("ts"), py::arg("interpolation_times"));
63-
64-
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
55+
pymio::bind_interpolate_result_methods(m);
6556

6657
pymio::iterable_enum<mio::oseirdb::InfectionState>(m, "InfectionState")
6758
.value("Susceptible", mio::oseirdb::InfectionState::Susceptible)

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

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
#include "compartments/compartmental_model.h"
2929
#include "epidemiology/age_group.h"
3030
#include "epidemiology/populations.h"
31+
#include "data/analyze_result.h"
3132

3233
//Includes from MEmilio
3334
#include "ode_sir/model.h"
@@ -51,17 +52,7 @@ inline std::string pretty_name<mio::osir::InfectionState>()
5152

5253
PYBIND11_MODULE(_simulation_osir, m)
5354
{
54-
m.def("interpolate_simulation_result",
55-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
56-
&mio::interpolate_simulation_result),
57-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
58-
59-
m.def("interpolate_simulation_result",
60-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
61-
&mio::interpolate_simulation_result),
62-
py::arg("ts"), py::arg("interpolation_times"));
63-
64-
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
55+
pymio::bind_interpolate_result_methods(m);
6556

6657
pymio::iterable_enum<mio::osir::InfectionState>(m, "InfectionState")
6758
.value("Susceptible", mio::osir::InfectionState::Susceptible)

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

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include "utils/custom_index_array.h"
2626
#include "utils/parameter_set.h"
2727
#include "epidemiology/populations.h"
28+
#include "data/analyze_result.h"
2829

2930
//Includes from MEmilio
3031
#include "sde_sir/model.h"
@@ -52,17 +53,7 @@ inline std::string pretty_name<mio::ssir::InfectionState>()
5253

5354
PYBIND11_MODULE(_simulation_ssir, m)
5455
{
55-
m.def("interpolate_simulation_result",
56-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const double)>(
57-
&mio::interpolate_simulation_result),
58-
py::arg("ts"), py::arg("abs_tol") = 1e-14);
59-
60-
m.def("interpolate_simulation_result",
61-
static_cast<mio::TimeSeries<double> (*)(const mio::TimeSeries<double>&, const std::vector<double>&)>(
62-
&mio::interpolate_simulation_result),
63-
py::arg("ts"), py::arg("interpolation_times"));
64-
65-
m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results<mio::TimeSeries<double>>);
56+
pymio::bind_interpolate_result_methods(m);
6657

6758
pymio::iterable_enum<mio::ssir::InfectionState>(m, "InfectionState")
6859
.value("Susceptible", mio::ssir::InfectionState::Susceptible)

0 commit comments

Comments
 (0)