Convert sympy Ausdrücke Funktion von numpy Arrays

Question

Ich habe ein System von ODEs in sympy geschrieben:

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] 

Ich möchte dies in eine Vektor Wertfunktion umwandeln, ein 1D numpy Array des x-Wertes zu akzeptieren, ein 1D-numpy-Array der k-Werte, das ein 1D-numpy-Array der an diesen Punkten ausgewerteten Gleichungen liefert. Die Unterschrift würde wie folgt aussehen:

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

Der Grund, warum ich möchte so etwas wie dieses ist diese Funktion scipy.integrate.odeint geben, so muss es schnell sein.

Versuch 1: subs

Natürlich habe ich einen Wrapper um subs schreiben:

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]) 

Aber das ist super langsam, vor allem für mein reales System, das viel größer ist und ausgeführt wird, viele Male. Ich muss diese Operation zu C-Code kompilieren.

Versuch 2: Theano

ich in der Nähe mit sympy's integration with theano bekommen:

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])))) 

Dies kompiliert jeden Ausdruck, aber alle diese Verpackung und der Ein- und Ausgänge Auspacken verlangsamt es zurück. Die von theano_function zurückgegebene Funktion akzeptiert numpige Arrays als Argumente, benötigt jedoch für jedes Symbol ein Array anstelle eines Elements für jedes Symbol. Dies ist das gleiche Verhalten für functify und ufunctify. Ich brauche das Broadcast-Verhalten nicht; Ich brauche es, um jedes Element des Arrays als ein anderes Symbol zu interpretieren.

Versuch 3: DeferredVector

Wenn ich DeferredVector verwenden kann ich eine Funktion machen die numpy Arrays akzeptiert, aber ich kann es nicht in C-Code kompilieren oder eine numpy Array zurück, ohne es selbst zu verpacken.

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]) 

Mit DeferredVector Ich brauche nicht die Eingänge entpacken, aber ich muss noch die Ausgänge packen. Auch kann ich lambdify verwenden, aber die ufuncify und theano_function Versionen versterben, so dass kein schneller C-Code generiert wird.

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 

Answer 1

Sie können die Sympy-Funktion lambdify verwenden. Zum Beispiel

from sympy import symbols, lambdify 
from sympy.parsing.sympy_parser import parse_expr 
import numpy as np 

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] 

# Convert each expression in syms to a function with signature f(x1, x2, k1, k2): 
funcs = [lambdify(xs + ks, f) for f in syms] 


# This is not exactly the same as the `my_odes` in the question. 
# `t` is included so this can be used with `scipy.integrate.odeint`. 
# The value returned by `sym.subs` is wrapped in a call to `float` 
# to ensure that the function returns python floats and not sympy Floats. 
def my_odes(x, t, k): 
    all_dict = dict(zip(xs, x)) 
    all_dict.update(dict(zip(ks, k))) 
    return np.array([float(sym.subs(all_dict)) for sym in syms]) 

def lambdified_odes(x, t, k): 
    x1, x2 = x 
    k1, k2 = k 
    xdot = [f(x1, x2, k1, k2) for f in funcs] 
    return xdot 


if __name__ == "__main__": 
    from scipy.integrate import odeint 

    k1 = 0.5 
    k2 = 1.0 
    init = [1.0, 0.0] 
    t = np.linspace(0, 1, 6) 
    sola = odeint(lambdified_odes, init, t, args=((k1, k2),)) 
    solb = odeint(my_odes, init, t, args=((k1, k2),)) 
    print(np.allclose(sola, solb)) 

True gedruckt wird, wenn das Skript ausgeführt wird.

Es ist viel schneller (beachten Sie die Änderung in Einheiten der Zeitergebnisse):

In [79]: t = np.linspace(0, 10, 1001) 

In [80]: %timeit sol = odeint(my_odes, init, t, args=((k1, k2),)) 
1 loops, best of 3: 239 ms per loop 

In [81]: %timeit sol = odeint(lambdified_odes, init, t, args=((k1, k2),)) 
1000 loops, best of 3: 610 µs per loop