{ "cells": [ { "cell_type": "markdown", "id": "b8b38e6c-6c1e-4875-b2e5-03bd488c9d29", "metadata": {}, "source": [ "# Interoperability with SymPy\n", "\n", "```{versionadded} 0.10.0\n", "\n", "```\n", "\n", "Starting with version 0.10, symbolic expressions in heyoka.py can be converted to/from [SymPy](https://www.sympy.org/en/index.html) expressions. This feature can be very useful, as it unlocks the possibility of using the wide arsenal of SymPy algorithms on heyoka.py's expressions. Possible use cases include:\n", "\n", "- the automatic simplification of heyoka.py's expressions,\n", "- the application of symbolic algorithms not available in heyoka.py (e.g., symbolic integration),\n", "- LaTeX pretty-printing.\n", "\n", "In this tutorial, we will show how SymPy interoperability can be used for expression simplification purposes, which ultimately leads to a measurable performance increase in the numerical integration of a system of ODEs.\n", "\n", "For this example, we will consider the relativistic dynamics of Mercury's orbit around the Sun, which, in the first post-Newtonian approximation, can be described by the following Hamiltonian:\n", "\n", "$$\n", "\\mathcal{H}_\\mathrm{1PN} \\left(v_x, v_y, v_z, x, y, z \\right) = \\frac{1}{2}v^2-\\frac{\\mu}{r} + \\varepsilon \\left(\\frac{\\mu^2}{2r^2} -\\frac{1}{8}v^4-\\frac{3}{2}\\frac{\\mu v^2}{r} \\right).\n", "$$\n", "\n", "where $\\left(v_x, v_y, v_z, x, y, z \\right)$ is Mercury's Cartesian state vector with respect to the Sun, $\\varepsilon = 1/c^2$, $v=\\sqrt{v_x^2+v_y^2+v_z^2}$, $r=\\sqrt{x^2+y^2+z^2}$ and $\\mu=GM_\\odot$ is the gravitational parameter of the system. See the [tutorial about Mercury's relativistic precession](<./mercury_precession.ipynb>) for more details about this dynamical system.\n", "\n", "Let us begin by defining the Hamiltonian as a symbolic expression in heyoka.py's expression system:" ] }, { "cell_type": "code", "execution_count": 1, "id": "508c0c08-1797-46bd-8c44-db71362baba5", "metadata": {}, "outputs": [], "source": [ "import heyoka as hy\n", "\n", "# Create the symbolic variables.\n", "vx, vy, vz, x, y, z = hy.make_vars(\"v_x\", \"v_y\", \"v_z\", \"x\", \"y\", \"z\")\n", "\n", "# mu and eps constants.\n", "mu = 0.01720209895 * 0.01720209895 * 365 * 365\n", "eps = 2.5037803127808595e-10\n", "\n", "# Auxiliary quantities.\n", "v2 = vx*vx + vy*vy + vz*vz\n", "r2 = x*x + y*y + z*z\n", "r = hy.sqrt(r2)\n", "\n", "# Initial conditions corresponding (roughly)\n", "# to Mercury's current orbit.\n", "ic = (-10.461611630114545,\n", " 6.667253667139184,\n", " 0.818635813947965,\n", " 0.16614243942411336,\n", " 0.2568228239702581,\n", " 0.0315338776710321)\n", "\n", "# Define the Hamiltonian.\n", "Ham = 1./2*v2 - mu/r + eps * (mu**2/(2*r2) - 1/8.*v2*v2 - 3./2.*mu*v2/r)" ] }, { "cell_type": "markdown", "id": "7a7d0a39-90b6-4cbe-9d79-dc7ac8b7f566", "metadata": {}, "source": [ "As usual, we are using Solar masses, astronomical units and years as units of measure.\n", "\n", "Let us now convert the heyoka.py Hamiltonian to a SymPy expression via the ``to_sympy()`` function and pretty-print it to screen:" ] }, { "cell_type": "code", "execution_count": 2, "id": "b197a3f3-b8f9-4961-b41d-f683c169695b", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle 0.5 v_{x}^{2} + 0.5 v_{y}^{2} + 0.5 v_{z}^{2} - 2.50378031278086 \\cdot 10^{-10} \\cdot \\left(0.125 v_{x}^{2} + 0.125 v_{y}^{2} + 0.125 v_{z}^{2}\\right) \\left(v_{x}^{2} + v_{y}^{2} + v_{z}^{2}\\right) - \\frac{2.50378031278086 \\cdot 10^{-10} \\cdot \\left(59.1343559232718 v_{x}^{2} + 59.1343559232718 v_{y}^{2} + 59.1343559232718 v_{z}^{2}\\right)}{\\sqrt{x^{2} + y^{2} + z^{2}}} + \\frac{3.89128862055816 \\cdot 10^{-7}}{2 x^{2} + 2 y^{2} + 2 z^{2}} - \\frac{39.4229039488479}{\\sqrt{x^{2} + y^{2} + z^{2}}}$" ], "text/plain": [ "0.5*v_x**2 + 0.5*v_y**2 + 0.5*v_z**2 - 2.50378031278086e-10*(0.125*v_x**2 + 0.125*v_y**2 + 0.125*v_z**2)*(v_x**2 + v_y**2 + v_z**2) - 2.50378031278086e-10*(59.1343559232718*v_x**2 + 59.1343559232718*v_y**2 + 59.1343559232718*v_z**2)/sqrt(x**2 + y**2 + z**2) + 3.89128862055816e-7/(2*x**2 + 2*y**2 + 2*z**2) - 39.4229039488479/sqrt(x**2 + y**2 + z**2)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hy.to_sympy(Ham)" ] }, { "cell_type": "markdown", "id": "28f61658-dcbc-40b0-817e-fe419cda9bb0", "metadata": {}, "source": [ "The SymPy expression can be converted back to heyoka.py using the ``from_sympy()`` function:" ] }, { "cell_type": "code", "execution_count": 3, "id": "38143cc5-b4b2-4d5f-a0ab-b10b4c1a58d8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((0.50000000000000000 * v_x**2.0000000000000000) + (0.50000000000000000 * v_y**2.0000000000000000) + (0.50000000000000000 * v_z**2.0000000000000000) + (3.8912886205581639e-07 / ((2.0000000000000000 * x**2.0000000000000000) + (2.0000000000000000 * y**2.0000000000000000) + (2.0000000000000000 * z**2.0000000000000000))) - (39.422903948847882 / (x**2.0000000000000000 + y**2.0000000000000000 + z**2.0000000000000000)**0.50000000000000000) - ((2.5037803127808595e-10 * ((59.134355923271826 * v_x**2.0000000000000000) + (59.134355923271826 * v_y**2.0000000000000000) + (59.134355923271826 * v_z**2.0000000000000000))) / (x**2.0000000000000000 + y**2.0000000000000000 + z**2.0000000000000000)**0.50000000000000000) - (2.5037803127808595e-10 * (v_x**2.0000000000000000 + v_y**2.0000000000000000 + v_z**2.0000000000000000) * ((0.12500000000000000 * v_x**2.0000000000000000) + (0.12500000000000000 * v_y**2.0000000000000000) + (0.12500000000000000 * v_z**2.0000000000000000))))" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hy.from_sympy(hy.to_sympy(Ham))" ] }, { "cell_type": "markdown", "id": "e13a4028-5d36-4c9a-a836-004f022c9b96", "metadata": {}, "source": [ "In order to formulate Hamilton's equations for this system we need, as usual, to differentiate the Hamiltonian with respect to the state variables. Let us do it with heyoka.py's ``diff()`` function:" ] }, { "cell_type": "code", "execution_count": 4, "id": "5cd6dd1b-5ca2-4ccf-8ecf-ef6f0fa80074", "metadata": {}, "outputs": [], "source": [ "# Formulate the equations of motion\n", "# using heyoka.py's expression system.\n", "ta = hy.taylor_adaptive(\n", " # Hamilton's equations.\n", " [(vx, -hy.diff(Ham, x)),\n", " (vy, -hy.diff(Ham, y)),\n", " (vz, -hy.diff(Ham, z)),\n", " (x, hy.diff(Ham, vx)),\n", " (y, hy.diff(Ham, vy)),\n", " (z, hy.diff(Ham, vz))],\n", " # Initial conditions.\n", " ic\n", ")" ] }, { "cell_type": "markdown", "id": "4d43a2c7-f598-4f1b-98f4-970e7bb8443a", "metadata": {}, "source": [ "heyoka.py's ``diff()`` function applies elementary calculus rules in order to differentiate a symbolic expression, performing only very basic simplifications in the process. We can get an estimate of the complexity of the resulting expressions by taking a look at the size of the Taylor decomposition for this dynamical system (this is the number of elementary subexpressions in which the equations of motion have been decomposed):" ] }, { "cell_type": "code", "execution_count": 5, "id": "6184e14a-5950-415a-aab5-9320f926673b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "92" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NOTE: subtract 12 to avoid counting the 6 definitions of the state variables\n", "# and the 6 definitions of their differential equations.\n", "len(ta.decomposition) - 12" ] }, { "cell_type": "markdown", "id": "27820d06-ae60-46f5-a0ba-a73e8431f9e4", "metadata": {}, "source": [ "Let us now try to simplify Hamilton's equations of motion using SymPy. We first define a small wrapper to invoke SymPy's ``simplify()`` function on a heyoka.py expression:" ] }, { "cell_type": "code", "execution_count": 6, "id": "a106c52f-0e13-4740-932d-09562edc7858", "metadata": {}, "outputs": [], "source": [ "import sympy\n", "\n", "def simplify(ex):\n", " return hy.from_sympy(sympy.simplify(hy.to_sympy(ex)))" ] }, { "cell_type": "markdown", "id": "c126f7b2-9819-497f-85e4-2674703cb8ab", "metadata": {}, "source": [ "Then, we create a new integrator with the simplified equations of motion:" ] }, { "cell_type": "code", "execution_count": 7, "id": "7196a47f-c1e8-48a6-b4a4-b29c175d215a", "metadata": {}, "outputs": [], "source": [ "ta_spy = hy.taylor_adaptive(\n", " # Hamilton's equations.\n", " [(vx, simplify(-hy.diff(Ham, x))),\n", " (vy, simplify(-hy.diff(Ham, y))),\n", " (vz, simplify(-hy.diff(Ham, z))),\n", " (x, simplify(hy.diff(Ham, vx))),\n", " (y, simplify(hy.diff(Ham, vy))),\n", " (z, simplify(hy.diff(Ham, vz)))],\n", " # Initial conditions.\n", " ic\n", ")" ] }, { "cell_type": "markdown", "id": "ba672d92-1ade-41ae-8930-cfcb791b3ebb", "metadata": {}, "source": [ "Let us take a look at the size of the Taylor decomposition for the new integrator:" ] }, { "cell_type": "code", "execution_count": 8, "id": "00aea412-94ae-4df3-a33c-d6fa60ffb3d3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "46" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(ta_spy.decomposition) - 12" ] }, { "cell_type": "markdown", "id": "9b941cd3-a85b-466d-af67-fb6d22a7e152", "metadata": {}, "source": [ "In other words, thanks to SymPy's expression simplification capabilities, we reduced by $\\sim 50\\%$ the symbolic complexity (i.e., the number of elementary subexpressions) of our ODE system.\n", "\n", "What is the performance impact of these simplifications? Let us find out:" ] }, { "cell_type": "code", "execution_count": 9, "id": "8994752f-e131-47d3-9fb8-bcafbc6a9532", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2.03 s, sys: 0 ns, total: 2.03 s\n", "Wall time: 2.03 s\n" ] }, { "data": { "text/plain": [ "(,\n", " 5.571678500915182e-06,\n", " 0.017389385665743787,\n", " 1030387,\n", " None,\n", " None)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "time ta_spy.propagate_until(10000.)" ] }, { "cell_type": "code", "execution_count": 10, "id": "085f7a99-3bdd-40bb-a258-629cc82248f3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2.23 s, sys: 57 µs, total: 2.23 s\n", "Wall time: 2.23 s\n" ] }, { "data": { "text/plain": [ "(,\n", " 0.006522499652827666,\n", " 0.01738986959999157,\n", " 1012008,\n", " None,\n", " None)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "time ta.propagate_until(10000.)" ] }, { "cell_type": "markdown", "id": "fd7375e4-ca0a-4090-99c0-f60afe2ba6d3", "metadata": {}, "source": [ "The performance increase is measurable but not as large as the reduction in symbolic complexity.\n", "\n", "As a word of warning, users should be aware that the general-purpose ``simplify()`` function from SymPy can be quite expensive from a computational point of view, especially for large expressions. SymPy provides various specialised [simplification primitives](https://docs.sympy.org/latest/tutorial/simplification.html) that can be considerably more efficient.\n", "\n", "## More details on the conversion process\n", "\n", "When converting a heyoka.py expression to SymPy, the following conventions are adopted:\n", "\n", "* variables are converted to SymPy ``Symbol``s,\n", "* constants are converted to SymPy numbers,\n", "* [runtime parameters](<./ODEs with parameters.ipynb>) are converted to SymPy ``Symbol``s conventionally called ``\"par[n]\"``, where ``n`` is the parameter index,\n", "* functions are mapped to the corresponding SymPy functions.\n", "\n", "Similarly, when converting a SymPy expression to heyoka.py, the following conventions are adopted:\n", "\n", "* SymPy ``Symbol``s are converted to heyoka.py variables, unless the name of the symbol is ``\"par[n]\"``, in which case the ``Symbol`` is converted to a runtime parameter,\n", "* SymPy numbers are converted to heyoka.py constants, if the conversion is lossless (otherwise an error is raised),\n", "* SymPy functions are mapped to heyoka.py functions.\n", "\n", "When converting a SymPy expression to heyoka.py, it is possible to customise the conversion process by passing a substitution dictionary to the ``from_sympy()`` function. Here's an example:" ] }, { "cell_type": "code", "execution_count": 11, "id": "205a39c6-cb0d-4b66-a843-43d1e8148b22", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(x * t)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create a few SymPy symbols.\n", "x, y, z = sympy.symbols(\"x y z\")\n", "\n", "# Convert the expression x*(y+z) to heyoka.py,\n", "# but provide a custom substitution rule for\n", "# the subexpression y+z.\n", "hy.from_sympy(x * (y+z), {y+z: hy.expression(\"t\")})" ] }, { "cell_type": "markdown", "id": "5ac12b1e-11f8-4de9-94bd-499d785a1710", "metadata": {}, "source": [ "Note that the same effect could be achieved by doing the substitution in SymPy *before* the conversion to heyoka.py. However, for large expressions, using the custom substitution dictionary in the ``from_sympy()`` function will perform better.\n", "\n", "## Handling rational numbers\n", "\n", "heyoka.py is very conservative when converting from SymPy expressions containing rational numbers. In particular, if a rational number cannot be converted *exactly* to a floating-point value, heyoka.py will raise an error. In order to avoid this, you can use the ``sympy.N`` function to forcibly convert all rational numbers to floating-point values:" ] }, { "cell_type": "code", "execution_count": 12, "id": "ccbcb4c4-8a3c-4b3e-a31d-3ac7eb2d1672", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.33333333333333331 * x)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hy.from_sympy(sympy.N(x / 3))" ] } ], "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 }