Add helpers in dtens to easily fetch gradient and Jacobian.

include/heyoka/expression.hpp

@@ -501,6 +501,8 @@ class HEYOKA_DLL_PUBLIC dtens
 
     [[nodiscard]] subrange get_derivatives(std::uint32_t, std::uint32_t) const;
     [[nodiscard]] subrange get_derivatives(std::uint32_t) const;
+    [[nodiscard]] std::vector<expression> get_gradient() const;
+    [[nodiscard]] std::vector<expression> get_jacobian() const;
 };
 
 HEYOKA_DLL_PUBLIC std::ostream &operator<<(std::ostream &, const dtens &);

src/expression_diff.cpp

@@ -12,6 +12,7 @@
 #include <cstdint>
 #include <deque>
 #include <functional>
+#include <iterator>
 #include <map>
 #include <memory>
 #include <ostream>
@@ -1423,6 +1424,44 @@ dtens::subrange dtens::get_derivatives(std::uint32_t component, std::uint32_t or
     return subrange{b, e};
 }
 
+std::vector<expression> dtens::get_gradient() const
+{
+    if (get_nouts() != 1u) {
+        throw std::invalid_argument(fmt::format("The gradient can be requested only for a function with a single "
+                                                "output, but the number of outputs is instead {}",
+                                                get_nouts()));
+    }
+
+    if (get_order() == 0u) {
+        throw std::invalid_argument("First-order derivatives are not available");
+    }
+
+    const auto sr = get_derivatives(0, 1);
+    std::vector<expression> retval;
+    retval.reserve(get_nvars());
+    std::transform(sr.begin(), sr.end(), std::back_inserter(retval), [](const auto &p) { return p.second; });
+
+    return retval;
+}
+
+std::vector<expression> dtens::get_jacobian() const
+{
+    if (get_nouts() == 0u) {
+        throw std::invalid_argument("Cannot return the Jacobian of a function with no outputs");
+    }
+
+    if (get_order() == 0u) {
+        throw std::invalid_argument("First-order derivatives are not available");
+    }
+
+    const auto sr = get_derivatives(1);
+    std::vector<expression> retval;
+    retval.reserve(get_nvars());
+    std::transform(sr.begin(), sr.end(), std::back_inserter(retval), [](const auto &p) { return p.second; });
+
+    return retval;
+}
+
 std::uint32_t dtens::get_nvars() const
 {
     if (p_impl->m_map.empty()) {

test/expression_diff_tensors.cpp

@@ -628,3 +628,74 @@ TEST_CASE("speelpenning complexity")
         s.compile();
     }
 }
+
+TEST_CASE("gradient")
+{
+    using Catch::Matchers::Message;
+
+    auto [x, y] = make_vars("x", "y");
+
+    // Error modes first.
+    {
+        dtens dt;
+        REQUIRE_THROWS_MATCHES(dt.get_gradient(), std::invalid_argument,
+                               Message("The gradient can be requested only for a function with a single "
+                                       "output, but the number of outputs is instead 0"));
+    }
+
+    {
+        auto dt = diff_tensors({x + y, x - y});
+        REQUIRE_THROWS_MATCHES(dt.get_gradient(), std::invalid_argument,
+                               Message("The gradient can be requested only for a function with a single "
+                                       "output, but the number of outputs is instead 2"));
+    }
+
+    {
+        auto dt = diff_tensors({x + y}, kw::diff_order = 0u);
+        REQUIRE_THROWS_MATCHES(dt.get_gradient(), std::invalid_argument,
+                               Message("First-order derivatives are not available"));
+    }
+
+    auto dt = diff_tensors({x + y}, kw::diff_order = 1u);
+    auto grad = dt.get_gradient();
+    REQUIRE(grad == std::vector{1_dbl, 1_dbl});
+}
+
+TEST_CASE("jacobian")
+{
+    using Catch::Matchers::Message;
+
+    auto [x, y] = make_vars("x", "y");
+
+    // Error modes first.
+    {
+        dtens dt;
+
+        REQUIRE_THROWS_MATCHES(dt.get_jacobian(), std::invalid_argument,
+                               Message("Cannot return the Jacobian of a function with no outputs"));
+    }
+
+    {
+        auto dt = diff_tensors({x + y}, kw::diff_order = 0u);
+        REQUIRE_THROWS_MATCHES(dt.get_jacobian(), std::invalid_argument,
+                               Message("First-order derivatives are not available"));
+    }
+
+    {
+        auto dt = diff_tensors({x + y}, kw::diff_order = 1u);
+        auto jac = dt.get_jacobian();
+        REQUIRE(jac == std::vector{1_dbl, 1_dbl});
+    }
+
+    {
+        auto dt = diff_tensors({x + y, x - y}, kw::diff_order = 1u);
+        auto jac = dt.get_jacobian();
+        REQUIRE(jac == std::vector{1_dbl, 1_dbl, 1_dbl, -1_dbl});
+    }
+
+    {
+        auto dt = diff_tensors({x + y, x - y, -x - y}, kw::diff_order = 1u);
+        auto jac = dt.get_jacobian();
+        REQUIRE(jac == std::vector{1_dbl, 1_dbl, 1_dbl, -1_dbl, -1_dbl, -1_dbl});
+    }
+}