Sympy autowrap (cython): defining “helpers” for sy

Question

I have a sympy.Matrix, called J_sym, that I would like to autowrap (preferably using the cython backend); the according symbols are stored in the list list_args.

However, the problem I run into is that apparently some sympy functions are not supported, in my case specifically sympy.Max and sympy.Heaviside.

Specifically, if I try

J_num_autowrap = autowrap(J_sym, backend="cython", args=list_args)

I get the following:

[...]
wrapped_code_10.c(4): warning C4013: 'Heaviside' undefined; assuming extern returning int
wrapped_code_10.c(4): warning C4013: 'Max' undefined; assuming extern returning int
[...]
: fatal error LNK1120: 2 unresolved externals
[...]
failed with exit status 1120

(a similar problem arises when one lambdifies using only numpy; but there one can easily pass a dictionary for these functions.)

Now, it seems that the "helpers" argument to autowrap should offer a solution; however, I can't work out how to properly use it. The docstring itself is absolutely cryptic, and I found only this link with an example:
https://github.com/sympy/sympy/issues/10572

I'm unsure how to properly use helpers here - can anyone help?

If I try something along these lines:

J_num_autowrap = autowrap(J_sym, backend="cython", args=list_args, helpers=("Max", max(x,y), [x,y]))

I get (a) an error if x, y are not declared - how can I do this with free variables? And (b) if I do the same after declaring x and y as sympy symbols, I just get the error "Truth value of relational can not be determined".

(The above link mentions that it is impossible to pass more than one helper anyway - even though it seems to me that this was to be fixed in 1.0. I still get the "ValueError: not enough values to unpack (expected 3, got 2)" when trying to pass 2 (gibberish) helper tuples. Not sure about the status - maybe I'm just out of luck here).

If there is any alternative way to generate C code from my matrix, I would appreciate any hints in that direction as well.

Answer 1

As a workaround for Max issue, you can use

helpers=("Max", (abs(x+y) + abs(x-y))/2, [x, y])

The expression (abs(x+y) + abs(x-y))/2 is mathematically equivalent to max(x, y) and does not generate any complaints from autowrap. A complete example:

from sympy import *
from sympy.utilities.autowrap import autowrap
x, y = symbols('x y')
f = autowrap(Max(2*x, y+1), args=[x, y], backend="cython", helpers=("Max", (abs(x+y) + abs(x-y))/2, [x, y]))
print([f(5, 6), f(1, 2)])    # outputs 10 and 3 

Similarly, Heaviside(x) is the same as (x + abs(x))/(2*x) -- except the latter expression is NaN when x=0. The uglier but safe version is (x + abs(x))/(2*abs(x) + 1e-300) where the added 1e-300 will almost never change the result.

Problem with multiple helpers

You want helpers for both Max and Heaviside. And this is where you hit an open issue: there is a bug in autowrap which makes it impossible to use more than one helper. Documentation suggests that the format would be like

helpers=[("Max", ..., [x, y]), ("Heaviside", ..., [x])]

But on this line autowrap does something that is convenient for one helper (we don't have to put it in a list), but fatal for more than one (extra layer of wrapping):

helpers = [helpers] if helpers else ()

Naturally, subsequent unpacking for name_h, expr_h, args_h in helpers fails.

Workaround for multiple helpers

The ufuncify handles helpers argument correctly. Unfortunately, ufuncify eventually calls autowrap for all backends other than NumPy, so the error still occurs in autowrap. But if you are willing to use NumPy backend, this is a solution:

from sympy.utilities.autowrap import ufuncify
x, y = symbols('x y', real=True)
my_helpers = [("Max", abs(x+y)/2 + abs(x-y)/2, [x, y]), ("Heaviside", (x + abs(x)) / (2*abs(x) + 1e-300), [x])]
f = ufuncify([x,y], Max(2*x, y+1) + Heaviside(x-4), backend="numpy", helpers=my_helpers)
print([f(5, 6), f(1, 2)])   $ outputs 11  and  3

Reducing to one helper

If you can change the expression to eliminate one of Max and Heaviside (or even both), Cython backend can be used. For example, maybe you only need Max(x, 0) and so you define a Python function "positive part":

pos = lambda x: (x+abs(x))/2

Then the autowrap has no problem with pos, and you only need to help with Heaviside:

f = autowrap(pos(9-x-y) + Heaviside(x-4), args = [x, y], backend = "cython", helpers=("Heaviside", (x + abs(x)) / (2*abs(x) + 1e-300), [x]))
print([f(5, 6), f(1, 2)])   #  outputs 1 and 6

How practical such replacements are depends on how your symbolic expression is obtained.