Convert sympy expressions to function of numpy arrays
<p>I have a system of ODEs written in sympy:</p>

<pre><code>from sympy.parsing.sympy_parser import parse_expr

xs = symbols('x1 x2')
ks = symbols('k1 k2')
strs = ['-k1 * x1**2 + k2 * x2', 'k1 * x1**2 - k2 * x2']
syms = [parse_expr(item) for item in strs]
</code></pre>

<p>I would like to convert this into a vector valued function, accepting a 1D numpy array of the x value, a 1D numpy array of the k values, returning a 1D numpy array of the equations evaluated at those points. The signature would look something like this:</p>

<pre><code>import numpy as np
x = np.array([3.5, 1.5])
k = np.array([4, 2])
xdot = my_odes(x, k)
</code></pre>

<p>The reason I want something like this is to give this function to <code>scipy.integrate.odeint</code>, so it needs to be fast.</p>

<p><strong>Attempt 1: subs</strong></p>

<p>Of course, I can write a wrapper around <code>subs</code>:</p>

<pre><code>def my_odes(x, k):
    all_dict = dict(zip(xs, x))
    all_dict.update(dict(zip(ks, k)))
    return np.array([sym.subs(all_dict) for sym in syms])
</code></pre>

<p>But this is super slow, especially for my real system which is much bigger and is run many times. I need to compile this operation to C code.</p>

<p><strong>Attempt 2: theano</strong></p>

<p>I can get close with <a href="http://docs.sympy.org/latest/modules/numeric-computation.html" rel="noreferrer">sympy's integration with theano</a>:</p>

<pre><code>from sympy.printing.theanocode import theano_function

f = theano_function(xs + ks, syms)

def my_odes(x, k):
    return np.array(f(*np.concatenate([x,k]))))
</code></pre>

<p>This compiles each expression, but all this packing and unpacking of the inputs and outputs slows it back down. The function returned by <code>theano_function</code> accepts numpy arrays as arguments, but it needs one array for each symbol rather than one element for each symbol. This is the same behavior for <code>functify</code> and <code>ufunctify</code> as well. I don't need the broadcast behavior; I need it to interpret each element of the array as a different symbol.</p>

<p><strong>Attempt 3: DeferredVector</strong></p>

<p>If I use <code>DeferredVector</code> I can make a function that accepts numpy arrays, but I can't compile it to C code or return a numpy array without packaging it myself.</p>

<pre><code>import numpy as np
import sympy as sp
from sympy import DeferredVector

x = sp.DeferredVector('x')
k =  sp.DeferredVector('k')
deferred_syms = [s.subs({'x1':x[0], 'x2':x[1], 'k1':k[0], 'k2':k[1]}) for s in syms]
f = [lambdify([x,k], s) for s in deferred_syms]

def my_odes(x, k):
    return np.array([f_i(x, k) for f_i in f])
</code></pre>

<p>Using <code>DeferredVector</code> I do not need to unpack the inputs, but I still need to pack the outputs. Also, I can use <code>lambdify</code>, but the <code>ufuncify</code> and <code>theano_function</code> versions perish, so no fast C code is generated.</p>

<pre><code>from sympy.utilities.autowrap import ufuncify
f = [ufuncify([x,k], s) for s in deferred_syms] # error

from sympy.printing.theanocode import theano_function
f = theano_function([x,k], deferred_syms) # error
</code></pre>