Merge pull request #379 from bluescarni/pr/lagham

Analytical mechanics
bluescarni · Dec 31, 2023 · 2de3a52 · 2de3a52
2 parents d46c2ac + f5589ef
commit 2de3a52
Show file tree

Hide file tree

Showing 14 changed files with 817 additions and 8 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -193,6 +193,8 @@ set(HEYOKA_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/expression_ops.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/expression_cfunc.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/expression_fix.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/lagrangian.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/hamiltonian.cpp"
     # "${CMAKE_CURRENT_SOURCE_DIR}/src/gp.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/logging.cpp"
     # VSOP2013 details.

diff --git a/README.md b/README.md
@@ -42,10 +42,13 @@ compilation via [LLVM](https://llvm.org/). Notable features include:
 
 * support for single-precision, double-precision, extended-precision (80-bit and 128-bit),
   and arbitrary-precision floating-point types,
-* the ability to maintain machine precision accuracy over
-  tens of billions of timesteps,
 * high-precision zero-cost dense output,
 * accurate and reliable event detection,
+* builtin support for analytical mechanics - bring your own Lagrangians/Hamiltonians
+  and let heyoka formulate and solve the equations of motion,
+* builtin support for machine learning applications via neural network models,
+* the ability to maintain machine precision accuracy over
+  tens of billions of timesteps,
 * batch mode integration to harness the power of modern
   [SIMD](https://en.wikipedia.org/wiki/SIMD) instruction sets
   (including AVX/AVX2/AVX-512/Neon/VSX),

diff --git a/doc/bibliography.rst b/doc/bibliography.rst
@@ -1,4 +1,4 @@
-.. only:: html or text
+.. only:: html or text or linkcheck
 
    Bibliography
    ============

diff --git a/doc/changelog.rst b/doc/changelog.rst
@@ -7,13 +7,17 @@ Changelog
 New
 ~~~
 
+- Add support for Lagrangian and Hamiltonian mechanics
+  (`#379 <https://github.com/bluescarni/heyoka/pull/379>`__).
 - It is now possible to pass a range of step callbacks to the
   ``propagate_*()`` functions. The individual callbacks will be
   automatically composed into a callback set
   (`#376 <https://github.com/bluescarni/heyoka/pull/376>`__).
 - New ``angle_reducer`` step callback to automatically reduce
   angular state variables to the :math:`\left[0, 2\pi\right)` range
   (`#376 <https://github.com/bluescarni/heyoka/pull/376>`__).
+- New ``callback`` module containing ready-made step and event callbacks
+  (`#376 <https://github.com/bluescarni/heyoka/pull/376>`__).
 
 Changes
 ~~~~~~~

diff --git a/doc/index.rst b/doc/index.rst
@@ -19,10 +19,13 @@ compilation via `LLVM <https://llvm.org/>`__. Notable features include:
 
 * support for single-precision, double-precision, extended-precision (80-bit and 128-bit),
   and arbitrary-precision floating-point types,
-* the ability to maintain machine precision accuracy over
-  tens of billions of timesteps,
 * high-precision zero-cost dense output,
 * accurate and reliable event detection,
+* builtin support for analytical mechanics - bring your own Lagrangians/Hamiltonians
+  and let heyoka formulate and solve the equations of motion,
+* builtin support for machine learning applications via neural network models,
+* the ability to maintain machine precision accuracy over
+  tens of billions of timesteps,
 * batch mode integration to harness the power of modern
   `SIMD <https://en.wikipedia.org/wiki/SIMD>`__ instruction sets
   (including AVX/AVX2/AVX-512/Neon/VSX),

diff --git a/include/heyoka/expression.hpp b/include/heyoka/expression.hpp
@@ -547,7 +547,7 @@ HEYOKA_BEGIN_NAMESPACE
 // returned by this function are optimised for evaluation. The users
 // can always unfix() and normalise() these expressions if needed.
 template <typename... KwArgs>
-dtens diff_tensors(const std::vector<expression> &v_ex, KwArgs &&...kw_args)
+dtens diff_tensors(const std::vector<expression> &v_ex, const KwArgs &...kw_args)
 {
     igor::parser p{kw_args...};
 
@@ -746,7 +746,7 @@ HEYOKA_CFUNC_EXTERN_INST(mppp::real)
 
 template <typename T, typename... KwArgs>
 inline std::vector<expression> add_cfunc(llvm_state &s, const std::string &name, const std::vector<expression> &fn,
-                                         KwArgs &&...kw_args)
+                                         const KwArgs &...kw_args)
 {
     igor::parser p{kw_args...};
 

diff --git a/include/heyoka/hamiltonian.hpp b/include/heyoka/hamiltonian.hpp
@@ -0,0 +1,26 @@
+// Copyright 2020, 2021, 2022, 2023 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_HAMILTONIAN_HPP
+#define HEYOKA_HAMILTONIAN_HPP
+
+#include <utility>
+#include <vector>
+
+#include <heyoka/config.hpp>
+#include <heyoka/detail/visibility.hpp>
+#include <heyoka/expression.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+HEYOKA_DLL_PUBLIC std::vector<std::pair<expression, expression>>
+hamiltonian(const expression &, const std::vector<expression> &, const std::vector<expression> &);
+
+HEYOKA_END_NAMESPACE
+
+#endif
diff --git a/include/heyoka/lagrangian.hpp b/include/heyoka/lagrangian.hpp
@@ -0,0 +1,28 @@
+// Copyright 2020, 2021, 2022, 2023 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_LAGRANGIAN_HPP
+#define HEYOKA_LAGRANGIAN_HPP
+
+#include <utility>
+#include <vector>
+
+#include <heyoka/config.hpp>
+#include <heyoka/detail/visibility.hpp>
+#include <heyoka/expression.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+HEYOKA_DLL_PUBLIC std::vector<std::pair<expression, expression>> lagrangian(const expression &,
+                                                                            const std::vector<expression> &,
+                                                                            const std::vector<expression> &,
+                                                                            const expression & = 0_dbl);
+
+HEYOKA_END_NAMESPACE
+
+#endif
diff --git a/src/hamiltonian.cpp b/src/hamiltonian.cpp
@@ -0,0 +1,147 @@
+// Copyright 2020, 2021, 2022, 2023 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/.
+
+#include <stdexcept>
+#include <unordered_set>
+#include <utility>
+#include <variant>
+#include <vector>
+
+#include <boost/safe_numerics/safe_integer.hpp>
+
+#include <fmt/core.h>
+
+#include <heyoka/config.hpp>
+#include <heyoka/expression.hpp>
+#include <heyoka/hamiltonian.hpp>
+#include <heyoka/kw.hpp>
+#include <heyoka/variable.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+namespace detail
+{
+
+namespace
+{
+
+void hamiltonian_impl_sanity_checks(const expression &H, const std::vector<expression> &qs,
+                                    const std::vector<expression> &ps)
+{
+    // Sanity checks on qs and ps.
+    if (qs.size() != ps.size()) {
+        throw std::invalid_argument(fmt::format(
+            "The number of generalised coordinates ({}) must be equal to the number of generalised momenta ({})",
+            qs.size(), ps.size()));
+    }
+
+    if (qs.empty()) {
+        throw std::invalid_argument("Cannot define a Hamiltonian without state variables");
+    }
+
+    for (const auto &q : qs) {
+        if (!std::holds_alternative<variable>(q.value())) {
+            throw std::invalid_argument(fmt::format(
+                "The list of generalised coordinates contains the expression '{}' which is not a variable", q));
+        }
+
+        if (std::get<variable>(q.value()).name().starts_with("__")) {
+            throw std::invalid_argument(
+                fmt::format("The list of generalised coordinates contains a variable with the invalid name '{}': names "
+                            "starting with '__' are reserved for internal use",
+                            std::get<variable>(q.value()).name()));
+        }
+    }
+
+    for (const auto &p : ps) {
+        if (!std::holds_alternative<variable>(p.value())) {
+            throw std::invalid_argument(
+                fmt::format("The list of generalised momenta contains the expression '{}' which is not a variable", p));
+        }
+
+        if (std::get<variable>(p.value()).name().starts_with("__")) {
+            throw std::invalid_argument(
+                fmt::format("The list of generalised momenta contains a variable with the invalid name '{}': names "
+                            "starting with '__' are reserved for internal use",
+                            std::get<variable>(p.value()).name()));
+        }
+    }
+
+    // Check for duplicates.
+    const std::unordered_set<expression> qs_set{qs.begin(), qs.end()};
+    const std::unordered_set<expression> ps_set{ps.begin(), ps.end()};
+
+    if (qs_set.size() != qs.size()) {
+        throw std::invalid_argument("The list of generalised coordinates contains duplicates");
+    }
+
+    if (ps_set.size() != ps.size()) {
+        throw std::invalid_argument("The list of generalised momenta contains duplicates");
+    }
+
+    for (const auto &q : qs) {
+        if (ps_set.contains(q)) {
+            throw std::invalid_argument(fmt::format("The list of generalised coordinates contains the expression '{}' "
+                                                    "which also appears as a generalised momentum",
+                                                    q));
+        }
+    }
+
+    // Sanity checks on H.
+    for (const auto &v : get_variables(H)) {
+        if (!qs_set.contains(expression{v}) && !ps_set.contains(expression{v})) {
+            throw std::invalid_argument(fmt::format(
+                "The Hamiltonian contains the variable '{}' which is not a generalised position or momentum", v));
+        }
+    }
+}
+
+} // namespace
+
+} // namespace detail
+
+std::vector<std::pair<expression, expression>> hamiltonian(const expression &H, const std::vector<expression> &qs,
+                                                           const std::vector<expression> &ps)
+{
+    using size_type = boost::safe_numerics::safe<decltype(qs.size())>;
+
+    // Sanity checks.
+    detail::hamiltonian_impl_sanity_checks(H, qs, ps);
+
+    // Cache the number of generalised coordinates/momenta.
+    const auto n_qs = size_type(qs.size());
+
+    // Assemble the diff arguments.
+    auto diff_args = qs;
+    diff_args.insert(diff_args.end(), ps.begin(), ps.end());
+
+    // Compute the tensor of derivatives of H up to order 1 wrt
+    // qs and ps.
+    const auto H_dt = diff_tensors({H}, kw::diff_args = diff_args, kw::diff_order = 1);
+
+    // Fetch the gradient.
+    auto grad = H_dt.get_gradient();
+
+    // Assemble the return value.
+    std::vector<std::pair<expression, expression>> ret;
+    ret.reserve(n_qs * 2);
+
+    // dq/dt.
+    for (size_type i = 0; i < n_qs; ++i) {
+        ret.emplace_back(qs[i], std::move(grad[n_qs + i]));
+    }
+
+    // dp/dt.
+    for (size_type i = 0; i < n_qs; ++i) {
+        ret.emplace_back(ps[i], -grad[i]);
+    }
+
+    return ret;
+}
+
+HEYOKA_END_NAMESPACE