Chuyển đổi sympy biểu đến chức năng của mảng NumPy

Question

Tôi có một hệ thống ODEs viết bằng sympy:

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] 

Tôi muốn chuyển đổi này vào một hàm vector có giá trị, chấp nhận một 1D NumPy mảng các giá trị x, một mảng có giá trị 1D của các giá trị k, trả về mảng numpy 1D của các phương trình được đánh giá tại các điểm đó. Chữ ký sẽ giống như thế này:

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

Lý do tôi muốn một cái gì đó như thế này là để cung cấp chức năng này để scipy.integrate.odeint, vì vậy nó cần phải được nhanh chóng.

Cố gắng 1: tàu ngầm

Tất nhiên, tôi có thể viết một wrapper xung quanh subs:

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

Nhưng điều này là siêu chậm, đặc biệt đối với hệ thống thật của tôi mà là lớn hơn nhiều và được điều hành nhiều lần. Tôi cần phải biên dịch hoạt động này sang mã C.

Cố gắng 2: theano

tôi có thể nhận được gần gũi với sympy's integration with theano:

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

này biên dịch mỗi biểu thức, nhưng tất cả bao bì này và giải nén của các đầu vào và đầu ra chậm nó xuống. Hàm trả về bởi theano_function chấp nhận các mảng numpy làm đối số, nhưng nó cần một mảng cho mỗi ký hiệu thay vì một phần tử cho mỗi ký hiệu. Đây cũng là hành vi tương tự đối với functify và ufunctify. Tôi không cần hành vi phát sóng; Tôi cần nó để giải thích từng phần tử của mảng như một biểu tượng khác.

Nỗ lực 3: DeferredVector

Nếu tôi sử dụng DeferredVector tôi có thể làm cho một chức năng chấp nhận mảng NumPy, nhưng tôi không thể biên dịch nó thành mã C hoặc trả lại một mảng NumPy mà không đóng gói nó bản thân mình.

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

Sử dụng DeferredVector Tôi không cần phải giải nén đầu vào, nhưng tôi vẫn cần đóng gói các đầu ra. Ngoài ra, tôi có thể sử dụng lambdify, nhưng các phiên bản ufuncify và theano_function bị hư hỏng, vì vậy không có mã C nhanh nào được tạo.

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

Bạn có thể sử dụng chức năng sympy lambdify. Ví dụ:

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 được in khi tập lệnh được chạy.

Nó là nhanh hơn nhiều (lưu ý sự thay đổi trong các đơn vị của kết quả thời gian):

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