Merge pull request #373 from bluescarni/pr/propagate_grid_api_changes

Change propagate_grid() to require a grid beginning with current time

.github/workflows/gha_ci.yml

@@ -34,7 +34,7 @@ jobs:
       - name: Build
         shell: pwsh
         run: |
-          conda install -y cmake 'llvmdev=13.*' tbb-devel tbb boost-cpp xtensor xtensor-blas blas blas-devel fmt spdlog sleef zlib libzlib 'mppp=1.*'
+          conda install -y cmake 'llvmdev=13.*' tbb-devel tbb libboost-devel xtensor xtensor-blas blas blas-devel fmt spdlog sleef zlib libzlib 'mppp=1.*'
           mkdir build
           cd build
           cmake ../ -G "Visual Studio 17 2022" -A x64 -DHEYOKA_BUILD_TESTS=yes -DHEYOKA_WITH_MPPP=yes -DHEYOKA_BUILD_TUTORIALS=ON -DHEYOKA_ENABLE_IPO=yes -DBoost_NO_BOOST_CMAKE=ON -DHEYOKA_WITH_SLEEF=yes -DMPPP_GMP_INCLUDE_DIR=C:\Miniconda\envs\test\Library\include -DMPPP_GMP_LIBRARY=C:\Miniconda\envs\test\Library\lib\mpir.lib
@@ -56,7 +56,7 @@ jobs:
       - name: Build
         shell: pwsh
         run: |
-          conda install -y cmake 'llvmdev=14.*' tbb-devel tbb boost-cpp xtensor xtensor-blas blas blas-devel fmt spdlog sleef zlib libzlib 'mppp=1.*'
+          conda install -y cmake 'llvmdev=14.*' tbb-devel tbb libboost-devel xtensor xtensor-blas blas blas-devel fmt spdlog sleef zlib libzlib 'mppp=1.*'
           mkdir build
           cd build
           cmake ../ -G "Visual Studio 17 2022" -A x64 -DHEYOKA_BUILD_TESTS=yes -DHEYOKA_WITH_MPPP=yes -DHEYOKA_BUILD_TUTORIALS=ON -DHEYOKA_ENABLE_IPO=yes -DBoost_NO_BOOST_CMAKE=ON -DHEYOKA_WITH_SLEEF=yes -DMPPP_GMP_INCLUDE_DIR=C:\Miniconda\envs\test\Library\include -DMPPP_GMP_LIBRARY=C:\Miniconda\envs\test\Library\lib\mpir.lib

doc/breaking_changes.rst

@@ -10,6 +10,20 @@ Breaking changes
 
 heyoka 4 includes several backwards-incompatible changes.
 
+API/behaviour changes
+~~~~~~~~~~~~~~~~~~~~~
+
+``propagate_grid()``
+^^^^^^^^^^^^^^^^^^^^
+
+The ``propagate_grid()`` functions of the adaptive integrators now require the first element of the
+time grid to be equal to the current integrator time. Previously, in case of a difference between the
+integrator time and the first grid point, heyoka would propagate the state of the system up to the
+first grid point with ``propagate_until()``.
+
+If you want to recover the previous behaviour, you will have to invoke manually ``propagate_until(grid[0])``
+before invoking ``propagate_grid()``.
+
 General
 ~~~~~~~
 

doc/changelog.rst

@@ -7,6 +7,11 @@ Changelog
 Changes
 ~~~~~~~
 
+- **BREAKING**: the ``propagate_grid()`` functions of the
+  adaptive integrators now require the first element of the
+  time grid to be equal to the current integrator time
+  (`#373 <https://github.com/bluescarni/heyoka/pull/373>`__).
+  This is a :ref:`breaking change <bchanges_4_0_0>`.
 - Move the declarations of all :ref:`keyword arguments <kwargs>`
   into the ``kw.hpp`` header
   (`#372 <https://github.com/bluescarni/heyoka/pull/372>`__).

doc/install.rst

@@ -7,9 +7,9 @@ Introduction
 ------------
 
 heyoka is written in modern C++, and it requires a compiler able to understand
