NumPy故障矢量化方法

Question

我有三个数组，a、b和c。形状分别为(N, 2)、(N, 3)、(N, 3)。

我需要比较b中的每行元素，并在a上更新同一行中的索引。我以为我知道如何向量化这个方法，但我认为我的维度是不正确的。

我目前拥有的是：

def to_cube(points):
    cube = np.empty((len(points), 3), dtype=np.half)
    delta = np.empty_like(cube)

    q = ((2 / 3) * points[:, 0]) / 0.1
    r = (((-1 / 3) * points[:, 0]) + ((np.sqrt(3) / 3) * points[:, 1])) / 0.1

    cube[:, 0] = np.round(q)
    cube[:, 1] = np.round(-q-r)
    cube[:, 2] = np.round(r)
    delta[:, 0] = np.abs(cube[:, 0] - q)
    delta[:, 1] = np.abs(cube[:, 1] - (-q-r))
    delta[:, 2] = np.abs(cube[:, 2] - r)

    if delta[:, 0] > delta[:, 1] and delta[:, 1] > delta[:, 2]:
        cube[:, 0] = -cube[:, 1] - cube[:, 2]
    elif delta[:, 1] > delta[:, 2]:
        cube[:, 1] = -cube[:, 0] - cube[:, 2]
    else:
        cube[:, 2] = -cube[:, 0] - cube[:, 1]

    return cube

这会抛出一个ValueError: The truth value of an array with more than one element is ambiguous。

查看条件后，可以清楚地看到，对delta[:, 0] > delta[:, 1]的第一次检查将返回形状(N, 1)的数组。如何将其更改为a中的每一行，获取该行上的适当索引，然后根据条件更新b中的同一行的特定索引？

编辑:示例

此示例假定q和r已经完成。这些矩阵表示cube和delta

>>> cube
array([[275.0, -400.0, 124.0]], dtype=float16) # so this is a (1, 3) but could be (N, 3)
>>> cube[0]
array([275.0, -400.0, 124.0], dtype=float16)
>>> delta
array([[5., 10., 3.]], dtype=float16)
>>> delta[0]
array([5., 10., 3.], dtype=float16)

现在，通过条件执行(值在以下内容中)：

if 5.0 > 10.0 and 10.0 > 3.0:
    cube[0] = -(-400.0) - 124.0
elif 10.0 > 3.0:
    cube[1] = -(275.0) - 124.0
else:
    cube[2] = -(275.0) - (-400.0)

return cube # array([275.0, -(275.0) - 124.0, 124.0], dtype=float16)

这显示了每行发生了什么，现在我需要对所有行都这样做。

编辑:潜在的解决方案(它是矢量化的吗？)

有一种方法可以确保使用for-range独立地访问行。

def to_cube(points):
    cube = np.empty((len(points), 3), dtype=np.half)
    delta = np.empty_like(cube)

    q = ((2 / 3) * points[:, 0]) / 0.1
    r = (((-1 / 3) * points[:, 0]) + ((np.sqrt(3) / 3) * points[:, 1])) / 0.1

    cube[:, 0] = np.round(q)
    cube[:, 1] = np.round(-q-r)
    cube[:, 2] = np.round(r)
    delta[:, 0] = np.abs(cube[:, 0] - q)
    delta[:, 1] = np.abs(cube[:, 1] - (-q-r))
    delta[:, 2] = np.abs(cube[:, 2] - r)

    for i in range(len(cube)):
        if delta[i, 0] > delta[i, 1] and delta[i, 1] > delta[i, 2]:
            cube[i, 0] = -cube[i, 1] - cube[i, 2]
        elif delta[i, 1] > delta[i, 2]:
            cube[i, 1] = -cube[i, 0] - cube[i, 2]
        else:
            cube[i, 2] = -cube[i, 0] - cube[i, 1]

    return cube

然而，我现在“循环”数组，似乎没有矢量化或广播。

Answer 1

对任何感兴趣的人来说，这就是我解决问题的方法：

def to_cube(points):
    cube = np.empty((len(points), 3), dtype=np.half)
    delta = np.empty_like(cube)

    q = ((2 / 3) * points[:, 0]) / 0.1
    r = (((-1 / 3) * points[:, 0]) + ((np.sqrt(3) / 3) * points[:, 1])) / 0.1

    cube[:, 0] = np.round(q)
    cube[:, 1] = np.round(-q-r)
    cube[:, 2] = np.round(r)
    delta[:, 0] = np.abs(cube[:, 0] - q)
    delta[:, 1] = np.abs(cube[:, 1] - (-q-r))
    delta[:, 2] = np.abs(cube[:, 2] - r)
    
    # define boolean arrays for where conditions exist
    rxc = ((delta[:, 0] > delta[:, 1]) & (delta[:, 1] > delta[:, 2]))
    ryc = (delta[:, 1] > delta[:, 2])
    rzc = ~(rxc + ryc)

    # update just those indices by condition
    cube[rxc, 0] = -cube[rxc, 1] - cube[rxc, 2]
    cube[ryc, 1] = -cube[ryc, 0] - cube[ryc, 2]
    cube[rzc, 2] = -cube[rzc, 0] - cube[rzc, 1]

    return cube

如果有人看到了改进的空间，我很想知道！

我的系统上的基准：

import numpy as np
from timeit import timeit

u = np.random.uniform

points = np.array([[u(0, 50), u(0, 50)] for _ in range(37000000)], dtype=np.half)

p = 'from __main__ import points, to_cube; to_cube(points)'

timeit(p, number=1)

# output: 17.94858811999