Merge pull request #171 from bluescarni/pr/heyoka_updates

hessian() helper

doc/changelog.rst

@@ -9,6 +9,9 @@ Changelog
 New
 ~~~
 
+- New convenience :func:`~heyoka.dtens.hessian()` method to fetch the Hessian
+  from a :class:`~heyoka.dtens` object
+  (`#171 <https://github.com/bluescarni/heyoka.py/pull/171>`__).
 - Compiled functions now support multithreaded parallelisation
   for batched evaluations
   (`#168 <https://github.com/bluescarni/heyoka.py/pull/168>`__).

doc/notebooks/computing_derivatives.ipynb

@@ -52,15 +52,7 @@
    "execution_count": 2,
    "id": "3d51a390-e9fa-4cad-8bdc-358b5ea3cf53",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[2024-02-12 07:49:15.252] [heyoka] [info] heyoka logger initialised\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "dt = hy.diff_tensors([x + hy.cos(y*z), hy.exp(x-z) + hy.log(y)],\n",
     "                     diff_args=[x, y, z], diff_order=1)"
@@ -384,6 +376,50 @@
     "dt.jacobian"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "419b5513-e301-4ef3-bccd-711ae0249497",
+   "metadata": {},
+   "source": [
+    "```{versionadded} 4.0.0\n",
+    "\n",
+    "```\n",
+    "\n",
+    "The convenience function {func}`~heyoka.dtens.hessian()` can be used to access the Hessian of a specific function component from a {class}`~heyoka.dtens` object. For instance:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "3a397e66-0b6f-406e-958c-d3d6e525f54f",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.0000000000000000,\n",
+       "        {({(2.0000000000000000 * y)} * {z**3.0000000000000000})},\n",
+       "        {({(3.0000000000000000 * z**2.0000000000000000)} * {y**2.0000000000000000})}],\n",
+       "       [{({(2.0000000000000000 * y)} * {z**3.0000000000000000})},\n",
+       "        {(2.0000000000000000 * {(x * z**3.0000000000000000)})},\n",
+       "        {({(3.0000000000000000 * z**2.0000000000000000)} * {({x} * {(2.0000000000000000 * y)})})}],\n",
+       "       [{({(3.0000000000000000 * z**2.0000000000000000)} * {y**2.0000000000000000})},\n",
+       "        {({(3.0000000000000000 * z**2.0000000000000000)} * {({x} * {(2.0000000000000000 * y)})})},\n",
+       "        {({(x * y**2.0000000000000000)} * {(3.0000000000000000 * {(2.0000000000000000 * z)})})}]],\n",
+       "      dtype=object)"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "dt = hy.diff_tensors([x * y**2 * z**3], [x, y, z], diff_order=2)\n",
+    "# Fetch the Hessian for the first (and only) component of the function.\n",
+    "dt.hessian(0)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "6090fe57-827e-435b-bfba-5c71267aad47",
@@ -410,7 +446,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 14,
    "id": "d2a60808-1a5c-44f0-9f60-a0dbd629006a",
    "metadata": {},
    "outputs": [
@@ -420,7 +456,7 @@
        "(x_1 * x_2 * x_3 * x_4 * x_5 * x_6 * x_7 * x_8)"
       ]
      },
-     "execution_count": 13,
+     "execution_count": 14,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -441,7 +477,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 15,
    "id": "118679d6-645f-49b5-af39-4ee15e926ada",
    "metadata": {},
    "outputs": [
@@ -458,7 +494,7 @@
        " (x_1 * x_2 * x_3 * x_4 * x_5 * x_6 * x_7)]"
       ]
      },
-     "execution_count": 14,
+     "execution_count": 15,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -680,7 +716,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,

doc/notebooks/lagrangian_propagator.ipynb

@@ -324,7 +324,7 @@
     "\n",
     "We can now proceed to the construction of an analytical expression for the state transition matrix (STM).\n",
     "\n",
-    "The STM is nothing but the Jacobian of the Lagrange propagator with respect to the initial conditions. In order to compute the derivatives, we employ the :func:`~heyoka.diff_tensors()` function as explained {ref}`here <computing_derivatives>`:"
+    "The STM is nothing but the Jacobian of the Lagrange propagator with respect to the initial conditions. In order to compute the derivatives, we employ the {func}`~heyoka.diff_tensors()` function as explained {ref}`here <computing_derivatives>`:"
    ]
   },
   {
@@ -457,7 +457,7 @@
    "source": [
     "we can indeed see that the the initial perturbation of $10^{-5}$ has been amplified by a factor of $3$ in the final state of $x$ as predicted by the STM.\n",
     "\n",
-    "Like the Lagrange propagator, the compiled function for the STM is also fully vectorised and SIMD-enabled. Note also that it is also possible via ``diff_tensors()`` to compute not only the Jacobian, but also higher-order tensors of derivatives (e.g., such as the Hessians)."
+    "Like the Lagrange propagator, the compiled function for the STM is also fully vectorised and SIMD-enabled. Note also that it is also possible via {func}`~heyoka.diff_tensors()` to compute not only the Jacobian, but also higher-order tensors of derivatives (e.g., such as the Hessians)."
    ]
   }
  ],

heyoka/_test_batch_integrator.py

@@ -355,9 +355,9 @@ def test_update_d_output(self):
 
             # Vector overload.
             with self.assertRaises(ValueError) as cm:
