python - 如何确定是否需要 numpy.vectorize()？

Question

有没有办法在运行时确定一个函数是否需要 numpy.vectorize() 才能按预期运行？

作为背景，我问这个是因为我在一个程序中使用 Numpy 从文献中可用的热力学函数(基于 CALPHAD)计算相图。对于给定的温度，人们评估自由能函数并确定接触上凹(二阶导数 > 0)的公共(public)切线曲线以定义相共存的组成范围。为此，最好直接定义二阶导数函数。在我尝试使用具有常数二阶导数的简单抛物线自由能函数进行测试之前，使用真正的自由能函数(不难获得其导数)一切都进行得很顺利。这破坏了我的算法，因为我没想到 numpy 广播会查看函数内部并决定它不需要广播。

困难归结为这种行为:

import numpy as np
def f(x):
   return(x*x)
def g(x):
   return(3.0)
def h(x):
   return(0*x+3.0)
def i(x):
   return(x-x+3.0)

x = np.linspace(1.0, 5.0, 5)

在 IPython 3.3.2 中运行会产生以下输出:

f(x) -> array([ 1., 4., 9., 16., 25.]) -- what I expected

g(x) -> 3.0 (note only 1 element, and a float, not ndarray) -- not naively expected

h(x) -> array([ 3., 3., 3., 3., 3.]) -- OK, fooled the broadcasting by having x do something

i(x) -> array([ 3., 3., 3., 3., 3.]) -- Same as h(x) but avoiding the multiply but with roundoff issues

现在我可以使用了

gv = np.vectorize(g)

得到

gv(x) -> array([ 3., 3., 3., 3., 3.]) -- expected behavior

如果我的程序(最终)接受任意用户输入的自由能函数，这将导致问题，除非所有用户都了解 numpy 内部广播魔法。或者，我可以反身 np.vectorize 一切来防止这种情况。问题在于函数在 numpy 中“正常工作”的成本。

即在 IPython 中使用 %timeit，

h(x) -> 100000 loops, best of 3: 3.45 µs per loop

如果我不必要地向量化 h(x)(即 hv = np.vectorize(h))，我会得到

hv(x) -> 10000 loops, best of 3: 43.2 µs per loop

因此，不必要的矢量化是一个巨大的惩罚(5 次函数评估需要 40 微秒)。

我想我可以对在小型 ndarray 上评估的函数的返回进行初步测试，看看返回类型是数组还是 float ，然后定义一个新函数(如果它是 float )，例如:

def gv(x):
   return(g(x)+0.0*x)

这看起来像是一场可怕的混战。

那么 - 在这种情况下，是否有更好的方法来“欺骗”numpy 使其有效广播？

Answer 1

解决所示问题。如果你想要一个新数组:

def g(x):
    return np.ones_like(x)*3

或者如果你想将数组中的所有元素就地设置为 3:

def g(x):
    x[:] = 3

请注意，这里没有 return 语句，因为您只是更新数组 x 以便所有元素均为 3。

如图所示，def g(x): return(3) 的问题是函数内部没有对 numpy 的引用。您声明 any 给定输入返回 3。声明 x=3 将遇到类似的问题，因为您正在更新指针 x 以指向 3 而不是 numpy 数组。虽然语句 x[:]=3 从 numpy.ndarray 类访问称为 View 的内部函数，而不是通常使用 = 简单更新指针的语句。