Make local copies of the step callbacks in ensemble propagations, and…

… docs updates.

doc/changelog.rst

@@ -7,6 +7,10 @@ Changelog
 New
 ~~~
 
+- The step callbacks can now optionally implement a ``pre_hook()``
+  member function that will be called before the first step
+  is taken by a ``propagate_*()`` function
+  (`#334 <https://github.com/bluescarni/heyoka/pull/334>`__).
 - The heyoka library now passes all ``clang-tidy`` checks
   (`#315 <https://github.com/bluescarni/heyoka/pull/315>`__).
 - Introduce several vectorised overloads in the expression
@@ -39,6 +43,10 @@ New
 Changes
 ~~~~~~~
 
+- The step callbacks are now copied in :ref:`ensemble propagations <tut_ensemble>`
+  rather then being shared among threads. The aim of this change
+  is to reduce the likelihood of data races
+  (`#334 <https://github.com/bluescarni/heyoka/pull/334>`__).
 - Comprehensive overhaul of the expression system, including:
   enhanced automatic simplification capabilities for sums,
   products and powers, removal of several specialised primitives

doc/tut_adaptive.rst

@@ -234,7 +234,7 @@ integration will be ``taylor_outcome::err_nf_state``.
 .. versionadded:: 0.7.0
 
 The ``propagate_for()`` and ``propagate_until()`` functions
-can be invoked with two additional optional keyword arguments:
+can be invoked with additional optional keyword arguments:
 
 - ``max_delta_t``: similarly to the ``step()`` function, this value
   represents the maximum timestep size in absolute value;
@@ -249,6 +249,21 @@ can be invoked with two additional optional keyword arguments:
   will continue after the invocation of the callback, otherwise the integration
   will be interrupted.
 
+  .. versionadded:: 1.0.0
+
+  Optionally, a function object callback can implement a ``pre_hook()`` member
+  function with signature
+
+  .. code-block:: c++
+
+     void pre_hook(taylor_adaptive<double> &);
+
+
+  that will be invoked once *before* the first step is taken
+  by the ``propagate_for()`` and ``propagate_until()`` functions.
+- ``c_output``: a boolean flag that enables :ref:`continuous output <tut_c_output>`.
+
+
 Propagation over a time grid
 ----------------------------
 
@@ -297,7 +312,7 @@ fact that they must be finite and ordered monotonically).
 .. versionadded:: 0.7.0
 
 The ``propagate_grid()`` function
-can be invoked with two additional optional keyword arguments:
+can be invoked with additional optional keyword arguments:
 
 - ``max_delta_t``: similarly to the ``step()`` function, this value
   represents the maximum timestep size in absolute value;
@@ -312,6 +327,19 @@ can be invoked with two additional optional keyword arguments:
   will continue after the invocation of the callback, otherwise the integration
   will be interrupted.
 
+  .. versionadded:: 1.0.0
+
+  Optionally, a function object callback can implement a ``pre_hook()`` member
+  function with signature
+
+  .. code-block:: c++
+
+     void pre_hook(taylor_adaptive<double> &);
+
+
+  that will be invoked once *before* the first step is taken
+  by the ``propagate_grid()`` function.
+
 Full code listing
 -----------------
 

doc/tut_ensemble.rst

@@ -66,7 +66,7 @@ The signature of the generator ``gen`` reads:
 
 That is, the generator takes in input a *copy* of the template integrator ``ta``
 and an iteration index ``idx`` in the ``[0, n_iter)`` range. ``gen`` is then expected
-to modify the copy of ``ta`` (e.g., by setting its initial conditions to specific
+to modify ``ta`` (e.g., by setting its initial conditions to specific
 values) and return it.
 
 The ``ensemble_propagate_until()`` function iterates over
@@ -151,22 +151,17 @@ In an ensemble propagation, it is thus important to keep in mind that
 the following actions may be performed
 concurrently by separate threads of execution:
 
-* invocation of the generator's call operator and of the call operator
-  of the callback that can (optionally) be passed to the ``propagate_*()``
-  functions. In other words, both the generator and the ``propagate_*()``
-  callback are shared among several threads of execution and used
-  concurrently;
-* copy construction of the events callbacks and invocation of the
-  call operator on the copies. That is, each thread of execution
+* invocation of the generator's call operator (that is, the generator is
+  shared among multiple threads and must support concurrent invocations
+  of its call operator);
+* copy construction of the callback that can (optionally) be passed to
+  the ``propagate_*()`` functions (that is, step callbacks are **not**
+  shared among multiple threads of execution, and a new copy is created for each
+  invocation of the ``propagate_*()`` functions);
+* copy construction of the events callbacks (that is, each thread of execution
   gets its own copy of the event callbacks thanks to the creation
-  of a new integrator object via the generator.
+  of a new integrator object via the generator).
 
 For instance, an event callback which performs write operations
 on a global variable without using some form of synchronisation
 will result in undefined behaviour when used in an ensemble propagation.
-
-Another example of unsafe usage is a ``propagate_*()`` callback that
-performs write operations into its own data member(s) without
-synchronisation: because the
-``propagate_*()`` callback is shared among several threads, such usage results
-in a data race.

src/ensemble_propagate.cpp

@@ -12,14 +12,14 @@
 #include <cstddef>
 #include <cstdint>
 #include <functional>
-#include <limits>
 #include <optional>
 #include <stdexcept>
 #include <tuple>
 #include <utility>
 #include <vector>
 
 #include <boost/numeric/conversion/cast.hpp>
