我有一套用交响乐写成的颂歌:
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值的1D numpy数组,k值的1D numpy数组,返回在这些点上计算的方程的1D numy数组。签名应该是这样的:
^{pr2}$我想要这样做的原因是把这个函数赋给scipy.integrate.odeint
,所以它需要很快。在
尝试1:subs
当然,我可以为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代码。在
尝试2:t否
我可以接近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]))))
这将编译每个表达式,但是所有这些输入和输出的打包和解包会使它慢下来。theano_function
返回的函数接受numpy数组作为参数,但它需要为每个符号使用一个数组,而不是为每个符号指定一个元素。这对functify
和ufunctify
也是相同的行为。我不需要广播行为;我需要它将数组的每个元素解释为不同的符号。在
尝试3:DeferredVector
如果我使用DeferredVector
我可以创建一个接受numpy数组的函数,但是我不能将它编译成C代码或返回一个numpy数组而不自己打包它。在
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
和{
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
您可以使用sypy函数^{} 。例如
运行脚本时打印
True
。在它要快得多(注意计时结果的单位变化):
^{pr2}$我写了a module named JiTCODE,它是为像你这样的问题量身定做的。 它接受符号表达式,将它们转换为C代码,在其周围包装一个Python扩展,对其进行编译,并加载它以与}一起使用。在
scipy.integrate.ode
或{您的示例如下所示:
然后您可以使用
ODE
,这与scipy.integrate.ode
的实例非常相似。在虽然您的应用程序不需要此函数,但您也可以提取并使用编译后的函数:
^{pr2}$注意,与您的规范不同,
k
不是作为NumPy数组传递的。对于大多数ODE应用程序,这不应该是相关的,因为硬编码控制参数更有效。在最后,请注意,对于这个小示例,由于}的开销,您可能无法获得最佳性能(另请参见SciPy Issue #8257或{a3})。
对于大型微分方程(如您所见),这种开销变得无关紧要。在
scipy.integrate.ode
或{相关问题 更多 >
编程相关推荐