From cae4953e9e1f17956b0678ad35885fe8ac9b27a0 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 8 Sep 2024 11:47:56 +0200 Subject: [PATCH] Initial work towards allowing empty initial state vectors in the integrators. --- include/heyoka/taylor.hpp | 47 +++++++++----- src/taylor_adaptive.cpp | 113 +++++++++++++++++++++------------- src/taylor_adaptive_batch.cpp | 24 ++++++-- test/taylor_adaptive.cpp | 12 ---- test/taylor_adaptive_mp.cpp | 16 ++++- 5 files changed, 133 insertions(+), 79 deletions(-) diff --git a/include/heyoka/taylor.hpp b/include/heyoka/taylor.hpp index fd7f65003..3e9575298 100644 --- a/include/heyoka/taylor.hpp +++ b/include/heyoka/taylor.hpp @@ -50,6 +50,7 @@ #include #include #include +#include #include #include #include @@ -533,6 +534,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada taylor_adaptive(); template + requires(!igor::has_unnamed_arguments()) explicit taylor_adaptive(std::vector> sys, std::vector state, const KwArgs &...kw_args) : taylor_adaptive(private_ctor_t{}, llvm_state(kw_args...)) @@ -540,20 +542,22 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada finalise_ctor(std::move(sys), std::move(state), kw_args...); } template - explicit taylor_adaptive(std::vector> sys, std::initializer_list state, - const KwArgs &...kw_args) - : taylor_adaptive(std::move(sys), std::vector(state), kw_args...) + requires(!igor::has_unnamed_arguments()) + explicit taylor_adaptive(std::vector> sys, const KwArgs &...kw_args) + : taylor_adaptive(std::move(sys), std::vector{}, kw_args...) { } template + requires(!igor::has_unnamed_arguments()) explicit taylor_adaptive(var_ode_sys sys, std::vector 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 - explicit taylor_adaptive(var_ode_sys sys, std::initializer_list state, const KwArgs &...kw_args) - : taylor_adaptive(std::move(sys), std::vector(state), kw_args...) + requires(!igor::has_unnamed_arguments()) + explicit taylor_adaptive(var_ode_sys sys, const KwArgs &...kw_args) + : taylor_adaptive(std::move(sys), std::vector{}, kw_args...) { } @@ -694,6 +698,15 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada } }; +template + requires(!igor::has_unnamed_arguments()) +explicit taylor_adaptive(std::vector>, std::initializer_list, + const KwArgs &...) -> taylor_adaptive; + +template + requires(!igor::has_unnamed_arguments()) +explicit taylor_adaptive(var_ode_sys, std::initializer_list, const KwArgs &...) -> taylor_adaptive; + // Prevent implicit instantiations. // NOLINTBEGIN #define HEYOKA_TAYLOR_ADAPTIVE_EXTERN_INST(F) \ @@ -932,6 +945,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch taylor_adaptive_batch(); template + requires(!igor::has_unnamed_arguments()) explicit taylor_adaptive_batch(std::vector> sys, std::vector state, std::uint32_t batch_size, const KwArgs &...kw_args) : taylor_adaptive_batch(private_ctor_t{}, llvm_state(kw_args...)) @@ -939,24 +953,13 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch finalise_ctor(std::move(sys), std::move(state), batch_size, kw_args...); } template - explicit taylor_adaptive_batch(std::vector> sys, std::initializer_list state, - std::uint32_t batch_size, const KwArgs &...kw_args) - : taylor_adaptive_batch(std::move(sys), std::vector(state), batch_size, kw_args...) - { - } - template + requires(!igor::has_unnamed_arguments()) explicit taylor_adaptive_batch(var_ode_sys sys, std::vector 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 - explicit taylor_adaptive_batch(var_ode_sys sys, std::initializer_list state, std::uint32_t batch_size, - const KwArgs &...kw_args) - : taylor_adaptive_batch(std::move(sys), std::vector(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> &get_propagate_res() const; }; +template + requires(!igor::has_unnamed_arguments()) +explicit taylor_adaptive_batch(std::vector>, std::initializer_list, std::uint32_t, + const KwArgs &...) -> taylor_adaptive_batch; + +template + requires(!igor::has_unnamed_arguments()) +explicit taylor_adaptive_batch(var_ode_sys, std::initializer_list, std::uint32_t, + const KwArgs &...) -> taylor_adaptive_batch; + // Prevent implicit instantiations. #define HEYOKA_TAYLOR_ADAPTIVE_BATCH_EXTERN_INST(F) extern template class taylor_adaptive_batch; diff --git a/src/taylor_adaptive.cpp b/src/taylor_adaptive.cpp index 0707e69a9..51e5f3880 100644 --- a/src/taylor_adaptive.cpp +++ b/src/taylor_adaptive.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -223,55 +224,84 @@ void taylor_adaptive::finalise_ctor_impl(sys_t vsys, std::vector 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) { + 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(*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::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) { - this->m_prec = prec ? boost::numeric_cast(*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) { + state.resize(boost::numeric_cast(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(n_orig_sv), static_cast(0)); +#if defined(HEYOKA_HAVE_REAL) + } +#endif + } else { +#if defined(HEYOKA_HAVE_REAL) + + if constexpr (std::same_as) { + // 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::finalise_ctor_impl(sys_t vsys, std::vector 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. diff --git a/src/taylor_adaptive_batch.cpp b/src/taylor_adaptive_batch.cpp index 6bff47ca1..5784bff6f 100644 --- a/src/taylor_adaptive_batch.cpp +++ b/src/taylor_adaptive_batch.cpp @@ -145,7 +145,6 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector 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,13 +156,29 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector 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::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(n_orig_sv), static_cast(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 @@ -171,9 +186,6 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector sta 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. diff --git a/test/taylor_adaptive.cpp b/test/taylor_adaptive.cpp index 9aa195674..663f8dbeb 100644 --- a/test/taylor_adaptive.cpp +++ b/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{}), 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") diff --git a/test/taylor_adaptive_mp.cpp b/test/taylor_adaptive_mp.cpp index aadf9990a..b7b139204 100644 --- a/test/taylor_adaptive_mp.cpp +++ b/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({{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({{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({{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({{x, x + par[0]}, {y, y + par[1]}}, {mppp::real{1, prec}, mppp::real{1, prec - 1}}, kw::compact_mode = cm,