-at least C++20. The library is regularly tested on
-a continuous integration pipeline which currently includes several
-compilers (GCC, Clang, MSVC) on several operating systems (Linux, OSX, Windows)
+at least C++20. Specifically, heyoka currently targets GCC>=10, clang>=14 and MSVC>=2022.
+The library is regularly tested on a continuous integration pipeline
+which includes several operating systems (Linux, OSX, Windows)
 and several CPU architectures (x86-64, 64-bit ARM and 64-bit PowerPC).
 
 heyoka has the following **mandatory** dependencies:

include/heyoka/detail/dfloat.hpp

@@ -103,7 +103,7 @@ inline std::pair<F, F> eft_add_dekker(F a, F b)
     auto x = a + b;
     auto y = (a - x) + b;
 
-    return {x, y};
+    return std::make_pair(std::move(x), std::move(y));
 }
 
 // Error-free transformation of the sum of two floating point numbers.
@@ -114,17 +114,17 @@ inline std::pair<F, F> eft_add_knuth(F a, F b)
 {
     auto x = a + b;
     auto z = x - a;
-    auto y = (a - (x - z)) + (b - z);
+    auto y = (a - (x - z)) + (std::move(b) - z);
 
-    return {x, y};
+    return std::make_pair(std::move(x), std::move(y));
 }
 
 // Normalise a double-length float.
 // Taken from:
 // https://github.com/fhajji/ntl/blob/6918e6b80336cee34f2131fcf71a58c72b931174/src/quad_float.cpp#L125
 // NOTE: this is based on the error-free trasformation requiring abs(x.hi) >= abs(x.lo).
 template <typename F>
-inline dfloat<F> normalise(const dfloat<F> &x)
+inline dfloat<F> normalise(dfloat<F> x)
 {
     // LCOV_EXCL_START
 #if !defined(NDEBUG)
@@ -136,9 +136,9 @@ inline dfloat<F> normalise(const dfloat<F> &x)
 #endif
     // LCOV_EXCL_STOP
 
-    auto [u, v] = eft_add_dekker(x.hi, x.lo);
+    auto [u, v] = eft_add_dekker(std::move(x.hi), std::move(x.lo));
 
-    return dfloat<F>(u, v);
+    return dfloat<F>(std::move(u), std::move(v));
 }
 
 // NOTE: taken with minimal adaptations from NTL.
@@ -161,7 +161,7 @@ inline dfloat<F> operator+(const dfloat<F> &a, const dfloat<F> &b)
     auto [u, v] = eft_add_dekker(x_hi, y_hi + x_lo);
     std::tie(u, v) = eft_add_dekker(u, v + y_lo);
 
-    return dfloat<F>(u, v);
+    return dfloat<F>(std::move(u), std::move(v));
 }
 
 // Subtraction.

include/heyoka/step_callback.hpp

