Introduction to the expression system¶

As we saw in the previous section, heyoka.py needs to be able to represent the right-hand side of an ODE system in symbolic form in order to be able to compute its high-order derivatives via automatic differentiation. heyoka.py represents generic mathematical expressions via a simple abstract syntax tree (AST) in which the internal nodes are n-ary functions and the leaf nodes can be:

symbolic variables,
numerical constants,
runtime parameters.

Both constants and parameters can be used to represent mathematical constants, the difference being that the value of a constant is determined when the expression is created, whereas the value of a parameter is loaded from a user-supplied data array at a later stage. Additionally, derivatives can be taken with respect to parameters.

The construction of the AST of an expression in heyoka.py can be accomplished via natural mathematical notation:

In [2]:

import heyoka as hy

# Define the symbolic variables x and y.
x, y = hy.make_vars("x", "y")

# Another way of creating a symbolic variable.
z = hy.expression("z")

# Create and print an expression.
print("The euclidean distance is: {}".format(hy.sqrt(x**2 + y**2)))

The euclidean distance is: (x**2.0000000000000000 + y**2.0000000000000000)**0.50000000000000000

Numerical constants can be created using any of the floating-point types supported by heyoka.py. For instance, on a typical Linux installation of heyoka.py on an x86 processor, one may write:

In [2]:

print("Double-precision 1.1: {}".format(hy.expression(1.1)))

import numpy as np
print("Single-precision 1.1: {}".format(hy.expression(np.float32("1.1"))))

print("Extended-precision 1.1: {}".format(hy.expression(np.longdouble("1.1"))))

print("Quadruple-precision 1.1: {}".format(hy.expression(hy.real128("1.1"))))

# NOTE: octuple precision has a
# 237-bit significand.
print("Octuple-precision 1.1: {}".format(hy.expression(hy.real("1.1", 237))))

Double-precision 1.1: 1.1000000000000001
Single-precision 1.1: 1.10000002
Extended-precision 1.1: 1.10000000000000000002
Quadruple-precision 1.1: 1.10000000000000000000000000000000008
Octuple-precision 1.1: 1.100000000000000000000000000000000000000000000000000000000000000000000004

Note that, while single and double precision are always supported in heyoka.py, support for extended-precision floating-point types varies depending on the software/hardware platform. Specifically:

on x86 processors, the NumPy {py:class}~numpy.longdouble type corresponds to 80-bit extended precision on most platforms (the exception being MSVC on Windows, where longdouble == float);
on some platforms (e.g., Linux ARM 64), the {py:class}~numpy.longdouble type implements the IEEE quadruple-precision floating-point format;
on some platforms where {py:class}~numpy.longdouble does not have quadruple precision, a nonstandard quadruple-precision type is instead available in C/C++ (this is the case, for instance, on x86-64 Linux and on some PowerPC platforms). On such platforms, and if the heyoka C++ library was compiled with support for the mp++ library, quadruple precision is supported via the real128 type (as shown above).

Note that the non-IEEE {py:class}~numpy.longdouble type available on some PowerPC platforms (which implements a double-length floating-point representation with 106 significant bits) is not supported by heyoka.py at this time.

Arbitrary-precision computations are supported by heyoka.py on all platforms via the real type, provided that the heyoka C++ library was compiled with support for the mp++ library. The real type implements a floating-point type whose precision can be set at runtime.

In addition to the standard mathematical operators, heyoka.py's expression system also supports several elementary and special functions, such as:

the square root,
exponentiation,
the basic trigonometric and hyperbolic functions, and their inverse counterparts,
the natural logarithm and exponential,
the standard logistic function (sigmoid),
the error function,
Kepler's elliptic anomaly and several other anomalies commonly used in astrodynamics.

In [3]:

# An expression involving a few elementary functions.
hy.cos(x + 2. * y) * hy.sqrt(z) - hy.exp(x)

Out[3]:

((cos((x + (2.0000000000000000 * y))) * sqrt(z)) - exp(x))

It must be emphasised that heyoka.py's expression system is not a full-fledged computer algebra system. In particular, its simplification capabilities are essentially non-existent. Because heyoka.py's performance is sensitive to the complexity of the ODEs, in order to achieve optimal performance it is important to ensure that the mathematical expressions supplied to heyoka.py are simplified as much as possible.

Starting form version 0.10, heyoka.py's expressions can be converted to/from SymPy expressions. It is thus possible to use SymPy for the automatic simplifcation of heyoka.py's expressions, and, more generally, to symbolically manipulate heyoka.py's expressions using the wide array of SymPy's capabilities. See the SymPy interoperability tutorial for a detailed example.