Merge pull request #449 from bluescarni/pr/ta_empty_init_state

Empty initial state vectors

doc/changelog.rst

@@ -7,6 +7,10 @@ Changelog
 New
 ~~~
 
+- It is now possible to initialise a Taylor integrator
+  with an empty initial state vector. This will result
+  in zero-initialization of the state vector
+  (`#449 <https://github.com/bluescarni/heyoka/pull/449 >`__).
 - Implement parallel compilation for Taylor integrators
   and compiled functions
   (`#446 <https://github.com/bluescarni/heyoka/pull/446>`__,

include/heyoka/taylor.hpp

@@ -50,6 +50,7 @@
 #include <heyoka/continuous_output.hpp>
 #include <heyoka/detail/dfloat.hpp>
 #include <heyoka/detail/fwd_decl.hpp>
+#include <heyoka/detail/igor.hpp>
 #include <heyoka/detail/llvm_fwd.hpp>
 #include <heyoka/detail/llvm_helpers.hpp>
 #include <heyoka/detail/type_traits.hpp>
@@ -533,27 +534,30 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada
     taylor_adaptive();
 
     template <typename... KwArgs>
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
     explicit taylor_adaptive(std::vector<std::pair<expression, expression>> sys, std::vector<T> state,
                              const KwArgs &...kw_args)
         : taylor_adaptive(private_ctor_t{}, llvm_state(kw_args...))
     {
         finalise_ctor(std::move(sys), std::move(state), kw_args...);
     }
     template <typename... KwArgs>
-    explicit taylor_adaptive(std::vector<std::pair<expression, expression>> sys, std::initializer_list<T> state,
-                             const KwArgs &...kw_args)
-        : taylor_adaptive(std::move(sys), std::vector<T>(state), kw_args...)
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
+    explicit taylor_adaptive(std::vector<std::pair<expression, expression>> sys, const KwArgs &...kw_args)
+        : taylor_adaptive(std::move(sys), std::vector<T>{}, kw_args...)
     {
     }
     template <typename... KwArgs>
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
     explicit taylor_adaptive(var_ode_sys sys, std::vector<T> state, const KwArgs &...kw_args)
         : taylor_adaptive(private_ctor_t{}, llvm_state(kw_args...))
     {
         finalise_ctor(std::move(sys), std::move(state), kw_args...);
     }
     template <typename... KwArgs>
-    explicit taylor_adaptive(var_ode_sys sys, std::initializer_list<T> state, const KwArgs &...kw_args)
-        : taylor_adaptive(std::move(sys), std::vector<T>(state), kw_args...)
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
+    explicit taylor_adaptive(var_ode_sys sys, const KwArgs &...kw_args)
+        : taylor_adaptive(std::move(sys), std::vector<T>{}, kw_args...)
     {
     }
 
@@ -694,6 +698,16 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada
     }
 };
 
+// Deduction guides to enable CTAD when the initial state is passed via std::initializer_list.
+template <typename T, typename... KwArgs>
+    requires(!igor::has_unnamed_arguments<KwArgs...>())
+explicit taylor_adaptive(std::vector<std::pair<expression, expression>>, std::initializer_list<T>,
+                         const KwArgs &...) -> taylor_adaptive<T>;
+
+template <typename T, typename... KwArgs>
+    requires(!igor::has_unnamed_arguments<KwArgs...>())
+explicit taylor_adaptive(var_ode_sys, std::initializer_list<T>, const KwArgs &...) -> taylor_adaptive<T>;
+
 // Prevent implicit instantiations.
 // NOLINTBEGIN
 #define HEYOKA_TAYLOR_ADAPTIVE_EXTERN_INST(F)                                                                          \
@@ -932,29 +946,32 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch
     taylor_adaptive_batch();
 
     template <typename... KwArgs>
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
     explicit taylor_adaptive_batch(std::vector<std::pair<expression, expression>> sys, std::vector<T> state,
                                    std::uint32_t batch_size, const KwArgs &...kw_args)
         : taylor_adaptive_batch(private_ctor_t{}, llvm_state(kw_args...))
     {
         finalise_ctor(std::move(sys), std::move(state), batch_size, kw_args...);
     }
     template <typename... KwArgs>
