Initial work towards allowing empty initial state vectors in the inte…

…grators.

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,15 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada
     }
 };
 
+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,31 +945,21 @@ 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...)
-    {
-    }
-    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...)
-    {
-    }
 
     taylor_adaptive_batch(const taylor_adaptive_batch &);
     taylor_adaptive_batch(taylor_adaptive_batch &&) noexcept;
@@ -1120,6 +1123,16 @@ 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;
 };
 
+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,84 @@ 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.
+            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 +311,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,7 +145,6 @@ 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);
@@ -157,23 +156,36 @@ void taylor_adaptive_batch<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> sta
         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) {
         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));
+    }
+
+    // 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::numeric_cast<decltype(state.size())>(n_orig_sv), 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

@@ -2191,18 +2191,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")

test/taylor_adaptive_mp.cpp

@@ -112,7 +112,7 @@ TEST_CASE("ctors prec")
                                         kw::opt_level = 0u),
             std::invalid_argument,
             Message("The precision deduced automatically from the initial state vector in a multiprecision "
-                    "adaptive Taylor integrator is 30, but values with different precision(s) were "
+                    "adaptive Taylor integrator is 30, but one or more values with a different precision were "
                     "detected in the state vector"));
 
         // Check that wrong precision in the pars if automatically corrected.
@@ -126,6 +126,14 @@ TEST_CASE("ctors prec")
         ta = taylor_adaptive<mppp::real>({{x, x + par[0]}}, {mppp::real{1, prec}}, kw::compact_mode = cm,
                                          kw::opt_level = 0u, kw::prec = prec + 1);
 
+        // Check that attempting to init with an empty state vector and no prec throws.
+        REQUIRE_THROWS_MATCHES(taylor_adaptive<mppp::real>({{x, x + par[0]}, {y, y + par[1]}}, {},
+                                                           kw::compact_mode = cm, kw::opt_level = 0u),
+                               std::invalid_argument,
+                               Message("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)"));
+
         // Check that time and state are automatically adjusted.
         REQUIRE(ta.get_dtime().first.get_prec() == prec + 1);
         REQUIRE(ta.get_dtime().second.get_prec() == prec + 1);
@@ -134,6 +142,12 @@ TEST_CASE("ctors prec")
         REQUIRE(std::all_of(ta.get_pars().begin(), ta.get_pars().end(),
                             [prec](const auto &val) { return val.get_prec() == prec + 1; }));
 
+        // Check init with an empty state vector and explicit precision.
+        ta = taylor_adaptive<mppp::real>({{x, x + par[0]}, {y, y + par[1]}}, {}, kw::compact_mode = cm,
+                                         kw::opt_level = 0u, kw::prec = 23);
+        REQUIRE(std::ranges::all_of(ta.get_state(), [](const auto &x) { return x.get_prec() == 23; }));
+        REQUIRE(std::ranges::all_of(ta.get_pars(), [](const auto &x) { return x.get_prec() == 23; }));
+
         // Check that it does not matter if the state vector has different precisions.
         ta = taylor_adaptive<mppp::real>({{x, x + par[0]}, {y, y + par[1]}},
                                          {mppp::real{1, prec}, mppp::real{1, prec - 1}}, kw::compact_mode = cm,