2

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.

DavidP
  • 105
  • 3
  • 12
  • I have no clue about how to get cython working correctly with sympy, but I assume you try this for performance reasons. If so, you might have a look at symengine: (https://github.com/symengine/symengine). Symengine is a project that strives to replace the core part of sympy with the symengine which is written in c++ and supposed to significantly speed up sympy's basic operations. It is still beta, so I don't know if it is of use for you at this stage. – laolux May 30 '17 at 21:05

1 Answers1

2

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.

  • Thank you for the extensive answer. While it is very clear, I must admit I am confused that `("Max", Max(x, y), [x, y])` should be doing anything in the first place - would that not amount to explaining `Max` by reference to `Max`? (Ignoring the problems with the helpers you mention.) Do I understand correctly that ufuncify could handle a list of helpers?). The workaround suggestion seems sensible - the matrix is a large Jacobian, so replacement would be very time intensive, but I could avoid `Max` from the get go. In my case, this will be even easier, all my terms have the form `Max(x,0)`. – DavidP May 31 '17 at 01:08
  • I revised the answer. You are right, replacing Max by Max was not a good idea: it only created an infinite recursion. `ufuncify` can handle multiple helpers, but only with Numpy backend. If you want Cython, you are limited to one helper - the other stuff would have to happen earlier. –  May 31 '17 at 02:11
  • Thanks once more. I got it to work by avoiding `Max` from the start; this also got rid of the `Heavisides` in the matrix, which just arised as part of the derivatives of the former - thus, I circumvented the helpers problem completely. The cythonized version is magnitudes faster using lambdify with numpy, so I'm quite happy. One last question is remaining, and I'm not sure it warrants an own post, so I'll just ask here: Is there any way to get the function to accept a numpy array, rather than a list? (I know I could just use f(*array), but it's slower and seems clumsy to me.) – DavidP May 31 '17 at 17:21
  • That's what `ufuncify` is for, to allow broadcasting. –  May 31 '17 at 18:42
  • Oh, I meant to pass all arguments for a single call as a numpy vector - similar to what one can do with lamdify+deferred vector. (Or actually, one has to do once one reaches 255 arguments . Is there a limit to the number of arguments for an autowrapped function when using cython?) – DavidP May 31 '17 at 19:03
  • Okay, now it does look like you have a separate question there. –  May 31 '17 at 21:18
  • I guess you are right, thanks though. Followup: https://stackoverflow.com/questions/44305644/sympy-autowrap-cython-limit-of-of-arguments-arguments-in-array-form – DavidP Jun 01 '17 at 11:01