Reverse Polish notation and Derivatives

With this code, we try to calculate the symbolic derivative, say d(y*x)/dx = y, with as simple means as possible, without using external libraries.

Specifically, this project aims to automatically calculate the Jacobian of a vector valued function used as right hand side in ordinary differential equations. For this purpose we wrote ode.sh. The ode.sh script uses the derivative binary of this package by default, but can also try to use maxima or yacas instead.

There are many software packages that can do this and much more (computer algebra systems). Here is a list of alternatives.

Summary

Especially the derivative program requires the input to be math in reverse Polish notation (rpn), e.g.: 1 2 + (for 1+2); this notation does not need, nor allow parentheses.

The dc program works with this notation, so it may be familiar to the reader.

RPN math is easier (citation needed) to process (at least for what we do here) as there are no parentheses and evaluation follows a very simple algorithm.

The input expression must be balanced, unlike: x y z +, which has an unused operand, or x +, which lacks an operand (dc would not care much about the former and print a warning about the latter).

An unbalanced expression produces an empty line in the output of any program in this repository (or possibly a line with an error message).

All programs read from stdin. In any case, we try to keep the line numbers of input and output on stdout the same (the derivative of input line n will be in output line n). There may be output on stderr, so redirection is necessary: 2>error.log.

Although it is of course possible to use derivative and simplify by themselves, you can also use to_rpn to convert more conventional math to rpn notation and to_infix to translate the result back. All standard math functions must have an @ to mark them as functions: -1 x * @exp is equivalent to exp(-x) in normal notation. This may change in the future.

Compiling

Both gcc and tcc will work and have been used to compile these sources. In the root directory of this repository (parent of src):

$ mkdir bin
$ make
$ make test

These commands will create the binaries derivative, simplify, and ll_test; make test will run test.sh.

The symbol component takes a macro that fixes the maximum length of variable names: MAX_NAME_SIZE. It can be defined on the command line, when compiling: gcc -DMAX_NAME_SIZE=10 [...], see the Makefile.

Installation is optional, see this note.

Installation

The default target directory for installation is /usr/local/bin:

$ sudo make install
$ # or perhaps
$ doas make install

Afterwards, these commands should work:

$ which derivative
/usr/local/bin/derivative
$ man derivative

If you prefer a different location, edit the PREFIX variable in Makefile.

Usage

The programs in this repository are meant to be used via pipes:

[...] | to_rpn | derivative x | simplify

They all read n lines from standard input (stdin) and write n lines to standard output (stdout) - unless something goes wrong.

Conversion to RPN

The program to_rpn reads stdin for math expressions in infix notation (the normal kind) and tries to convert them into reverse Polish notation:

$ bin/to_rpn < math.txt

This program takes no command line arguments. If the math expressions come in larger numbers (or from a pipe), from a source that just uses function names without any @ signs, it is probably convenient to do sth like this on the input:

sed -r  's/(exp|log|sin|cos|tan)/@\1/g' original_input.txt | to_rpn

to replace exp with @exp (and similar).

Caveat: Ambiguity and Spacing

This program tries to convert an item it encounters to a number first, only if that fails does it consider other possibilities: the item could then be an operator, function, or variable name. This creates an ambiguity with the unary minus (and plus). Consequently, some strings are not interpreted right: 1-2 is read as +1 and -2 (a sequence of two numbers), but 1 - 2 is understood correctly as the operation minus(1,2), a-b works perfectly fine because b is not a valid number. The bad: unary minus does not work at all for variables -a+1.

The dc program resolves this ambiguity by using the underscore to denote negative numbers: dc -e '_1 1 + p' prints 0; but, we don't do that underscore thing.

