Using Sympy's lambdify with IPython.parallel

For a numerical computation I am working on, I need to define a large Liouvillian matrix. Rather than program it in element by element, which is tedious and error-prone, I have used Sympy to construct it algebraically, and then use lambdify to make numpy matrices for use in numerical work.

Here I construct a dumb sympy matrix that doesn't mean anything:

In [1]:
import sympy as s
from sympy.abc import x,y
s.init_printing()

element = lambda n, m : m * x**n if (n+m) % 3 else y
L = s.Matrix([[element(n,m) for m in range(9)] for n in range(9)])
L
Out[1]:
$$\left[\begin{matrix}y & 1 & 2 & y & 4 & 5 & y & 7 & 8\\0 & x & y & 3 x & 4 x & y & 6 x & 7 x & y\\0 & y & 2 x^{2} & 3 x^{2} & y & 5 x^{2} & 6 x^{2} & y & 8 x^{2}\\y & x^{3} & 2 x^{3} & y & 4 x^{3} & 5 x^{3} & y & 7 x^{3} & 8 x^{3}\\0 & x^{4} & y & 3 x^{4} & 4 x^{4} & y & 6 x^{4} & 7 x^{4} & y\\0 & y & 2 x^{5} & 3 x^{5} & y & 5 x^{5} & 6 x^{5} & y & 8 x^{5}\\y & x^{6} & 2 x^{6} & y & 4 x^{6} & 5 x^{6} & y & 7 x^{6} & 8 x^{6}\\0 & x^{7} & y & 3 x^{7} & 4 x^{7} & y & 6 x^{7} & 7 x^{7} & y\\0 & y & 2 x^{8} & 3 x^{8} & y & 5 x^{8} & 6 x^{8} & y & 8 x^{8}\end{matrix}\right]$$

In the case of this example, I could have constructed a numpy matrix directly using the same nested loop, but this isn't true of the matrix in my actual problem. Anyway it is nice to see it written out in algebraic notation before plugging in numbers.

I use lambdify to get a Numpy matrix for my numerical work:

In [2]:
numer_L = s.lambdify((x,y), L, 'numpy')
numer_L(3,4)
Out[2]:
matrix([[    4,     1,     2,     4,     4,     5,     4,     7,     8],
        [    0,     3,     4,     9,    12,     4,    18,    21,     4],
        [    0,     4,    18,    27,     4,    45,    54,     4,    72],
        [    4,    27,    54,     4,   108,   135,     4,   189,   216],
        [    0,    81,     4,   243,   324,     4,   486,   567,     4],
        [    0,     4,   486,   729,     4,  1215,  1458,     4,  1944],
        [    4,   729,  1458,     4,  2916,  3645,     4,  5103,  5832],
        [    0,  2187,     4,  6561,  8748,     4, 13122, 15309,     4],
        [    0,     4, 13122, 19683,     4, 32805, 39366,     4, 52488]])

Say I want to do a computation involving this matrix (determinant, say) evaluated at multiple values of $y$:

In [3]:
# in series
import numpy
s_result = list(map(lambda y: numpy.linalg.det(numer_L(3,y)), range(30)))
s_result
Out[3]:
$$\begin{bmatrix}0.0, & 9.9284346199e-51, & 1.17528889921e-49, & 0.0, & 0.0, & 0.0, & 0.0, & 0.0, & 0.0, & 1.32486805867e-62, & -1.46681392565e-30, & 0.0, & 0.0, & -7.24249763767e-49, & 0.0, & 2.48513164374e-45, & 0.0, & 0.0, & 6.83960894448e-29, & 0.0, & 0.0, & 2.98238783422e-29, & -2.62584900178e-45, & 0.0, & 0.0, & 0.0, & 0.0, & 1.22565792285e-25, & 0.0, & 0.0\end{bmatrix}$$
In [4]:
# set up parallel environment. 2 engines started with `ipcluster start -n 2`
from IPython.parallel import Client
rc = Client()
dview = rc[:]
In [5]:
# in parallel
# do imports and push our lambda function over
dview.execute('import numpy')
dview.push(dict(numer_L=numer_L))
p_result = dview.map_sync(lambda y: numpy.linalg.det(numer_L(3,y)), range(30))
p_result
[0:apply]: 
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)<string> in <module>()
<ipython-input-5-1f431230550c> in <lambda>(y)
/Users/tkb/.virtualenvs/sympy/lib/python2.7/site-packages/numpy/__init__.pyc in <lambda>(x, y)
NameError: global name 'ImmutableMatrix' is not defined

