From c3b60c588f85d11830011d27e6071bce244b0fbf Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 4 Mar 2024 18:17:30 +0100 Subject: [PATCH 1/7] Minor internal doc bits. --- include/heyoka/func.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/include/heyoka/func.hpp b/include/heyoka/func.hpp index 2bc78c2a8..687efe309 100644 --- a/include/heyoka/func.hpp +++ b/include/heyoka/func.hpp @@ -121,9 +121,9 @@ template struct HEYOKA_DLL_PUBLIC_INLINE_CLASS func_iface_impl : public Base, tanuki::iface_impl_helper { [[nodiscard]] const std::string &get_name() const final { - // NOTE: make sure we are invoking the member functions - // from func_base (these functions could have been overriden - // in the derived class). + // NOTE: make sure we are invoking the member function from func_base, + // as in principle there could be a get_name() function in the derived + // function class that hides it. return static_cast(this->value()).get_name(); } @@ -153,6 +153,9 @@ struct HEYOKA_DLL_PUBLIC_INLINE_CLASS func_iface_impl : public Base, tanuki::ifa [[nodiscard]] const std::vector &args() const final { + // NOTE: make sure we are invoking the member function from func_base, + // as in principle there could be an args() function in the derived + // function class that hides it. return static_cast(this->value()).args(); } std::pair get_mutable_args_range() final From 8a30de3570ebce3496b4590b3dedfae30b54ce51 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 6 Mar 2024 08:59:57 +0100 Subject: [PATCH 2/7] Some initial code for tseries conversion. --- CMakeLists.txt | 1 + include/heyoka/func.hpp | 29 ++++++ include/heyoka/math/sum.hpp | 4 + include/heyoka/tseries.hpp | 76 ++++++++++++++ src/func.cpp | 19 ++++ src/math/sum.cpp | 31 ++++++ src/tseries.cpp | 193 ++++++++++++++++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/tseries.cpp | 49 +++++++++ 9 files changed, 403 insertions(+) create mode 100644 include/heyoka/tseries.hpp create mode 100644 src/tseries.cpp create mode 100644 test/tseries.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 37c2b73c0..e9f34f8ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -174,6 +174,7 @@ set(HEYOKA_SRC_FILES "${CMAKE_CURRENT_SOURCE_DIR}/src/nt_event.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/t_event.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/dtens.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/src/tseries.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/cfunc_class.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/continuous_output.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/detail/cm_utils.cpp" diff --git a/include/heyoka/func.hpp b/include/heyoka/func.hpp index 687efe309..c5437832e 100644 --- a/include/heyoka/func.hpp +++ b/include/heyoka/func.hpp @@ -31,6 +31,7 @@ #include #include #include +#include // Current archive version is 2. // Changelog: @@ -115,6 +116,13 @@ concept func_has_taylor_decompose = requires(T &&x, taylor_dc_t &dc) { } -> std::same_as; }; +template +concept func_has_to_tseries = requires(const T &x, funcptr_map &m, std::uint32_t order) { + { + x.to_tseries(m, order) + } -> std::same_as; +}; + // Function interface implementation. template requires is_udf @@ -271,6 +279,22 @@ struct HEYOKA_DLL_PUBLIC_INLINE_CLASS func_iface_impl : public Base, tanuki::ifa fmt::format("Taylor diff in compact mode is not implemented for the function '{}'", get_name())); } } + + [[nodiscard]] bool has_to_tseries() const final + { + return func_has_to_tseries; + } + [[nodiscard]] tseries to_tseries(funcptr_map &m, std::uint32_t order) const final + { + if constexpr (func_has_to_tseries) { + return this->value().to_tseries(m, order); + } + + // LCOV_EXCL_START + assert(false); + throw; + // LCOV_EXCL_STOP + } }; // The function interface. @@ -318,6 +342,9 @@ struct HEYOKA_DLL_PUBLIC func_iface { virtual llvm::Function *taylor_c_diff_func(llvm_state &, llvm::Type *, std::uint32_t, std::uint32_t, bool) const = 0; + [[nodiscard]] virtual bool has_to_tseries() const = 0; + [[nodiscard]] virtual tseries to_tseries(funcptr_map &, std::uint32_t) const = 0; + template using impl = func_iface_impl; }; @@ -422,6 +449,8 @@ class HEYOKA_DLL_PUBLIC func const std::vector &, llvm::Value *, llvm::Value *, std::uint32_t, std::uint32_t, std::uint32_t, std::uint32_t, bool) const; llvm::Function *taylor_c_diff_func(llvm_state &, llvm::Type *, std::uint32_t, std::uint32_t, bool) const; + + [[nodiscard]] tseries to_tseries(detail::funcptr_map &, std::uint32_t order) const; }; namespace detail diff --git a/include/heyoka/math/sum.hpp b/include/heyoka/math/sum.hpp index 526870b63..9112ad28f 100644 --- a/include/heyoka/math/sum.hpp +++ b/include/heyoka/math/sum.hpp @@ -15,11 +15,13 @@ #include #include +#include #include #include #include #include #include +#include HEYOKA_BEGIN_NAMESPACE @@ -55,6 +57,8 @@ class HEYOKA_DLL_PUBLIC sum_impl : public func_base std::uint32_t, std::uint32_t, std::uint32_t, bool) const; llvm::Function *taylor_c_diff_func(llvm_state &, llvm::Type *, std::uint32_t, std::uint32_t, bool) const; + + tseries to_tseries(detail::funcptr_map &, std::uint32_t) const; }; HEYOKA_DLL_PUBLIC expression sum_split(const expression &, std::uint32_t); diff --git a/include/heyoka/tseries.hpp b/include/heyoka/tseries.hpp new file mode 100644 index 000000000..a3de4c359 --- /dev/null +++ b/include/heyoka/tseries.hpp @@ -0,0 +1,76 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef HEYOKA_TSERIES_HPP +#define HEYOKA_TSERIES_HPP + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +HEYOKA_BEGIN_NAMESPACE + +class HEYOKA_DLL_PUBLIC tseries +{ + struct impl; + std::unique_ptr m_impl; + + struct ptag; + template + HEYOKA_DLL_LOCAL explicit tseries(ptag, T, std::uint32_t); + +public: + tseries() = delete; + explicit tseries(std::vector); + explicit tseries(const variable &, std::uint32_t); + explicit tseries(number, std::uint32_t); + explicit tseries(param, std::uint32_t); + tseries(const tseries &); + tseries(tseries &&) noexcept; + tseries &operator=(const tseries &); + tseries &operator=(tseries &&) noexcept; + ~tseries(); + + [[nodiscard]] std::uint32_t get_order() const; + [[nodiscard]] const std::vector &get_cfs() const; +}; + +HEYOKA_DLL_PUBLIC std::ostream &operator<<(std::ostream &, const tseries &); + +namespace detail +{ + +HEYOKA_DLL_PUBLIC tseries to_tseries_impl(funcptr_map &, const expression &, std::uint32_t); + +} // namespace detail + +HEYOKA_DLL_PUBLIC std::vector to_tseries(const std::vector &, std::uint32_t); + +HEYOKA_END_NAMESPACE + +// fmt formatter for tseries, implemented +// on top of the streaming operator. +namespace fmt +{ + +template <> +struct formatter : fmt::ostream_formatter { +}; + +} // namespace fmt + +#endif diff --git a/src/func.cpp b/src/func.cpp index 0f6ef85ac..1779c150a 100644 --- a/src/func.cpp +++ b/src/func.cpp @@ -52,6 +52,7 @@ #include #include #include +#include #include HEYOKA_BEGIN_NAMESPACE @@ -859,6 +860,24 @@ llvm::Function *llvm_c_eval_func_helper(const std::string &name, } // namespace detail +tseries func::to_tseries(detail::funcptr_map &m, std::uint32_t order) const +{ + if (!m_func->has_to_tseries()) [[unlikely]] { + throw not_implemented_error( + fmt::format("The function '{}' does not support conversion to a Taylor series", get_name())); + } + + auto ret = m_func->to_tseries(m, order); + + if (ret.get_order() != order) [[unlikely]] { + throw std::invalid_argument(fmt::format("The conversion of the function '{}' to a Taylor series produced a " + "series of order {}, but a series of order {} was expected instead", + get_name(), ret.get_order(), order)); + } + + return ret; +} + HEYOKA_END_NAMESPACE // s11n implementation for null_func. diff --git a/src/math/sum.cpp b/src/math/sum.cpp index 7301be9bc..976cc4479 100644 --- a/src/math/sum.cpp +++ b/src/math/sum.cpp @@ -19,6 +19,8 @@ #include #include +#include + #include #include @@ -49,6 +51,7 @@ #include #include #include +#include #include #include @@ -412,6 +415,34 @@ llvm::Function *sum_impl::taylor_c_diff_func(llvm_state &s, llvm::Type *fp_t, st return sum_taylor_c_diff_func_impl(s, fp_t, *this, n_uvars, batch_size); } +tseries sum_impl::to_tseries(detail::funcptr_map &m, std::uint32_t order) const +{ + using su32_t = boost::safe_numerics::safe; + + // Convert the arguments to tseries. + std::vector tseries_args; + tseries_args.reserve(args().size()); + + for (const auto &arg : args()) { + tseries_args.push_back(to_tseries_impl(m, arg, order)); + } + + // Compute the coefficients of the output tseries. + std::vector tseries_cfs, tmp; + tseries_cfs.reserve(su32_t(order) + 1); + for (su32_t i = 0; i <= order; ++i) { + tmp.clear(); + + for (decltype(args().size()) arg_idx = 0; arg_idx < args().size(); ++arg_idx) { + tmp.push_back(tseries_args[arg_idx].get_cfs()[i]); + } + + tseries_cfs.push_back(sum(tmp)); + } + + return tseries(std::move(tseries_cfs)); +} + // Helper to split the input sum 'e' into nested sums, each // of which will have at most 'split' arguments. // If 'e' is not a sum, or if it is a sum with no more than diff --git a/src/tseries.cpp b/src/tseries.cpp new file mode 100644 index 000000000..246066a42 --- /dev/null +++ b/src/tseries.cpp @@ -0,0 +1,193 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +HEYOKA_BEGIN_NAMESPACE + +struct tseries::impl { + std::vector m_cf; +}; + +struct tseries::ptag { +}; + +tseries::tseries(std::vector vex) : m_impl(std::make_unique(impl{std::move(vex)})) +{ + if (m_impl->m_cf.empty()) [[unlikely]] { + throw std::invalid_argument("Cannot construct a tseries from an empty list of coefficients"); + } + + // LCOV_EXCL_START + if (m_impl->m_cf.size() > std::numeric_limits::max()) [[unlikely]] { + throw std::overflow_error("Overflow detected during the construction of a tseries: the order is too large"); + } + // LCOV_EXCL_STOP +} + +tseries::tseries(const variable &v, std::uint32_t order) +{ + using su32_t = boost::safe_numerics::safe; + + std::vector cf_vec; + cf_vec.reserve(su32_t(order) + 1); + for (su32_t i = 0; i <= order; ++i) { + cf_vec.emplace_back(fmt::format("cf_{}_{}", static_cast(i), v.name())); + } + + m_impl = std::make_unique(impl{std::move(cf_vec)}); +} + +template +tseries::tseries(ptag, T x, std::uint32_t order) +{ + using su32_t = boost::safe_numerics::safe; + + std::vector cf_vec; + cf_vec.reserve(su32_t(order) + 1); + cf_vec.emplace_back(std::move(x)); + for (std::uint32_t i = 0; i < order; ++i) { + cf_vec.emplace_back(0.); + } + + m_impl = std::make_unique(impl{std::move(cf_vec)}); +} + +tseries::tseries(number n, std::uint32_t order) : tseries(ptag{}, std::move(n), order) {} + +tseries::tseries(param p, std::uint32_t order) : tseries(ptag{}, std::move(p), order) {} + +tseries::tseries(const tseries &other) : m_impl(std::make_unique(*other.m_impl)) {} + +tseries::tseries(tseries &&) noexcept = default; + +tseries &tseries::operator=(const tseries &other) +{ + if (this != &other) { + // NOTE: this will correctly revive a moved-from object. + *this = tseries(other); + } + + return *this; +} + +tseries &tseries::operator=(tseries &&) noexcept = default; + +tseries::~tseries() = default; + +std::uint32_t tseries::get_order() const +{ + assert(!m_impl->m_cf.empty()); + + return static_cast(m_impl->m_cf.size() - 1u); +} + +const std::vector &tseries::get_cfs() const +{ + return m_impl->m_cf; +} + +std::ostream &operator<<(std::ostream &os, const tseries &ts) +{ + os << fmt::format("{}", ts.get_cfs()); + return os; +} + +namespace detail +{ + +tseries to_tseries_impl(funcptr_map &func_map, const expression &ex, std::uint32_t order) +{ + return std::visit( + [order, &func_map](const auto &arg) { + using type = std::remove_cvref_t; + + if constexpr (std::same_as) { + const auto *f_id = arg.get_ptr(); + + // Check if we already converted ex to tseries. + if (auto it = func_map.find(f_id); it != func_map.end()) { + return it->second; + } + + auto ret = arg.to_tseries(func_map, order); + + // Put the return value in the cache. + [[maybe_unused]] const auto [_, flag] = func_map.emplace(f_id, ret); + // NOTE: an expression cannot contain itself. + assert(flag); + + return ret; + } else { + return tseries(arg, order); + } + }, + ex.value()); +} + +} // namespace detail + +std::vector to_tseries(const std::vector &v_ex_, std::uint32_t order) +{ + // Pre-process v_ex with several transformations. + + // Transform sums into subs. + auto v_ex = detail::sum_to_sub(v_ex_); + + // Split sums. + v_ex = detail::split_sums_for_decompose(v_ex); + + // Transform prods into divs. + v_ex = detail::prod_to_div_taylor_diff(v_ex); + + // Split prods. + // NOTE: split must be 2 here as we implement only + // binary multiplication for tseries. + v_ex = detail::split_prods_for_decompose(v_ex, 2); + + // Unfix. + // NOTE: unfix is the last step, as we want to keep expressions + // fixed in the previous preprocessing steps. + v_ex = unfix(v_ex); + + detail::funcptr_map func_map; + + std::vector ret; + ret.reserve(v_ex.size()); + + for (const auto &ex : v_ex) { + ret.push_back(detail::to_tseries_impl(func_map, ex, order)); + } + + return ret; +} + +HEYOKA_END_NAMESPACE diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 534067424..a1066905e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -152,6 +152,7 @@ ADD_HEYOKA_TESTCASE(lagrangian) ADD_HEYOKA_TESTCASE(hamiltonian) ADD_HEYOKA_TESTCASE(cfunc) ADD_HEYOKA_TESTCASE(cfunc_multieval) +ADD_HEYOKA_TESTCASE(tseries) if(HEYOKA_WITH_MPPP AND mp++_WITH_MPFR) ADD_HEYOKA_TESTCASE(event_detection_mp) diff --git a/test/tseries.cpp b/test/tseries.cpp new file mode 100644 index 000000000..6a2b32d7f --- /dev/null +++ b/test/tseries.cpp @@ -0,0 +1,49 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#include +#include + +#include "catch.hpp" + +using namespace heyoka; + +TEST_CASE("tseries") {} + +TEST_CASE("to_tseries") +{ + auto ts = to_tseries({42_dbl}, 12)[0]; + + REQUIRE(ts.get_order() == 12u); + REQUIRE(ts.get_cfs()[0] == 42_dbl); + + for (auto i = 1u; i <= 12u; ++i) { + REQUIRE(ts.get_cfs()[i] == 0_dbl); + } + + ts = to_tseries({par[42]}, 12)[0]; + + REQUIRE(ts.get_order() == 12u); + REQUIRE(ts.get_cfs()[0] == par[42]); + + for (auto i = 1u; i <= 12u; ++i) { + REQUIRE(ts.get_cfs()[i] == 0_dbl); + } + + ts = to_tseries({"x"_var}, 12)[0]; + + REQUIRE(ts.get_order() == 12u); + + for (auto i = 0u; i <= 12u; ++i) { + REQUIRE(ts.get_cfs()[i] == expression(fmt::format("cf_{}_{}", i, "x"))); + } + + fmt::println("{}", to_tseries({"x"_var + "y"_var}, 12)[0]); +} From 46f437fd77b85c31cf6f985a2e18f95103a726f5 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 6 Mar 2024 09:19:37 +0100 Subject: [PATCH 3/7] Minor cleanup. --- src/taylor_00.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/taylor_00.cpp b/src/taylor_00.cpp index 7c1378e19..41fb9c2c1 100644 --- a/src/taylor_00.cpp +++ b/src/taylor_00.cpp @@ -10,7 +10,6 @@ #include #include -#include #include #include #include @@ -120,8 +119,6 @@ taylor_determine_h(llvm_state &s, llvm::Type *fp_t, } #endif - using std::exp; - auto &builder = s.builder(); llvm::Value *max_abs_state = nullptr, *max_abs_diff_o = nullptr, *max_abs_diff_om1 = nullptr; From 16ae6c8a80f405899f3d0f995c717060ff9cfc32 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Fri, 8 Mar 2024 08:26:43 +0100 Subject: [PATCH 4/7] Initial plumbing for handling the ad_mode option in the Taylor integrators. --- include/heyoka/detail/i_data.hpp | 4 +++ include/heyoka/detail/taylor_common.hpp | 23 ++++++++++-- include/heyoka/kw.hpp | 1 + include/heyoka/taylor.hpp | 47 +++++++++++++++++++------ src/detail/i_data.cpp | 12 ++++--- src/taylor_adaptive.cpp | 27 ++++++++++---- src/taylor_adaptive_batch.cpp | 27 ++++++++++---- src/taylor_stream_ops.cpp | 16 +++++++++ 8 files changed, 128 insertions(+), 29 deletions(-) diff --git a/include/heyoka/detail/i_data.hpp b/include/heyoka/detail/i_data.hpp index 088b3c316..8ca0318a8 100644 --- a/include/heyoka/detail/i_data.hpp +++ b/include/heyoka/detail/i_data.hpp @@ -43,6 +43,8 @@ struct taylor_adaptive::i_data { std::uint32_t m_order{}; // Tolerance. T m_tol{}; + // AD mode. + taylor_ad_mode m_ad_mode = taylor_ad_mode::classic; // High accuracy. bool m_high_accuracy{}; // Compact mode. @@ -105,6 +107,8 @@ struct taylor_adaptive_batch::i_data { std::uint32_t m_order{}; // Tolerance. T m_tol{}; + // AD mode. + taylor_ad_mode m_ad_mode = taylor_ad_mode::classic; // High accuracy. bool m_high_accuracy{}; // Compact mode. diff --git a/include/heyoka/detail/taylor_common.hpp b/include/heyoka/detail/taylor_common.hpp index ff9bef87f..c007a231a 100644 --- a/include/heyoka/detail/taylor_common.hpp +++ b/include/heyoka/detail/taylor_common.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -165,14 +166,32 @@ inline llvm::Function *taylor_c_diff_func_numpar(llvm_state &s, llvm::Type *fp_t // precision due to the way precision propagation works in mp++. I don't // think there's any negative consequence here. template -std::uint32_t taylor_order_from_tol(T tol) +std::uint32_t taylor_order_from_tol(T tol, taylor_ad_mode ad_mode) { using std::ceil; using std::isfinite; using std::log; + using std::log2; + + assert(ad_mode == taylor_ad_mode::classic || ad_mode == taylor_ad_mode::tseries); // Determine the order from the tolerance. - auto order_f = ceil(-log(tol) / 2 + 1); + auto order_f = [&]() { + if (ad_mode == taylor_ad_mode::classic) { + return ceil(-log(tol) / 2 + 1); + } else { +#if defined(HEYOKA_HAVE_REAL) + if constexpr (std::same_as) { + return ceil(-log(tol) / log2(mppp::real(3, tol.get_prec())) + 1); + } else { +#endif + return ceil(-log(tol) / log2(static_cast(3)) + 1); +#if defined(HEYOKA_HAVE_REAL) + } +#endif + } + }(); + // LCOV_EXCL_START if (!isfinite(order_f)) { throw std::invalid_argument( diff --git a/include/heyoka/kw.hpp b/include/heyoka/kw.hpp index f46a485a3..ef9854737 100644 --- a/include/heyoka/kw.hpp +++ b/include/heyoka/kw.hpp @@ -43,6 +43,7 @@ IGOR_MAKE_NAMED_ARGUMENT(check_prec); IGOR_MAKE_NAMED_ARGUMENT(tol); IGOR_MAKE_NAMED_ARGUMENT(t_events); IGOR_MAKE_NAMED_ARGUMENT(nt_events); +IGOR_MAKE_NAMED_ARGUMENT(ad_mode); // NOTE: these are used for constructing events. IGOR_MAKE_NAMED_ARGUMENT(callback); IGOR_MAKE_NAMED_ARGUMENT(cooldown); diff --git a/include/heyoka/taylor.hpp b/include/heyoka/taylor.hpp index 911cadb5c..bbbcf22e1 100644 --- a/include/heyoka/taylor.hpp +++ b/include/heyoka/taylor.hpp @@ -184,9 +184,14 @@ enum class taylor_outcome : std::int64_t { HEYOKA_DLL_PUBLIC std::ostream &operator<<(std::ostream &, taylor_outcome); +// Enum to represent the automatic differentiation mode used by Taylor integrators. +enum class taylor_ad_mode { classic, tseries }; + +HEYOKA_DLL_PUBLIC std::ostream &operator<<(std::ostream &, taylor_ad_mode); + HEYOKA_END_NAMESPACE -// fmt formatter for taylor_outcome, implemented on top of the streaming operator. +// fmt formatters for taylor_outcome and taylor_ad_mode, implemented on top of the streaming operators. namespace fmt { @@ -194,6 +199,10 @@ template <> struct formatter : fmt::ostream_formatter { }; +template <> +struct formatter : fmt::ostream_formatter { +}; + } // namespace fmt HEYOKA_BEGIN_NAMESPACE @@ -266,7 +275,20 @@ auto taylor_adaptive_common_ops(const KwArgs &...kw_args) } }(); - return std::tuple{high_accuracy, std::move(tol), compact_mode, std::move(pars), parallel_mode}; + // AD mode (defaults to classic). + auto ad_mode = [&p]() -> taylor_ad_mode { + if constexpr (p.has(kw::ad_mode)) { + if constexpr (std::same_as, taylor_ad_mode>) { + return p(kw::ad_mode); + } else { + static_assert(always_false_v, "Invalid type for the 'ad_mode' keyword argument."); + } + } else { + return taylor_ad_mode::classic; + } + }(); + + return std::tuple{high_accuracy, std::move(tol), compact_mode, std::move(pars), parallel_mode, ad_mode}; } // Small helper to construct a default value for the max_delta_t @@ -470,7 +492,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada // Private implementation-detail constructor machinery. void finalise_ctor_impl(const std::vector> &, std::vector, std::optional, std::optional, bool, bool, std::vector, std::vector, - std::vector, bool, std::optional); + std::vector, bool, std::optional, taylor_ad_mode); template void finalise_ctor(const std::vector> &sys, std::vector state, const KwArgs &...kw_args) @@ -491,7 +513,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada } }(); - auto [high_accuracy, tol, compact_mode, pars, parallel_mode] + auto [high_accuracy, tol, compact_mode, pars, parallel_mode, ad_mode] = detail::taylor_adaptive_common_ops(kw_args...); // Extract the terminal events, if any. @@ -526,7 +548,8 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada }(); finalise_ctor_impl(sys, std::move(state), std::move(tm), std::move(tol), high_accuracy, compact_mode, - std::move(pars), std::move(tes), std::move(ntes), parallel_mode, std::move(prec)); + std::move(pars), std::move(tes), std::move(ntes), parallel_mode, std::move(prec), + ad_mode); } } @@ -570,6 +593,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive : public detail::taylor_ada [[nodiscard]] bool get_high_accuracy() const; [[nodiscard]] bool get_compact_mode() const; [[nodiscard]] std::uint32_t get_dim() const; + [[nodiscard]] taylor_ad_mode get_ad_mode() const; [[nodiscard]] T get_time() const; void set_time(T); @@ -839,7 +863,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch // Private implementation-detail constructor machinery. void finalise_ctor_impl(const std::vector> &, std::vector, std::uint32_t, std::vector, std::optional, bool, bool, std::vector, std::vector, - std::vector, bool); + std::vector, bool, taylor_ad_mode); template void finalise_ctor(const std::vector> &sys, std::vector state, std::uint32_t batch_size, const KwArgs &...kw_args) @@ -860,7 +884,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch } }(); - auto [high_accuracy, tol, compact_mode, pars, parallel_mode] + auto [high_accuracy, tol, compact_mode, pars, parallel_mode, ad_mode] = detail::taylor_adaptive_common_ops(kw_args...); // Extract the terminal events, if any. @@ -882,7 +906,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch }(); finalise_ctor_impl(sys, std::move(state), batch_size, std::move(tm), std::move(tol), high_accuracy, - compact_mode, std::move(pars), std::move(tes), std::move(ntes), parallel_mode); + compact_mode, std::move(pars), std::move(tes), std::move(ntes), parallel_mode, ad_mode); } } @@ -927,6 +951,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch [[nodiscard]] bool get_high_accuracy() const; [[nodiscard]] bool get_compact_mode() const; [[nodiscard]] std::uint32_t get_dim() const; + [[nodiscard]] taylor_ad_mode get_ad_mode() const; [[nodiscard]] const std::vector &get_time() const; [[nodiscard]] const T *get_time_data() const; @@ -1132,14 +1157,16 @@ namespace detail // - 3: removed the mr flag from the terminal event callback siganture, // which resulted also in changes in the event detection data structure. // - 4: switched to pimpl implementation for i_data. -inline constexpr int taylor_adaptive_s11n_version = 4; +// - 5: added the m_ad_mode member to i_data. +inline constexpr int taylor_adaptive_s11n_version = 5; // Boost s11n class version history for taylor_adaptive_batch: // - 1: added the m_state_vars and m_rhs members. // - 2: removed the mr flag from the terminal event callback siganture, // which resulted also in changes in the event detection data structure. // - 3: switched to pimpl implementation for i_data. -inline constexpr int taylor_adaptive_batch_s11n_version = 3; +// - 4: added the m_ad_mode member to i_data. +inline constexpr int taylor_adaptive_batch_s11n_version = 4; } // namespace detail diff --git a/src/detail/i_data.cpp b/src/detail/i_data.cpp index 58874fe71..71aec2d75 100644 --- a/src/detail/i_data.cpp +++ b/src/detail/i_data.cpp @@ -97,6 +97,7 @@ void taylor_adaptive::i_data::save(boost::archive::binary_oarchive &ar, unsig ar << m_dc; ar << m_order; ar << m_tol; + ar << m_ad_mode; ar << m_high_accuracy; ar << m_compact_mode; ar << m_pars; @@ -117,6 +118,7 @@ void taylor_adaptive::i_data::load(boost::archive::binary_iarchive &ar, unsig ar >> m_dc; ar >> m_order; ar >> m_tol; + ar >> m_ad_mode; ar >> m_high_accuracy; ar >> m_compact_mode; ar >> m_pars; @@ -138,7 +140,7 @@ taylor_adaptive::i_data::i_data(llvm_state s) : m_llvm(std::move(s)) template taylor_adaptive::i_data::i_data(const i_data &other) : m_state(other.m_state), m_time(other.m_time), m_llvm(other.m_llvm), m_dim(other.m_dim), m_dc(other.m_dc), - m_order(other.m_order), m_tol(other.m_tol), m_high_accuracy(other.m_high_accuracy), + m_order(other.m_order), m_tol(other.m_tol), m_ad_mode(other.m_ad_mode), m_high_accuracy(other.m_high_accuracy), m_compact_mode(other.m_compact_mode), m_pars(other.m_pars), m_tc(other.m_tc), m_last_h(other.m_last_h), m_d_out(other.m_d_out), m_state_vars(other.m_state_vars), m_rhs(other.m_rhs) { @@ -184,6 +186,7 @@ void taylor_adaptive_batch::i_data::save(boost::archive::binary_oarchive &ar, ar << m_dc; ar << m_order; ar << m_tol; + ar << m_ad_mode; ar << m_high_accuracy; ar << m_compact_mode; ar << m_pars; @@ -222,6 +225,7 @@ void taylor_adaptive_batch::i_data::load(boost::archive::binary_iarchive &ar, ar >> m_dc; ar >> m_order; ar >> m_tol; + ar >> m_ad_mode; ar >> m_high_accuracy; ar >> m_compact_mode; ar >> m_pars; @@ -260,9 +264,9 @@ template taylor_adaptive_batch::i_data::i_data(const i_data &other) : m_batch_size(other.m_batch_size), m_state(other.m_state), m_time_hi(other.m_time_hi), m_time_lo(other.m_time_lo), m_llvm(other.m_llvm), m_dim(other.m_dim), m_dc(other.m_dc), m_order(other.m_order), m_tol(other.m_tol), - m_high_accuracy(other.m_high_accuracy), m_compact_mode(other.m_compact_mode), m_pars(other.m_pars), - m_tc(other.m_tc), m_last_h(other.m_last_h), m_d_out(other.m_d_out), m_pinf(other.m_pinf), m_minf(other.m_minf), - m_delta_ts(other.m_delta_ts), m_step_res(other.m_step_res), m_prop_res(other.m_prop_res), + m_ad_mode(other.m_ad_mode), m_high_accuracy(other.m_high_accuracy), m_compact_mode(other.m_compact_mode), + m_pars(other.m_pars), m_tc(other.m_tc), m_last_h(other.m_last_h), m_d_out(other.m_d_out), m_pinf(other.m_pinf), + m_minf(other.m_minf), m_delta_ts(other.m_delta_ts), m_step_res(other.m_step_res), m_prop_res(other.m_prop_res), m_ts_count(other.m_ts_count), m_min_abs_h(other.m_min_abs_h), m_max_abs_h(other.m_max_abs_h), m_cur_max_delta_ts(other.m_cur_max_delta_ts), m_pfor_ts(other.m_pfor_ts), m_t_dir(other.m_t_dir), m_rem_time(other.m_rem_time), m_time_copy_hi(other.m_time_copy_hi), m_time_copy_lo(other.m_time_copy_lo), diff --git a/src/taylor_adaptive.cpp b/src/taylor_adaptive.cpp index 73e2fbce8..972837939 100644 --- a/src/taylor_adaptive.cpp +++ b/src/taylor_adaptive.cpp @@ -163,7 +163,7 @@ void taylor_adaptive::finalise_ctor_impl(const std::vector time, std::optional tol, bool high_accuracy, bool compact_mode, std::vector pars, std::vector tes, std::vector ntes, bool parallel_mode, - [[maybe_unused]] std::optional prec) + [[maybe_unused]] std::optional prec, taylor_ad_mode ad_mode) { HEYOKA_TAYLOR_REF_FROM_I_DATA(m_step_f); HEYOKA_TAYLOR_REF_FROM_I_DATA(m_state); @@ -175,6 +175,7 @@ void taylor_adaptive::finalise_ctor_impl(const std::vector::finalise_ctor_impl(const std::vector::finalise_ctor_impl(const std::vector::finalise_ctor_impl(const std::vector taylor_ad_mode::tseries) [[unlikely]] { + throw std::invalid_argument("An invalid enumerator was specified for the 'ad_mode' keyword argument"); + } + m_ad_mode = ad_mode; + // Store the dimension of the system. m_dim = boost::numeric_cast(sys.size()); @@ -330,7 +337,7 @@ void taylor_adaptive::finalise_ctor_impl(const std::vector(m_llvm.context()); @@ -1615,6 +1622,12 @@ std::uint32_t taylor_adaptive::get_dim() const return m_i_data->m_dim; } +template +taylor_ad_mode taylor_adaptive::get_ad_mode() const +{ + return m_i_data->m_ad_mode; +} + template T taylor_adaptive::get_time() const { diff --git a/src/taylor_adaptive_batch.cpp b/src/taylor_adaptive_batch.cpp index db463ef42..f2d4aa2c5 100644 --- a/src/taylor_adaptive_batch.cpp +++ b/src/taylor_adaptive_batch.cpp @@ -74,7 +74,8 @@ void taylor_adaptive_batch::finalise_ctor_impl(const std::vector state, std::uint32_t batch_size, std::vector time, std::optional tol, bool high_accuracy, bool compact_mode, std::vector pars, std::vector tes, - std::vector ntes, bool parallel_mode) + std::vector ntes, bool parallel_mode, + taylor_ad_mode ad_mode) { // NOTE: this must hold because tol == 0 is interpreted // as undefined in finalise_ctor(). @@ -96,6 +97,7 @@ void taylor_adaptive_batch::finalise_ctor_impl(const std::vector::finalise_ctor_impl(const std::vector::finalise_ctor_impl(const std::vector::epsilon(); } + // Validate and store the AD mode. + if (ad_mode < taylor_ad_mode::classic || ad_mode > taylor_ad_mode::tseries) [[unlikely]] { + throw std::invalid_argument("An invalid enumerator was specified for the 'ad_mode' keyword argument"); + } + m_ad_mode = ad_mode; + // Store the dimension of the system. m_dim = boost::numeric_cast(sys.size()); @@ -182,7 +191,7 @@ void taylor_adaptive_batch::finalise_ctor_impl(const std::vector(m_llvm.context()); @@ -1959,6 +1968,12 @@ std::uint32_t taylor_adaptive_batch::get_dim() const return m_i_data->m_dim; } +template +taylor_ad_mode taylor_adaptive_batch::get_ad_mode() const +{ + return m_i_data->m_ad_mode; +} + template const std::vector &taylor_adaptive_batch::get_time() const { diff --git a/src/taylor_stream_ops.cpp b/src/taylor_stream_ops.cpp index c3038122d..4761a15b6 100644 --- a/src/taylor_stream_ops.cpp +++ b/src/taylor_stream_ops.cpp @@ -259,4 +259,20 @@ std::ostream &operator<<(std::ostream &os, taylor_outcome oc) #undef HEYOKA_TAYLOR_ENUM_STREAM_CASE +std::ostream &operator<<(std::ostream &os, taylor_ad_mode m) +{ + switch (m) { + case taylor_ad_mode::classic: + os << "classic"; + break; + case taylor_ad_mode::tseries: + os << "tseries"; + break; + default: + os << "invalid"; + } + + return os; +} + HEYOKA_END_NAMESPACE From 1470496dc63cd44782d9087b8957295fe2e2fa8c Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 17 Mar 2024 10:28:50 +0100 Subject: [PATCH 5/7] Small doc addition. --- doc/tut_cfunc.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/tut_cfunc.rst b/doc/tut_cfunc.rst index 00d398001..55b3cac8c 100644 --- a/doc/tut_cfunc.rst +++ b/doc/tut_cfunc.rst @@ -3,6 +3,8 @@ Compiled functions ================== +.. versionadded:: 4.0.0 + .. cpp:namespace-push:: heyoka heyoka can compile just-in-time (JIT) multivariate vector functions defined From 43b4623b82c629cc2a969325d34dbf91eb87a908 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 2 May 2024 19:10:32 +0200 Subject: [PATCH 6/7] [skip ci] From c1dbbfbde23bebd8a47af09709e72cfc22496ae6 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 2 May 2024 19:40:37 +0200 Subject: [PATCH 7/7] Fixes for recent changes. [skip ci] --- include/heyoka/func.hpp | 2 +- src/tseries.cpp | 5 ----- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/include/heyoka/func.hpp b/include/heyoka/func.hpp index 8060169ac..ef2c5e55f 100644 --- a/include/heyoka/func.hpp +++ b/include/heyoka/func.hpp @@ -263,7 +263,7 @@ struct HEYOKA_DLL_PUBLIC_INLINE_CLASS func_iface_impl : public Base { [[nodiscard]] tseries to_tseries(funcptr_map &m, std::uint32_t order) const final { if constexpr (func_has_to_tseries) { - return this->value().to_tseries(m, order); + return getval(this).to_tseries(m, order); } // LCOV_EXCL_START diff --git a/src/tseries.cpp b/src/tseries.cpp index 246066a42..822cfd9da 100644 --- a/src/tseries.cpp +++ b/src/tseries.cpp @@ -173,11 +173,6 @@ std::vector to_tseries(const std::vector &v_ex_, std::uint3 // binary multiplication for tseries. v_ex = detail::split_prods_for_decompose(v_ex, 2); - // Unfix. - // NOTE: unfix is the last step, as we want to keep expressions - // fixed in the previous preprocessing steps. - v_ex = unfix(v_ex); - detail::funcptr_map func_map; std::vector ret;