+#include <boost/safe_numerics/safe_integer.hpp>
 
 #include <tbb/blocked_range.h>
 #include <tbb/parallel_for.h>
@@ -43,9 +43,9 @@
 
 // NOTE: these actions will be performed concurrently from
 // multiple threads of execution:
-// - invocation of the generator's call operator and of the propagate callback,
-// - copy construction of the events' callbacks and invocation of the call operator
-//   on the copies.
+// - invocation of the generator's call operator,
+// - copy construction of the events' and step callbacks and
+//   invocation of the call operator on the copies.
 
 HEYOKA_BEGIN_NAMESPACE
 
@@ -78,10 +78,13 @@ ensemble_propagate_until_impl(const taylor_adaptive<T> &ta, T t, std::size_t n_i
             // Generate the integrator for the current iteration.
             auto local_ta = gen(ta, i);
 
+            // Make a local copy of the callback.
+            auto local_cb = cb;
+
             // Do the propagation.
             auto loc_ret
                 = local_ta.propagate_until(t, kw::max_steps = max_steps, kw::max_delta_t = max_delta_t,
-                                           kw::callback = cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
+                                           kw::callback = local_cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
 
             // Assign the results.
             opt_retval[i].emplace(std::tuple_cat(std::make_tuple(std::move(local_ta)), std::move(loc_ret)));
@@ -123,10 +126,13 @@ ensemble_propagate_for_impl(const taylor_adaptive<T> &ta, T delta_t, std::size_t
             // Generate the integrator for the current iteration.
             auto local_ta = gen(ta, i);
 
+            // Make a local copy of the callback.
+            auto local_cb = cb;
+
             // Do the propagation.
             auto loc_ret
                 = local_ta.propagate_for(delta_t, kw::max_steps = max_steps, kw::max_delta_t = max_delta_t,
-                                         kw::callback = cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
+                                         kw::callback = local_cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
 
             // Assign the results.
             opt_retval[i].emplace(std::tuple_cat(std::make_tuple(std::move(local_ta)), std::move(loc_ret)));
@@ -166,9 +172,12 @@ ensemble_propagate_grid_impl(const taylor_adaptive<T> &ta, std::vector<T> grid,
             // Generate the integrator for the current iteration.
             auto local_ta = gen(ta, i);
 
+            // Make a local copy of the callback.
+            auto local_cb = cb;
+
             // Do the propagation.
             auto loc_ret = local_ta.propagate_grid(grid, kw::max_steps = max_steps, kw::max_delta_t = max_delta_t,
-                                                   kw::callback = cb);
+                                                   kw::callback = local_cb);
 
             // Assign the results.
             opt_retval[i].emplace(std::tuple_cat(std::make_tuple(std::move(local_ta)), std::move(loc_ret)));
@@ -246,10 +255,13 @@ ensemble_propagate_until_batch_impl(
             // Generate the integrator for the current iteration.
             auto local_ta = gen(ta, i);
 
+            // Make a local copy of the callback.
+            auto local_cb = cb;
+
             // Do the propagation.
             auto loc_ret
                 = local_ta.propagate_until(t, kw::max_steps = max_steps, kw::max_delta_t = max_delta_ts,
-                                           kw::callback = cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
+                                           kw::callback = local_cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
 
             // Assign the results.
             opt_retval[i].emplace(std::move(local_ta), std::move(loc_ret));
@@ -290,10 +302,13 @@ ensemble_propagate_for_batch_impl(
             // Generate the integrator for the current iteration.
             auto local_ta = gen(ta, i);
 
+            // Make a local copy of the callback.
+            auto local_cb = cb;
+
             // Do the propagation.
             auto loc_ret
                 = local_ta.propagate_for(delta_t, kw::max_steps = max_steps, kw::max_delta_t = max_delta_ts,
-                                         kw::callback = cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
+                                         kw::callback = local_cb, kw::write_tc = write_tc, kw::c_output = with_c_out);
 
             // Assign the results.
             opt_retval[i].emplace(std::move(local_ta), std::move(loc_ret));
@@ -320,13 +335,8 @@ std::vector<std::tuple<taylor_adaptive_batch<T>, std::vector<T>>> ensemble_propa
     assert(batch_size != 0u); // LCOV_EXCL_LINE
 
     // Splat out the time grid.
-    // LCOV_EXCL_START
-    if (grid_.size() > std::numeric_limits<decltype(grid_.size())>::max() / batch_size) {
-        throw std::overflow_error("Overflow detected in an ensemble propagation");
-    }
-    // LCOV_EXCL_STOP
     std::vector<T> grid;
-    grid.reserve(grid_.size() * batch_size);
+    grid.reserve(boost::safe_numerics::safe<decltype(grid.size())>(grid_.size()) * batch_size);
     for (auto gval : grid_) {
         for (std::uint32_t i = 0; i < batch_size; ++i) {
             grid.push_back(gval);
@@ -350,9 +360,12 @@ std::vector<std::tuple<taylor_adaptive_batch<T>, std::vector<T>>> ensemble_propa
             // Generate the integrator for the current iteration.
             auto local_ta = gen(ta, i);
 
+            // Make a local copy of the callback.
+            auto local_cb = cb;
+
             // Do the propagation.
             auto loc_ret = local_ta.propagate_grid(grid, kw::max_steps = max_steps, kw::max_delta_t = max_delta_ts,
-                                                   kw::callback = cb);
+                                                   kw::callback = local_cb);
 
             // Assign the results.
             opt_retval[i].emplace(std::move(local_ta), std::move(loc_ret));