diff --git a/README.md b/README.md index 050fbfb1..c453c727 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,8 @@ compilation via [LLVM](https://llvm.org/). Notable features include: * accurate and reliable event detection, * builtin support for analytical mechanics - bring your own Lagrangians/Hamiltonians and let heyoka.py formulate and solve the equations of motion, +* builtin support for high-order variational equations - compute not only the solution, + but also its partial derivatives, * builtin support for machine learning applications via neural network models, * the ability to maintain machine precision accuracy over tens of billions of timesteps, diff --git a/doc/api_exsys.rst b/doc/api_exsys.rst index faed6fc3..6cc31324 100644 --- a/doc/api_exsys.rst +++ b/doc/api_exsys.rst @@ -34,6 +34,7 @@ Attributes :toctree: autosummary_generated par + time Enums ----- diff --git a/doc/api_var_ode_sys.rst b/doc/api_var_ode_sys.rst new file mode 100644 index 00000000..0aa44d2e --- /dev/null +++ b/doc/api_var_ode_sys.rst @@ -0,0 +1,24 @@ +.. _api_var_ode_sys: + +Variational ODE systems +======================= + +.. currentmodule:: heyoka + +Classes +------- + +.. autosummary:: + :toctree: autosummary_generated + :template: custom-class-template.rst + + var_ode_sys + +Enums +----- + +.. autosummary:: + :toctree: autosummary_generated + :template: custom-enum-template.rst + + var_args diff --git a/doc/basic_tutorials.rst b/doc/basic_tutorials.rst index 8f9a82ff..7f598ee5 100644 --- a/doc/basic_tutorials.rst +++ b/doc/basic_tutorials.rst @@ -34,3 +34,4 @@ Basic notebooks/Non-autonomous systems.ipynb notebooks/Dense output.ipynb notebooks/Event detection + notebooks/var_ode_sys.ipynb diff --git a/doc/changelog.rst b/doc/changelog.rst index b1ebf5b1..a86a84f0 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -9,6 +9,11 @@ Changelog New ~~~ +- Add support for variational ODE systems and Taylor map computation + (`#177 `__). +- Add thermonets: neural, differentiable, high-performance + models for the Earth's thermosphere density + (`#176 `__). - Add a vectorised implementation of ``diff()`` (`#173 `__). diff --git a/doc/index.md b/doc/index.md index c1d1cc90..acb9c2dc 100644 --- a/doc/index.md +++ b/doc/index.md @@ -16,6 +16,8 @@ compilation via [LLVM](https://llvm.org/). Notable features include: - accurate and reliable event detection, - builtin support for analytical mechanics - bring your own Lagrangians/Hamiltonians and let heyoka.py formulate and solve the equations of motion, +- builtin support for high-order variational equations - compute not only the solution, + but also its partial derivatives, - builtin support for machine learning applications via neural network models, - the ability to maintain machine precision accuracy over tens of billions of timesteps, @@ -107,6 +109,7 @@ examples_others api_exsys api_integrators +api_var_ode_sys api_lagham api_model ``` diff --git a/doc/notebooks/Batch mode overview.ipynb b/doc/notebooks/Batch mode overview.ipynb index e80155d3..3773df22 100644 --- a/doc/notebooks/Batch mode overview.ipynb +++ b/doc/notebooks/Batch mode overview.ipynb @@ -22,10 +22,10 @@ "It is important to emphasise that batch mode does not reduce\n", "the CPU time required to integrate a system of ODEs. Rather, as a fine-grained\n", "form of data parallelism, batch mode allows to integrate multiple ODE systems in parallel\n", - "at no additional cost, and it is thus most useful when the need arise\n", + "at (almost) no additional cost, and it is thus most useful when the need arise\n", "to integrate the same ODE system with different initial conditions and parameters.\n", "\n", - "Although batch mode can in principle be used with all floating-point types supported\n", + "Although batch mode can in principle be used with all the fundamental C++ floating-point types supported\n", "by heyoka.py, in practice at this time no CPU provides SIMD instructions for extended-precision\n", "datatypes. Thus, heyoka.py currently enables batch mode only for single and double precision\n", "computations. Specifically, in this tutorial we will consider the application of batch mode to\n", diff --git a/doc/notebooks/Customising the adaptive integrator.ipynb b/doc/notebooks/Customising the adaptive integrator.ipynb index 178cf959..e6a0a51e 100644 --- a/doc/notebooks/Customising the adaptive integrator.ipynb +++ b/doc/notebooks/Customising the adaptive integrator.ipynb @@ -128,6 +128,8 @@ "id": "national-investing", "metadata": {}, "source": [ + "(ta_compact_mode)=\n", + "\n", "Compact mode\n", "------------\n", "\n", @@ -290,7 +292,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/doc/notebooks/The variational equations.ipynb b/doc/notebooks/The variational equations.ipynb index 63b0966e..51f6c1b3 100644 --- a/doc/notebooks/The variational equations.ipynb +++ b/doc/notebooks/The variational equations.ipynb @@ -5,12 +5,17 @@ "id": "excited-uganda", "metadata": {}, "source": [ + "```{note}\n", + "\n", + "Starting with version 5.0.0, heyoka.py features builtin support for the formulation of the [variational equations](./var_ode_sys.ipynb). This tutorial is now deprecated and it will be removed in a future version.\n", + "\n", + "```\n", + "\n", "The variational equations\n", "=========================\n", "\n", "In this tutorial, we will show how it is possible to exploit heyoka.py's [expression system](<./The expression system.ipynb>) to implement first-order variational equations for a system of ODEs.\n", "\n", - "\n", "First, let's start with a brief recap on the variational equations. For simplicity, we will consider an autonomous ODE system in the variables $x$ and $y$ (the extension to more variables just complicates the notation):\n", "\n", "$$\n", @@ -292,7 +297,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.13" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/doc/notebooks/ensemble_mode.ipynb b/doc/notebooks/ensemble_mode.ipynb index 0caffb2b..ae868f47 100644 --- a/doc/notebooks/ensemble_mode.ipynb +++ b/doc/notebooks/ensemble_mode.ipynb @@ -9,6 +9,10 @@ "\n", "# Ensemble propagations\n", "\n", + "```{versionadded} 0.17.0\n", + "\n", + "```\n", + "\n", "Starting with version 0.17.0, heyoka.py offers support for\n", "*ensemble propagations*. In ensemble mode, multiple distinct\n", "instances of the same ODE system are integrated in parallel,\n", diff --git a/doc/notebooks/var_ode_sys.ipynb b/doc/notebooks/var_ode_sys.ipynb new file mode 100644 index 00000000..347a92c3 --- /dev/null +++ b/doc/notebooks/var_ode_sys.ipynb @@ -0,0 +1,785 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e916542d-c443-4403-bd35-e56122b2823b", + "metadata": {}, + "source": [ + "(var_ode_sys)=\n", + "\n", + "# Variational ODEs\n", + "\n", + "```{versionadded} 5.0.0\n", + "\n", + "```\n", + "\n", + "Consider a system of differential equations in the standard form\n", + "\n", + "$$\n", + "\\frac{d\\boldsymbol{x}}{dt} = f\\left(\\boldsymbol{x}, \\boldsymbol{\\alpha}, t\\right),\n", + "$$\n", + "\n", + "where $\\boldsymbol{x}$ is the vector of state variables and $\\boldsymbol{\\alpha}$ a vector of parameters. For a given set of initial conditions $\\boldsymbol{x}_0$ at time $t_0$, the solution of this system will be\n", + "\n", + "$$\n", + "\\boldsymbol{x} = \\boldsymbol{x}\\left(t, \\boldsymbol{x}_0, t_0, \\boldsymbol{\\alpha} \\right).\n", + "$$\n", + "\n", + "When solving numerically initial-value problems, it is often useful to compute not only the solution, but also its partial derivatives with respect to the initial conditions and/or the parameters. The derivatives with respect to the initial conditions, for instance, are needed for the computation of [chaos indicators](https://en.wikipedia.org/wiki/Lyapunov_exponent) and for [uncertainty propagation](https://en.wikipedia.org/wiki/Propagation_of_uncertainty), and they can also be used to propagate a small neighborhood in phase space around the initial conditions. The derivatives with respect to the parameters of the system are required when formulating optimisation and inversion problems such as orbit determination, trajectory optimisation and training of neural networks in [neural ODEs](./NeuralODEs.ipynb).\n", + "\n", + "There are two main methods for the computation of the partial derivatives. The first one is based on the application of automatic differentiation (AD) techniques directly to the numerical integration algorithm. This can be done either by replacing the algebra of floating-point numbers with the algebra of (generalised) [dual numbers](https://en.wikipedia.org/wiki/Dual_number) (aka truncated Taylor polynomials), or via [differentiable programming](https://en.wikipedia.org/wiki/Differentiable_programming) techniques. The former approach is used by libraries such as [pyaudi](https://github.com/darioizzo/audi), [desolver](https://github.com/Microno95/desolver) and [TaylorIntegration.jl](https://docs.sciml.ai/TaylorIntegration/stable/jet_transport/), while differentiable programming is popular in the machine learning community with projects such as [PyTorch](https://pytorch.org/), [JAX](https://jax.readthedocs.io/en/latest/) and [TensorFlow](https://www.tensorflow.org/). Differentiable programming is also popular in the [Julia programming language](https://en.wikipedia.org/wiki/Julia_(programming_language)) community.\n", + "\n", + "The second method is based on the formulation of the *variational equations*, that is, differential equations satisfied by the partial derivatives which are added to and solved together with the original ODEs. For instance, we can formulate differential equations for the first-order derivatives with respect to the initial conditions via elementary calculus:\n", + "\n", + "$$\n", + "\\frac{d}{dt}\\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0} = \\frac{\\partial }{\\partial \\boldsymbol{x}_0} \\frac{d \\boldsymbol{x}}{dt} = \\frac{\\partial f}{\\partial \\boldsymbol{x}_0} = \\frac{\\partial f}{\\partial \\boldsymbol{x}} \\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0}.\n", + "$$\n", + "\n", + "The variational ODE system then reads\n", + "\n", + "$$\n", + "\\begin{cases}\n", + "\\frac{d\\boldsymbol{x}}{dt} = f\\left(\\boldsymbol{x}, \\boldsymbol{\\alpha}, t\\right) \\\\\n", + "\\frac{d}{dt}\\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0} = \\frac{\\partial f}{\\partial \\boldsymbol{x}} \\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0}\n", + "\\end{cases},\n", + "$$\n", + "\n", + "and the original state vector $\\boldsymbol{x}$ has been extended to include the variational state variables $\\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0}$.\n", + "\n", + "heyoka.py adopts the variational approach for the computation of the partial derivatives, supporting the formulation of variational ODEs at arbitrary differentiation orders and with respect to any combination of initial conditions, parameters and initial time. In this tutorial, we will explore this feature and show a couple of interesting use cases.\n", + "\n", + "Before beginning, however, let us point out for clarity (and for the benefit of the search engines indexing this page) that in the scientific literature there is a bewildering variety of different names and monikers used when discussing partial derivatives of ODEs and their applications. Here is a (partial) list:\n", + "\n", + "- in the astrodynamics community, the term *differential algebra* is often used to refer to the computation of partial derivatives via truncated Taylor polynomials (e.g., see [this paper](https://link.springer.com/article/10.1007/s10569-010-9283-5)). The term actually originates from the community of beam physics, where it has been used in the context of the theoretical modelling of particle accelerators since the 90s (e.g., see [this review](https://www.bmtdynamics.org/pub/papers/DAHAPE12/DAHAPE12.pdf));\n", + "- in the mathematical community, the term *jet transport* is sometimes used to refer to the propagation of a small neighborhood in phase space around the initial conditions via the Taylor series constructed form the partial derivatives (e.g., see [this paper](http://www.maia.ub.es/~angel/varis/granada09.pdf)). In heyoka.py, we refer to a similar idea as {ref}`Taylor map evaluation `;\n", + "- in the Julia programming language community, the term *local sensitivity analysis* refers to the computation of the partial derivatives via the variational equations, while *discrete sensitivity analysis* refers to the computation of the partial derivatives by directly differentiating the numerical method's steps (e.g., see [this review](https://arxiv.org/abs/1812.01892));\n", + "- in the space engineering community, the term *state transition tensors* is sometimes used to indicate the generalisations of the [state transition matrix](https://en.wikipedia.org/wiki/State-transition_matrix) (which in turn is built from the first-order partial derivatives) to higher differentiation orders.\n", + "\n", + "## Constructing a variational ODE system\n", + "\n", + "Let us begin with the definition of a simple ODE system:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "2c186e41-5526-4b7d-a11a-0ce38332c9ea", + "metadata": {}, + "outputs": [], + "source": [ + "import heyoka as hy\n", + "\n", + "# Create the symbolic variables x and v.\n", + "x, v = hy.make_vars(\"x\", \"v\")\n", + "\n", + "# Create an ODE system.\n", + "sys = [(x, v), (v, hy.cos(hy.time) - hy.par[0] * v - hy.sin(x))]" + ] + }, + { + "cell_type": "markdown", + "id": "1d9126e2-93b3-49a2-8fdb-33e31a161e35", + "metadata": {}, + "source": [ + "This is the forced damped pendulum system already considered in [another tutorial](<./Non-autonomous systems.ipynb>), where we have introduced the air friction coefficient as the [runtime parameter](<./ODEs with parameters.ipynb>) ``par[0]``.\n", + "\n", + "We then proceed to create a {class}`~heyoka.var_ode_sys`:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "ba980238-26c2-417c-817f-92fb2dbe3ec1", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the variational ODE system.\n", + "vsys = hy.var_ode_sys(sys, hy.var_args.vars, order=2)" + ] + }, + { + "cell_type": "markdown", + "id": "bab348cb-ca8b-49f0-bc40-0c44ff091521", + "metadata": {}, + "source": [ + "Upon construction, {class}`~heyoka.var_ode_sys` formulates the variational equations for the input ODE system ``sys`` up to the specified differentiation ``order``. The second argument specifies with respect to which quantities the variational equations are formulated. In this case, we used the ``vars`` enumerator of the {class}`~heyoka.var_args` enum: this means that the variational equations will be formulated with respect to the initial conditions of the state variables. In a completely equivalent manner, we could have written instead:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "220429e4-fd82-4495-873e-a08ae8c09a9a", + "metadata": {}, + "outputs": [], + "source": [ + "vsys = hy.var_ode_sys(sys, [x, v], order=2)" + ] + }, + { + "cell_type": "markdown", + "id": "6c86a6e8-20d2-4d78-8846-46ed6628e823", + "metadata": {}, + "source": [ + "In this case, instead of a {class}`~heyoka.var_args` enumerator, we passed an explicit list of state variables with respect to whose initial conditions we want to formulate the variational equations. Please refer to the documentation of {class}`~heyoka.var_ode_sys` for an exhaustive explanation of what can be passed as second argument to the constructor.\n", + "\n", + "Let us explore a bit the {class}`~heyoka.var_ode_sys` class. First of all, we can access the variational system of equations:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "16e3bc0a-e7d5-499c-88bd-e17aaa6c1657", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(x, v),\n", + " (v, ((cos(t) - (p0 * v)) - sin(x))),\n", + " (∂[(0, 1)]x, ∂[(0, 1)]v),\n", + " (∂[(1, 1)]x, ∂[(1, 1)]v),\n", + " (∂[(0, 1)]v, (-(∂[(0, 1)]x * cos(x)) - (∂[(0, 1)]v * p0))),\n", + " (∂[(1, 1)]v, (-(∂[(1, 1)]x * cos(x)) - (∂[(1, 1)]v * p0))),\n", + " (∂[(0, 2)]x, ∂[(0, 2)]v),\n", + " (∂[(0, 1), (1, 1)]x, ∂[(0, 1), (1, 1)]v),\n", + " (∂[(1, 2)]x, ∂[(1, 2)]v),\n", + " (∂[(0, 2)]v,\n", + " (-(∂[(0, 2)]v * p0) - (((∂[(0, 1)]x * -sin(x)) * ∂[(0, 1)]x) + (∂[(0, 2)]x * cos(x))))),\n", + " (∂[(0, 1), (1, 1)]v,\n", + " (-(∂[(0, 1), (1, 1)]v * p0) - (((∂[(0, 1)]x * -sin(x)) * ∂[(1, 1)]x) + (∂[(0, 1), (1, 1)]x * cos(x))))),\n", + " (∂[(1, 2)]v,\n", + " (-(∂[(1, 2)]v * p0) - (((∂[(1, 1)]x * -sin(x)) * ∂[(1, 1)]x) + (∂[(1, 2)]x * cos(x)))))]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.sys" + ] + }, + { + "cell_type": "markdown", + "id": "c969185b-7df4-4f95-9ab1-04af1205abdd", + "metadata": {}, + "source": [ + "The first two equations are from the original system of ODEs, while the remaining ones are the variational equations. The names of the variational variables begin with the $\\partial$ symbol, followed by a sparse multiindex encoding of the differentiation indices. For instance, the variational variable ``∂[(0, 1)]x`` is the first-order derivative of $x$ with respect to the first variational argument, that is, $\\frac{\\partial x}{\\partial x_0}$. Similarly, ``∂[(0, 1), (1, 1)]x`` is the second order derivative of $x$ with respect to both variational arguments, that is, $\\frac{\\partial^2 x}{\\partial x_0 \\partial y_0}$. The ordering of the variational equations follows the same scheme explained in the tutorial about [computing derivatives](<./computing_derivatives.ipynb>).\n", + "\n", + "We can also query other properties of ``vsys``. For instance, the differentiation order:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "e4715f11-f4c3-4953-8d6e-7d01cd14d32a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.order" + ] + }, + { + "cell_type": "markdown", + "id": "27e8ac8e-1d81-4309-84d7-f18f16fc2d06", + "metadata": {}, + "source": [ + "The number of state variables in the original ODE system:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "2fad08b8-23d0-4572-a74c-edc42ff93c3b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.n_orig_sv" + ] + }, + { + "cell_type": "markdown", + "id": "f5d76bce-6e58-4e87-b5ba-98fdb953dcc1", + "metadata": {}, + "source": [ + "And the list of arguments with respect to which the variational equations are formulated:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9369b6f1-9b2c-4e7e-b70e-21679dd62fda", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[x, v]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.vargs" + ] + }, + { + "cell_type": "markdown", + "id": "276d5cc7-615c-44fe-b08d-48a9cee19fcb", + "metadata": {}, + "source": [ + "## Constructing a variational integrator\n", + "\n", + "After the construction of a variational ODE system, we are now ready to construct a variational integrator. We can do this by simply passing the variational ODE system (instead of the original, non-variational ODE system) as first input argument to the constructor:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8c51dcf8-a089-42c2-9b84-5e2d5c22da24", + "metadata": {}, + "outputs": [], + "source": [ + "# Construct a variational integrator.\n", + "ta_var = hy.taylor_adaptive(vsys, [.2, .3], pars=[.4], compact_mode=True)" + ] + }, + { + "cell_type": "markdown", + "id": "5925b0ce-1b04-4174-86ef-593c7ffade74", + "metadata": {}, + "source": [ + "Note how we constructed the integrator with {ref}`compact mode ` enabled: the formulation of the variational equations, especially at high differentiation orders, greatly increases the size of the symbolic expressions that need to be just-in-time compiled during the creation of the integrator. By enabling compact mode, we keep the compilation time at manageable levels.\n", + "\n", + "Let us inspect the integrator:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "844c0ae8-9f32-4d09-a900-8098301863ff", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "C++ datatype : double\n", + "Tolerance : 2.220446049250313e-16\n", + "High accuracy : false\n", + "Compact mode : true\n", + "Taylor order : 20\n", + "Dimension : 12\n", + "Time : 0\n", + "State : [0.2, 0.3, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0]\n", + "Parameters : [0.4]\n", + "Variational order : 2" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var" + ] + }, + { + "cell_type": "markdown", + "id": "afc46400-dc3a-4e1d-bae6-aed5771d46a6", + "metadata": {}, + "source": [ + "The screen output informs us that ``ta_var`` is a variational integrator of order 2. We can also see that, although on construction we passed the initial conditions only for the ``x`` and ``v`` state variables, the integrator automatically set up appropriate initial conditions for the variational variables. Indeed, with respect to a non-variational integrator, the state vector has been augmented to store also the values of the variational variables:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "81cbba40-07fb-45b1-9622-29e9c966b273", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.2, 0.3, 1. , 0. , 0. , 1. , 0. , 0. , 0. , 0. , 0. , 0. ])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.state" + ] + }, + { + "cell_type": "markdown", + "id": "8e5b8bd0-bf91-41f5-baf6-c1cb2d2e2e7a", + "metadata": {}, + "source": [ + "Alternatively, instead of relying on the integrator to automatically set up the initial conditions of the variational variables, we could also pass a full 12-elements vector of initial conditions in input.\n", + "\n", + "We are now ready to proceed to our first variational integration:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "afa9f717-bc57-4f19-9a25-198bff3cdd14", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.11757215, -1.24940656, -0.425478 , 0.41983649, -0.19171818,\n", + " -0.51871994, 0.27771857, 0.22392433, 0.60414865, -0.1122785 ,\n", + " -0.12689065, 0.12443791])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Propagate until t=3.\n", + "ta_var.propagate_until(3.)\n", + "\n", + "# Print the full state vector.\n", + "ta_var.state" + ] + }, + { + "cell_type": "markdown", + "id": "a51db4d7-5c9f-446f-ae1f-eb9fc759c94d", + "metadata": {}, + "source": [ + "That is all fine and good, but how do we fetch the values of the derivatives we are interested in at the end of an integration? As mentioned earlier, the partial derivatives are ordered in the state vector following the same criterion explained in the tutorial about [computing derivatives](<./computing_derivatives.ipynb>): first by total order of differentiation, then by component (i.e., the derivatives of ``x`` precede the derivatives of ``v``) and finally by reverse lexicographic order with respect to the differentiation multiindices. However, navigating by hand this ordering scheme can be complicated, especially at high differentiation orders.\n", + "\n", + "Variational integrators provide a couple of methods that facilitate the task of locating specific derivatives in the state vector. The first helper is ``get_vslice()``. This method takes as input a differentiation order and, optionally, a component index, and returns a {class}`slice` into the state vector corresponding to the range of indices for the requested derivatives. Let us see a couple of examples:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "a78531a3-a1f9-4a60-bc21-5c37ab13f6ad", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "slice(6, 12, None)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Fetch the range of all order-2 derivatives.\n", + "sl = ta_var.get_vslice(order=2)\n", + "sl" + ] + }, + { + "cell_type": "markdown", + "id": "bfd3b6b5-8cd6-4c92-be86-d1125cee17bf", + "metadata": {}, + "source": [ + "That is, the order-2 derivatives are between indices 6 and 12 in the state vector. We can use ``sl`` to index directly into the state vector:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "3fed6708-3e5e-44c2-871a-0237d005937e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.27771857, 0.22392433, 0.60414865, -0.1122785 , -0.12689065,\n", + " 0.12443791])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Fetch all order-2 derivatives from the state vector.\n", + "ta_var.state[sl]" + ] + }, + { + "cell_type": "markdown", + "id": "455cbfba-964d-4cbe-b6c4-efdde7f85ad2", + "metadata": {}, + "source": [ + "If we are interested only in the order-2 derivatives of ``v``, we can pass the additional ``component`` keyword argument:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "ff3ae54f-455d-44d3-b166-a0a607434d6a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([-0.1122785 , -0.12689065, 0.12443791])" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Fetch the range of the order-2 derivatives for v.\n", + "# component=1 means the second original state variable,\n", + "# i.e., v (component=0 would fetch the derivatives for x).\n", + "sl = ta_var.get_vslice(order=2, component=1)\n", + "\n", + "# Fetch the order-2 derivatives for v.\n", + "ta_var.state[sl]" + ] + }, + { + "cell_type": "markdown", + "id": "6fa1d1cd-840e-491c-aefb-535518c7bd2e", + "metadata": {}, + "source": [ + "Often fetching the values of the derivatives is not enough, and we also need to access the differentiation multiindices associated to each derivative. In order to do this, we can use the ``get_mindex()`` method, which takes in input a single index into the state vector and returns the corresponding differentiation multiindex.\n", + "\n", + "Let us see a couple of examples:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "5ae39909-0554-4def-b04b-4ed69ade36b7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[0, 0, 0]" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.get_mindex(i=0)" + ] + }, + { + "cell_type": "markdown", + "id": "1018d0c7-e7f3-463a-8e88-618d1600be31", + "metadata": {}, + "source": [ + "At ``i=0`` in the state vector we have the order-0 derivative of the first state variable - that is, a complicated way of saying that we have the current value of ``x``. Let us see a more useful example:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "064adc95-4ed0-4c6e-b927-57f8a607483c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Multiindex: [0, 2, 0], derivative value: 0.27771856625922825\n", + "Multiindex: [0, 1, 1], derivative value: 0.2239243290399873\n", + "Multiindex: [0, 0, 2], derivative value: 0.6041486469599123\n", + "Multiindex: [1, 2, 0], derivative value: -0.1122785039896587\n", + "Multiindex: [1, 1, 1], derivative value: -0.12689065229074553\n", + "Multiindex: [1, 0, 2], derivative value: 0.12443790781972988\n" + ] + } + ], + "source": [ + "# Fetch the range of all order-2 derivatives.\n", + "sl = ta_var.get_vslice(order=2)\n", + "\n", + "# Print the multiindices and associated values.\n", + "for idx in range(sl.start, sl.stop):\n", + " print(f\"Multiindex: {ta_var.get_mindex(idx)}, derivative value: {ta_var.state[idx]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "d2ae2f2e-3b67-47be-a863-b17f44fb26c2", + "metadata": {}, + "source": [ + "Recall that in a multiindex the first number refers to the component index (i.e., 0 for ``x`` and 1 for ``v``), while the remaining indices refer to the differentiation orders with respect to the variational arguments.\n", + "\n", + "(taylor_map)=\n", + "## Taylor map evaluation\n", + "\n", + "One of the most useful applications of the variational equations is the ability to compute how a small perturbation on the initial conditions and/or parameters of the system affects the current state of the system, and to do it quickly (i.e., without having to repeat the numerical integration with the updated initial conditions/parameters). This is accomplished by using the values of the partial derivatives to construct and evaluate the multivariate Taylor series of the solution around the original initial conditions/parameters of the system. This approach, when applied to perturbations on the initial conditions, is sometimes called *jet transport* in the literature. Here, more generally, we will call it *evaluation of the Taylor map*.\n", + "\n", + "Variational integrators provide a specific method called ``eval_taylor_map()`` to construct and evaluate the Taylor map. Let us see a simple example. We begin by re-creating from scratch our variational integrator:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "229a9a4e-bf5a-4b57-a599-e4571ebd7507", + "metadata": {}, + "outputs": [], + "source": [ + "ta_var = hy.taylor_adaptive(vsys, [.2, .3], pars=[.4], compact_mode=True)" + ] + }, + { + "cell_type": "markdown", + "id": "b09de854-aa19-4f08-a36d-66b926bd1d6b", + "metadata": {}, + "source": [ + "We define two small displacements on the state variables ``x`` and ``v``:" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "7c2d30d2-17cd-4ab4-889a-8bcb72bb9b11", + "metadata": {}, + "outputs": [], + "source": [ + "dx = 1e-4\n", + "dv = -2e-4" + ] + }, + { + "cell_type": "markdown", + "id": "eb9e6427-fb2c-4f11-91b2-85d9d45577b7", + "metadata": {}, + "source": [ + "And we create a non-variational integrator with displaced initial conditions with respect to the variational one:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "1ab3b27e-3506-47e0-a0bb-35cd2076238f", + "metadata": {}, + "outputs": [], + "source": [ + "# Non-variational integrator with displaced\n", + "# initial conditions.\n", + "ta = hy.taylor_adaptive(sys, [.2 + dx, .3 + dv], pars=[.4], compact_mode=True)" + ] + }, + { + "cell_type": "markdown", + "id": "8fa77407-e47a-4338-b65e-f314a4a7afdd", + "metadata": {}, + "source": [ + "Next, we propagate both integrators up to $t=3$:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "42a8ca75-4eca-4da2-b33e-838085fd497f", + "metadata": {}, + "outputs": [], + "source": [ + "ta_var.propagate_until(3.)\n", + "ta.propagate_until(3.);" + ] + }, + { + "cell_type": "markdown", + "id": "80f7bd5e-606b-425c-bfe8-e3b409a1975b", + "metadata": {}, + "source": [ + "Clearly, the values of ``x`` and ``v`` will differ in the two integrators due to the different initial conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "fa2c2e21-81f3-4e2a-a219-296b940bb3b9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Non-variational state: [ 0.11744564 -1.24932198]\n", + "Variational state : [ 0.11757215 -1.24940656]\n" + ] + } + ], + "source": [ + "print(f\"Non-variational state: {ta.state}\")\n", + "print(f\"Variational state : {ta_var.state[:2]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "b102f4fd-63dc-45a6-81c1-1cf925aa8b15", + "metadata": {}, + "source": [ + "We can now use the ``eval_taylor_map()`` method on the variational integrator to compute the effect of the displacements ``dx`` and ``dv`` on the state of the system:" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "b7f778cc-db85-43aa-bf29-5445d6170958", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.11744564, -1.24932198])" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.eval_taylor_map([dx, dv])" + ] + }, + { + "cell_type": "markdown", + "id": "a3f8b0d2-0446-4dd5-a8fb-0539863ed161", + "metadata": {}, + "source": [ + "``eval_taylor_map()`` takes in input a vector of displacements (one for each variational argument), and computes and evaluates the Taylor map, returning a reference to an internal array in ``ta_var`` storing the result of the evaluation (i.e., the updated values of the state variables). We can see that, in this specific case, the evaluation of the Taylor map reproduces accurately the state vector computed by the non-variational integrator with displaced initial conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "d561e6f0-2500-4faf-9caa-f8bf29ff99e4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Taylor map error: [6.58681443e-13 6.68798350e-13]\n" + ] + } + ], + "source": [ + "print(f'Taylor map error: {ta_var.eval_taylor_map([dx, dv]) - ta.state}')" + ] + }, + { + "cell_type": "markdown", + "id": "fb423c90-3458-4939-82ed-97326baf8be8", + "metadata": {}, + "source": [ + "Note that the Taylor map state vector can also be fetched via the ``tstate`` property of the integrator:" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "d1dc660b-aa56-47d0-bd83-d2f323b12dfd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.11744564, -1.24932198])" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.tstate" + ] + }, + { + "cell_type": "markdown", + "id": "0152b819-815d-4082-9e26-4da55c731d57", + "metadata": {}, + "source": [ + "```{warning}\n", + "Accessing the Taylor map state vector via ``tstate`` will **NOT** trigger any Taylor map evaluation, it will merely return a reference to the internal array storing the result of the evaluation. It is your responsibility to ensure that you called ``eval_taylor_map()`` before accessing this array via ``tstate``.\n", + "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "fc4df22b-1edf-4700-987b-0b663999fb59", + "metadata": {}, + "source": [ + "## A note on computational efficiency\n", + "\n", + "{class}`~heyoka.var_ode_sys` uses internally the {func}`~heyoka.diff_tensors()` and {class}`~heyoka.dtens` API to formulate the variational equations. This means that the computation of the symbolic derivatives is performed in an efficient manner. For instance, reverse-mode symbolic automatic differentiation will be employed when computing the first-order variationals of ODE systems containing a large number of parameters (e.g., in [neural ODEs](./NeuralODEs.ipynb)).\n", + "\n", + "See the [computing derivatives](<./computing_derivatives.ipynb>) tutorial for a more in-depth discussion of how heyoka.py computes symbolic derivatives." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/heyoka/CMakeLists.txt b/heyoka/CMakeLists.txt index 24991218..699e7e49 100644 --- a/heyoka/CMakeLists.txt +++ b/heyoka/CMakeLists.txt @@ -30,6 +30,8 @@ set(HEYOKA_PY_PYTHON_FILES _test_vsop2013.py _test_elp2000.py _test_lagham.py + _test_var_ode_sys.py + _test_var_integrator.py model/__init__.py callback/__init__.py ) @@ -58,6 +60,7 @@ set(_HEYOKA_PY_CORE_SOURCES dtypes.cpp custom_casters.cpp expose_expression.cpp + expose_var_ode_sys.cpp expose_batch_integrators.cpp numpy_memory.cpp expose_models.cpp diff --git a/heyoka/__init__.py b/heyoka/__init__.py index 4c724a15..a358d79e 100644 --- a/heyoka/__init__.py +++ b/heyoka/__init__.py @@ -382,3 +382,20 @@ def _create_par(): >>> p0 = par[0] # p0 will represent the parameter value at index 0 """ + + +# Machinery for the time attribute. +def _create_time(): + from . import core + + return core._time + + +time = _create_time() +""" +Time placeholder. + +This global object is an :py:class:`~heyoka.expression` which is used to represent +time (i.e., the independent variable) in differential equations. + +""" diff --git a/heyoka/_sympy_utils.py b/heyoka/_sympy_utils.py index baf33f6a..c3a20903 100644 --- a/heyoka/_sympy_utils.py +++ b/heyoka/_sympy_utils.py @@ -142,7 +142,7 @@ def _build_fmap(): if not _with_sympy: return None - from . import core, pi + from . import core, pi, time as htime retval = {} @@ -179,7 +179,7 @@ def mul_wrapper(*args): retval[_spy.Function("heyoka_kepE")] = core.kepE retval[_spy.Function("heyoka_kepF")] = core.kepF retval[_spy.Function("heyoka_kepDE")] = core.kepDE - retval[_spy.Function("heyoka_time")] = lambda: core.time + retval[_spy.Function("heyoka_time")] = lambda: htime return retval diff --git a/heyoka/_test_expression.py b/heyoka/_test_expression.py index 90b99cb4..7d677e1c 100644 --- a/heyoka/_test_expression.py +++ b/heyoka/_test_expression.py @@ -342,3 +342,14 @@ def test_relu_wrappers(self): self.assertEqual(leaky_relup(0.0)(x * y), relup(x * y)) self.assertEqual(leaky_relu(0.1)(x * y + y), relu(x * y + y, 0.1)) self.assertEqual(leaky_relup(0.1)(x * y + y), relup(x * y + y, 0.1)) + + def test_dfun(self): + from . import make_vars, dfun + + x, y = make_vars("x", "y") + + self.assertEqual(str(dfun("f", [x, y])), "(∂^0 f)") + self.assertEqual( + str(dfun(name="f", args=[x, y], didx=[(0, 3), (1, 5)])), + "(∂^8 f)/(∂a0^3 ∂a1^5)", + ) diff --git a/heyoka/_test_model.py b/heyoka/_test_model.py index 8e6383be..cd4cb3bb 100644 --- a/heyoka/_test_model.py +++ b/heyoka/_test_model.py @@ -327,7 +327,11 @@ def test_nrlmsise00(self): "h", "lat", "lon", "f107", "f107a", "ap" ) nrlmsise00 = model.nrlmsise00_tn( - geodetic=[h, lat, lon], f107=f107, f107a=f107a, ap=ap, time_expr=time / 86400 + geodetic=[h, lat, lon], + f107=f107, + f107a=f107a, + ap=ap, + time_expr=time / 86400, ) nrlmsise00_cf = cfunc([nrlmsise00], vars=[h, lat, lon, f107, f107a, ap]) @@ -353,21 +357,32 @@ def test_nrlmsise00(self): def test_jb08(self): from . import model, make_vars, cfunc, time - h, lat, lon, f107a, f107, s107a, s107, m107a, m107, y107a, y107, dDstdT = ( - make_vars( - "h", - "lat", - "lon", - "f107a", - "f107", - "s107a", - "s107", - "m107a", - "m107", - "y107a", - "y107", - "dDstdT", - ) + ( + h, + lat, + lon, + f107a, + f107, + s107a, + s107, + m107a, + m107, + y107a, + y107, + dDstdT, + ) = make_vars( + "h", + "lat", + "lon", + "f107a", + "f107", + "s107a", + "s107", + "m107a", + "m107", + "y107a", + "y107", + "dDstdT", ) jb08 = model.jb08_tn( geodetic=[h, lat, lon], diff --git a/heyoka/_test_sympy.py b/heyoka/_test_sympy.py index 3c3d5863..fc7ae17b 100644 --- a/heyoka/_test_sympy.py +++ b/heyoka/_test_sympy.py @@ -192,7 +192,16 @@ def test_func_conversion(self): import sympy as spy - from . import core, make_vars, from_sympy, to_sympy, pi, sum as hsum, prod + from . import ( + core, + make_vars, + from_sympy, + to_sympy, + pi, + sum as hsum, + prod, + time as htime, + ) from .model import nbody @@ -320,8 +329,8 @@ def test_func_conversion(self): self.assertEqual(to_sympy(core.sigmoid(hx + hy)), 1.0 / (1.0 + spy.exp(-x - y))) - self.assertEqual(core.time, from_sympy(spy.Function("heyoka_time")())) - self.assertEqual(to_sympy(core.time), spy.Function("heyoka_time")()) + self.assertEqual(htime, from_sympy(spy.Function("heyoka_time")())) + self.assertEqual(to_sympy(htime), spy.Function("heyoka_time")()) with self.assertRaises(TypeError) as cm: from_sympy(abs(x)) diff --git a/heyoka/_test_var_integrator.py b/heyoka/_test_var_integrator.py new file mode 100644 index 00000000..d74bbdc0 --- /dev/null +++ b/heyoka/_test_var_integrator.py @@ -0,0 +1,236 @@ +# Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +# +# This file is part of the heyoka.py 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/. + + +import unittest as _ut + + +class var_integrator_test_case(_ut.TestCase): + def test_scalar(self): + from . import ( + make_vars, + var_ode_sys, + var_args, + cos, + sin, + par, + time, + taylor_adaptive, + core, + ) + from sys import getrefcount + import numpy as np + + x, v = make_vars("x", "v") + + orig_sys = [(x, v), (v, cos(time) - par[0] * v - sin(x))] + + vsys = var_ode_sys(orig_sys, var_args.vars, order=2) + + ta = taylor_adaptive(vsys, [0.2, 0.3], pars=[0.4], time=0.5, compact_mode=True) + + self.assertTrue(ta.dim > 2) + self.assertEqual(ta.n_orig_sv, 2) + self.assertTrue(ta.is_variational) + self.assertEqual(ta.vorder, 2) + self.assertEqual(ta.vargs, [x, v]) + + # Check that the refcount increases by when accessing tstate. + rc = getrefcount(ta) + ts = ta.tstate + self.assertEqual(getrefcount(ta), rc + 1) + self.assertEqual(ts.shape, (2,)) + self.assertTrue(np.all(ts == [0.0, 0.0])) + + # Test that tstate cannot be written to. + with self.assertRaises(ValueError) as cm: + ta.tstate[0] = 0.5 + + self.assertEqual(ta.get_vslice(order=0), slice(0, 2, None)) + self.assertEqual(ta.get_vslice(order=0, component=1), slice(1, 2, None)) + + self.assertEqual(ta.get_vslice(order=1), slice(2, 6, None)) + self.assertEqual(ta.get_vslice(order=1, component=1), slice(4, 6, None)) + + self.assertEqual(ta.get_mindex(0), [0, 0, 0]) + self.assertEqual(ta.get_mindex(1), [1, 0, 0]) + self.assertEqual(ta.get_mindex(2), [0, 1, 0]) + self.assertEqual(ta.get_mindex(3), [0, 0, 1]) + self.assertEqual(ta.get_mindex(4), [1, 1, 0]) + self.assertEqual(ta.get_mindex(i=5), [1, 0, 1]) + + ta.propagate_until(3.0) + + ts2 = ta.eval_taylor_map([0.0, 0.0]) + self.assertEqual(getrefcount(ta), rc + 2) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + ts2 = ta.eval_taylor_map(np.array([0.0, 0.0])) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + with self.assertRaises(TypeError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0], dtype=np.int32)) + self.assertTrue( + "Invalid dtype detected for the inputs of a Taylor map evaluation:" + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0, 0.0, 0.0])[::2]) + self.assertTrue( + "Invalid inputs array detected in a Taylor map evaluation: the array is not C-style contiguous, please " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([[0.0, 0.0], [0.0, 0.0]])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 2 dimensions, " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0, 0.0, 0.0])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 4 elements, " + in str(cm.exception) + ) + + if not hasattr(core, "real"): + return + + from . import real + + prec = 14 + + ta = taylor_adaptive( + vsys, + [real(0.2, prec), real(0.3, prec)], + pars=[real(0.4, prec)], + time=real(0.5, prec), + compact_mode=True, + fp_type=real, + ) + + ta.propagate_until(real(3.0, prec)) + + with self.assertRaises(ValueError): + ta.eval_taylor_map(np.empty((2,), dtype=real)) + + with self.assertRaises(ValueError): + ta.eval_taylor_map(np.array([real(0.0, prec), real(0.0, prec - 1)])) + + ts2 = ta.eval_taylor_map(np.array([real(0.0, prec), real(0.0, prec)])) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + def test_batch(self): + from . import ( + make_vars, + var_ode_sys, + var_args, + cos, + sin, + par, + time, + taylor_adaptive_batch, + core, + ) + from sys import getrefcount + import numpy as np + + x, v = make_vars("x", "v") + + orig_sys = [(x, v), (v, cos(time) - par[0] * v - sin(x))] + + vsys = var_ode_sys(orig_sys, var_args.vars, order=2) + + ta = taylor_adaptive_batch( + vsys, + [[0.2, 0.21], [0.3, 0.31]], + pars=[[0.4, 0.41]], + time=[0.5, 0.51], + compact_mode=True, + ) + + self.assertTrue(ta.dim > 2) + self.assertEqual(ta.n_orig_sv, 2) + self.assertTrue(ta.is_variational) + self.assertEqual(ta.vorder, 2) + self.assertEqual(ta.vargs, [x, v]) + + # Check that the refcount increases by when accessing tstate. + rc = getrefcount(ta) + ts = ta.tstate + self.assertEqual(getrefcount(ta), rc + 1) + self.assertEqual(ts.shape, (2, 2)) + self.assertTrue(np.all(ts == [[0.0, 0.0], [0.0, 0.0]])) + + # Test that tstate cannot be written to. + with self.assertRaises(ValueError) as cm: + ta.tstate[0] = 0.5 + + self.assertEqual(ta.get_vslice(order=0), slice(0, 2, None)) + self.assertEqual(ta.get_vslice(order=0, component=1), slice(1, 2, None)) + + self.assertEqual(ta.get_vslice(order=1), slice(2, 6, None)) + self.assertEqual(ta.get_vslice(order=1, component=1), slice(4, 6, None)) + + self.assertEqual(ta.get_mindex(0), [0, 0, 0]) + self.assertEqual(ta.get_mindex(1), [1, 0, 0]) + self.assertEqual(ta.get_mindex(2), [0, 1, 0]) + self.assertEqual(ta.get_mindex(3), [0, 0, 1]) + self.assertEqual(ta.get_mindex(4), [1, 1, 0]) + self.assertEqual(ta.get_mindex(i=5), [1, 0, 1]) + + ta.propagate_until(3.0) + + ts2 = ta.eval_taylor_map([[0.0, 0.0], [0.0, 0.0]]) + self.assertEqual(getrefcount(ta), rc + 2) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + ts2 = ta.eval_taylor_map(np.array([[0.0, 0.0], [0.0, 0.0]])) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + with self.assertRaises(TypeError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0], dtype=np.int32)) + self.assertTrue( + "Invalid dtype detected for the inputs of a Taylor map evaluation:" + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0, 0.0, 0.0])[::2]) + self.assertTrue( + "Invalid inputs array detected in a Taylor map evaluation: the array is not C-style contiguous, please " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 1 dimension(s), " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([[0.0, 0.0]])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 1 row(s), but it must have 2 row(s) instead" + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([[0.0], [0.0]])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 1 column(s), but it must have 2 column(s) instead" + in str(cm.exception) + ) diff --git a/heyoka/_test_var_ode_sys.py b/heyoka/_test_var_ode_sys.py new file mode 100644 index 00000000..c44cee2e --- /dev/null +++ b/heyoka/_test_var_ode_sys.py @@ -0,0 +1,69 @@ +# Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +# +# This file is part of the heyoka.py 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/. + + +import unittest as _ut + + +class var_ode_sys_test_case(_ut.TestCase): + def test_enum(self): + from . import var_args + + self.assertEqual(var_args.vars | var_args.time | var_args.params, var_args.all) + self.assertTrue((var_args.vars | var_args.time) & var_args.time) + + def test_basic(self): + from . import make_vars, var_ode_sys, var_args, sin, par, time + from copy import copy, deepcopy + from pickle import dumps, loads + + x, v = make_vars("x", "v") + + orig_sys = [(x, v), (v, -par[0] * sin(x) + time)] + + vsys = var_ode_sys(orig_sys, var_args.vars) + + self.assertEqual(orig_sys, vsys.sys[:2]) + self.assertEqual(vsys.vargs, [x, v]) + self.assertEqual(vsys.n_orig_sv, 2) + self.assertEqual(vsys.order, 1) + + vsys = var_ode_sys(sys=orig_sys, args=[v, time, x], order=2) + + self.assertEqual(orig_sys, vsys.sys[:2]) + self.assertEqual(vsys.vargs, [v, time, x]) + self.assertEqual(vsys.n_orig_sv, 2) + self.assertEqual(vsys.order, 2) + + vsys = var_ode_sys(orig_sys, var_args.vars | var_args.params) + + self.assertEqual(orig_sys, vsys.sys[:2]) + self.assertEqual(vsys.vargs, [x, v, par[0]]) + self.assertEqual(vsys.n_orig_sv, 2) + self.assertEqual(vsys.order, 1) + + vsys2 = copy(vsys) + + self.assertEqual(orig_sys, vsys2.sys[:2]) + self.assertEqual(vsys2.vargs, [x, v, par[0]]) + self.assertEqual(vsys2.n_orig_sv, 2) + self.assertEqual(vsys2.order, 1) + + vsys2 = deepcopy(vsys) + + self.assertEqual(orig_sys, vsys2.sys[:2]) + self.assertEqual(vsys2.vargs, [x, v, par[0]]) + self.assertEqual(vsys2.n_orig_sv, 2) + self.assertEqual(vsys2.order, 1) + + vsys2 = loads(dumps(vsys)) + + self.assertEqual(orig_sys, vsys2.sys[:2]) + self.assertEqual(vsys2.vargs, [x, v, par[0]]) + self.assertEqual(vsys2.n_orig_sv, 2) + self.assertEqual(vsys2.order, 1) diff --git a/heyoka/cfunc.cpp b/heyoka/cfunc.cpp index 5b55ea67..8dd564bc 100644 --- a/heyoka/cfunc.cpp +++ b/heyoka/cfunc.cpp @@ -123,19 +123,6 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) // Fetch the dtype corresponding to T. const auto dt = get_dtype(); - // Small helper to facilitate the conversion of an iterable into - // a contiguous NumPy array of type dt. - const auto as_carray = [dt](const py::iterable &v) { - using namespace pybind11::literals; - - py::array ret = py::module_::import("numpy").attr("ascontiguousarray")(v, "dtype"_a = py::dtype(dt)); - - assert(ret.dtype().num() == dt); - assert(is_npy_array_carray(ret)); - - return ret; - }; - // Fetch the inputs array. // NOTE: this will either fetch the existing array, or convert // the input iterable to a py::array on the fly. @@ -162,7 +149,7 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) return std::move(v); } else { - return as_carray(v); + return as_carray(v, dt); } }, inputs_ob); @@ -261,7 +248,7 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) return std::move(v); } else { - return as_carray(v); + return as_carray(v, dt); } }, *pars_ob); @@ -328,7 +315,7 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) return std::move(v); } else { - return as_carray(v); + return as_carray(v, dt); } }, std::get(tm)); diff --git a/heyoka/common_utils.cpp b/heyoka/common_utils.cpp index 4524e98e..f89e9493 100644 --- a/heyoka/common_utils.cpp +++ b/heyoka/common_utils.cpp @@ -9,8 +9,13 @@ #include #include +#include #include #include +#include + +#include +#include #include #include @@ -33,6 +38,7 @@ #endif +#include #include #include "common_utils.hpp" @@ -108,4 +114,62 @@ bool with_pybind11_eh_impl() } // namespace detail +std::pair +dtens_t_it::operator()(const std::pair &p) const +{ + const auto &[sv_idx, ex] = p; + + return std::make_pair( + sparse_to_dense(sv_idx, boost::numeric_cast(dt->get_nargs())), ex); +} + +heyoka::dtens::v_idx_t dtens_t_it::sparse_to_dense(const heyoka::dtens::sv_idx_t &sv_idx, + heyoka::dtens::v_idx_t::size_type nargs) +{ + // Init the dense vector from the component index. + heyoka::dtens::v_idx_t ret{sv_idx.first}; + + // Transform the sparse index/order pairs into dense format. + // NOTE: no overflow check needed on ++idx because dtens ensures that + // the number of variables can be represented by std::uint32_t. + std::uint32_t idx = 0; + for (auto it = sv_idx.second.begin(); it != sv_idx.second.end(); ++idx) { + if (it->first == idx) { + // The current index shows up in the sparse vector, + // fetch the corresponding order and move to the next + // element of the sparse vector. + ret.push_back(it->second); + assert(it->second != 0u); + ++it; + } else { + // The current index does not show up in the sparse + // vector, set the order to zero. + ret.push_back(0); + } + } + + // Sanity check on the number of diff variables + // inferred from the sparse vector. + assert(ret.size() - 1u <= nargs); + + // Pad missing values at the end of ret. + ret.resize(boost::safe_numerics::safe(nargs) + 1); + + return ret; +} + +// Small helper to facilitate the conversion of an iterable into +// a contiguous NumPy array of type dt. +py::array as_carray(const py::iterable &v, int dt) +{ + using namespace pybind11::literals; + + py::array ret = py::module_::import("numpy").attr("ascontiguousarray")(v, "dtype"_a = py::dtype(dt)); + + assert(ret.dtype().num() == dt); + assert(is_npy_array_carray(ret)); + + return ret; +} + } // namespace heyoka_py diff --git a/heyoka/common_utils.hpp b/heyoka/common_utils.hpp index c5b4afa8..83f5e0bb 100644 --- a/heyoka/common_utils.hpp +++ b/heyoka/common_utils.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #if defined(__GLIBCXX__) @@ -24,6 +25,7 @@ #include +#include #include // NOTE: implementation of Py_SET_TYPE() for Python < 3.9. See: @@ -198,6 +200,20 @@ bool with_pybind11_eh(const F &f) } } +// Functor to transform on-the-fly the content of a dtens +// from sparse format into dense format. +struct dtens_t_it { + const heyoka::dtens *dt = nullptr; + + std::pair + operator()(const std::pair &) const; + + // Helper to implement the conversion from sparse to dense format. + static heyoka::dtens::v_idx_t sparse_to_dense(const heyoka::dtens::sv_idx_t &, heyoka::dtens::v_idx_t::size_type); +}; + +py::array as_carray(const py::iterable &, int); + } // namespace heyoka_py #endif diff --git a/heyoka/core.cpp b/heyoka/core.cpp index 2145e5e9..93858165 100644 --- a/heyoka/core.cpp +++ b/heyoka/core.cpp @@ -61,6 +61,7 @@ #include "expose_models.hpp" #include "expose_real.hpp" #include "expose_real128.hpp" +#include "expose_var_ode_sys.hpp" #include "logging.hpp" #include "pickle_wrappers.hpp" #include "setup_sympy.hpp" @@ -197,6 +198,9 @@ PYBIND11_MODULE(core, m) // Expression. heypy::expose_expression(m); + // var_ode_sys. + heypy::expose_var_ode_sys(m); + // Models. heypy::expose_models(m); diff --git a/heyoka/docstrings.cpp b/heyoka/docstrings.cpp index 43bf1a8c..128de087 100644 --- a/heyoka/docstrings.cpp +++ b/heyoka/docstrings.cpp @@ -535,4 +535,147 @@ A few checks are run on the input arguments. Specifically, the number of geodesi )"; } +std::string var_ode_sys() +{ + return R"(Class to represent variational ODE systems. + +.. note:: + + A :ref:`tutorial ` explaining the use of this + class is available. + +)"; +} + +std::string var_ode_sys_init() +{ + return R"(__init__(self, sys: list[tuple[expression, expression]], args: var_args | list[expression], order: int = 1) + +Constructor. + +A variational ODE system is constructed from two mandatory arguments: the original ODE system +*sys* and an *args* parameter representing the quantities with respect to which the variational +equations will be formulated. + +If *args* is of type :class:`~heyoka.var_args`, then the variational equations will be constructed +with respect to arguments deduced from *sys*. E.g., if *sys* contains the two state variables +:math:`x` and :math:`y` and *args* is the *vars* enumerator of :class:`~heyoka.var_args`, then +the variational equations will be formulated with respect to the initial conditions of +:math:`x` and :math:`y`. Similarly, if *sys* contains two parameters ``par[0]`` and ``par[1]`` +and *args* is the *params* enumerator of :class:`~heyoka.var_args`, then +the variational equations will be formulated with respect to the two parameters. + +If *args* is a list of :class:`~heyoka.expression`, then the variational equations will be formulated +with respect to the quantities contained in the list. Specifically: + +- variable expression are used to request derivatives with respect to the intial conditions + for that state variable; +- parameter expressions are used to request derivatives with respect to those parameters; +- :attr:`heyoka.time` is used to request the derivative with respect to the initial integration time. + +Several checks are run on the input arguments. Specifically: + +- *sys* must be a valid ODE system; +- if *args* is of type :class:`~heyoka.var_args`, it must be either one of the valid enumerators + or a combination of the valid enumerators; +- if *args* is a list of expressions: + + - it must not be empty, + - it must consists only of variables, parameters or the :attr:`heyoka.time` placeholder, + - it must not contain duplicates, + - any variable expression must refer to a state variable in *sys*. + +Additionally, the differentiation order *order* must be at least 1. + +:param sys: the input ODE system. +:param args: the variational arguments. +:param order: the differentiation order. + +:raises ValueError: if one or more input arguments are malformed, as explained above. + +)"; +} + +std::string var_ode_sys_sys() +{ + return R"(The full system of equations (including partials). + +:rtype: list[tuple[expression, expression]] + +)"; +} + +std::string var_ode_sys_vargs() +{ + return R"(The list of variational arguments. + +:rtype: list[expression] + +)"; +} + +std::string var_ode_sys_n_orig_sv() +{ + return R"(The number of original state variables. + +:rtype: int + +)"; +} + +std::string var_ode_sys_order() +{ + return R"(The differentitation order. + +:rtype: int + +)"; +} + +std::string var_args() +{ + return R"(Enum for selecting variational arguments. + +Values of this enum can be used in the constructor of :class:`~heyoka.var_ode_sys()` to select +the arguments with respect to which the variational equations will be formulated. + +The enumerators can be combined with the logical OR ``|`` operator. + +Examples: + >>> from heyoka import var_args + >>> va = var_args.vars | var_args.params # Select differentiation wrt all initial conditions and parameters + >>> var_args.vars | var_args.params | var_args.time == var_args.all + True + +)"; +} + +std::string var_args_vars() +{ + return R"(Differentiate with respect to the initial conditions of all state variables. + +)"; +} + +std::string var_args_params() +{ + return R"(Differentiate with respect to all runtime parameters. + +)"; +} + +std::string var_args_time() +{ + return R"(Differentiate with respect to the initial integration time. + +)"; +} + +std::string var_args_all() +{ + return R"(Differentiate with respect to the initial conditions of all state variables, all runtime parameters and the initial integration time. + +)"; +} + } // namespace heyoka_py::docstrings diff --git a/heyoka/docstrings.hpp b/heyoka/docstrings.hpp index c455226b..67cb6c17 100644 --- a/heyoka/docstrings.hpp +++ b/heyoka/docstrings.hpp @@ -44,11 +44,24 @@ std::string diff_args_all(); std::string lagrangian(); std::string hamiltonian(); -// Models +// Models. std::string cart2geo(); std::string nrlmsise00_tn(); std::string jb08_tn(); +// var_ode_sys() and related. +std::string var_args(); +std::string var_args_vars(); +std::string var_args_params(); +std::string var_args_time(); +std::string var_args_all(); +std::string var_ode_sys(); +std::string var_ode_sys_init(); +std::string var_ode_sys_sys(); +std::string var_ode_sys_vargs(); +std::string var_ode_sys_n_orig_sv(); +std::string var_ode_sys_order(); + } // namespace heyoka_py::docstrings #endif diff --git a/heyoka/expose_batch_integrators.cpp b/heyoka/expose_batch_integrators.cpp index 42e87c15..b4c28770 100644 --- a/heyoka/expose_batch_integrators.cpp +++ b/heyoka/expose_batch_integrators.cpp @@ -6,9 +6,12 @@ // 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 +#include #include #include #include +#include #include #include #include @@ -31,6 +34,7 @@ #include #include #include +#include #include "common_utils.hpp" #include "custom_casters.hpp" @@ -43,6 +47,7 @@ namespace heyoka_py { namespace py = pybind11; +namespace hey = heyoka; namespace detail { @@ -50,10 +55,29 @@ namespace detail namespace { +// Helper to fetch the tstate from a variational integrator. +// Extracted here for re-use. +template +auto fetch_tstate(const py::object &o) +{ + const auto *ta = py::cast *>(o); + + assert(ta->get_tstate().size() % ta->get_batch_size() == 0u); + const auto n_orig_sv = boost::numeric_cast(ta->get_n_orig_sv()); + const auto bs = boost::numeric_cast(ta->get_batch_size()); + + auto ret + = py::array(py::dtype(get_dtype()), py::array::ShapeContainer{n_orig_sv, bs}, ta->get_tstate().data(), o); + + // Ensure the returned array is read-only. + ret.attr("flags").attr("writeable") = false; + + return ret; +} + template void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) { - namespace hey = heyoka; namespace kw = hey::kw; using namespace pybind11::literals; @@ -64,10 +88,11 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) using sys_t = std::vector>; // Implementation of the ctor. - auto tab_ctor_impl = [](const sys_t &sys, const py::iterable &state_ob, std::optional time_ob, - std::optional pars_ob, T tol, bool high_accuracy, bool compact_mode, - std::vector tes, std::vector ntes, bool parallel_mode, unsigned opt_level, - bool force_avx512, bool slp_vectorize, bool fast_math) { + auto tab_ctor_impl = [](std::variant vsys, const py::iterable &state_ob, + std::optional time_ob, std::optional pars_ob, T tol, + bool high_accuracy, bool compact_mode, std::vector tes, std::vector ntes, + bool parallel_mode, unsigned opt_level, bool force_avx512, bool slp_vectorize, + bool fast_math) { // Fetch the dtype corresponding to T. const auto dt = get_dtype(); @@ -115,77 +140,83 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) pars = py::cast>(pars_arr.attr("flatten")()); } - if (time_ob) { - // Times provided. - py::array time_arr = *time_ob; - if (time_arr.ndim() != 1 || boost::numeric_cast(time_arr.shape(0)) != batch_size) { - py_throw(PyExc_ValueError, - fmt::format("Invalid time vector passed to the constructor of a batch integrator: " - "the expected array shape is ({}), but the input array has either the wrong " - "number of dimensions or the wrong shape", - batch_size) - .c_str()); - } - auto time = py::cast>(time_arr); - - // NOTE: GIL release is fine here even if the events contain - // Python objects, as the event vectors are moved in - // upon construction and thus we should never end up calling - // into the interpreter. - py::gil_scoped_release release; - - return hey::taylor_adaptive_batch{sys, - std::move(state), - batch_size, - kw::time = std::move(time), - kw::tol = tol, - kw::high_accuracy = high_accuracy, - kw::compact_mode = compact_mode, - kw::pars = std::move(pars), - kw::t_events = std::move(tes), - kw::nt_events = std::move(ntes), - kw::parallel_mode = parallel_mode, - kw::opt_level = opt_level, - kw::force_avx512 = force_avx512, - kw::slp_vectorize = slp_vectorize, - kw::fast_math = fast_math}; - } else { - // Times not provided. - - // NOTE: GIL release is fine here even if the events contain - // Python objects, as the event vectors are moved in - // upon construction and thus we should never end up calling - // into the interpreter. - py::gil_scoped_release release; - - return hey::taylor_adaptive_batch{sys, - std::move(state), - batch_size, - kw::tol = tol, - kw::high_accuracy = high_accuracy, - kw::compact_mode = compact_mode, - kw::pars = std::move(pars), - kw::t_events = std::move(tes), - kw::nt_events = std::move(ntes), - kw::parallel_mode = parallel_mode, - kw::opt_level = opt_level, - kw::force_avx512 = force_avx512, - kw::slp_vectorize = slp_vectorize, - kw::fast_math = fast_math}; - } + return std::visit( + [&](auto &sys) { + if (time_ob) { + // Times provided. + py::array time_arr = *time_ob; + if (time_arr.ndim() != 1 || boost::numeric_cast(time_arr.shape(0)) != batch_size) { + py_throw( + PyExc_ValueError, + fmt::format("Invalid time vector passed to the constructor of a batch integrator: " + "the expected array shape is ({}), but the input array has either the wrong " + "number of dimensions or the wrong shape", + batch_size) + .c_str()); + } + auto time = py::cast>(time_arr); + + // NOTE: GIL release is fine here even if the events contain + // Python objects, as the event vectors are moved in + // upon construction and thus we should never end up calling + // into the interpreter. + py::gil_scoped_release release; + + return hey::taylor_adaptive_batch{std::move(sys), + std::move(state), + batch_size, + kw::time = std::move(time), + kw::tol = tol, + kw::high_accuracy = high_accuracy, + kw::compact_mode = compact_mode, + kw::pars = std::move(pars), + kw::t_events = std::move(tes), + kw::nt_events = std::move(ntes), + kw::parallel_mode = parallel_mode, + kw::opt_level = opt_level, + kw::force_avx512 = force_avx512, + kw::slp_vectorize = slp_vectorize, + kw::fast_math = fast_math}; + } else { + // Times not provided. + + // NOTE: GIL release is fine here even if the events contain + // Python objects, as the event vectors are moved in + // upon construction and thus we should never end up calling + // into the interpreter. + py::gil_scoped_release release; + + return hey::taylor_adaptive_batch{std::move(sys), + std::move(state), + batch_size, + kw::tol = tol, + kw::high_accuracy = high_accuracy, + kw::compact_mode = compact_mode, + kw::pars = std::move(pars), + kw::t_events = std::move(tes), + kw::nt_events = std::move(ntes), + kw::parallel_mode = parallel_mode, + kw::opt_level = opt_level, + kw::force_avx512 = force_avx512, + kw::slp_vectorize = slp_vectorize, + kw::fast_math = fast_math}; + } + }, + vsys); }; py::class_> tab_c(m, fmt::format("taylor_adaptive_batch_{}", suffix).c_str(), py::dynamic_attr{}); tab_c - .def(py::init([tab_ctor_impl](const sys_t &sys, const py::iterable &state, std::optional time, - std::optional pars, T tol, bool high_accuracy, bool compact_mode, - std::vector tes, std::vector ntes, bool parallel_mode, - unsigned opt_level, bool force_avx512, bool slp_vectorize, bool fast_math) { - return tab_ctor_impl(sys, state, std::move(time), std::move(pars), tol, high_accuracy, compact_mode, - std::move(tes), std::move(ntes), parallel_mode, opt_level, force_avx512, - slp_vectorize, fast_math); + .def(py::init([tab_ctor_impl](std::variant vsys, const py::iterable &state, + std::optional time, std::optional pars, T tol, + bool high_accuracy, bool compact_mode, std::vector tes, + std::vector ntes, bool parallel_mode, unsigned opt_level, + bool force_avx512, bool slp_vectorize, bool fast_math) { + return tab_ctor_impl(std::move(vsys), state, std::move(time), std::move(pars), tol, high_accuracy, + compact_mode, std::move(tes), std::move(ntes), parallel_mode, opt_level, + force_avx512, slp_vectorize, fast_math); }), "sys"_a, "state"_a, "time"_a = py::none{}, "pars"_a = py::none{}, "tol"_a.noconvert() = static_cast(0), "high_accuracy"_a = false, "compact_mode"_a = false, "t_events"_a = py::list{}, "nt_events"_a = py::list{}, @@ -509,6 +540,103 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) .def_property_readonly("compact_mode", &hey::taylor_adaptive_batch::get_compact_mode) .def_property_readonly("high_accuracy", &hey::taylor_adaptive_batch::get_high_accuracy) .def_property_readonly("with_events", &hey::taylor_adaptive_batch::with_events) + // Variational-specific bits. + .def_property_readonly("n_orig_sv", &hey::taylor_adaptive_batch::get_n_orig_sv) + .def_property_readonly("is_variational", &hey::taylor_adaptive_batch::is_variational) + .def_property_readonly("vargs", &hey::taylor_adaptive_batch::get_vargs) + .def_property_readonly("vorder", &hey::taylor_adaptive_batch::get_vorder) + .def_property_readonly("tstate", &fetch_tstate) + .def( + "get_vslice", + [](const hey::taylor_adaptive_batch &ta, std::uint32_t order, std::optional component) { + const auto ret = component ? ta.get_vslice(*component, order) : ta.get_vslice(order); + + return py::slice(boost::numeric_cast(ret.first), + boost::numeric_cast(ret.second), {}); + }, + "order"_a, "component"_a = py::none{}) + .def( + "get_mindex", + [](const hey::taylor_adaptive_batch &ta, std::uint32_t i) { + const auto &ret = ta.get_mindex(i); + + return dtens_t_it::sparse_to_dense( + ret, boost::numeric_cast(ta.get_vargs().size())); + }, + "i"_a) + .def("eval_taylor_map", + [](py::object &o, std::variant in) { + auto *ta = py::cast *>(o); + + // Fetch the dtype corresponding to T. + const auto dt = get_dtype(); + + // Fetch the inputs array. + // NOTE: this will either fetch the existing array, or convert + // the input iterable to a py::array on the fly. + const auto inputs = std::visit( + [&](U &v) { + if constexpr (std::same_as) { + // Check the dtype. + if (v.dtype().num() != dt) [[unlikely]] { + py_throw( + PyExc_TypeError, + fmt::format( + "Invalid dtype detected for the inputs of a Taylor map evaluation: the " + "expected dtype is '{}', but the dtype of the inputs array is '{}' instead", + str(py::dtype(dt)), str(v.dtype())) + .c_str()); + } + + // Check that the inputs array is a C-style contiguous array. + if (!is_npy_array_carray(v)) [[unlikely]] { + py_throw(PyExc_ValueError, + "Invalid inputs array detected in a Taylor map evaluation: the array is not " + "C-style contiguous, please " + "consider using numpy.ascontiguousarray() to turn it into one"); + } + + return std::move(v); + } else { + return as_carray(v, dt); + } + }, + in); + + // Validate the number of dimensions for the inputs. + // NOTE: this needs to be done regardless of the original type of in. + if (inputs.ndim() != 2) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} dimension(s), " + "but it must have 2 dimensions instead", + inputs.ndim()) + .c_str()); + } + + // Validate the shape for the inputs. + if (boost::numeric_cast(inputs.shape(0)) != ta->get_n_orig_sv()) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} row(s), " + "but it must have {} row(s) instead", + inputs.shape(0), ta->get_n_orig_sv()) + .c_str()); + } + + if (boost::numeric_cast(inputs.shape(1)) != ta->get_batch_size()) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} column(s), " + "but it must have {} column(s) instead", + inputs.shape(1), ta->get_batch_size()) + .c_str()); + } + + // Run the evaluation. + const auto *data_ptr = static_cast(inputs.data()); + ta->eval_taylor_map(std::ranges::subrange(data_ptr, data_ptr + inputs.size())); + + // Return the tstate. + return fetch_tstate(o); + }) // Event detection. .def_property_readonly("with_events", &hey::taylor_adaptive_batch::with_events) .def_property_readonly("te_cooldowns", &hey::taylor_adaptive_batch::get_te_cooldowns) diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index 8d5fb8ab..049c8649 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -24,7 +24,6 @@ #include #include -#include #include #include @@ -52,67 +51,12 @@ #include "common_utils.hpp" #include "custom_casters.hpp" #include "docstrings.hpp" +#include "expose_expression.hpp" #include "pickle_wrappers.hpp" namespace heyoka_py { -namespace detail -{ - -namespace -{ - -template -using uncvref_t = std::remove_cv_t>; - -// Functor to transform on-the-fly the content of a dtens -// from sparse format into dense format. -struct dtens_t_it { - const heyoka::dtens *dt = nullptr; - - std::pair - operator()(const std::pair &p) const - { - const auto &[sv_idx, ex] = p; - - // Init the dense vector from the component index. - heyoka::dtens::v_idx_t ret{sv_idx.first}; - - // Transform the sparse index/order pairs into dense format. - // NOTE: no overflow check needed on ++idx because dtens ensures that - // the number of variables can be represented by std::uint32_t. - std::uint32_t idx = 0; - for (auto it = sv_idx.second.begin(); it != sv_idx.second.end(); ++idx) { - if (it->first == idx) { - // The current index shows up in the sparse vector, - // fetch the corresponding order and move to the next - // element of the sparse vector. - ret.push_back(it->second); - assert(it->second != 0u); - ++it; - } else { - // The current index does not show up in the sparse - // vector, set the order to zero. - ret.push_back(0); - } - } - - // Sanity check on the number of diff variables - // inferred from the sparse vector. - assert(ret.size() - 1u <= dt->get_nargs()); - - // Pad missing values at the end of ret. - ret.resize(boost::safe_numerics::safe(dt->get_nargs()) + 1); - - return std::make_pair(std::move(ret), ex); - } -}; - -} // namespace - -} // namespace detail - namespace py = pybind11; // NOTE: regarding single-precision support: we only expose the expression ctor, @@ -381,8 +325,8 @@ void expose_expression(py::module_ &m) [](const mvf_arg &e, const mvf_arg &M) { return std::visit( [](const auto &a, const auto &b) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; if constexpr (!std::is_same_v && !std::is_same_v) { py_throw(PyExc_TypeError, "At least one of the arguments of kepE() must be an expression"); @@ -400,9 +344,9 @@ void expose_expression(py::module_ &m) [](const mvf_arg &h, const mvf_arg &k, const mvf_arg &lam) { return std::visit( [](const auto &a, const auto &b, const auto &c) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; - using tp3 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; + using tp3 = std::remove_cvref_t; constexpr auto tp1_num = static_cast(!std::is_same_v); constexpr auto tp2_num = static_cast(!std::is_same_v); @@ -438,9 +382,9 @@ void expose_expression(py::module_ &m) [](const mvf_arg &s0, const mvf_arg &c0, const mvf_arg &DM) { return std::visit( [](const auto &a, const auto &b, const auto &c) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; - using tp3 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; + using tp3 = std::remove_cvref_t; constexpr auto tp1_num = static_cast(!std::is_same_v); constexpr auto tp2_num = static_cast(!std::is_same_v); @@ -477,8 +421,8 @@ void expose_expression(py::module_ &m) [](const mvf_arg &y, const mvf_arg &x) { return std::visit( [](const auto &a, const auto &b) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; if constexpr (!std::is_same_v && !std::is_same_v) { py_throw(PyExc_TypeError, "At least one of the arguments of atan2() must be an expression"); @@ -490,8 +434,21 @@ void expose_expression(py::module_ &m) }, "y"_a.noconvert(), "x"_a.noconvert()); + // dfun(). + m.def( + "dfun", + [](std::string name, std::vector args, + std::optional>> didx) { + if (didx) { + return hey::dfun(std::move(name), std::move(args), std::move(*didx)); + } else { + return hey::dfun(std::move(name), std::move(args)); + } + }, + "name"_a, "args"_a, "didx"_a = py::none{}); + // Time. - m.attr("time") = hey::time; + m.attr("_time") = hey::time; // pi. m.attr("pi") = hey::pi; @@ -553,7 +510,7 @@ void expose_expression(py::module_ &m) const auto s_idx = boost::numeric_cast::difference_type>(idx); - return detail::dtens_t_it{&dt}(dt.begin()[s_idx]); + return dtens_t_it{&dt}(dt.begin()[s_idx]); }); dtens_cl.def("__contains__", [](const hey::dtens &dt, const std::variant &v_idx_) { @@ -563,8 +520,8 @@ void expose_expression(py::module_ &m) dtens_cl.def( "__iter__", [](const hey::dtens &dt) { - auto t_begin = boost::iterators::make_transform_iterator(dt.begin(), detail::dtens_t_it{&dt}); - auto t_end = boost::iterators::make_transform_iterator(dt.end(), detail::dtens_t_it{&dt}); + auto t_begin = boost::iterators::make_transform_iterator(dt.begin(), dtens_t_it{&dt}); + auto t_end = boost::iterators::make_transform_iterator(dt.end(), dtens_t_it{&dt}); return py::make_key_iterator(t_begin, t_end); }, @@ -586,8 +543,8 @@ void expose_expression(py::module_ &m) [](const hey::dtens &dt, std::uint32_t order, std::optional component) { const auto sr = component ? dt.get_derivatives(*component, order) : dt.get_derivatives(order); - auto t_begin = boost::iterators::make_transform_iterator(sr.begin(), detail::dtens_t_it{&dt}); - auto t_end = boost::iterators::make_transform_iterator(sr.end(), detail::dtens_t_it{&dt}); + auto t_begin = boost::iterators::make_transform_iterator(sr.begin(), dtens_t_it{&dt}); + auto t_end = boost::iterators::make_transform_iterator(sr.end(), dtens_t_it{&dt}); return std::vector(t_begin, t_end); }, diff --git a/heyoka/expose_var_ode_sys.cpp b/heyoka/expose_var_ode_sys.cpp new file mode 100644 index 00000000..60f823ba --- /dev/null +++ b/heyoka/expose_var_ode_sys.cpp @@ -0,0 +1,61 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka.py 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 +#include +#include + +#include +#include + +#include + +#include "common_utils.hpp" +#include "docstrings.hpp" +#include "expose_var_ode_sys.hpp" +#include "pickle_wrappers.hpp" + +namespace heyoka_py +{ + +namespace py = pybind11; + +void expose_var_ode_sys(py::module_ &m) +{ + namespace hey = heyoka; + // NOLINTNEXTLINE(google-build-using-namespace) + using namespace pybind11::literals; + + // var_args enum. + py::enum_(m, "var_args", docstrings::var_args().c_str()) + .value("vars", hey::var_args::vars, docstrings::var_args_vars().c_str()) + .value("params", hey::var_args::params, docstrings::var_args_params().c_str()) + .value("time", hey::var_args::time, docstrings::var_args_time().c_str()) + .value("all", hey::var_args::all, docstrings::var_args_all().c_str()) + .def("__or__", [](hey::var_args a, hey::var_args b) { return a | b; }) + .def("__and__", [](hey::var_args a, hey::var_args b) { return a & b; }); + + // var_ode_sys class. + py::class_ v_cl(m, "var_ode_sys", py::dynamic_attr{}, docstrings::var_ode_sys().c_str()); + v_cl.def(py::init([](const std::vector> &sys, + const std::variant> &args, + std::uint32_t order) { return hey::var_ode_sys(sys, args, order); }), + "sys"_a, "args"_a, "order"_a = static_cast(1), docstrings::var_ode_sys_init().c_str()); + v_cl.def_property_readonly("sys", &hey::var_ode_sys::get_sys, docstrings::var_ode_sys_sys().c_str()); + v_cl.def_property_readonly("vargs", &hey::var_ode_sys::get_vargs, docstrings::var_ode_sys_vargs().c_str()); + v_cl.def_property_readonly("n_orig_sv", &hey::var_ode_sys::get_n_orig_sv, + docstrings::var_ode_sys_n_orig_sv().c_str()); + v_cl.def_property_readonly("order", &hey::var_ode_sys::get_order, docstrings::var_ode_sys_order().c_str()); + // Copy/deepcopy. + v_cl.def("__copy__", copy_wrapper); + v_cl.def("__deepcopy__", deepcopy_wrapper, "memo"_a); + // Pickle support. + v_cl.def(py::pickle(&pickle_getstate_wrapper, &pickle_setstate_wrapper)); +} + +} // namespace heyoka_py diff --git a/heyoka/expose_var_ode_sys.hpp b/heyoka/expose_var_ode_sys.hpp new file mode 100644 index 00000000..1b452567 --- /dev/null +++ b/heyoka/expose_var_ode_sys.hpp @@ -0,0 +1,23 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka.py 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_PY_EXPOSE_VAR_ODE_SYS_HPP +#define HEYOKA_PY_EXPOSE_VAR_ODE_SYS_HPP + +#include + +namespace heyoka_py +{ + +namespace py = pybind11; + +void expose_var_ode_sys(py::module_ &); + +} // namespace heyoka_py + +#endif diff --git a/heyoka/taylor_expose_integrator.cpp b/heyoka/taylor_expose_integrator.cpp index d90cf4c7..a92af47d 100644 --- a/heyoka/taylor_expose_integrator.cpp +++ b/heyoka/taylor_expose_integrator.cpp @@ -10,9 +10,11 @@ #include #include +#include #include -#include +#include #include +#include #include #include #include @@ -49,6 +51,7 @@ #include #include #include +#include #include "common_utils.hpp" #include "custom_casters.hpp" @@ -57,6 +60,12 @@ #include "step_cb_utils.hpp" #include "taylor_expose_integrator.hpp" +#if defined(HEYOKA_HAVE_REAL) + +#include "expose_real.hpp" + +#endif + namespace heyoka_py { @@ -78,6 +87,22 @@ constexpr bool default_cm = #endif ; +// Helper to fetch the tstate from a variational integrator. +// Extracted here for re-use. +template +auto fetch_tstate(const py::object &o) +{ + const auto *ta = py::cast *>(o); + auto ret = py::array(py::dtype(get_dtype()), + py::array::ShapeContainer{boost::numeric_cast(ta->get_tstate().size())}, + ta->get_tstate().data(), o); + + // Ensure the returned array is read-only. + ret.attr("flags").attr("writeable") = false; + + return ret; +} + // Implementation of the exposition of the scalar integrators. template void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) @@ -91,30 +116,35 @@ void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) using sys_t = std::vector>; py::class_> cl(m, (fmt::format("taylor_adaptive_{}", suffix)).c_str(), py::dynamic_attr{}); - cl.def(py::init([](const sys_t &sys, std::vector state, T time, std::vector pars, T tol, bool high_accuracy, - bool compact_mode, std::vector tes, std::vector ntes, bool parallel_mode, - unsigned opt_level, bool force_avx512, bool slp_vectorize, bool fast_math, long long prec) { + cl.def(py::init([](std::variant vsys, std::vector state, T time, std::vector pars, + T tol, bool high_accuracy, bool compact_mode, std::vector tes, std::vector ntes, + bool parallel_mode, unsigned opt_level, bool force_avx512, bool slp_vectorize, bool fast_math, + long long prec) { // NOTE: GIL release is fine here even if the events contain // Python objects, as the event vectors are moved in // upon construction and thus we should never end up calling // into the interpreter. py::gil_scoped_release release; - return hey::taylor_adaptive{sys, - std::move(state), - kw::time = time, - kw::tol = tol, - kw::high_accuracy = high_accuracy, - kw::compact_mode = compact_mode, - kw::pars = std::move(pars), - kw::t_events = std::move(tes), - kw::nt_events = std::move(ntes), - kw::parallel_mode = parallel_mode, - kw::opt_level = opt_level, - kw::force_avx512 = force_avx512, - kw::slp_vectorize = slp_vectorize, - kw::fast_math = fast_math, - kw::prec = prec}; + return std::visit( + [&](auto &sys) { + return hey::taylor_adaptive{std::move(sys), + std::move(state), + kw::time = time, + kw::tol = tol, + kw::high_accuracy = high_accuracy, + kw::compact_mode = compact_mode, + kw::pars = std::move(pars), + kw::t_events = std::move(tes), + kw::nt_events = std::move(ntes), + kw::parallel_mode = parallel_mode, + kw::opt_level = opt_level, + kw::force_avx512 = force_avx512, + kw::slp_vectorize = slp_vectorize, + kw::fast_math = fast_math, + kw::prec = prec}; + }, + vsys); }), "sys"_a, "state"_a.noconvert(), "time"_a.noconvert() = static_cast(0), "pars"_a.noconvert() = py::list{}, "tol"_a.noconvert() = static_cast(0), "high_accuracy"_a = false, "compact_mode"_a = default_cm, @@ -344,6 +374,105 @@ void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) .def_property_readonly("compact_mode", &hey::taylor_adaptive::get_compact_mode) .def_property_readonly("high_accuracy", &hey::taylor_adaptive::get_high_accuracy) .def_property_readonly("with_events", &hey::taylor_adaptive::with_events) + // Variational-specific bits. + .def_property_readonly("n_orig_sv", &hey::taylor_adaptive::get_n_orig_sv) + .def_property_readonly("is_variational", &hey::taylor_adaptive::is_variational) + .def_property_readonly("vargs", &hey::taylor_adaptive::get_vargs) + .def_property_readonly("vorder", &hey::taylor_adaptive::get_vorder) + .def_property_readonly("tstate", &fetch_tstate) + .def( + "get_vslice", + [](const hey::taylor_adaptive &ta, std::uint32_t order, std::optional component) { + const auto ret = component ? ta.get_vslice(*component, order) : ta.get_vslice(order); + + return py::slice(boost::numeric_cast(ret.first), + boost::numeric_cast(ret.second), {}); + }, + "order"_a, "component"_a = py::none{}) + .def( + "get_mindex", + [](const hey::taylor_adaptive &ta, std::uint32_t i) { + const auto &ret = ta.get_mindex(i); + + return dtens_t_it::sparse_to_dense( + ret, boost::numeric_cast(ta.get_vargs().size())); + }, + "i"_a) + .def("eval_taylor_map", + [](py::object &o, std::variant in) { + auto *ta = py::cast *>(o); + + // Fetch the dtype corresponding to T. + const auto dt = get_dtype(); + + // Fetch the inputs array. + // NOTE: this will either fetch the existing array, or convert + // the input iterable to a py::array on the fly. + const auto inputs = std::visit( + [&](U &v) { + if constexpr (std::same_as) { + // Check the dtype. + if (v.dtype().num() != dt) [[unlikely]] { + py_throw( + PyExc_TypeError, + fmt::format( + "Invalid dtype detected for the inputs of a Taylor map evaluation: the " + "expected dtype is '{}', but the dtype of the inputs array is '{}' instead", + str(py::dtype(dt)), str(v.dtype())) + .c_str()); + } + + // Check that the inputs array is a C-style contiguous array. + if (!is_npy_array_carray(v)) [[unlikely]] { + py_throw(PyExc_ValueError, + "Invalid inputs array detected in a Taylor map evaluation: the array is not " + "C-style contiguous, please " + "consider using numpy.ascontiguousarray() to turn it into one"); + } + + return std::move(v); + } else { + return as_carray(v, dt); + } + }, + in); + + // Validate the number of dimensions for the inputs. + // NOTE: this needs to be done regardless of the original type of in. + if (inputs.ndim() != 1) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} dimensions, " + "but it must have 1 dimension instead", + inputs.ndim()) + .c_str()); + } + + // Validate the shape for the inputs. + if (boost::numeric_cast(inputs.shape(0)) != ta->get_n_orig_sv()) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} elements, " + "but it must have {} elements instead", + inputs.shape(0), ta->get_n_orig_sv()) + .c_str()); + } + +#if defined(HEYOKA_HAVE_REAL) + + // Run the checks specific for mppp::real. + if constexpr (std::same_as) { + // Check that the inputs array contains values with the correct precision. + pyreal_check_array(inputs, ta->get_prec()); + } + +#endif + + // Run the evaluation. + const auto *data_ptr = static_cast(inputs.data()); + ta->eval_taylor_map(std::ranges::subrange(data_ptr, data_ptr + inputs.shape(0))); + + // Return the tstate. + return fetch_tstate(o); + }) // Event detection. .def_property_readonly("with_events", &hey::taylor_adaptive::with_events) .def_property_readonly("te_cooldowns", &hey::taylor_adaptive::get_te_cooldowns) diff --git a/heyoka/test.py b/heyoka/test.py index 7f186fff..c263dfc5 100644 --- a/heyoka/test.py +++ b/heyoka/test.py @@ -2079,6 +2079,8 @@ def run_test_suite(): _test_vsop2013, _test_elp2000, _test_lagham, + _test_var_ode_sys, + _test_var_integrator, ) import numpy as np from .model import nbody @@ -2098,6 +2100,10 @@ def run_test_suite(): suite.addTest( tl.loadTestsFromTestCase(_test_batch_integrator.batch_integrator_test_case) ) + suite.addTest(tl.loadTestsFromTestCase(_test_var_ode_sys.var_ode_sys_test_case)) + suite.addTest( + tl.loadTestsFromTestCase(_test_var_integrator.var_integrator_test_case) + ) suite.addTest(tl.loadTestsFromTestCase(_test_lagham.lagham_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_vsop2013.vsop2013_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_elp2000.elp2000_test_case)) diff --git a/tools/circleci_conda_heyoka_head_310.sh b/tools/circleci_conda_heyoka_head_310.sh index a92d60f0..499fde14 100644 --- a/tools/circleci_conda_heyoka_head_310.sh +++ b/tools/circleci_conda_heyoka_head_310.sh @@ -45,7 +45,8 @@ cd $HEYOKA_PY_PROJECT_DIR cd doc -make html linkcheck doctest +# make html linkcheck doctest +make html doctest set +e set +x diff --git a/tools/circleci_conda_heyoka_head_39.sh b/tools/circleci_conda_heyoka_head_39.sh index 8fa01aa9..c353beeb 100644 --- a/tools/circleci_conda_heyoka_head_39.sh +++ b/tools/circleci_conda_heyoka_head_39.sh @@ -45,7 +45,8 @@ cd $HEYOKA_PY_PROJECT_DIR cd doc -make html linkcheck doctest +#make html linkcheck doctest +make html doctest set +e set +x diff --git a/tools/gha_conda_docs.sh b/tools/gha_conda_docs.sh index ff24a2df..8d410ecb 100644 --- a/tools/gha_conda_docs.sh +++ b/tools/gha_conda_docs.sh @@ -58,12 +58,7 @@ cd $HEYOKA_PY_PROJECT_DIR cd doc -# NOTE: run linkcheck only if we are on a pull request. -if [[ -n "${GITHUB_HEAD_REF}" ]]; then - make html linkcheck doctest -else - make html doctest -fi +make html doctest set +e set +x