Sympy autowrap (cython)：为sympy.Max、sympy.Heaviside定义“助手”

Question

我有一个名为J_sym的sympy.Matrix，我希望自动捕获它(最好使用cython后端)；相应的符号存储在列表list_args中。

然而，我遇到的问题是，显然不支持某些症状函数，在我的例子中特别是sympy.Max和sympy.Heaviside。

具体来说，如果我尝试

J_num_autowrap = autowrap(J_sym, backend="cython", args=list_args)

我得到了以下信息：

[...]
wrapped_code_10.c(4): warning C4013: 'Heaviside' undefined; assuming extern returning int
wrapped_code_10.c(4): warning C4013: 'Max' undefined; assuming extern returning int
[...]
: fatal error LNK1120: 2 unresolved externals
[...]
failed with exit status 1120

(当只使用numpy进行lambdifies时，也会出现类似的问题；但是可以很容易地传递这些函数的字典。)

现在，似乎autowrap的"helpers“参数应该提供一个解决方案；然而，我不知道如何正确地使用它。文档字符串本身绝对是隐秘的，我只找到了这个链接和一个例子：

https://github.com/sympy/sympy/issues/10572

我不确定如何在这里正确使用助手--有人能帮上忙吗？

如果我尝试这样的方法：

J_num_autowrap = autowrap(J_sym, backend="cython", args=list_args, helpers=("Max", max(x,y), [x,y]))

如果x，y没有声明，我会得到(a)一个错误--我如何用自由变量来做这件事？以及(b)如果我在将x和y声明为渐近符号后执行相同的操作，则只会得到错误消息“关系的真值无法确定”。

(上面的链接提到，无论如何都不可能传递多个帮助器-尽管在我看来，这似乎应该在1.0中修复。当我尝试传递2个(乱七八糟的)帮助器元组时，我仍然得到"ValueError:没有足够的值来解包(期望3，得到2)“。不确定状态-也许我在这里只是运气不佳)。

如果有任何从我的矩阵生成C代码的替代方法，我也希望在这个方向上有任何提示。

Answer 1

作为Max问题的一种解决方法，您可以使用

helpers=("Max", (abs(x+y) + abs(x-y))/2, [x, y])

表达式(abs(x+y) + abs(x-y))/2在数学上等同于max(x, y)，并且不会从autowrap中产生任何抱怨。一个完整的例子：

from sympy import *
from sympy.utilities.autowrap import autowrap
x, y = symbols('x y')
f = autowrap(Max(2*x, y+1), args=[x, y], backend="cython", helpers=("Max", (abs(x+y) + abs(x-y))/2, [x, y]))
print([f(5, 6), f(1, 2)])    # outputs 10 and 3 

类似地，Heaviside(x)与(x + abs(x))/(2*x)相同--除了后一个表达式在x=0时为NaN。更丑陋但安全的版本是(x + abs(x))/(2*abs(x) + 1e-300)，其中添加的1e-300几乎永远不会改变结果。

多个帮助器的问题

Max和Heaviside都需要的帮助器。这就是你遇到open issue的地方:在自动抓取中有一个错误，这使得它不可能使用多个帮助器。Documentation建议格式如下所示

helpers=[("Max", ..., [x, y]), ("Heaviside", ..., [x])]

但是on this line自动抓取对于一个帮助器来说是方便的(我们不需要把它放在列表中)，但是对于多个帮助器来说却是致命的(额外的包装层)：

helpers = [helpers] if helpers else ()

当然，随后的for name_h, expr_h, args_h in helpers解包会失败。

多个辅助对象的解决方法

ufuncify正确处理helpers参数。不幸的是，ufuncify最终会为除NumPy之外的所有后端调用autowrap，因此该错误仍然发生在自动抓取中。但如果你愿意使用NumPy后端，这是一个解决方案：

from sympy.utilities.autowrap import ufuncify
x, y = symbols('x y', real=True)
my_helpers = [("Max", abs(x+y)/2 + abs(x-y)/2, [x, y]), ("Heaviside", (x + abs(x)) / (2*abs(x) + 1e-300), [x])]
f = ufuncify([x,y], Max(2*x, y+1) + Heaviside(x-4), backend="numpy", helpers=my_helpers)
print([f(5, 6), f(1, 2)])   $ outputs 11  and  3

减少到一个帮手

如果您可以更改表达式以消除Max和Heaviside中的一个(甚至两者)，则可以使用Cython后端。例如，您可能只需要Max(x，0)，因此您定义了一个Python函数"positive part"：

pos = lambda x: (x+abs(x))/2

那么pos的自动抓取就没有问题了，你只需要帮助Heaviside就行了：

f = autowrap(pos(9-x-y) + Heaviside(x-4), args = [x, y], backend = "cython", helpers=("Heaviside", (x + abs(x)) / (2*abs(x) + 1e-300), [x]))
print([f(5, 6), f(1, 2)])   #  outputs 1 and 6

这种替换的实用性取决于符号表达式的获取方式。