-                ta.update_d_output(np.array([0.3, 0.4, 0.45, 0.46], dtype=fp_t))[0] = (
-                    fp_t(0.5)
-                )
+                ta.update_d_output(np.array([0.3, 0.4, 0.45, 0.46], dtype=fp_t))[
+                    0
+                ] = fp_t(0.5)
 
             d_out2 = ta.update_d_output(np.array([0.3, 0.4, 0.45, 0.46], dtype=fp_t))
             self.assertEqual(d_out2.shape, (2, 4))

heyoka/_test_dtens.py

@@ -18,6 +18,32 @@ def test_gradient(self):
         dt = diff_tensors([x - y], [x, y])
         self.assertEqual(dt.gradient, [ex(1.0), ex(-1.0)])
 
+    def test_hessian(self):
+        from . import diff_tensors, make_vars, expression as ex, cfunc
+        import numpy as np
+
+        x, y, z = make_vars("x", "y", "z")
+
+        dt = diff_tensors([x * y**2 * z**3], [x, y, z], diff_order=2)
+
+        h = dt.hessian(component=0)
+
+        self.assertTrue(np.all(h - h.T) == ex(0.0))
+
+        cf = cfunc(h[np.triu_indices(3)].flatten(), [x, y, z])
+
+        xval, yval, zval = 1.0, 2.0, 3.0
+
+        out = cf([xval, yval, zval])
+
+        self.assertEqual(len(out), 6)
+        self.assertEqual(out[0], 0.0)
+        self.assertEqual(out[1], 2 * yval * zval * zval * zval)
+        self.assertEqual(out[2], 3 * yval * yval * zval * zval)
+        self.assertEqual(out[3], 2 * xval * zval * zval * zval)
+        self.assertEqual(out[4], 6 * xval * yval * zval * zval)
+        self.assertEqual(out[5], 6 * xval * yval * yval * zval)
+
     def test_jacobian(self):
         from . import diff_tensors, make_vars, expression as ex, diff_args
         import numpy as np

heyoka/_test_sympy.py

@@ -126,7 +126,9 @@ def test_number_conversion(self):
                     )
                     self.assertEqual(
                         from_sympy(Rational(2**expo + 1, 2**128)),
-                        expression(np.longdouble(2**expo + 1) / np.longdouble(2**128)),
+                        expression(
+                            np.longdouble(2**expo + 1) / np.longdouble(2**128)
+                        ),
                     )
 
         # Too high precision.

heyoka/docstrings.cpp

@@ -140,7 +140,7 @@ directly by the user.
 
 std::string dtens_get_derivatives()
 {
-    return R"(get_derivatives(order: int, component: Optional[int] = None) -> list[tuple[list[int], expression]]
+    return R"(get_derivatives(self, order: int, component: Optional[int] = None) -> list[tuple[list[int], expression]]
 
 Get the derivatives for the specified order and component.
 
@@ -196,7 +196,7 @@ std::string dtens_nargs()
 
 std::string dtens_index_of()
 {
-    return R"(index_of(vidx: list[int] | tuple[int, list[tuple[int, int]]]) -> int
+    return R"(index_of(self, vidx: list[int] | tuple[int, list[tuple[int, int]]]) -> int
 
 Get the position corresponding to the input indices vector.
 
@@ -243,6 +243,25 @@ std::string dtens_jacobian()
 )";
 }
 
+std::string dtens_hessian()
+{
+    return R"(hessian(self, component: int) -> numpy.ndarray[expression]
+
+Hessian of a component.
+
+.. versionadded:: 4.0.0
+
+This method will return the Hessian of the selected function component as a 2D array.
+
+:param component: the index of the function component whose Hessian will be returned.
+
+:return: the Hessian of the selected component.
+
+:raises ValueError: if *component* is invalid or if the derivative order is not at least 2.
+
+)";
+}
+
 std::string diff_tensors()
 {
     return R"(diff_tensors(func: list[expression], diff_args: list[expression] | diff_args, diff_order: int = 1) -> dtens

heyoka/docstrings.hpp

@@ -28,6 +28,7 @@ std::string dtens_index_of();
 std::string dtens_args();
 std::string dtens_gradient();
 std::string dtens_jacobian();
+std::string dtens_hessian();
 std::string dtens_nouts();
 std::string dtens_nargs();
 

heyoka/expose_expression.cpp

@@ -621,6 +621,30 @@ void expose_expression(py::module_ &m)
                                                          boost::numeric_cast<py::ssize_t>(dt.get_nargs())});
         },
         docstrings::dtens_jacobian().c_str());
+    // Hessian.
+    dtens_cl.def(
+        "hessian",
+        [](const hey::dtens &dt, std::uint32_t component) {
+            py::list h = py::cast(dt.get_hessian(component));
+
+            // Reconstruct the Hessian from its representation as a list
+            // of component of the upper triangular part.
+            const auto nargs = dt.get_nargs();
+
+            auto np = py::module_::import("numpy");
+            auto arr1 = np.attr("full")(py::make_tuple(nargs, nargs), 0., "dtype"_a = "object");
+
+            auto ui = np.attr("triu_indices")(nargs);
+            arr1[ui] = h;
+
+            auto arr2 = arr1.attr("T").attr("copy")();
+            auto di = np.attr("diag_indices")(nargs);
+
+            arr2[di] = 0.;
+
+            return arr1 + arr2;
+        },
+        "component"_a, docstrings::dtens_hessian().c_str());
     // Copy/deepcopy.
     dtens_cl.def("__copy__", copy_wrapper<hey::dtens>);
     dtens_cl.def("__deepcopy__", deepcopy_wrapper<hey::dtens>, "memo"_a);