tseries

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"

include/heyoka/detail/i_data.hpp

@@ -43,6 +43,8 @@ struct taylor_adaptive<T>::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<T>::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.

include/heyoka/detail/taylor_common.hpp

@@ -20,6 +20,7 @@
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <concepts>
 #include <cstdint>
 #include <functional>
 #include <limits>
@@ -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 <typename T>
-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<T, mppp::real>) {
+                return ceil(-log(tol) / log2(mppp::real(3, tol.get_prec())) + 1);
+            } else {
+#endif
+                return ceil(-log(tol) / log2(static_cast<T>(3)) + 1);
+#if defined(HEYOKA_HAVE_REAL)
+            }
+#endif
+        }
+    }();
+
     // LCOV_EXCL_START
     if (!isfinite(order_f)) {
         throw std::invalid_argument(

include/heyoka/func.hpp

@@ -31,6 +31,7 @@
 #include <heyoka/detail/visibility.hpp>
 #include <heyoka/exceptions.hpp>
 #include <heyoka/s11n.hpp>
+#include <heyoka/tseries.hpp>
 
 // Current archive version is 2.
 // Changelog:
@@ -100,6 +101,11 @@ concept func_has_taylor_decompose = requires(T &&x, taylor_dc_t &dc) {
     { std::move(x).taylor_decompose(dc) } -> std::same_as<taylor_dc_t::size_type>;
 };
 
+template <typename T>
+concept func_has_to_tseries = requires(const T &x, funcptr_map<tseries> &m, std::uint32_t order) {
+    { x.to_tseries(m, order) } -> std::same_as<tseries>;
+};
+
 // Function interface implementation.
 template <typename Base, typename Holder, typename T>
     requires is_udf<T>
@@ -249,6 +255,22 @@ struct HEYOKA_DLL_PUBLIC_INLINE_CLASS func_iface_impl : public Base {
                 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<T>;
+    }
+    [[nodiscard]] tseries to_tseries(funcptr_map<tseries> &m, std::uint32_t order) const final
+    {
+        if constexpr (func_has_to_tseries<T>) {
+            return getval<Holder>(this).to_tseries(m, order);
+        }
+
+        // LCOV_EXCL_START
+        assert(false);
+        throw;
+        // LCOV_EXCL_STOP
+    }
 };
 
 // The function interface.
@@ -293,6 +315,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<tseries> &, std::uint32_t) const = 0;
+
     template <typename Base, typename Holder, typename T>
     using impl = func_iface_impl<Base, Holder, T>;
 };
@@ -395,6 +420,8 @@ class HEYOKA_DLL_PUBLIC func
                              const std::vector<llvm::Value *> &, 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<tseries> &, std::uint32_t order) const;
 };
 
 namespace detail

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);

include/heyoka/math/sum.hpp

@@ -14,11 +14,13 @@
 #include <vector>
 
 #include <heyoka/config.hpp>
+#include <heyoka/detail/func_cache.hpp>
 #include <heyoka/detail/fwd_decl.hpp>
 #include <heyoka/detail/llvm_fwd.hpp>
 #include <heyoka/detail/visibility.hpp>
 #include <heyoka/func.hpp>
 #include <heyoka/s11n.hpp>
+#include <heyoka/tseries.hpp>
 
 HEYOKA_BEGIN_NAMESPACE
 
@@ -52,6 +54,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<tseries> &, std::uint32_t) const;
 };
 
 HEYOKA_DLL_PUBLIC expression sum_split(const expression &, std::uint32_t);

include/heyoka/taylor.hpp

@@ -145,16 +145,25 @@ 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
 {
 
 template <>
 struct formatter<heyoka::taylor_outcome> : fmt::ostream_formatter {
 };
 
+template <>
+struct formatter<heyoka::taylor_ad_mode> : fmt::ostream_formatter {
+};
+
 } // namespace fmt
 
 HEYOKA_BEGIN_NAMESPACE
@@ -227,7 +236,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<std::remove_cvref_t<decltype(p(kw::ad_mode))>, taylor_ad_mode>) {
+                return p(kw::ad_mode);
+            } else {
+                static_assert(always_false_v<T>, "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
@@ -433,7 +455,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::pair<expression, expression>> &, std::vector<T>, std::optional<T>,
                             std::optional<T>, bool, bool, std::vector<T>, std::vector<t_event_t>,
-                            std::vector<nt_event_t>, bool, std::optional<long long>);
+                            std::vector<nt_event_t>, bool, std::optional<long long>, taylor_ad_mode);
     template <typename... KwArgs>
     void finalise_ctor(const std::vector<std::pair<expression, expression>> &sys, std::vector<T> state,
                        const KwArgs &...kw_args)
@@ -454,7 +476,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<T>(kw_args...);
 
             // Extract the terminal events, if any.
@@ -489,7 +511,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);
         }
     }
 
@@ -533,6 +556,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);
@@ -806,7 +830,7 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS taylor_adaptive_batch
     // Private implementation-detail constructor machinery.
     void finalise_ctor_impl(const std::vector<std::pair<expression, expression>> &, std::vector<T>, std::uint32_t,
                             std::vector<T>, std::optional<T>, bool, bool, std::vector<T>, std::vector<t_event_t>,
-                            std::vector<nt_event_t>, bool);
+                            std::vector<nt_event_t>, bool, taylor_ad_mode);
     template <typename... KwArgs>
     void finalise_ctor(const std::vector<std::pair<expression, expression>> &sys, std::vector<T> state,
                        std::uint32_t batch_size, const KwArgs &...kw_args)
@@ -829,7 +853,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<T>(kw_args...);
 
             // Extract the terminal events, if any.
@@ -851,7 +875,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);
         }
     }
 
@@ -896,6 +920,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<T> &get_time() const;
     [[nodiscard]] const T *get_time_data() const;
@@ -1103,14 +1128,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
 

include/heyoka/tseries.hpp

@@ -0,0 +1,76 @@
+// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani ([email protected]), Dario Izzo ([email protected])
+//
+// 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 <cstdint>
+#include <memory>
+#include <ostream>
+#include <vector>
+
+#include <fmt/core.h>
+#include <fmt/ostream.h>
+
+#include <heyoka/config.hpp>
+#include <heyoka/detail/func_cache.hpp>
+#include <heyoka/detail/fwd_decl.hpp>
+#include <heyoka/detail/visibility.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+class HEYOKA_DLL_PUBLIC tseries
+{
+    struct impl;
+    std::unique_ptr<impl> m_impl;
+
+    struct ptag;
+    template <typename T>
+    HEYOKA_DLL_LOCAL explicit tseries(ptag, T, std::uint32_t);
+
+public:
+    tseries() = delete;
+    explicit tseries(std::vector<expression>);
+    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<expression> &get_cfs() const;
+};
+
+HEYOKA_DLL_PUBLIC std::ostream &operator<<(std::ostream &, const tseries &);
+
+namespace detail
+{
+
+HEYOKA_DLL_PUBLIC tseries to_tseries_impl(funcptr_map<tseries> &, const expression &, std::uint32_t);
+
+} // namespace detail
+
+HEYOKA_DLL_PUBLIC std::vector<tseries> to_tseries(const std::vector<expression> &, std::uint32_t);
+
+HEYOKA_END_NAMESPACE
+
+// fmt formatter for tseries, implemented
+// on top of the streaming operator.
+namespace fmt
+{
+
+template <>
+struct formatter<heyoka::tseries> : fmt::ostream_formatter {
+};
+
+} // namespace fmt
+
+#endif