array-broadcasting in sp.lambdify, Lambdify and LambdifyCSE

[raphaelquast]

Hi,
i tried to implement symengine's LambdifyCSE to speedup my code..
however, while investigating this possibility, i found that LambdifyCSE
seems to have different array-broadcasting than Lambdify (and also sympy's lambdify)

x = var('x')
funcs = [sp.cos(x)**i for i in range(5)]

#ordinary sympy
f_sympy = sp.lambdify((x), funcs, modules = 'numpy')

#ordinary symengine
f_seng = Lambdify([x], funcs)

#symengine with cse
f_seng_cse = LambdifyCSE([x], funcs)

ideally, the functions should give the same output for the same input...
however, in order to get the same output i have to do the following:
(aside of the fact that sympy's lambdify treats constants rather strange...)

a = np.array([1, 2, 3])

s1 = f_sympy(a)
s2 = f_seng(a).T
s3 = f_seng_cse(a[:,np.newaxis]).T

furthermore calling

f_seng_cse(a)

gives

"ValueError: all the input arrays must have same number of dimensions"

which is rather strange since there is only one input-array....

for higher dimensional inputs this generalizes to

a = np.array([[1,2,3],[1,2,3],[1,2,3]])
s1 = f_sympy(a)
s2 = f_seng(a).T
s3 = f_seng_cse(a[:,:,np.newaxis]).T

In order to get the same outputs with functions taking more than one input-arguments,
i have to use

x, y = var('x, y')

funcs = [sp.cos(x)**i + sp.cos(y)**i for i in range(5)]

f_sympy = sp.lambdify((x, y),funcs, modules = 'numpy')
f_seng = Lambdify([x, y], funcs)
f_seng_cse = LambdifyCSE([x, y], funcs)

a = np.array([[1,2,3],[1,2,3],[1,2,3]])

s1 = f_sympy(a, a)
s2 = np.reshape(f_seng(np.reshape(np.ravel([a,a], order = 'F'), (-1, 2), order = 'C')).T, (-1, *a.shape), order = 'F')
s3 = np.reshape(f_seng_cse(np.reshape(np.ravel([a,a], order = 'F'), (-1, 2), order = 'C')).T, (-1, *a.shape), order = 'F')

which is certainly not very readable any more...
is there any obvious reason for this or is this some kind of a bug?
or even better, could you eventually tell me a clear way on how to perform correct
array-broadcasting of multidimensional functions in LambdifyCSE ?

Thanks!

[bjodah]

Hi @raphaelquast,

Thank you for reporting this. I didn't see this issue until now -- feel free to ping me about Lambdify (you can see authors under git blame).

I would consider this a bug. LambifyCSE should work just like Lambdify (the user should not need to pass np.newaxis).

To call f_seng in your second example I would do it this way:

>>> inp = np.empty(a.shape + (2,))
>>> inp[..., 0] = a
>>> inp[..., 1] = a
>>> f_seng(inp)

(since both x & y are "a" it's not obvious to me if this gives you what you want?)

The assumptions made for broadcasting were only dictated by cache locality (and C ordered
output array)

Actually, the way Lambdify is implemented now is that a list of length 3 is treated as matrix of shape (1, 3) in some sense, e.g.:

In [42]: symengine.Lambdify([x], [x, x**2, x**3])([4, 5])
Out[42]: 
array([[   4.,   16.,   64.],
       [   5.,   25.,  125.]])

But if we would have passed the same list to Matrix it would have been interpreted as a (3, 1) shaped matrix:

In [44]: symengine.Matrix([x, x**2, x**3])
Out[44]: 
[x]
[x**2]
[x**3]

In [45]: symengine.Matrix([x, x**2, x**3]).shape
Out[45]: (3, 1)

That might be inconsistent enough that we should change what Lambdify does here.
That is: right now the extra input shape is appended to the "left end" of the shape tuple, perhaps it should be the "right end" (then we should change the ordering of the output array to column major).

[raphaelquast]

hi,
thanks for the answer (and sorry for taking so long to respond)
the examples i provided have just been the most minimal examples i could
think of that would reproduce the problem. (which resulted in using only
a single array 'a' for all inputs out of lazyness)

my actual task can be described as follows:

i need to evaluate a list of ~10-20 rather lengthy functions (>100 000 characters as string)
that are all structured in such a way that CSE would definitely make a lot of sense...
the functions have at least 4 input-arguments that will be either scalars, 1d- or 2d arrays
(depending on user-input) and thus i need reliable array-broadcasting to actually evaluate them.

for the time being i got back to using symengine just to perform the symbolic evaluation
of the functions, and then sympify the results to use sympy's cse and lambdify to generate
my own version of cse-based functions which ended up being a more robust and fast
approach than using symengines Lambdify-possibilities to take over the evaluation.

the major drawback here is that sympifying takes up to some minutes to convert the
derived symengine-based formulas back to sympy-based formulas, but symengines power in
performing fast symbolic equation-expansions still pays off...

i will keep an eye on the development of symengines Lambdify-possibilities, and as soon
as the broadcasting-rules can be considered "save" i'll give it another try...

[raphaelquast]

@bjodah
Hey, i've now finally found the time to install the dev-version and look at how
the new Lambdify's behave, and it seems that the ordinary Lambdify function
now works as expected.... Thanks!

Unfortunately i found that LambdifyCSE is still not working correctly.

...this can already be seen using the initial examples:
(in #188 isuruf mentioned that fixing LambdifyCSE is not a problem once symengines
cse is ready. does this mean that it will anyway be depreciated soon?)

import sympy as sp
from symengine import Lambdify, LambdifyCSE, var

x = var('x')
funcs = [sp.cos(x)**i for i in range(5)]

#ordinary sympy
f_sympy = sp.lambdify((x), funcs, modules = 'numpy')

#ordinary symengine
f_seng = Lambdify([x], funcs, order='F')

#symengine with cse
f_seng_cse = LambdifyCSE([x], funcs, order='F')

a = np.array([1, 2, 3])

s1 = f_sympy(a)
# OK
s2 = f_seng(a)
# OK
s3 = f_seng_cse(a)
# yields:
# *ValueError: all the input array dimensions except for the concatenation axis must match exactly*

[isuruf]

LambdifyCSE and Lambdify are merged now, so this issue is gone now.