-    explicit taylor_adaptive_batch(std::vector<std::pair<expression, expression>> sys, std::initializer_list<T> state,
-                                   std::uint32_t batch_size, const KwArgs &...kw_args)
-        : taylor_adaptive_batch(std::move(sys), std::vector<T>(state), batch_size, kw_args...)
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
+    explicit taylor_adaptive_batch(std::vector<std::pair<expression, expression>> sys, std::uint32_t batch_size,
+                                   const KwArgs &...kw_args)
+        : taylor_adaptive_batch(std::move(sys), std::vector<T>{}, batch_size, kw_args...)
     {
     }
     template <typename... KwArgs>
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
     explicit taylor_adaptive_batch(var_ode_sys sys, std::vector<T> state, std::uint32_t batch_size,
                                    const KwArgs &...kw_args)
         : taylor_adaptive_batch(private_ctor_t{}, llvm_state(kw_args...))
     {
         finalise_ctor(std::move(sys), std::move(state), batch_size, kw_args...);
     }
     template <typename... KwArgs>
-    explicit taylor_adaptive_batch(var_ode_sys sys, std::initializer_list<T> state, std::uint32_t batch_size,
-                                   const KwArgs &...kw_args)
-        : taylor_adaptive_batch(std::move(sys), std::vector<T>(state), batch_size, kw_args...)
+        requires(!igor::has_unnamed_arguments<KwArgs...>())
+    explicit taylor_adaptive_batch(var_ode_sys sys, std::uint32_t batch_size, const KwArgs &...kw_args)
+        : taylor_adaptive_batch(std::move(sys), std::vector<T>{}, batch_size, kw_args...)
     {
     }
 
@@ -1120,6 +1137,17 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch
     [[nodiscard]] const std::vector<std::tuple<taylor_outcome, T, T, std::size_t>> &get_propagate_res() const;
 };
 
+// Deduction guides to enable CTAD when the initial state is passed via std::initializer_list.
+template <typename T, typename... KwArgs>
+    requires(!igor::has_unnamed_arguments<KwArgs...>())
+explicit taylor_adaptive_batch(std::vector<std::pair<expression, expression>>, std::initializer_list<T>, std::uint32_t,
+                               const KwArgs &...) -> taylor_adaptive_batch<T>;
+
+template <typename T, typename... KwArgs>
+    requires(!igor::has_unnamed_arguments<KwArgs...>())
+explicit taylor_adaptive_batch(var_ode_sys, std::initializer_list<T>, std::uint32_t,
+                               const KwArgs &...) -> taylor_adaptive_batch<T>;
+
 // Prevent implicit instantiations.
 #define HEYOKA_TAYLOR_ADAPTIVE_BATCH_EXTERN_INST(F) extern template class taylor_adaptive_batch<F>;
 

src/taylor_adaptive.cpp

@@ -11,6 +11,7 @@
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <concepts>
 #include <cstddef>
 #include <cstdint>
 #include <initializer_list>
@@ -223,55 +224,86 @@ void taylor_adaptive<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> state,
     // just re-validate for peace of mind.
     validate_ode_sys(sys, tes, ntes);
 
