diff --git a/doc/breaking_changes.rst b/doc/breaking_changes.rst index 94f7451c5..056104aeb 100644 --- a/doc/breaking_changes.rst +++ b/doc/breaking_changes.rst @@ -8,6 +8,19 @@ Breaking changes 6.0.0 ----- +API/behaviour changes +~~~~~~~~~~~~~~~~~~~~~ + +In heyoka 6.0.0, the array of parameter values passed to the constructor +of a Taylor integrator must either be empty (in which case the parameter +values will be zero-inited), or it must have the correct size. + +In previous versions, heyoka would pad with zeroes the array of parameter values +if its size was less than the correct one. + +General +~~~~~~~ + - Support for LLVM<15 has been dropped. .. _bchanges_4_0_0: diff --git a/doc/changelog.rst b/doc/changelog.rst index 98534d8e1..316d27b2a 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -10,7 +10,7 @@ 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 `__). + (`#449 `__). - Implement parallel compilation for Taylor integrators and compiled functions (`#446 `__, @@ -25,6 +25,12 @@ New Changes ~~~~~~~ +- **BREAKING**: the array of parameter values passed to the + constructor of a Taylor integrator must now either be empty + (in which case the parameter values will be zero-inited), + or have the correct size + (`#451 `__). + This is a :ref:`breaking change `. - **BREAKING**: the minimum required LLVM version is now 15 (`#444 `__). This is a :ref:`breaking change `. diff --git a/src/taylor_adaptive.cpp b/src/taylor_adaptive.cpp index 093c49c2e..02e9ec15e 100644 --- a/src/taylor_adaptive.cpp +++ b/src/taylor_adaptive.cpp @@ -257,7 +257,7 @@ void taylor_adaptive::finalise_ctor_impl(sys_t vsys, std::vector state, mppp::real{mppp::real_kind::zero, this->get_prec()}); } else { #endif - state.resize(boost::numeric_cast(n_orig_sv), static_cast(0)); + state.resize(boost::numeric_cast(n_orig_sv)); #if defined(HEYOKA_HAVE_REAL) } #endif @@ -358,13 +358,40 @@ void taylor_adaptive::finalise_ctor_impl(sys_t vsys, std::vector state, #endif } - // Parameter values. + // Compute the total number of params, including the rhs and the event functions. + const auto tot_n_pars = detail::tot_n_pars_in_ode_sys(sys, tes, ntes); + + // Check/setup pars. + if (pars.empty()) { +#if defined(HEYOKA_HAVE_REAL) + if constexpr (std::same_as) { + // For mppp::real, ensure that the parameter + // values all have the inferred precision. + pars.resize(boost::numeric_cast(tot_n_pars), + mppp::real{mppp::real_kind::zero, this->get_prec()}); + } else { +#endif + pars.resize(boost::numeric_cast(tot_n_pars)); +#if defined(HEYOKA_HAVE_REAL) + } +#endif + } else if (pars.size() != tot_n_pars) [[unlikely]] { + throw std::invalid_argument(fmt::format( + "Invalid number of parameter values passed to the constructor of an adaptive " + "Taylor integrator: {} parameter value(s) were passed, but the ODE system contains {} parameter(s)", + pars.size(), tot_n_pars)); + } + + // Assign the parameter values. m_pars = std::move(pars); #if defined(HEYOKA_HAVE_REAL) // Enforce the correct precision for mppp::real. - if constexpr (std::is_same_v) { + // NOTE: this is redundant in case pars was originally empty, + // but let us keep it like this in order not to complicate further + // the logic. + if constexpr (std::same_as) { for (auto &val : m_pars) { val.prec_round(this->get_prec()); } @@ -409,9 +436,6 @@ void taylor_adaptive::finalise_ctor_impl(sys_t vsys, std::vector state, // Store the dimension of the system. m_dim = boost::numeric_cast(sys.size()); - // Compute the total number of params, including the rhs and the event functions. - const auto tot_n_pars = detail::tot_n_pars_in_ode_sys(sys, tes, ntes); - // Do we have events? const auto with_events = !tes.empty() || !ntes.empty(); @@ -449,27 +473,6 @@ void taylor_adaptive::finalise_ctor_impl(sys_t vsys, std::vector state, compact_mode, parallel_mode, m_order); } - // Fix m_pars' size, if necessary. - if (m_pars.size() < tot_n_pars) { -#if defined(HEYOKA_HAVE_REAL) - if constexpr (std::is_same_v) { - // For mppp::real, ensure that the appended parameter - // values all have the inferred precision. - m_pars.resize(boost::numeric_cast(tot_n_pars), - mppp::real{mppp::real_kind::zero, this->get_prec()}); - } else { -#endif - m_pars.resize(boost::numeric_cast(tot_n_pars)); -#if defined(HEYOKA_HAVE_REAL) - } -#endif - } else if (m_pars.size() > tot_n_pars) { - throw std::invalid_argument(fmt::format( - "Excessive number of parameter values passed to the constructor of an adaptive " - "Taylor integrator: {} parameter value(s) were passed, but the ODE system contains only {} parameter(s)", - m_pars.size(), tot_n_pars)); - } - // Log runtimes in trace mode. spdlog::stopwatch sw; diff --git a/src/taylor_adaptive_batch.cpp b/src/taylor_adaptive_batch.cpp index 55fc67501..eac130da4 100644 --- a/src/taylor_adaptive_batch.cpp +++ b/src/taylor_adaptive_batch.cpp @@ -147,7 +147,6 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector sta m_batch_size = batch_size; 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; @@ -173,7 +172,7 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector sta 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(n_orig_sv) * m_batch_size, static_cast(0)); + state.resize(boost::safe_numerics::safe(n_orig_sv) * m_batch_size); } // Assign the state. @@ -252,6 +251,22 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector sta } // LCOV_EXCL_STOP + // Check/setup pars. + using su32_t = boost::safe_numerics::safe; + const auto pars_req_size = su32_t(tot_n_pars) * m_batch_size; + if (pars.empty()) { + pars.resize(pars_req_size); + } else if (pars.size() != pars_req_size) [[unlikely]] { + throw std::invalid_argument( + fmt::format("Invalid number of parameter values passed to the constructor of an adaptive " + "Taylor integrator in batch mode: {} parameter value(s) were passed, but the ODE " + "system contains {} parameter(s) (in batches of {})", + pars.size(), tot_n_pars, m_batch_size)); + } + + // Assign pars. + m_pars = std::move(pars); + // Do we have events? const auto with_events = !tes.empty() || !ntes.empty(); @@ -288,19 +303,6 @@ void taylor_adaptive_batch::finalise_ctor_impl(sys_t vsys, std::vector sta high_accuracy, compact_mode, parallel_mode, m_order); } - // Fix m_pars' size, if necessary. - using su32_t = boost::safe_numerics::safe; - const auto pars_req_size = su32_t(tot_n_pars) * m_batch_size; - if (m_pars.size() < pars_req_size) { - m_pars.resize(pars_req_size); - } else if (m_pars.size() > pars_req_size) { - throw std::invalid_argument( - fmt::format("Excessive number of parameter values passed to the constructor of an adaptive " - "Taylor integrator in batch mode: {} parameter value(s) were passed, but the ODE " - "system contains only {} parameter(s) (in batches of {})", - m_pars.size(), tot_n_pars, m_batch_size)); - } - // Log runtimes in trace mode. spdlog::stopwatch sw; diff --git a/test/taylor_adaptive.cpp b/test/taylor_adaptive.cpp index 94cd1b89a..5dad867aa 100644 --- a/test/taylor_adaptive.cpp +++ b/test/taylor_adaptive.cpp @@ -1562,9 +1562,9 @@ TEST_CASE("cb interrupt") } } -// Test excessive number of params provided upon construction +// Test wrong number of params provided upon construction // of the integrator. -TEST_CASE("param too many") +TEST_CASE("param wrong number") { using Catch::Matchers::Message; @@ -1576,9 +1576,36 @@ TEST_CASE("param too many") kw::pars = std::vector{1., 2.}, kw::t_events = {t_event(v - par[0])}}), std::invalid_argument, - Message( - "Excessive number of parameter values passed to the constructor of an adaptive " - "Taylor integrator: 2 parameter value(s) were passed, but the ODE system contains only 1 parameter(s)")); + Message("Invalid number of parameter values passed to the constructor of an adaptive " + "Taylor integrator: 2 parameter value(s) were passed, but the ODE system contains 1 parameter(s)")); + + REQUIRE_THROWS_MATCHES( + (void)(taylor_adaptive{{prime(x) = v, prime(v) = -9.8 * sin(x)}, + {0.05, 0.025}, + kw::pars = std::vector{1.}, + kw::t_events = {t_event(v - par[1])}}), + std::invalid_argument, + Message("Invalid number of parameter values passed to the constructor of an adaptive " + "Taylor integrator: 1 parameter value(s) were passed, but the ODE system contains 2 parameter(s)")); +} + +// Test to check that the parameter values are zeroed out +// when an empty pars array is passed on construction. +TEST_CASE("param auto setup") +{ + using Catch::Matchers::Message; + + auto [x, v] = make_vars("x", "v"); + + auto ta = taylor_adaptive{ + {prime(x) = v, prime(v) = -9.8 * sin(x)}, {0.05, 0.025}, kw::t_events = {t_event(v - par[3])}}; + REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0.}); + + ta = taylor_adaptive{{prime(x) = v, prime(v) = -9.8 * sin(x)}, + {0.05, 0.025}, + kw::pars = std::vector{}, + kw::t_events = {t_event(v - par[3])}}; + REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0.}); } // Test case for bug: parameters in event equations are ignored diff --git a/test/taylor_adaptive_batch.cpp b/test/taylor_adaptive_batch.cpp index c74b907a9..84e918b5a 100644 --- a/test/taylor_adaptive_batch.cpp +++ b/test/taylor_adaptive_batch.cpp @@ -950,9 +950,9 @@ TEST_CASE("cb interrupt") } } -// Test excessive number of params provided upon construction +// Test wrong number of params provided upon construction // of the integrator. -TEST_CASE("param too many") +TEST_CASE("param wrong number") { using Catch::Matchers::Message; @@ -963,9 +963,9 @@ TEST_CASE("param too many") 2u, kw::pars = std::vector{1., 2., 3.}}), std::invalid_argument, - Message("Excessive number of parameter values passed to the constructor of an adaptive " + Message("Invalid number of parameter values passed to the constructor of an adaptive " "Taylor integrator in batch mode: 3 parameter value(s) were passed, but the ODE " - "system contains only 1 parameter(s) " + "system contains 1 parameter(s) " "(in batches of 2)")); REQUIRE_THROWS_MATCHES((void)(taylor_adaptive_batch{{prime(x) = v, prime(v) = -9.8 * sin(x)}, @@ -974,10 +974,42 @@ TEST_CASE("param too many") kw::pars = std::vector{1., 2., 3.}, kw::t_events = {t_event_batch(v - par[0])}}), std::invalid_argument, - Message("Excessive number of parameter values passed to the constructor of an adaptive " + Message("Invalid number of parameter values passed to the constructor of an adaptive " "Taylor integrator in batch mode: 3 parameter value(s) were passed, but the ODE " - "system contains only 1 parameter(s) " + "system contains 1 parameter(s) " "(in batches of 2)")); + + REQUIRE_THROWS_MATCHES((void)(taylor_adaptive_batch{{prime(x) = v + par[3], prime(v) = -9.8 * sin(x)}, + {0.05, 0.06, 0.025, 0.026}, + 2u, + kw::pars = std::vector{1., 2., 3.}}), + std::invalid_argument, + Message("Invalid number of parameter values passed to the constructor of an adaptive " + "Taylor integrator in batch mode: 3 parameter value(s) were passed, but the ODE " + "system contains 4 parameter(s) " + "(in batches of 2)")); +} + +// Test to check that the parameter values are zeroed out +// when an empty pars array is passed on construction. +TEST_CASE("param auto setup") +{ + using Catch::Matchers::Message; + + auto [x, v] = make_vars("x", "v"); + + auto ta = taylor_adaptive_batch{{prime(x) = v, prime(v) = -9.8 * sin(x)}, + {0.05, 0.06, 0.025, 0.026}, + 2u, + kw::t_events = {t_event_batch(v - par[3])}}; + REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0., 0., 0., 0., 0.}); + + ta = taylor_adaptive_batch{{prime(x) = v, prime(v) = -9.8 * sin(x)}, + {0.05, 0.06, 0.025, 0.026}, + 2u, + kw::pars = std::vector{}, + kw::t_events = {t_event_batch(v - par[3])}}; + REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0., 0., 0., 0., 0.}); } // Test case for bug: parameters in event equations are ignored @@ -1251,7 +1283,7 @@ TEST_CASE("stream output") auto ta = taylor_adaptive_batch{{prime(x) = v - par[1], prime(v) = -9.8 * sin(x + par[0])}, {fp_t(0.), fp_t(0.01), fp_t(0.5), fp_t(0.51)}, 2u, - kw::pars = std::vector{fp_t(-1e-4), fp_t(-1.1e-4)}}; + kw::pars = std::vector{fp_t(-1e-4), fp_t(-1.1e-4), fp_t(0), fp_t(0)}}; std::ostringstream oss; diff --git a/test/taylor_adaptive_mp.cpp b/test/taylor_adaptive_mp.cpp index 655a3962c..f9171d43a 100644 --- a/test/taylor_adaptive_mp.cpp +++ b/test/taylor_adaptive_mp.cpp @@ -161,10 +161,10 @@ TEST_CASE("ctors prec") [prec](const auto &val) { return val.get_prec() == prec + 1; })); // Try with several different precisions for the data. - ta = taylor_adaptive({{x, x + par[0]}, {y, y + par[1]}}, - {mppp::real{1, prec}, mppp::real{1, prec - 1}}, kw::compact_mode = cm, - kw::opt_level = 0u, kw::prec = prec + 1, kw::pars = {mppp::real{-1, prec - 2}}, - kw::time = mppp::real{0, prec + 3}, kw::tol = mppp::real{1e-1, prec + 4}); + ta = taylor_adaptive( + {{x, x + par[0]}, {y, y + par[1]}}, {mppp::real{1, prec}, mppp::real{1, prec - 1}}, kw::compact_mode = cm, + kw::opt_level = 0u, kw::prec = prec + 1, kw::pars = {mppp::real{-1, prec - 2}, mppp::real{-1, prec - 3}}, + kw::time = mppp::real{0, prec + 3}, kw::tol = mppp::real{1e-1, prec + 4}); REQUIRE(ta.get_tol().get_prec() == prec + 1); REQUIRE(ta.get_dtime().first.get_prec() == prec + 1); REQUIRE(ta.get_dtime().second.get_prec() == prec + 1); @@ -1178,3 +1178,30 @@ TEST_CASE("empty init state") REQUIRE(std::ranges::all_of(ta.get_state(), [&](const auto &val) { return val.get_prec() == prec; })); } } + +TEST_CASE("pars handling") +{ + const auto dyn = model::pendulum(kw::gconst = par[3]); + + const auto prec = 23; + + // Empty pars array. + { + auto ta = taylor_adaptive{dyn, kw::prec = prec}; + REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0.}); + REQUIRE(std::ranges::all_of(ta.get_pars(), [&](const auto &val) { return val.get_prec() == prec; })); + } + { + auto ta = taylor_adaptive{dyn, kw::prec = prec, kw::pars = std::vector{}}; + REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0.}); + REQUIRE(std::ranges::all_of(ta.get_pars(), [&](const auto &val) { return val.get_prec() == prec; })); + } + + // Pars array with correct size but wrong precision. + { + auto ta = taylor_adaptive{var_ode_sys(dyn, var_args::vars), kw::prec = prec, + kw::pars = std::vector(4, mppp::real{})}; + REQUIRE(ta.get_state() == std::vector{0.0, 0.0, 1.0, 0.0, 0.0, 1.0}); + REQUIRE(std::ranges::all_of(ta.get_state(), [&](const auto &val) { return val.get_prec() == prec; })); + } +}