为函数中变量的固定值优化Python数学代码

Question

我有一个很长的数学公式（只是为了把你放在上下文中：它有293095个字符），实际上它是一个python函数的主体。此函数有15个输入参数，如下所示：

def math_func(t,X,P,n1,n2,R,r):
    x,y,z = X
    a,b,c = P
    u1,v1,w1 = n1
    u2,v2,w2 = n2
    return <long math formula>

公式使用简单的数学运算+ - * ** /和一个对arctan的函数调用。以下是它的摘录：

^{pr2}$

现在的重点是，在实践中，这个函数的批量求值将针对P,n1,n2,R和r的固定值进行，这将把自由变量集减少到只有4个，并且“理论上”参数较少的公式应该更快。在

所以问题是：如何在Python中实现这种优化？

我知道我可以把所有的东西放在一个字符串中，然后做一些replace，compile和{}之类的

formula = formula.replace('r','1').replace('R','2')....
code = compile(formula,'formula-name','eval')
math_func = lambda t,x,y,z: eval(code)

如果某些操作（比如power）被它们的值所替代，那就更好了，例如18*r**6*t*u1**2*u2**2*v1**2应该变成{}而不是{}。我认为compile应该这样做，但无论如何我不确定。是否compile实际执行此优化？

我的解决方案加快了计算速度，但如果我能把它压缩得更多，那就更好了。注意：最好使用标准Python（稍后我可以试试Cython）。在

一般来说，我对pythonic实现我的目标的方法很感兴趣，也许需要一些额外的库：什么是一种比较好的方法？我的解决方案是一个好的方法吗？在

编辑：（提供更多上下文）

这个巨大的表达式是在一个圆弧上的符号线积分的输出。在空间中，弧由半径r、两个正交法向量（如二维版本中的x和y轴）n1=(u1,v1,w1)、n2=(u2,v2,w2)和中心P=(a,b,c)。剩下的部分是我执行积分X=(x,y,z)和我正在集成的函数的参数R。在

Sympy和Maple只需花很长时间来计算，实际输出来自Mathematica。在

如果您对这里的公式很好奇，它是（伪伪代码）：

G(u) = P + r*(1-u**2)/(1+u**2)*n1 + r*2*u/(1+u**2)*n2
integral of (1-|X-G(t)|^2/R^2)^3 over t

Answer 1

你可以使用Sympy：

>>> from sympy import symbols
>>> x,y,z,a,b,c,u1,v1,w1,u2,v2,w2,t,r = symbols("x,y,z,a,b,c,u1,v1,w1,u2,v2,w2,t,r")
>>> r=u1=u2=v1=1
>>> a = 18*r**6*t*u1**2*u2**2*v1**2
>>> a
18*t

然后可以创建如下Python函数：

^{pr2}$

这个f函数实际上就是18*t：

>>> import dis
>>> dis.dis(f)
  1           0 LOAD_CONST               1 (18)
              3 LOAD_FAST                0 (_Dummy_18)
              6 BINARY_MULTIPLY
              7 RETURN_VALUE

如果要将生成的代码编译成机器代码，可以尝试使用JIT编译器，例如Numba、Theano或{a3}。在

Answer 2

我将如何解决这个问题：

1. compile()将函数转换为AST（抽象语法树），而不是普通的字节码函数-有关详细信息，请参阅标准ast模块。在
2. 遍历AST，将对固定参数的所有引用替换为其固定值。有一些库，比如macropy，可能对这个有用，我没有任何具体的建议。在
3. 再次遍历AST，执行任何可能启用的优化，例如Mult（1，X）=>；X。您不必担心两个常量之间的操作，因为Python（从2.6开始）已经对此进行了优化。在
4. compile()将AST转换为正常函数。就这么说吧，希望速度增加了足够的数量来证明所有的预优化是正确的。在