The output consists of space separated rpn expressions, so no mixups of this sort are likely to happen downstream (apart from negated variables, the doesn't work in rpn either).

derivative

$ cd RPN-Derivative
$ to_rpn < math.txt | derivative t

The derivative function takes one mandatory command line argument, the name of a variable t. The derivative of the input expressions will be calculated with respect to the named variable.

simplify

simplify tries to reduce the length of an rpn expression by eliminating obviously unnecessary operations like x 0 +.

$ echo "x 0 0 + +" | simplify [n]

This program has one optional parameter n (an integer): simplification will be repeated n times. This is equivalent to repeated calls to simplify:

$ echo "x 0 0 + +" | simplify
x 0 +
$ echo "x 0 0 + +" | simplify | simplify
x
$ echo "x 0 0 + +" | simplify 2
x

Simplify tries to reduce fractions by finding common factors. If a common factor is found within a sum, it will be factored out. Not all obvious common factors are detected.

Note

To calculate the derivative of x*y/(a+b) with respect to y, we would type:

$ echo "x y * a b + /" | derivative y | simplify 3

this produces the output:

x a b + * a b + a b + * /

which is x*(a+b)/((a+b)*(a+b)) in infix notation (so x/(a+b)). This is almost right (technically correct, but unnecessary). This fraction is reduced correctly, now Currently, no amount of simplifying will reduce the fraction unless it reduces to 1; This is still true: simplify doesn't really know math, it merely tries some simple pattern recognition.

If you encounter a specific case where you automate code creation for a computational task and it always produces a difficult type of fraction, you can use sed and regular expressions to further manipulate the result. The file reduce_fraction.sed contains some examples:

$ echo '(pow(x,3)*(4/x))' | sed -E -f sh/reduce_fraction.sed
(pow(x,(3)-1)*4)
$ echo '(pow(x,3)/x)' | sed -E -f sh/reduce_fraction.sed
(pow(x,(3)-1))

a simple case

In this example d[x*y*(y+2)]/dy, the derivative is just not problematic at all and does not require any special treatment:

$ echo "x y * y 2 + *" | derivative y | simplify 2
x y 2 + * x y * +

which is x*(y+2) + x*y (OK).

rpn to infix

The output of the programs lilsted in previous Sections can be converted back into infix notation like this:

$ bin/to_infix < rpn_math.txt

Optional parameters

-s (safe fractions) --safe-fractions=1e-10

This option is useful when denominators are known/required/assumed to be non-negative. With -s, to_infix will add a small safety constant to denominators to prevent division by zero. This could of course be terribly wrong, because it does change the printed expression, we assume that the user is aware of the risks.

In part, this is a safety measure for cases where reduction of fractions didn't work:

x*x/x

The above is a perfectly fine fraction, but it cannot be naïvely evaluated at x=0 (numerically); it has to be reduced to x. In contrast, a similar fraction

x*x/(1e-16 + x)

would work numerically without smart fraction reduction for all arguments x ≥ 0. And finally the expression:

x*x/(1e-16 + fabs(x))

works for all x (but is not equal to the original x everywhere, or indeed anywhere). So, a procedure now can call the resulting math more carelessly and not crash while doing so.

By default the safety constant is the smallest positive double precision floating point value (rounded by printing):

echo 'x x * x /' | to_infix -s
((x*x)/(2.22507e-308 + fabs(x)))

Obviously, this should only be used when stability is more important than accuracy: i.e. the function must never crash. This is sometimes the case with optimization (a model is called thousands of times in a row, with sketchy trial arguments). And the above example is for illustration only, because:

$ echo 'x x * x /' | simplify 2 | to_infix
x

Example

Like all other programs in this repository, to_infix reads from standard input and writes to standard output.

$ echo "@pow(t,3)" | to_rpn | derivative t | simplify 4 | to_infix
(pow(t,3)*(3/t))

There will be many superfluous parentheses to produce safe expressions (they can be used as sub-expressions). The @ symbol will be stripped from functions to make it easier to use the expressions elsewhere.

ODE code generator

sh/ode.sh is a more general c-code generating script. It assumes that the working directory contains text files (tsv files), with specific names and content:

File	Mandatory	Content
Constant.txt	no	constant names and values
Parameters.txt	yes	parameter names and values
Variables.txt	yes	state variable names, initial values, and a unit or measurement
Expression.txt	no	math expression, names and formulae
Function.txt	no	terms that define a vector valued output function

See sh/README.md

Mathematical Functions

To make it easy convenient to parse math expressions, standard math functions shall be prefixed with an @. This is to avoid confusion with variable names. It would have been possible to distinguish functions and variable names by string matching and reserving a list of names for functions.

A function may be not implemented (yet) and thus its name be understood as a variable name (e.g. sinh). But with a leading @ an error will occur if the function is not known.

Currently: @exp, @log, @sin, @cos, @pow are known.

Limitations: Many

0. There is no unary minus for variables: -a+b doesn't work, b-a does obviously as does 0-a+b
1. Very few operator symbols and functions are and will ever be supported. All logical opertors are missing, bitwise operators and integer arithmetic (e.g. remainder/modulus) as they are difficult to differentiate.
2. There are very few checks to see if the input is a valid RPN expression.
3. All variables must be one letter the maximum length of variables can be set at compile time.
4. Conversion back into infix notation writes many unnecessary parentheses; this is probably fine.
5. Because simplify doesn't reliably reduce all fractions, some perfectly finite derivatives cannot be evaluated at, e.g.: x=0.

This tool-set will never become a scripting language like dc, with multiple registers and user defined macros (or other things that are too hard for me to implement).

Documentation

The doc folder contains developer documentation in html form, created with codedoc.

The src folder contains further notes on the code organization.

Small Tests

No formal testing framework is used here, but make test will run a series of tests from a shell script.

It is possible to manually test simple cases numerically, without external software. We use dc to do this, as it is probably installed on every unix like/related system. The dc program wants reverse polish notation as input, so it is perfectly well suited to check output from bin/derivative before it has been converted to R or C code (see the next Section).

There are two shell scripts in the tests folder that help with the testing:

1. numerical.sh
2. eval.sh

The first script calculates the finite difference approximation of a derivative a given an rpn expression for f: (f(x+h)-f(x-h))/(2*h); the second script evaluates any rpn expression using dc (it does all necessary substitutionsso that dc will accept the expression), e.g.:

$ echo "x a @pow" | tests/numerical.sh x 2 0.0001 | tests/eval.sh a=3
12.0000005000

The above instruction will calculate the finite difference at x=2 and h=1e-4 and pass the expression on to eval.sh, which in turn will substitute all occurences of a with 3.

Finally, eval.sh calls dc with the following instruction:

2.0001 3 ^ 1.9999 3 ^ - 2 0.0001 * / p

The last dc command item, p, prints the result (12):

$ dc -e '2.0001 3 ^ 1.9999 3 ^ - 2 0.0001 * / p'
12

which is correct (but was rounded). We can compare this to the analytical derivative:

$ echo 'x a @pow' | bin/derivative x | bin/simplify 3 | tests/eval.sh a=3 x=2
12.0000000000
$ echo "x a @pow" | tests/numerical.sh x 2 0.0001 | tests/eval.sh a=3
12.0000005000

note normally dc will floor fractions. To avoid heavy losses in accuracy, we set dc's precision (k) to 10 digits: dc -e '3 2 / p' prints 1 (because floor(3/2) is 1); dc -e '3 k 3 2 / p' prints 1.500 with precision k=3.