-    // Run an immediate check on state. This is a bit redundant with other checks
-    // later (e.g., state.size() must be consistent with the ODE definition, which in
-    // turn cannot consist of zero equations), but it's handy to do it here so that,
-    // e.g., we can immediately infer the precision if T == mppp::real.
-    if (state.empty()) {
-        throw std::invalid_argument("Cannot initialise an adaptive integrator with an empty state vector");
+#if defined(HEYOKA_HAVE_REAL)
+
+    // Setup the m_prec data member for multiprecision integrations.
+    if constexpr (std::same_as<T, mppp::real>) {
+        if (!prec && state.empty()) [[unlikely]] {
+            throw std::invalid_argument("A multiprecision integrator can be initialised with an empty state vector "
+                                        "only if the desired precision is explicitly passed to the constructor (we "
+                                        "cannot deduce the desired precision from an empty state vector)");
+        }
+
+        // The precision is either passed by the user or automatically inferred from the (nonempty) state vector.
+        this->m_prec = prec ? boost::numeric_cast<mpfr_prec_t>(*prec) : state[0].get_prec();
     }
 
-    // Assign the state.
-    m_state = std::move(state);
+#endif
 
-#if defined(HEYOKA_HAVE_REAL)
+    // Fetch the original number of equations/state variables.
+    const auto n_orig_sv = is_variational ? std::get<1>(vsys).get_n_orig_sv()
+                                          : boost::numeric_cast<std::uint32_t>(std::get<0>(vsys).size());
+    // NOTE: this is ensured by validate_ode_sys().
+    assert(n_orig_sv != 0u);
 
-    // Setup the precision: it is either passed by the user
-    // or automatically inferred from the state vector.
-    // NOTE: this must be done early so that the precision of the integrator
-    // is available for other checks later.
-    if constexpr (std::is_same_v<T, mppp::real>) {
-        this->m_prec = prec ? boost::numeric_cast<mpfr_prec_t>(*prec) : m_state[0].get_prec();
-
-        if (prec) {
-            // If the user explicitly specifies a precision, enforce that precision
-            // on the state vector.
-            // NOTE: if the user specifies an invalid precision, we are sure
-            // here that an exception will be thrown: m_state is not empty
-            // and prec_round() will check the input precision value.
-            for (auto &val : m_state) {
-                // NOTE: use directly this->m_prec in order to avoid
-                // triggering an assertion in get_prec() if a bogus
-                // prec value was provided by the user.
-                val.prec_round(this->m_prec);
-            }
+    // Zero init the state vector, if empty.
+    if (state.empty()) {
+        // NOTE: we will perform further initialisation for the variational quantities
+        // at a later stage, if needed.
+
+#if defined(HEYOKA_HAVE_REAL)
+        if constexpr (std::same_as<T, mppp::real>) {
+            state.resize(boost::numeric_cast<decltype(state.size())>(n_orig_sv),
+                         mppp::real{mppp::real_kind::zero, this->get_prec()});
         } else {
-            // If the precision is automatically deduced, ensure that
-            // the same precision is used for all initial values.
-            // NOTE: the automatically-deduced precision will be a valid one,
-            // as it is taken from a valid mppp::real (i.e., m_state[0]).
-            if (std::any_of(m_state.begin() + 1, m_state.end(),
-                            [this](const auto &val) { return val.get_prec() != this->get_prec(); })) {
-                throw std::invalid_argument(
-                    fmt::format("The precision deduced automatically from the initial state vector in a multiprecision "
-                                "adaptive Taylor integrator is {}, but values with different precision(s) were "
-                                "detected in the state vector",
-                                this->get_prec()));
+#endif
+            state.resize(boost::numeric_cast<decltype(state.size())>(n_orig_sv), static_cast<T>(0));
+#if defined(HEYOKA_HAVE_REAL)
+        }
+#endif
+    } else {
+#if defined(HEYOKA_HAVE_REAL)
+
+        if constexpr (std::same_as<T, mppp::real>) {
+            // The user passed in a non-empty multiprecision state. We need
+            // to either enforce the desired precision, or to check that the
+            // values in the state have all the same precision.
+            if (prec) {
+                // If the user explicitly specifies a precision, enforce that precision
+                // on the state vector.
+                // NOTE: if the user specifies an invalid precision, we are sure
+                // here that an exception will be thrown: state is not empty
+                // and prec_round() will check the input precision value.
+                for (auto &val : state) {
+                    // NOTE: use directly this->m_prec in order to avoid
+                    // triggering an assertion in get_prec() if a bogus
+                    // prec value was provided by the user.
+                    val.prec_round(this->m_prec);
+                }
+            } else {
+                // If the precision is automatically deduced, ensure that
+                // the same precision is used for all initial values.
+                // NOTE: the automatically-deduced precision will be a valid one,
+                // as it is taken from a valid mppp::real (i.e., state[0]).
+                if (std::any_of(state.begin() + 1, state.end(),
+                                [this](const auto &val) { return val.get_prec() != this->get_prec(); })) {
+                    throw std::invalid_argument(fmt::format(
+                        "The precision deduced automatically from the initial state vector in a multiprecision "
+                        "adaptive Taylor integrator is {}, but one or more values with a different precision were "
+                        "detected in the state vector",
+                        this->get_prec()));
+                }
             }
         }
-    }
 
 #endif
+    }
+
+    // NOTE: at this point, state must be a non-empty vector.
+    assert(!state.empty());
+
+    // Assign the state.
+    m_state = std::move(state);
 
     // Check the input state size.
     // NOTE: keep track of whether or not we need to automatically setup the initial
@@ -281,9 +313,6 @@ void taylor_adaptive<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> state,
     bool auto_ic_setup = false;
     if (m_state.size() != sys.size()) {
         if (is_variational) {
-            // Fetch the original number of equations/state variables.
-            const auto n_orig_sv = std::get<1>(vsys).get_n_orig_sv();
-
             if (m_state.size() == n_orig_sv) [[likely]] {
                 // Automatic setup of the initial conditions for the derivatives wrt
                 // variables and parameters.

src/taylor_adaptive_batch.cpp

@@ -145,35 +145,47 @@ void taylor_adaptive_batch<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> sta
 
     // Init the data members.
     m_batch_size = batch_size;
-    m_state = std::move(state);
     m_time_hi = std::move(time);
     m_time_lo.resize(m_time_hi.size());
     m_pars = std::move(pars);
     m_high_accuracy = high_accuracy;
     m_compact_mode = compact_mode;
 
     // Check several input params.
-    if (m_batch_size == 0u) {
+    if (m_batch_size == 0u) [[unlikely]] {
         throw std::invalid_argument("The batch size in an adaptive Taylor integrator cannot be zero");
     }
 
-    if (m_state.size() % m_batch_size != 0u) {
+    if (state.size() % m_batch_size != 0u) [[unlikely]] {
         throw std::invalid_argument(
             fmt::format("Invalid size detected in the initialization of an adaptive Taylor "
                         "integrator: the state vector has a size of {}, which is not a multiple of the batch size ({})",
-                        m_state.size(), m_batch_size));
+                        state.size(), m_batch_size)); // LCOV_EXCL_LINE
+    }
+
+    // Fetch the original number of equations/state variables.
+    const auto n_orig_sv = is_variational ? std::get<1>(vsys).get_n_orig_sv()
+                                          : boost::numeric_cast<std::uint32_t>(std::get<0>(vsys).size());
+    // NOTE: this is ensured by validate_ode_sys().
+    assert(n_orig_sv != 0u);
+
+    // Zero init the state vector, if empty.
+    if (state.empty()) {
+        // NOTE: we will perform further initialisation for the variational quantities
+        // at a later stage, if needed.
+        state.resize(boost::safe_numerics::safe<decltype(state.size())>(n_orig_sv) * m_batch_size, static_cast<T>(0));
     }
 
+    // Assign the state.
+    m_state = std::move(state);
+
     // NOTE: keep track of whether or not we need to automatically setup the initial
     // conditions in a variational integrator. This is needed because we need
     // to delay the automatic ic setup for the derivatives wrt the initial time until
     // after we have correctly set up state, pars and time in the integrator.
     bool auto_ic_setup = false;
     if (m_state.size() / m_batch_size != sys.size()) {
         if (is_variational) {
-            // Fetch the original number of equations/state variables.
-            const auto n_orig_sv = std::get<1>(vsys).get_n_orig_sv();
-
             if (m_state.size() / m_batch_size == n_orig_sv) [[likely]] {
                 // Automatic setup of the initial conditions for the derivatives wrt
                 // variables and parameters.

test/taylor_adaptive.cpp

@@ -51,6 +51,7 @@
 #include <heyoka/math/sum.hpp>
 #include <heyoka/math/time.hpp>
 #include <heyoka/model/nbody.hpp>
+#include <heyoka/model/pendulum.hpp>
 #include <heyoka/number.hpp>
 #include <heyoka/s11n.hpp>
 #include <heyoka/step_callback.hpp>
@@ -2191,18 +2192,6 @@ TEST_CASE("ctad")
 #endif
 }
 
-// Test to trigger the empty state check early in the
-// construction process.
-TEST_CASE("early empty state check")
-{
-    using Catch::Matchers::Message;
-
-    auto [x, v] = make_vars("x", "v");
-
-    REQUIRE_THROWS_MATCHES(taylor_adaptive({prime(x) = v, prime(v) = -x}, std::vector<double>{}), std::invalid_argument,
-                           Message("Cannot initialise an adaptive integrator with an empty state vector"));
-}
-
 // Test to check that event callbacks which alter the time coordinate
 // result in an exception being thrown.
 TEST_CASE("event cb time")
@@ -2537,3 +2526,18 @@ TEST_CASE("invalid initial state")
         Message("Inconsistent sizes detected in the initialization of an adaptive Taylor "
                 "integrator: the state vector has a dimension of 1, while the number of equations is 2"));
 }
+
+TEST_CASE("empty init state")
+{
+    const auto dyn = model::pendulum();
+
+    {
+        auto ta = taylor_adaptive<double>{dyn};
+        REQUIRE(ta.get_state() == std::vector{0., 0.});
+    }
+
+    {
+        auto ta = taylor_adaptive<double>{var_ode_sys(dyn, var_args::vars)};
+        REQUIRE(ta.get_state() == std::vector{0.0, 0.0, 1.0, 0.0, 0.0, 1.0});
+    }
+}