请注意，Python永远不会单独优化像1*X这样的东西，因为它无法知道运行时X是什么类型的-它可能是以任意方式实现乘法运算的类的实例，因此结果不一定是X。只有你知道所有变量都是普通数，遵守算法，使此优化有效。在

Answer 3

解决这样问题的“正确方法”是以下一种或多种：

1. 找到更有效的配方
2. 象征性地简化和减少术语
3. 使用矢量化（例如NumPy）
4. 将双关键转到已经优化的低级库（例如，在隐式执行强表达式优化的C或Fortran语言中，而不是Python，后者确实nada）。

不过，让我们暂时假设，方法1、3和4不可用，您必须使用Python来实现这一点。然后简化和“提升”公共子表达式是您的主要工具。

好消息是，机会有很多。例如，表达式^{cd1>}重复26次。您可以只分配^{{cd2>}一次，然后每次替换^{cd1>}即可节省25个计算。

当您开始在这里查找公共表达式时，您会发现它们在任何地方都。把这个过程机械化会很好的，对吧？通常，这需要一个完整的表达式解析器（例如来自^{cd4>}模块），这是一个指数时间优化问题。但你的表情有点特殊。虽然长而多变，但并不特别复杂。它没有什么内部的假肢分组，所以我们可以用更快更脏的方法来摆脱。

在“如何”之前，结果代码为：

sa = r**6                      # 26 occurrences
sb = u1**2                     # 5 occurrences
sc = u2**2                     # 5 occurrences
sd = v1**2                     # 5 occurrences
se = u1**4                     # 4 occurrences
sf = u2**3                     # 3 occurrences
sg = u1**3                     # 3 occurrences
sh = v1**4                     # 3 occurrences
si = u2**4                     # 3 occurrences
sj = v1**3                     # 3 occurrences
sk = v2**2                     # 1 occurrence
sl = v1**6                     # 1 occurrence
sm = v1**5                     # 1 occurrence
sn = u1**6                     # 1 occurrence
so = u1**5                     # 1 occurrence
sp = u2**6                     # 1 occurrence
sq = u2**5                     # 1 occurrence
sr = 6*sa                      # 6 occurrences
ss = 3*sa                      # 5 occurrences
st = ss*t                      # 5 occurrences
su = 12*sa                     # 4 occurrences
sv = sa*t                      # 3 occurrences
sw = v1*v2                     # 5 occurrences
sx = sj*v2                     # 3 occurrences
sy = 24*sv                     # 3 occurrences
sz = 15*sv                     # 2 occurrences
sA = sr*u1                     # 2 occurrences
sB = sy*u1                     # 2 occurrences
sC = sb*sc                     # 2 occurrences
sD = st*se                     # 2 occurrences

# revised formula
sv*sn - sr*so*u2 - sz*se*sc +
20*sa*sg*sf + sz*sb*si - sA*sq -
sv*sp + sD*sd - su*sg*u2*sd -
18*sv*sC*sd + su*u1*sf*sd +
st*si*sd + st*sb*sh - sA*u2*sh -
st*sc*sh + sv*sl - sr*se*sw -
sy*sg*u2*sw + 36*sa*sC*sw +
sB*sf*sw - sr*si*sw -
su*sb*sx - sB*u2*sx +
su*sc*sx - sr*sm*v2 - sD*sk

这样就避免了81次计算。只是个粗略的伤口。甚至可以进一步改善结果。例如，可以预先计算子表达式^{cd5}和^{cd6>}。但我们会把下一层留一天。

注意，这不包括起始^{cd7>}。大多数简化可以（而且需要）在表达的核心上进行，而不是从外部术语上。所以我现在把它们去掉了，一旦计算出公共核心，它们就可以被添加回去。

你怎么做到的？

^{pr2}$

使用正则表达式解析是一项棘手的工作。那次旅行容易出错、悲伤和遗憾。我通过提升一些不需要严格的指数来防止不良结果，并将随机值插入公式的前后，以确保它们都给出相同的结果。如果这是生产代码，我建议使用“双关到C”策略。但如果你不能。。。