[1:apply]: 
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)<string> in <module>()
<ipython-input-5-1f431230550c> in <lambda>(y)
/Users/tkb/.virtualenvs/sympy/lib/python2.7/site-packages/numpy/__init__.pyc in <lambda>(x, y)
NameError: global name 'ImmutableMatrix' is not defined

That didn't work because apparently the lambda function needs ImmutableMatrix defined, which I hadn't ever heard of, and is not even the type of the matrix we lambdified:

In [6]:
type(L)
Out[6]:
sympy.matrices.dense.MutableDenseMatrix

In any case, I don't want my engines running any Sympy code. The task I want to distribute is numerical, not algebraic, and hopefully lambdify has generated numpy code that can run on its own.

My current workaround

I inspected the lambdify source code and found a lambdastr function that ought to show me the code being generated:

In [7]:
from sympy.utilities.lambdify import lambdastr
lstr = lambdastr((x,y), L, dummify=True)
lstr
Out[7]:
'lambda x,y: (ImmutableMatrix([[y, 1, 2, y, 4, 5, y, 7, 8], [0, x, y, 3*x, 4*x, y, 6*x, 7*x, y], [0, y, 2*x**2, 3*x**2, y, 5*x**2, 6*x**2, y, 8*x**2], [y, x**3, 2*x**3, y, 4*x**3, 5*x**3, y, 7*x**3, 8*x**3], [0, x**4, y, 3*x**4, 4*x**4, y, 6*x**4, 7*x**4, y], [0, y, 2*x**5, 3*x**5, y, 5*x**5, 6*x**5, y, 8*x**5], [y, x**6, 2*x**6, y, 4*x**6, 5*x**6, y, 7*x**6, 8*x**6], [0, x**7, y, 3*x**7, 4*x**7, y, 6*x**7, 7*x**7, y], [0, y, 2*x**8, 3*x**8, y, 5*x**8, 6*x**8, y, 8*x**8]]))'

It looks like this ought to work if I just import numpy.matrix as ImmutableMatrix, but no dice:

In [8]:
# in parallel
# do imports and push our lambda function over
dview.execute('import numpy')
dview.execute('from numpy import matrix as ImmutableMatrix')
dview.push(dict(numer_L=numer_L))
p_result = dview.map_sync(lambda y: numpy.linalg.det(numer_L(3,y)), range(30))
p_result
[0:apply]: 
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)<string> in <module>()
<ipython-input-8-f4a5fb6c80d1> in <lambda>(y)
/Users/tkb/.virtualenvs/sympy/lib/python2.7/site-packages/numpy/__init__.pyc in <lambda>(x, y)
NameError: global name 'ImmutableMatrix' is not defined

[1:apply]: 
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)<string> in <module>()
<ipython-input-8-f4a5fb6c80d1> in <lambda>(y)
/Users/tkb/.virtualenvs/sympy/lib/python2.7/site-packages/numpy/__init__.pyc in <lambda>(x, y)
NameError: global name 'ImmutableMatrix' is not defined

I got it to work by pushing the generated code as a string, and doing the eval myself:

In [9]:
# in parallel
# do imports and push our lambda function over
dview.execute('import numpy')
dview.execute('from numpy import matrix as ImmutableMatrix')
dview.push(dict(lstr=lstr))
p_result = dview.map_sync(lambda y: numpy.linalg.det(eval(lstr)(3,y)), range(30))
p_result
Out[9]:
$$\begin{bmatrix}0.0, & 9.9284346199e-51, & 1.17528889921e-49, & 0.0, & 0.0, & 0.0, & 0.0, & 0.0, & 0.0, & 1.32486805867e-62, & -1.46681392565e-30, & 0.0, & 0.0, & -7.24249763767e-49, & 0.0, & 2.48513164374e-45, & 0.0, & 0.0, & 6.83960894448e-29, & 0.0, & 0.0, & 2.98238783422e-29, & -2.62584900178e-45, & 0.0, & 0.0, & 0.0, & 0.0, & 1.22565792285e-25, & 0.0, & 0.0\end{bmatrix}$$
In [10]:
p_result == s_result
Out[10]:
True

So, I kind of got things working, but it doesn't feel like the right way to do it, probably due to the way lambdify works. Is there a way to generate more parallelizable code from Sympy?

Versions

In [11]:
import sys
sys.version
Out[11]:
'2.7.3 (default, Nov 29 2012, 17:18:24) \n[GCC 4.2.1 Compatible Apple Clang 4.1 ((tags/Apple/clang-421.11.66))]'
In [12]:
import IPython
IPython.version_info
Out[12]:
(1, 1, 0, '')
In [13]:
s.__version__ # sympy
Out[13]:
'0.7.4.1'
In [14]:
numpy.__version__
Out[14]:
'1.8.0'