diff --git a/cpp/memilio/data/analyze_result.h b/cpp/memilio/data/analyze_result.h index 430cc0e28d..e78e0c6508 100644 --- a/cpp/memilio/data/analyze_result.h +++ b/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 +#include #include 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 -TimeSeries interpolate_simulation_result(const TimeSeries& simulation_result, const FP abs_tol = 1e-14); +TimeSeries interpolate_simulation_result(const TimeSeries& simulation_result, + const FP abs_tol = FP{100.} * Limits::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 TimeSeries interpolate_simulation_result(const TimeSeries& simulation_result, - const std::vector& interpolation_times); /** @@ -180,14 +185,15 @@ FP result_distance_2norm(const std::vector>& result1, template TimeSeries interpolate_simulation_result(const TimeSeries& 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 tps(static_cast(day_max) - static_cast(day0) + 1); diff --git a/cpp/memilio/math/integrator.h b/cpp/memilio/math/integrator.h index 1d9aa83332..a5f6ad8a48 100755 --- a/cpp/memilio/math/integrator.h +++ b/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::zero_tolerance(); assert(tmax > t0); assert(dt > 0); - const size_t num_steps = + const auto num_steps = static_cast(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(dt, dt_copy, Limits::zero_tolerance()); + m_is_adaptive |= !floating_point_equal(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 { diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/data/analyze_result.h b/pycode/memilio-simulation/memilio/simulation/bindings/data/analyze_result.h new file mode 100644 index 0000000000..35a8fab37c --- /dev/null +++ b/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 +* +* 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& 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 (*)(const mio::TimeSeries&, 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 (*)(const mio::TimeSeries&, const std::vector&)>( + &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>, + "Interpolates results of all runs with evenly spaced, integer time points that represent whole days."); +} + +} // namespace pymio + +#endif // PYMIO_DATA_ANALYZE_RESULT_H diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp index 13af660254..0cd77f8223 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/omseirs4.cpp +++ b/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() PYBIND11_MODULE(_simulation_omseirs4, m) { // interpolation helpers - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); // InfectionState enum pymio::iterable_enum(m, "InfectionState") diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp index 799477ba9d..46facf0d03 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecir.cpp +++ b/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); PYBIND11_MODULE(_simulation_osecir, m) { // https://github.com/pybind/pybind11/issues/1153 - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); m.def("ensemble_mean", &mio::ensemble_mean); m.def("ensemble_percentile", &mio::ensemble_percentile); diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp index 3ef274963e..288c95fd95 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osecirvvs.cpp +++ b/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); PYBIND11_MODULE(_simulation_osecirvvs, m) { - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); + pymio::bind_interpolate_result_methods(m); pymio::iterable_enum(m, "InfectionState") .value("SusceptibleNaive", mio::osecirvvs::InfectionState::SusceptibleNaive) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/oseir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/oseir.cpp index d62c8c9493..42bd8edc4c 100755 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/oseir.cpp +++ b/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() PYBIND11_MODULE(_simulation_oseir, m) { - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); pymio::iterable_enum(m, "InfectionState") .value("Susceptible", mio::oseir::InfectionState::Susceptible) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/oseirdb.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/oseirdb.cpp index 1a28476729..68c1a991a1 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/oseirdb.cpp +++ b/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() PYBIND11_MODULE(_simulation_oseirdb, m) { - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); pymio::iterable_enum(m, "InfectionState") .value("Susceptible", mio::oseirdb::InfectionState::Susceptible) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/osir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/osir.cpp index 4f386e4eb0..6fe0cc20c7 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/osir.cpp +++ b/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() PYBIND11_MODULE(_simulation_osir, m) { - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); pymio::iterable_enum(m, "InfectionState") .value("Susceptible", mio::osir::InfectionState::Susceptible) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/ssir.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/ssir.cpp index 9503153e46..56d510ccd0 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/ssir.cpp +++ b/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() PYBIND11_MODULE(_simulation_ssir, m) { - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); pymio::iterable_enum(m, "InfectionState") .value("Susceptible", mio::ssir::InfectionState::Susceptible) diff --git a/pycode/memilio-simulation/memilio/simulation/bindings/models/ssirs.cpp b/pycode/memilio-simulation/memilio/simulation/bindings/models/ssirs.cpp index 77904208ef..e85d0372ed 100644 --- a/pycode/memilio-simulation/memilio/simulation/bindings/models/ssirs.cpp +++ b/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() PYBIND11_MODULE(_simulation_ssirs, m) { - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const double)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("abs_tol") = 1e-14); - - m.def("interpolate_simulation_result", - static_cast (*)(const mio::TimeSeries&, const std::vector&)>( - &mio::interpolate_simulation_result), - py::arg("ts"), py::arg("interpolation_times")); - - m.def("interpolate_ensemble_results", &mio::interpolate_ensemble_results>); + pymio::bind_interpolate_result_methods(m); pymio::iterable_enum(m, "InfectionState") .value("Susceptible", mio::ssirs::InfectionState::Susceptible)