Merge pull request #451 from bluescarni/pr/new_par_init

Forbid integrator init with a pars array of the wrong size
bluescarni · Sep 10, 2024 · a3c2022 · a3c2022
2 parents 822b0cd + ccbf41f
commit a3c2022
Show file tree

Hide file tree

Showing 7 changed files with 169 additions and 59 deletions.
diff --git 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
@@ -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 <https://github.com/bluescarni/heyoka/pull/449 >`__).
+  (`#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>`__,
@@ -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 <https://github.com/bluescarni/heyoka/pull/451>`__).
+  This is a :ref:`breaking change <bchanges_6_0_0>`.
 - **BREAKING**: the minimum required LLVM version is now 15
   (`#444 <https://github.com/bluescarni/heyoka/pull/444>`__).
   This is a :ref:`breaking change <bchanges_6_0_0>`.

diff --git a/src/taylor_adaptive.cpp b/src/taylor_adaptive.cpp
@@ -257,7 +257,7 @@ void taylor_adaptive<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> state,
                          mppp::real{mppp::real_kind::zero, this->get_prec()});
         } else {
 #endif
-            state.resize(boost::numeric_cast<decltype(state.size())>(n_orig_sv), static_cast<T>(0));
+            state.resize(boost::numeric_cast<decltype(state.size())>(n_orig_sv));
 #if defined(HEYOKA_HAVE_REAL)
         }
 #endif
@@ -358,13 +358,40 @@ void taylor_adaptive<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> 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<T, mppp::real>) {
+            // For mppp::real, ensure that the parameter
+            // values all have the inferred precision.
+            pars.resize(boost::numeric_cast<decltype(pars.size())>(tot_n_pars),
+                        mppp::real{mppp::real_kind::zero, this->get_prec()});
+        } else {
+#endif
+            pars.resize(boost::numeric_cast<decltype(pars.size())>(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<T, mppp::real>) {
+    // 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<T, mppp::real>) {
         for (auto &val : m_pars) {
             val.prec_round(this->get_prec());
         }
@@ -409,9 +436,6 @@ void taylor_adaptive<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> state,
     // Store the dimension of the system.
     m_dim = boost::numeric_cast<std::uint32_t>(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<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> 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<T, mppp::real>) {
-            // For mppp::real, ensure that the appended parameter
-            // values all have the inferred precision.
-            m_pars.resize(boost::numeric_cast<decltype(m_pars.size())>(tot_n_pars),
-                          mppp::real{mppp::real_kind::zero, this->get_prec()});
-        } else {
-#endif
-            m_pars.resize(boost::numeric_cast<decltype(m_pars.size())>(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
@@ -147,7 +147,6 @@ void taylor_adaptive_batch<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> 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<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> 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<decltype(state.size())>(n_orig_sv) * m_batch_size, static_cast<T>(0));
+        state.resize(boost::safe_numerics::safe<decltype(state.size())>(n_orig_sv) * m_batch_size);
     }
 
     // Assign the state.
@@ -252,6 +251,22 @@ void taylor_adaptive_batch<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> sta
     }
     // LCOV_EXCL_STOP
 
+    // Check/setup pars.
+    using su32_t = boost::safe_numerics::safe<std::uint32_t>;
+    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<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> sta
                                                high_accuracy, compact_mode, parallel_mode, m_order);
     }
 
-    // Fix m_pars' size, if necessary.
-    using su32_t = boost::safe_numerics::safe<std::uint32_t>;
-    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
@@ -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<double>(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<double>{{prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                       {0.05, 0.025},
+                                       kw::pars = std::vector{1.},
+                                       kw::t_events = {t_event<double>(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<double>{
+        {prime(x) = v, prime(v) = -9.8 * sin(x)}, {0.05, 0.025}, kw::t_events = {t_event<double>(v - par[3])}};
+    REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0.});
+
+    ta = taylor_adaptive<double>{{prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                 {0.05, 0.025},
+                                 kw::pars = std::vector<double>{},
+                                 kw::t_events = {t_event<double>(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
@@ -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<double>{{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<double>(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<double>{{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<double>{{prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                            {0.05, 0.06, 0.025, 0.026},
+                                            2u,
+                                            kw::t_events = {t_event_batch<double>(v - par[3])}};
+    REQUIRE(ta.get_pars() == std::vector{0., 0., 0., 0., 0., 0., 0., 0.});
+
+    ta = taylor_adaptive_batch<double>{{prime(x) = v, prime(v) = -9.8 * sin(x)},
+                                       {0.05, 0.06, 0.025, 0.026},
+                                       2u,
+                                       kw::pars = std::vector<double>{},
+                                       kw::t_events = {t_event_batch<double>(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<fp_t>{{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
@@ -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<mppp::real>({{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<mppp::real>(
+            {{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<mppp::real>{dyn, kw::prec = prec};
+        REQUIRE(ta.get_pars() == std::vector<mppp::real>{0., 0., 0., 0.});
+        REQUIRE(std::ranges::all_of(ta.get_pars(), [&](const auto &val) { return val.get_prec() == prec; }));
+    }
+    {
+        auto ta = taylor_adaptive<mppp::real>{dyn, kw::prec = prec, kw::pars = std::vector<mppp::real>{}};
+        REQUIRE(ta.get_pars() == std::vector<mppp::real>{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<mppp::real>{var_ode_sys(dyn, var_args::vars), kw::prec = prec,
+                                              kw::pars = std::vector<mppp::real>(4, mppp::real{})};
+        REQUIRE(ta.get_state() == std::vector<mppp::real>{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; }));
+    }
+}