将 sympy 表达式转换为 numpy 数组的函数 [英] Convert sympy expressions to function of numpy arrays

问题描述

我有一个用 sympy 编写的 ODE 系统:

I have a system of ODEs written in 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]

我想将其转换为向量值函数，接受 x 值的一维 numpy 数组、k 值的一维 numpy 数组，返回在这些点计算的方程的一维 numpy 数组.签名看起来像这样:

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:

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

我想要这样的原因是为了把这个函数交给scipy.integrate.odeint，所以它需要很快.

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

尝试 1:订阅

当然，我可以为 subs 编写一个包装器:

Of course, I can write a wrapper around 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])

但这太慢了，特别是对于我的真实系统来说，它要大得多并且运行了很多次.我需要将此操作编译为 C 代码.

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.

尝试 2:theano

我可以接近 sympy 与 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]))))

这会编译每个表达式，但是所有这些输入和输出的打包和解包都会减慢它的速度.theano_function 返回的函数接受 numpy 数组作为参数，但它需要每个符号一个数组，而不是每个符号一个元素.这与 functify 和 ufunctify 的行为相同.我不需要广播行为；我需要它将数组的每个元素解释为不同的符号.

This compiles each expression, but all this packing and unpacking of the inputs and outputs slows it back down. The function returned by theano_function 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 functify and ufunctify as well. I don't need the broadcast behavior; I need it to interpret each element of the array as a different symbol.

尝试 3:DeferredVector

如果我使用 DeferredVector 我可以制作一个接受 numpy 数组的函数，但是我不能将它编译成 C 代码或返回一个 numpy 数组而不自己打包.

If I use DeferredVector 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.

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

使用 DeferredVector 我不需要解包输入，但我仍然需要打包输出.另外，我可以使用 lambdify，但是 ufuncify 和 theano_function 版本已经消失，所以没有生成快速的 C 代码.

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

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

推荐答案

你可以使用sympy函数<代码>lambdify.例如，

You can use the sympy function lambdify. For example,

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 在脚本运行时打印.

要快得多(注意计时结果单位的变化):

It is much faster (note the change in units of the timing results):

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

这篇关于将 sympy 表达式转换为 numpy 数组的函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！