ilmorez ilmorez - 5 months ago 23
Python Question

SymPy cannot lambdify Product

I'm using SymPy 1.0 and Python 2.7. I want to compute the sum of first 100 integer numbers:

This code runs succesfully

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Sum(x[i], (i, 0, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s, 'numpy')
s_lambda(np.arange(101))


And gives
5050
as expected. But when I try to do the same with a
Product
instead of a
Sum
:

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 0, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s, 'numpy')
s_lambda(np.arange(101))


I got a
NameError: global name 'Product' is not defined

What am I doing wrong? Is there a workaround to get what I want?

Edit 1:
And what if I don't know in advance the limit of the
Product
. Let's say something like

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
n = sy.symbols('n', integer=True, positive=True)
s = sy.Product(x[i], (i, 0, n))
s_lambda = sy.lambdify((sy.DeferredVector('x'), n) s.doit(), 'numpy')
s_lambda(np.arange(101), 5)


Edit 2:
I'm trying to find a workaround.
NameError: global name 'Product' is not defined
error is given because of this:


lambdastr((sy.DeferredVector('x'), n), p)


That gives:

lambda x,n: (Product(x[i], (i, 0, n)))


While for the
Sum
we got a working lambda function:

lambda x,n: ((builtins.sum(x[i] for i in range(0, n+1))))


At this point the problem revolves around the definition of the
Product
function. According to the manual I can inject via a
dict
my definition of a function


def my_prod(a, b):
# my implementation
pass

my_fun = {"Product" : my_prod}
f = sy.lambdify((sy.DeferredVector('x'), n), p, modules=['numpy', my_fun])
f([1,2,3,4,5], 2)


Problem is,
list indices must be integers, not Symbol
error is raised when I try to use the lambdified function
f
. I guess this is due to
i
that is a symbol while it is supposed to be an integer. I can't understand why it's not passed the actual
integer
before trying to call
my_prod
as it is for the
Sum
case.

Answer

When the limits of the Product are known in advance

You could work around the problem by calling .doit() to expand the Product into its component parts:

In [104]: s = sy.Product(x[i], (i, 1, 10)); s
Out[104]: Product(x[i], (i, 1, 10))

In [105]: s.doit()
Out[105]: x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9]*x[10]

For example,

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 1, 10))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(11)))

prints

3628800

However, if you use .doit() with sy.Product(x[i], (i, 1, 100)) then you'll get an arithmetic overflow since np.arange(101) has dtype int32 or int64 (depending on your OS) and the product 100!

In [109]: math.factorial(100)
Out[109]: 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

is too large be stored in either an int32 or int64 array value.

In [118]: np.iinfo('int64').max
Out[118]: 9223372036854775807

In [119]: np.iinfo('int64').max < math.factorial(100)
Out[119]: True

Thus,

s = sy.Product(x[i], (i, 1, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(101)))

raises a

RuntimeWarning: overflow encountered in long_scalars

and erroneously prints 0.


If you change the input from an array of dtype int64 to a list of Python ints then the product can be computed correctly:

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
s = sy.Product(x[i], (i, 1, 100))
s_lambda = sy.lambdify(sy.DeferredVector('x'), s.doit(), 'numpy')
print(s_lambda(np.arange(101).tolist()))

prints

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

When the limits of the Product are not known in advance

The work-around (AFAICS) becomes more complicated. If you use a debugger to trace the code path followed when Sum is used you'll find that LambdaPrinter._print_Sum is called to convert Sum(x[i], (i, 0, n)) to the expression builtins.sum(x[i] for i in range(0, n+1)).

If we add a _print_Product method to NumPyPrinter (a subclass of LambdaPrinter), then we can get lambdify to successfully convert Product into an expression that NumPy can evaluate:

import sympy as sy
from sympy.tensor import IndexedBase, Idx
import numpy as np
import sympy.printing.lambdarepr as SPL

def _print_Product(self, expr):
    loops = (
        'for {i} in range({a}, {b}+1)'.format(
            i=self._print(i),
            a=self._print(a),
            b=self._print(b))
        for i, a, b in expr.limits)
    return '(prod([{function} {loops}]))'.format(
        function=self._print(expr.function),
        loops=' '.join(loops))
SPL.NumPyPrinter._print_Product = _print_Product

x = sy.IndexedBase('x')
i = sy.symbols('i', cls=Idx)
n = sy.symbols('n', integer=True, positive=True)
s = sy.Product(x[i], (i, 1, n))
s_lambda = sy.lambdify((sy.DeferredVector('x'), n), s, 'numpy')
print(s_lambda(np.arange(101), 5))

prints

120
Comments