Notebook

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.

Constants and parameters are mathematically equivalent, the only 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.

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

In [1]:

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*x + y*y + z*z)))

The euclidean distance is: sqrt(((x**2 + y**2) + z**2))

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("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
Extended-precision 1.1: 1.10000000000000000002
Quadruple-precision 1.1: 1.10000000000000000000000000000000008
Octuple-precision 1.1: 1.100000000000000000000000000000000000000000000000000000000000000000000004

Note that, while double precision is always supported in heyoka.py via the standard float Python type, support for extended-precision floating-point types varies depending on the software/hardware platform. Specifically:

on x86 processors, the 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., Unix ARM 64), the longdouble type implements the IEEE quadruple-precision floating-point format;
on some platforms where 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 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 the following elementary functions (with more to come in the near future):

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,
the inverse of Kepler's elliptic equation.

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 very limited. 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.

Note that, 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.