当函数包含条件时,使用 Numpy 将函数应用于数组

Applying a function to an array using Numpy when the function contains a condition

当函数包含条件时,我很难将函数应用于数组。我有一个低效的解决方法,正在寻找一种有效(快速)的方法。在一个简单的例子中:

pts = np.linspace(0,1,11)
def fun(x, y):
    if x > y:
        return 0
    else:
        return 1

现在,如果我 运行:

result = fun(pts, pts)

然后我得到错误

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

if x > y 行提出。我的低效解决方法给出了正确的结果但速度太慢:

result = np.full([len(pts)]*2, np.nan)
for i in range(len(pts)):
    for j in range(len(pts)):
        result[i,j] = fun(pts[i], pts[j])

以更好(更重要的是,更快)的方式获得此信息的最佳方法是什么?

当函数包含条件时,我很难将函数应用于数组。我有一个低效的解决方法,正在寻找一种有效(快速)的方法。在一个简单的例子中:

pts = np.linspace(0,1,11)
def fun(x, y):
    if x > y:
        return 0
    else:
        return 1

现在,如果我 运行:

result = fun(pts, pts)

然后我得到错误

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

if x > y 行提出。我的低效解决方法给出了正确的结果但速度太慢:

result = np.full([len(pts)]*2, np.nan)
for i in range(len(pts)):
    for j in range(len(pts)):
        result[i,j] = fun(pts[i], pts[j])

以更好(更重要的是,更快)的方式获得此信息的最佳方法是什么?

编辑:使用

def fun(x, y):
    if x > y:
        return 0
    else:
        return 1
x = np.array(range(10))
y = np.array(range(10))
xv,yv = np.meshgrid(x,y)
result = fun(xv, yv)  

仍然提出相同的 ValueError

错误非常明显 - 假设您有

x = np.array([1,2])
y = np.array([2,1])

这样

(x>y) == np.array([0,1])

您的 if np.array([0,1]) 语句的结果应该是什么?是真的还是假的? numpy 告诉你这是模棱两可的。使用

(x>y).all()

(x>y).any()

是明确的,因此 numpy 为您提供解决方案 - 任何单元格对满足条件,或所有单元格对 - 都是明确的真值。您必须自己定义 向量 x 大于向量 y.

的确切含义

xy的所有对进行操作的numpy解决方案使得x[i]>y[j]是使用网格生成所有对:

>>> import numpy as np
>>> x=np.array(range(10))
>>> y=np.array(range(10))
>>> xv,yv=np.meshgrid(x,y)
>>> xv[xv>yv]
array([1, 2, 3, 4, 5, 6, 7, 8, 9, 2, 3, 4, 5, 6, 7, 8, 9, 3, 4, 5, 6, 7, 8,
       9, 4, 5, 6, 7, 8, 9, 5, 6, 7, 8, 9, 6, 7, 8, 9, 7, 8, 9, 8, 9, 9])
>>> yv[xv>yv]
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
       2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 7, 8])

要么将 xvyv 发送到 fun,要么在函数中创建网格,具体取决于哪个更有意义。这会生成所有 xi,yjxi>yj。如果您只需要 return xv>yv 的实际索引,其中每个单元格 ij 对应 x[i]y[j]。在你的情况下:

def fun(x, y):
    xv,yv=np.meshgrid(x,y)
    return xv>yv

将 return 一个矩阵,其中 fun(x,y)[i][j] 如果 x[i]>y[j] 为真,否则为假。或者

return  np.where(xv>yv)

将 return 两个索引对数组的元组,这样

for i,j in fun(x,y):

也会保证x[i]>y[j]

In [253]: x = np.random.randint(0,10,5)
In [254]: y = np.random.randint(0,10,5)
In [255]: x
Out[255]: array([3, 2, 2, 2, 5])
In [256]: y
Out[256]: array([2, 6, 7, 6, 5])
In [257]: x>y
Out[257]: array([ True, False, False, False, False])
In [258]: np.where(x>y,0,1)
Out[258]: array([0, 1, 1, 1, 1])

要对这两个一维数组进行笛卡尔比较,请重塑其中一个,使其可以使用 broadcasting:

In [259]: x[:,None]>y
Out[259]: 
array([[ True, False, False, False, False],
       [False, False, False, False, False],
       [False, False, False, False, False],
       [False, False, False, False, False],
       [ True, False, False, False, False]])
In [260]: np.where(x[:,None]>y,0,1)
Out[260]: 
array([[0, 1, 1, 1, 1],
       [1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1],
       [0, 1, 1, 1, 1]])

您的 if 函数仅适用于标量输入。如果给定数组,a>b 会生成一个布尔数组,不能在 if 语句中使用。您的迭代有效,因为它传递了标量值。对于一些复杂的函数,这是你能做的最好的(np.vectorize 可以使迭代更简单,但不会更快)。

我的答案是查看数组比较,并从中得出答案。在这种情况下,参数 3 where 可以很好地将布尔数组映射到所需的 1/0。还有其他方法可以进行此映射。

您的双循环需要额外的编码层,广播 None

对于更复杂的示例,或者如果您正在处理的数组有点大,或者如果您可以写入已经预先分配的数组,您可以考虑 Numba.

例子

import numba as nb
import numpy as np

@nb.njit()
def fun(x, y):
  if x > y:
    return 0
  else:
    return 1

@nb.njit(parallel=False)
#@nb.njit(parallel=True)
def loop(x,y):
  result=np.empty((x.shape[0],y.shape[0]),dtype=np.int32)
  for i in nb.prange(x.shape[0]):
    for j in range(y.shape[0]):
      result[i,j] = fun(x[i], y[j])
  return result

@nb.njit(parallel=False)
def loop_preallocated(x,y,result):
  for i in nb.prange(x.shape[0]):
    for j in range(y.shape[0]):
      result[i,j] = fun(x[i], y[j])
  return result

时间

x = np.array(range(1000))
y = np.array(range(1000))

#Compilation overhead of the first call is neglected

res=np.where(x[:,None]>y,0,1) -> 2.46ms
loop(single_threaded)         -> 1.23ms
loop(parallel)                -> 1.0ms
loop(single_threaded)*        -> 0.27ms
loop(parallel)*               -> 0.058ms

*可能受缓存影响。测试你自己的例子。