From 67d117ef74be063f3586e4097f1994c3197668b9 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 29 Feb 2024 11:43:44 +0100 Subject: [PATCH 1/4] Expose the hessian() helper. --- heyoka/expose_expression.cpp | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index 2ed9c7b7..bcd787a0 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -621,6 +621,27 @@ void expose_expression(py::module_ &m) boost::numeric_cast(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; + }); // Copy/deepcopy. dtens_cl.def("__copy__", copy_wrapper); dtens_cl.def("__deepcopy__", deepcopy_wrapper, "memo"_a); From 78042e0d5744911fdd8a94cb72e852b0e3b0e6e4 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 29 Feb 2024 11:48:11 +0100 Subject: [PATCH 2/4] Fix Lagrange propagator notebook with proper link to diff_tensors(). --- doc/notebooks/lagrangian_propagator.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/notebooks/lagrangian_propagator.ipynb b/doc/notebooks/lagrangian_propagator.ipynb index 74e89ad6..8650b515 100644 --- a/doc/notebooks/lagrangian_propagator.ipynb +++ b/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 `:" + "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 `:" ] }, { @@ -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)." ] } ], From 63de525f0016cfc3c66cf90867ecc15c2dcfdb51 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Fri, 1 Mar 2024 08:50:55 +0100 Subject: [PATCH 3/4] Docs and testing for the hessian() helper. --- doc/changelog.rst | 3 ++ doc/notebooks/computing_derivatives.ipynb | 64 ++++++++++++++++++----- heyoka/_test_batch_integrator.py | 6 +-- heyoka/_test_dtens.py | 26 +++++++++ heyoka/_test_sympy.py | 4 +- heyoka/docstrings.cpp | 19 +++++++ heyoka/docstrings.hpp | 1 + heyoka/expose_expression.cpp | 31 ++++++----- 8 files changed, 122 insertions(+), 32 deletions(-) diff --git a/doc/changelog.rst b/doc/changelog.rst index e5109789..9e6cbe18 100644 --- a/doc/changelog.rst +++ b/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 `__). - Compiled functions now support multithreaded parallelisation for batched evaluations (`#168 `__). diff --git a/doc/notebooks/computing_derivatives.ipynb b/doc/notebooks/computing_derivatives.ipynb index fd404756..4fb9253e 100644 --- a/doc/notebooks/computing_derivatives.ipynb +++ b/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, diff --git a/heyoka/_test_batch_integrator.py b/heyoka/_test_batch_integrator.py index 88ecb505..a8c1370c 100644 --- a/heyoka/_test_batch_integrator.py +++ b/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)) diff --git a/heyoka/_test_dtens.py b/heyoka/_test_dtens.py index 46f8bbcc..eb8bbda1 100644 --- a/heyoka/_test_dtens.py +++ b/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 diff --git a/heyoka/_test_sympy.py b/heyoka/_test_sympy.py index 492a2831..452d4461 100644 --- a/heyoka/_test_sympy.py +++ b/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. diff --git a/heyoka/docstrings.cpp b/heyoka/docstrings.cpp index 53de747c..ab587c75 100644 --- a/heyoka/docstrings.cpp +++ b/heyoka/docstrings.cpp @@ -243,6 +243,25 @@ std::string dtens_jacobian() )"; } +std::string dtens_hessian() +{ + return R"(hessian(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 diff --git a/heyoka/docstrings.hpp b/heyoka/docstrings.hpp index db3adefc..c3882b85 100644 --- a/heyoka/docstrings.hpp +++ b/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(); diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index bcd787a0..d6a7160c 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -622,26 +622,29 @@ void expose_expression(py::module_ &m) }, 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)); + 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(); + // 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 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 ui = np.attr("triu_indices")(nargs); + arr1[ui] = h; - auto arr2 = arr1.attr("T").attr("copy")(); - auto di = np.attr("diag_indices")(nargs); + auto arr2 = arr1.attr("T").attr("copy")(); + auto di = np.attr("diag_indices")(nargs); - arr2[di] = 0.; + arr2[di] = 0.; - return arr1 + arr2; - }); + return arr1 + arr2; + }, + "component"_a, docstrings::dtens_hessian().c_str()); // Copy/deepcopy. dtens_cl.def("__copy__", copy_wrapper); dtens_cl.def("__deepcopy__", deepcopy_wrapper, "memo"_a); From fa6f809044c4728ed305527f92d1601d0013402d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Fri, 1 Mar 2024 09:02:58 +0100 Subject: [PATCH 4/4] Fix missing self argument in the docstrings of methods. --- heyoka/docstrings.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/heyoka/docstrings.cpp b/heyoka/docstrings.cpp index ab587c75..f0c81e6b 100644 --- a/heyoka/docstrings.cpp +++ b/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. @@ -245,7 +245,7 @@ std::string dtens_jacobian() std::string dtens_hessian() { - return R"(hessian(component: int) -> numpy.ndarray[expression] + return R"(hessian(self, component: int) -> numpy.ndarray[expression] Hessian of a component.