@@ -95,13 +95,16 @@ struct HEYOKA_DLL_PUBLIC_INLINE_CLASS step_cb_ref_iface {
 template <typename TA>
 struct HEYOKA_DLL_PUBLIC_INLINE_CLASS step_cb_ifaceT {
     template <typename Holder, typename T>
-    using type = tanuki::composite_wrap_interfaceT<callable<bool(TA &)>, pre_hook_wrap_t<TA>>::template type<Holder, T>;
+    // NOTE: clang 14 requires typename here, hopefully it does
+    // not do any harm in other compilers.
+    using type =
+        typename tanuki::composite_wrap_interfaceT<callable<bool(TA &)>, pre_hook_wrap_t<TA>>::template type<Holder, T>;
 };
 
 // Configuration.
 template <typename TA>
 inline constexpr auto step_cb_wrap_config = tanuki::config<empty_callable, step_cb_ref_iface<TA>::template type>{
-    // Similarly to std::function, ensure that callable can store
+    // Similarly to std::function, ensure that step_callback can store
     // in static storage pointers and reference wrappers.
     // NOTE: reference wrappers are not guaranteed to have the size
     // of a pointer, but in practice that should always be the case.

include/heyoka/taylor.hpp

@@ -1292,6 +1292,15 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada
     // - continuous output, if requested (only for propagate_for/until()).
     // NOTE: the min/max timesteps are well-defined
     // only if at least 1-2 steps were taken successfully.
+    // NOTE: the propagate_*() functions are not guaranteed to bring
+    // the integrator time *exactly* to the requested final time. This
+    // ultimately stems from the fact that in floating-point arithmetics
+    // in general a + (b - a) != b, and this happens regardless of the
+    // use of a double-length time representation. This occurrence however
+    // seems to be pretty rare in practice, so for the time being we leave
+    // this as it is and just document the corner-case behaviour. Perhaps
+    // in the future we can offer a stronger guarantee, which however will
+    // result in a more complicated logic.
     template <typename... KwArgs>
     std::tuple<taylor_outcome, T, T, std::size_t, std::optional<continuous_output<T>>>
     propagate_until(T t, KwArgs &&...kw_args)

src/taylor_00.cpp

@@ -1645,6 +1645,7 @@ taylor_adaptive<T>::propagate_grid_impl(std::vector<T> grid, std::size_t max_ste
 
 #if defined(HEYOKA_HAVE_REAL)
 
+    // Ensure all grid points have the correct precision for real.
     if constexpr (std::is_same_v<T, mppp::real>) {
         for (auto &val : grid) {
             val.prec_round(this->get_prec());
@@ -1678,15 +1679,19 @@ taylor_adaptive<T>::propagate_grid_impl(std::vector<T> grid, std::size_t max_ste
         }
     }
 
+    // Require that the user provides a grid starting from the
+    // current integrator time, modulo the double-length correction.
+    if (m_time.hi != grid[0]) {
+        throw std::invalid_argument(
+            fmt::format("When invoking propagate_grid(), the first element of the time grid "
+                        "must match the current time coordinate - however, the first element of the time grid has a "
+                        "value of {}, while the current time coordinate is {}",
+                        grid[0], m_time.hi));
+    }
+
     // Pre-allocate the return value.
     std::vector<T> retval;
-    // LCOV_EXCL_START
-    if (get_dim() > std::numeric_limits<decltype(retval.size())>::max() / grid.size()) {
-        throw std::overflow_error("Overflow detected in the creation of the return value of propagate_grid() in an "
-                                  "adaptive Taylor integrator");
-    }
-    // LCOV_EXCL_STOP
-    retval.reserve(grid.size() * get_dim());
+    retval.reserve(boost::safe_numerics::safe<decltype(retval.size())>(grid.size()) * get_dim());
 
     // Initial values for the counters
     // and the min/max abs of the integration
@@ -1701,6 +1706,10 @@ taylor_adaptive<T>::propagate_grid_impl(std::vector<T> grid, std::size_t max_ste
     T min_h = detail::num_inf_like(max_delta_t), max_h = detail::num_zero_like(max_delta_t);
 
     // Propagate the system up to the first grid point.
+    // This is necessary in order to account for the fact
+    // that the user cannot pass as starting point in the grid
+    // a time coordinate which is *exactly* equal to m_time,
+    // due to the usage of a double-length representation.
     // NOTE: we pass write_tc = true because some grid
     // points after the first one might end up being
     // calculated via dense output *before*
@@ -3572,24 +3581,36 @@ std::vector<T> taylor_adaptive_batch<T>::propagate_grid_impl(const std::vector<T
         }
     }
 
+    // Require that the user provides a grid starting from the
+    // current integrator time, modulo the double-length correction.
+    for (std::uint32_t i = 0; i < m_batch_size; ++i) {
+        if (m_time_hi[i] != grid_ptr[i]) {
+            throw std::invalid_argument(
+                fmt::format("When invoking propagate_grid(), the first element of the time grid "
+                            "must match the current time coordinate - however, the first element of the time grid at "
+                            "batch index {} has a "
+                            "value of {}, while the current time coordinate is {}",
+                            i, grid_ptr[i], m_time_hi[i]));
+        }
+    }
+
     // Pre-allocate the return value.
     std::vector<T> retval;
-    // LCOV_EXCL_START
-    if (get_dim() > std::numeric_limits<decltype(retval.size())>::max() / grid.size()) {
-        throw std::overflow_error("Overflow detected in the creation of the return value of propagate_grid() in an "
-                                  "adaptive Taylor integrator in batch mode");
-    }
-    // LCOV_EXCL_STOP
     // NOTE: fill with NaNs, so that the missing entries
     // are signalled with NaN if we exit early.
-    retval.resize(grid.size() * get_dim(), std::numeric_limits<T>::quiet_NaN());
+    retval.resize(boost::safe_numerics::safe<decltype(retval.size())>(grid.size()) * get_dim(),
+                  std::numeric_limits<T>::quiet_NaN());
 
     // NOTE: this is a buffer of size m_batch_size
     // that is used in various places as temp storage.
     std::vector<T> pgrid_tmp;
     pgrid_tmp.resize(boost::numeric_cast<decltype(pgrid_tmp.size())>(m_batch_size));
 
     // Propagate the system up to the first batch of grid points.
+    // This is necessary in order to account for the fact
+    // that the user cannot pass as starting point in the grid
+    // a time coordinate which is *exactly* equal to m_time,
+    // due to the usage of a double-length representation.
     // NOTE: we pass write_tc = true because some grid
     // points after the first one might end up being
     // calculated via dense output *before*

test/dfloat_time.cpp

@@ -123,6 +123,7 @@ TEST_CASE("scalar test")
 
             const auto grid = std::vector{fp_t{1.}, fp_t{10.}, final_time};
 
+            ta.propagate_until(fp_t{1.});
             const auto out = std::get<4>(ta.propagate_grid(grid));
 
             for (auto j = 0u; j < 3u; ++j) {
@@ -147,6 +148,7 @@ TEST_CASE("scalar test")
 
             const auto grid = std::vector{fp_t{-1.}, fp_t{-10.}, -final_time};
 
+            ta.propagate_until(fp_t{-1.});
             const auto out = std::get<4>(ta.propagate_grid(grid));
 
             for (auto j = 0u; j < 3u; ++j) {
@@ -251,6 +253,7 @@ TEST_CASE("batch test")
 
             const auto grid = std::vector{fp_t{1.}, fp_t{2.}, fp_t{10.}, fp_t{20.}, fp_t(10000.), fp_t(11000.)};
 
+            ta.propagate_until({fp_t{1.}, fp_t{2}});
             const auto out = ta.propagate_grid(grid);
 
             for (auto j = 0u; j < 3u; ++j) {
@@ -281,6 +284,7 @@ TEST_CASE("batch test")
 
             const auto grid = std::vector{fp_t{-1.}, fp_t{-2.}, fp_t{-10.}, fp_t{-20.}, fp_t(-10000.), fp_t(-11000.)};
 
+            ta.propagate_until({fp_t{-1.}, fp_t{-2}});
             const auto out = ta.propagate_grid(grid);
 
             for (auto j = 0u; j < 3u; ++j) {

test/step_callback.cpp

@@ -438,7 +438,7 @@ TEST_CASE("step_callback pre_hook")
         REQUIRE(ta0.get_pars()[0] == 1.5);
 
         REQUIRE_THROWS_MATCHES(
-            ta0.propagate_grid({3., 4.}, kw::callback = tm_cb{}), std::runtime_error,
+            ta0.propagate_grid({ta0.get_time(), ta0.get_time() + 1}, kw::callback = tm_cb{}), std::runtime_error,
             Message("The invocation of the callback passed to propagate_grid() resulted in the alteration of the "
                     "time coordinate of the integrator - this is not supported"));
     }
@@ -479,7 +479,8 @@ TEST_CASE("step_callback pre_hook")
         REQUIRE(ta0.get_pars()[1] == 1.5);
 
         REQUIRE_THROWS_MATCHES(
-            ta0.propagate_grid({3., 3., 4., 4.}, kw::callback = tm_cb{}), std::runtime_error,
+            ta0.propagate_grid({ta0.get_time()[0], ta0.get_time()[1], 4., 4.}, kw::callback = tm_cb{}),
+            std::runtime_error,
             Message("The invocation of the callback passed to propagate_grid() resulted in the alteration of the "
                     "time coordinate of the integrator - this is not supported"));
     }