NEP 10 — 优化迭代器/UFunc 性能#

作者:

Mark Wiebe <mwwiebe@gmail.com>

内容类型:

text/x-rst

创建:

2010年11月25日

状态:

最终

目录#

摘要#

本 NEP 提出用一个新的迭代器替换 NumPy 迭代器和多迭代器，该迭代器设计得更灵活，并允许更缓存友好的数据访问。新迭代器还包含了大部分核心 ufunc 功能，从而可以轻松地在不完全符合 ufunc 模具的上下文中获得当前 ufunc 的优势。主要优势包括

  • 自动重新排序以查找缓存友好的访问模式

  • 标准和可自定义的广播

  • 自动类型/字节序/对齐转换

  • 可选缓冲以最大程度地减少转换内存使用

  • 可选输出数组，在未提供时自动分配

  • 自动输出或通用类型选择

此迭代器设计的大部分内容已实现，并取得了令人鼓舞的结果。构造开销略有增加（a.flat：0.5 us，nditer(a)：1.4 us 和 broadcast(a,b)：1.4 us，nditer([a,b])：2.2 us），但是，如示例所示，它已经可以改进纯 Python 代码与迭代器结合使用时内置 NumPy 机制的性能。一个示例重写了 np.add，对一些 Fortran 连续数组获得了四倍的改进，另一个示例将图像合成代码从 1.4 秒改进到 180 毫秒。

该实现试图考虑 NumPy 2.0 重构中做出的设计决策，使其将来集成到 libndarray 相对简单。

动机#

NumPy 默认从 UFunc 返回 C 连续数组。在处理结构不同的数据时，这可能导致极差的内存访问模式。一个简单的计时示例说明了这一点，从将 Fortran 连续数组加在一起产生了超过八倍的性能损失。所有计时都在 Athlon 64 X2 4200+ 上使用 NumPy 2.0dev（2010 年 11 月 22 日）和 64 位操作系统进行。

In [1]: import numpy as np
In [2]: a = np.arange(1000000,dtype=np.float32).reshape(10,10,10,10,10,10)
In [3]: b, c, d = a.copy(), a.copy(), a.copy()

In [4]: timeit a+b+c+d
10 loops, best of 3: 28.5 ms per loop

In [5]: timeit a.T+b.T+c.T+d.T
1 loops, best of 3: 237 ms per loop

In [6]: timeit a.T.ravel('A')+b.T.ravel('A')+c.T.ravel('A')+d.T.ravel('A')
10 loops, best of 3: 29.6 ms per loop

在这种情况下，通过切换到内存视图、添加，然后重新整形回来，可以轻松恢复性能。为了进一步检查问题并了解它并不总是那么容易解决，让我们考虑一下在 NumPy 中处理图像缓冲区的简单代码。

图像合成示例#

对于更现实的示例，请考虑图像缓冲区。图像通常以 Fortran 连续顺序存储，并且颜色通道可以被视为结构化的“RGB”类型或长度为 3 的额外维度。生成的内存布局既不是 C 连续也不是 Fortran 连续，但由于 ndarray 的灵活性，可以直接在 NumPy 中使用。这看起来很理想，因为它使内存布局与典型的 C 或 C++ 图像代码兼容，同时在 Python 中提供自然访问。获取像素 (x,y) 的颜色只需“image[x,y]”。

事实证明，此布局在 NumPy 中的性能非常差。以下代码创建了两个黑色图像，并在其上执行“覆盖”合成操作。

In [9]: image1 = np.zeros((1080,1920,3), dtype=np.float32).swapaxes(0,1)
In [10]: alpha1 = np.zeros((1080,1920,1), dtype=np.float32).swapaxes(0,1)
In [11]: image2 = np.zeros((1080,1920,3), dtype=np.float32).swapaxes(0,1)
In [12]: alpha2 = np.zeros((1080,1920,1), dtype=np.float32).swapaxes(0,1)
In [13]: def composite_over(im1, al1, im2, al2):
   ....:     return (im1 + (1-al1)*im2, al1 + (1-al1)*al2)

In [14]: timeit composite_over(image1,alpha1,image2,alpha2)
1 loops, best of 3: 3.51 s per loop

如果我们放弃方便的布局，并使用 C 连续默认值，则性能会提高大约七倍。

In [16]: image1 = np.zeros((1080,1920,3), dtype=np.float32)
In [17]: alpha1 = np.zeros((1080,1920,1), dtype=np.float32)
In [18]: image2 = np.zeros((1080,1920,3), dtype=np.float32)
In [19]: alpha2 = np.zeros((1080,1920,1), dtype=np.float32)

In [20]: timeit composite_over(image1,alpha1,image2,alpha2)
1 loops, best of 3: 581 ms per loop

但这并非全部，因为事实证明广播 alpha 通道也会造成性能损失。如果我们使用 3 个值而不是一个值的 alpha 通道，则会得到

In [21]: image1 = np.zeros((1080,1920,3), dtype=np.float32)
In [22]: alpha1 = np.zeros((1080,1920,3), dtype=np.float32)
In [23]: image2 = np.zeros((1080,1920,3), dtype=np.float32)
In [24]: alpha2 = np.zeros((1080,1920,3), dtype=np.float32)

In [25]: timeit composite_over(image1,alpha1,image2,alpha2)
1 loops, best of 3: 313 ms per loop

最后比较一下，让我们看看当我们使用一维数组以确保只有一个循环执行计算时的性能。

In [26]: image1 = np.zeros((1080*1920*3), dtype=np.float32)
In [27]: alpha1 = np.zeros((1080*1920*3), dtype=np.float32)
In [28]: image2 = np.zeros((1080*1920*3), dtype=np.float32)
In [29]: alpha2 = np.zeros((1080*1920*3), dtype=np.float32)

In [30]: timeit composite_over(image1,alpha1,image2,alpha2)
1 loops, best of 3: 312 ms per loop

为了获得参考性能数字，我用 C 直接实现了这个简单的操作（注意使用与 NumPy 相同的编译选项）。如果我模拟 Python 代码的内存分配和布局，性能大约为 0.3 秒，与 NumPy 的性能非常一致。将操作组合成一个传递将时间缩短至大约 0.15 秒。

此示例的一个细微变化是使用一个具有四个通道的内存块 (1920,1080,4) 而不是单独的图像和 alpha。这在图像处理应用程序中更为常见，以下是使用 C 连续布局的情况。

In [31]: image1 = np.zeros((1080,1920,4), dtype=np.float32)
In [32]: image2 = np.zeros((1080,1920,4), dtype=np.float32)
In [33]: def composite_over(im1, im2):
   ....:     ret = (1-im1[:,:,-1])[:,:,np.newaxis]*im2
   ....:     ret += im1
   ....:     return ret

In [34]: timeit composite_over(image1,image2)
1 loops, best of 3: 481 ms per loop

要查看所提议的新迭代器实现可以带来的改进，请转到文档底部靠近 API 提案之后的示例继续部分。

改善缓存一致性#

为了从 UFunc 调用中获得最佳性能，内存读取模式应尽可能规则。现代 CPU 尝试预测内存读/写模式并在提前填充缓存。最可预测的模式是所有输入和输出以相同的顺序进行顺序处理。

我建议默认情况下，UFunc 输出的内存布局应尽可能接近输入。只要存在歧义或不匹配，它默认为 C 连续布局。

要了解如何实现这一点，我们首先考虑广播后所有输入的步长。通过确定一组步长是否兼容和/或模棱两可，我们可以确定最大限度地提高一致性的输出内存布局。

在广播中，输入形状首先通过预先添加单一维度转换为广播形状，然后创建广播步长，其中任何单一维度的步长都设置为零。

步长也可能是负数，在某些情况下，可以将其标准化以适合以下讨论。如果特定轴的所有步长都是负数或零，则可以在适当地调整基本数据指针后对该维度的步长取反。

这是一个关于三个具有 C 连续布局的输入如何导致广播步长的示例。为了简化起见，示例使用 1 的元素大小。

输入形状

(5,3,7)

(5,3,1)

(1,7)

广播形状

(5,3,7)

(5,3,1)

(1,1,7)

广播步长

(21,7,1)

(3,1,0)

(0,0,1)

兼容步长 - 如果存在轴的排列，使得集合中每个步长（不包括为零的条目）的步长都递减，则一组步长是兼容的。

上面的示例使用恒等排列满足定义。在动机图像示例中，如果我们分离或不分离颜色和 alpha 信息，则步长会略有不同。此处证明兼容性的排列是转置 (0,1)。

输入/广播形状

图像 (1920, 1080, 3)

Alpha (1920, 1080, 1)

广播步长（分离）

(3,5760,1)

(1,1920,0)

广播步长（一起）

(4,7680,1)

(4,7680,0)

模棱两可的步长 - 如果存在多个轴的排列，使得集合中每个步长（不包括为零的条目）的步长都递减，则一组兼容步长是模棱两可的。

这通常发生在集合中的步长在每个轴上都有 0 步长时。最简单的示例是在二维中，如下所示。

广播形状

(1,3)

(5,1)

广播步长

(0,1)

(1,0)

但是，也可能存在没有单个输入强制整个布局的明确兼容步长，例如在此示例中

广播形状

(1,3,4)

(5,3,1)

广播步长

(0,4,1)

(3,1,0)

面对模棱两可，我们可以选择完全放弃步长兼容的事实，或者尝试通过添加额外的约束来解决模棱两可。我认为适当的选择是通过选择最接近 C 连续但仍与输入步长兼容的内存布局来解决它。

输出布局选择算法#

我们希望生成的输出 ndarray 内存布局如下

一致/明确的步长

唯一的、一致的布局

一致/模棱两可的步长

最接近 C 连续的一致布局

不一致的步长

C 连续

以下是计算输出布局排列的算法的伪代码。

perm = range(ndim) # Identity, i.e. C-contiguous
# Insertion sort, ignoring 0-strides
# Note that the sort must be stable, and 0-strides may
# be reordered if necessary, but should be moved as little
# as possible.
for i0 = 1 to ndim-1:
    # ipos is where perm[i0] will get inserted
    ipos = i0
    j0 = perm[i0]
    for i1 = i0-1 to 0:
        j1 = perm[i1]
        ambig, shouldswap = True, False
        # Check whether any strides are ordered wrong
        for strides in broadcast_strides:
            if strides[j0] != 0 and strides[j1] != 0:
                if strides[j0] > strides[j1]:
                    # Only set swap if it's still ambiguous.
                    if ambig:
                        shouldswap = True
                else:
                    # Set swap even if it's not ambiguous,
                    # because not swapping is the choice
                    # for conflicts as well.
                    shouldswap = False
                ambig = False
        # If there was an unambiguous comparison, either shift ipos
        # to i1 or stop looking for the comparison
        if not ambig:
            if shouldswap:
                ipos = i1
            else:
                break
    # Insert perm[i0] into the right place
    if ipos != i0:
       for i1 = i0-1 to ipos:
         perm[i1+1] = perm[i1]
       perm[ipos] = j0
# perm is now the closest consistent ordering to C-contiguous
return perm

合并维度#

在许多情况下，内存布局允许使用一维循环而不是跟踪迭代器中的多个坐标。现有代码在数据为 C 连续时已经利用了这一点，但由于我们正在重新排序轴，因此我们可以更普遍地应用此优化。

一旦迭代步长已排序为单调递减，任何可以合并的维度都并排放置。如果对于所有操作数，按步长[i+1] 步进 shape[i+1] 次与按步长[i] 步进相同，或者步长[i+1]*shape[i+1] == 步长[i]，则可以将维度 i 和 i+1 合并为一个维度。

以下是合并的伪代码。

# Figure out which pairs of dimensions can be coalesced
can_coalesce = [False]*ndim
for strides, shape in zip(broadcast_strides, broadcast_shape):
    for i = 0 to ndim-2:
        if strides[i+1]*shape[i+1] == strides[i]:
            can_coalesce[i] = True
# Coalesce the types
new_ndim = ndim - count_nonzero(can_coalesce)
for strides, shape in zip(broadcast_strides, broadcast_shape):
    j = 0
    for i = 0 to ndim-1:
        # Note that can_coalesce[ndim-1] is always False, so
        # there is no out-of-bounds access here.
        if can_coalesce[i]:
            shape[i+1] = shape[i]*shape[i+1]
        else:
            strides[j] = strides[i]
            shape[j] = shape[i]
            j += 1

内部循环专门化#

专门化完全由内部循环函数处理，因此此优化与其他优化无关。已经进行了一些专门化，例如减少操作。该想法在 https://projects.scipy.org.cn/numpy/wiki/ProjectIdeas 中提到，“使用 intrinsics（SSE 指令）加速 NumPy 中的低级循环。”

以下是涵盖加法/减法/乘法/除法重要情况的两个参数函数的一些可能性。

  • 第一个或第二个参数是单个值（即 0 步长值）并且不与输出发生别名。arr = arr + 1; arr = 1 + arr

    • 可以加载一次常量，而不是每次都从内存中重新加载它

  • 步长与数据类型的尺寸匹配。例如，C 连续或 Fortran 连续数据

    • 可以使用简单的循环而不使用步长

  • 步长与数据类型的尺寸匹配，并且它们都为 16 字节对齐（或与 16 字节对齐的偏移量相同）

    • 可以使用 SSE 同时处理多个值

  • 第一个输入和输出是相同的值（即减少操作）。

    • 这已经在现有代码中的许多 UFunc 中专门化了

以上情况通常不是互斥的，例如，当步长与数据类型大小匹配时，常量参数可以与 SSE 结合使用，并且减少也可以使用 SSE 进行优化。

实现细节#

除了内部循环专门化之外，讨论的优化会显着影响 ufunc_object.c 和用于执行广播的 PyArrayIterObject/PyArrayMultiIterObject。通常，应该能够在需要时模拟当前行为，但我认为默认情况下应该生成和操作将提供最佳性能的内存布局。

为了支持新的缓存友好型行为，我们为任何 order= 参数引入了一个新的选项“K”（表示“保持”）。

建议的‘order=’标志如下所示

‘C’

C 连续布局

‘F’

Fortran 连续布局

‘A’

如果输入具有 Fortran 连续布局，则为‘F’，否则为‘C’（“任意连续”）

‘K’

等效于‘C’的布局，然后是轴的一些排列，尽可能接近输入的布局（“保持布局”）

或者作为枚举类型

/* For specifying array memory layout or iteration order */
typedef enum {
        /* Fortran order if inputs are all Fortran, C otherwise */
        NPY_ANYORDER=-1,
        /* C order */
        NPY_CORDER=0,
        /* Fortran order */
        NPY_FORTRANORDER=1,
        /* An order as close to the inputs as possible */
        NPY_KEEPORDER=2
} NPY_ORDER;

也许一个好的策略是首先实现此处讨论的功能，而无需更改默认值。一旦它们被实现并经过充分测试，默认值就可以从 order='C' 更改为 order='K'，在所有合适的地方。此外，UFuncs 应该获得一个 order= 参数来控制其输出的布局。

迭代器可以进行自动类型转换，我已经创建了一系列逐渐更宽松的类型转换规则。也许对于 2.0 版本，NumPy 可以采用这种枚举类型作为其处理类型转换的首选方式。

/* For specifying allowed casting in operations which support it */
typedef enum {
        /* Only allow identical types */
        NPY_NO_CASTING=0,
        /* Allow identical and byte swapped types */
        NPY_EQUIV_CASTING=1,
        /* Only allow safe casts */
        NPY_SAFE_CASTING=2,
        /* Allow safe casts and casts within the same kind */
        NPY_SAME_KIND_CASTING=3,
        /* Allow any casts */
        NPY_UNSAFE_CASTING=4
} NPY_CASTING;

迭代器重写#

根据代码分析，对现有的迭代对象进行重构以实现这些优化似乎非常困难。此外，迭代器的一些用法需要修改内部值或标志，因此使用迭代器的代码无论如何都必须更改。因此，我们建议创建一个新的迭代器对象，该对象包含现有的迭代器功能并将其扩展以考虑优化。

替换迭代器的高级目标包括

  • 内存使用量小，内存分配次数少。

  • 简单情况（如扁平数组）应该几乎没有开销。

  • 将单迭代和多迭代合并到一个对象中。

应该提供给用户代码的功能

  • 以 C、Fortran 或“最快”（默认）顺序迭代。

  • 如果需要，跟踪 C 样式或 Fortran 样式的扁平索引（现有迭代器始终跟踪 C 样式索引）。这可以独立于迭代顺序进行。

  • 如果需要，跟踪坐标（现有迭代器需要手动更改内部迭代器标志才能保证这一点）。

  • 跳过最后一个内部维度的迭代，以便可以使用内部循环对其进行处理。

  • 跳转到数组中的特定坐标。

  • 迭代任意子集的轴（例如，支持一次减少多个轴）。

  • 如果提供了 NULL 输入，则能够自动分配输出参数。这些输出应该具有与迭代顺序匹配的内存布局，并且是 order='K' 支持的机制。

  • 自动复制和/或缓冲不满足类型/字节序/对齐要求的输入。无论执行了哪些缓冲或复制操作，调用者的迭代内部循环都应该相同。

实现说明

  • 用户代码绝不能触碰迭代器内部。这允许将来大幅更改内部内存布局，如果找到更高性能的实现策略。

  • 使用函数指针而不是宏进行迭代。这样，可以为常见情况创建专门化，例如当 ndim 很小时，对于不同的标志设置，以及当迭代的数组数量很少时。此外，可以规定一个迭代模式，该模式首先复制函数指针，以允许编译器将函数指针保存在寄存器中。

  • 动态创建内存布局，以最大程度地减少迭代器占用的缓存行数（对于 LP64，sizeof(PyArrayIterObject) 约为 2.5KB，而加法等二元运算需要三个这样的迭代器用于多迭代器）。

  • 将 C-API 对象与 Python 引用计数隔离，以便可以从 C 自然地使用它。然后，Python 对象成为 C 迭代器的包装器。这类似于 PEP 3118 中 Py_buffer 和 memoryview 的设计分离。

建议的迭代器内存布局#

以下结构描述了迭代器内存。所有项目都打包在一起，这意味着标志、ndim 和 niter 的不同值将产生略微不同的布局。

struct {
    /* Flags indicate what optimizations have been applied, and
     * affect the layout of this struct. */
    uint32 itflags;
    /* Number of iteration dimensions.  If FLAGS_HASCOORDS is set,
     * it matches the creation ndim, otherwise it may be smaller.  */
    uint16 ndim;
    /* Number of objects being iterated.  This is fixed at creation time. */
    uint16 niter;

    /* The number of times the iterator will iterate */
    intp itersize;

    /* The permutation is only used when FLAGS_HASCOORDS is set,
     * and is placed here so its position depends on neither ndim
     * nor niter. */
    intp perm[ndim];

    /* The data types of all the operands */
    PyArray_Descr *dtypes[niter];
    /* Backups of the starting axisdata 'ptr' values, to support Reset */
    char *resetdataptr[niter];
    /* Backup of the starting index value, to support Reset */
    npy_intp resetindex;

    /* When the iterator is destroyed, Py_XDECREF is called on all
       these objects */
    PyObject *objects[niter];

    /* Flags indicating read/write status and buffering
     * for each operand. */
    uint8 opitflags[niter];
    /* Padding to make things intp-aligned again */
    uint8 padding[];

    /* If some or all of the inputs are being buffered */
    #if (flags&FLAGS_BUFFERED)
    struct buffer_data {
        /* The size of the buffer, and which buffer we're on.
         * the i-th iteration has i = buffersize*bufferindex+pos
         */
        intp buffersize;
        /* For tracking position inside the buffer */
        intp size, pos;
        /* The strides for the pointers */
        intp stride[niter];
        /* Pointers to the data for the current iterator position.
         * The buffer_data.value ptr[i] equals either
         * axis_data[0].ptr[i] or buffer_data.buffers[i] depending
         * on whether copying to the buffer was necessary.
         */
        char* ptr[niter];
        /* Functions to do the copyswap and casting necessary */
        transferfn_t readtransferfn[niter];
        void *readtransferdata[niter];
        transferfn_t writetransferfn[niter];
        void *writetransferdata[niter];
        /* Pointers to the allocated buffers for operands
         * which the iterator determined needed buffering
         */
        char *buffers[niter];
    };
    #endif /* FLAGS_BUFFERED */

    /* Data per axis, starting with the most-frequently
     * updated, and in decreasing order after that. */
    struct axis_data {
        /* The shape of this axis */
        intp shape;
        /* The current coordinate along this axis */
        intp coord;
        /* The operand and index strides for this axis */
        intp stride[niter];
        #if (flags&FLAGS_HASINDEX)
            intp indexstride;
        #endif
        /* The operand pointers and index values for this axis */
        char* ptr[niter];
        #if (flags&FLAGS_HASINDEX)
            intp index;
        #endif
    }[ndim];
};

axis_data 结构数组的顺序是按增量更新的快速性递增排列的。如果 perm 是单位矩阵，这意味着它与 C 顺序相反。这样做是为了使最常访问的数据项最靠近结构的开头（其中包含常见属性），从而提高缓存一致性。它还简化了 iternext 调用，同时使 getcoord 和相关函数稍微复杂了一些。

建议的迭代器 API#

现有的迭代器 API 包括 PyArrayIter_Check、PyArray_Iter* 和 PyArray_ITER_* 等函数。多迭代器数组包括 PyArray_MultiIter*、PyArray_Broadcast 和 PyArray_RemoveSmallest。新的迭代器设计用单个对象和关联的 API 替换了所有这些功能。新 API 的一个目标是，现有迭代器的所有用途都应该能够在不付出太多努力的情况下替换为新的迭代器。

选择的 C-API 命名约定基于 numpy-refactor 分支中的约定，其中 libndarray 将数组命名为 NpyArray，并将函数命名为 NpyArray_*。迭代器命名为 NpyIter，函数命名为 NpyIter_*。

Python 暴露将迭代器命名为 np.nditer。此迭代器的一种可能的发布策略是在 1.X（1.6？）版本中添加迭代器，但 NumPy 代码不使用它。然后，2.0 版本可以发布并将其完全集成。如果选择此策略，则应在 1.X 版本发布之前尽可能最终确定命名约定和 API。名称 np.iter 不能使用，因为它与 Python 内置的 iter 冲突。我建议在 Python 中使用名称 np.nditer，因为它目前未使用。

除了为新迭代器设定的性能目标外，API 似乎还可以重构以更好地支持一些常见的 NumPy 编程习惯用法。

通过将 UFunc 代码中的一些功能移动到迭代器中，它应该可以简化扩展代码，这些代码希望在不完全符合 UFunc 范例的情况下模拟 UFunc 行为。特别是，模拟 UFunc 缓冲行为并非易事。

旧 -> 新迭代器 API 转换#

对于常规迭代器

PyArray_IterNew

NpyIter_New

PyArray_IterAllButAxis

NpyIter_New + axes 参数 **或** 迭代器标志 NPY_ITER_NO_INNER_ITERATION

PyArray_BroadcastToShape

**不支持**（但如果需要，可以支持）

PyArrayIter_Check

需要在 Python 暴露中添加此功能

PyArray_ITER_RESET

NpyIter_Reset

PyArray_ITER_NEXT

来自 NpyIter_GetIterNext 的函数指针

PyArray_ITER_DATA

NpyIter_GetDataPtrArray

PyArray_ITER_GOTO

NpyIter_GotoCoords

PyArray_ITER_GOTO1D

NpyIter_GotoIndex

PyArray_ITER_NOTDONE

函数指针 iternext 的返回值

对于多迭代器

PyArray_MultiIterNew

NpyIter_MultiNew

PyArray_MultiIter_RESET

NpyIter_Reset

PyArray_MultiIter_NEXT

来自 NpyIter_GetIterNext 的函数指针

PyArray_MultiIter_DATA

NpyIter_GetDataPtrArray

PyArray_MultiIter_NEXTi

**不支持**（始终是同步迭代）

PyArray_MultiIter_GOTO

NpyIter_GotoCoords

PyArray_MultiIter_GOTO1D

NpyIter_GotoIndex

PyArray_MultiIter_NOTDONE

函数指针 iternext 的返回值

PyArray_Broadcast

由 NpyIter_MultiNew 处理

PyArray_RemoveSmallest

迭代器标志 NPY_ITER_NO_INNER_ITERATION

对于其他 API 调用

PyArray_ConvertToCommonType

迭代器标志 NPY_ITER_COMMON_DTYPE

迭代器指针类型#

迭代器结构是内部生成的，但仍然需要一个类型来在将错误类型传递到 API 时提供警告和/或错误。我们使用不完整结构的 typedef 来实现这一点

typedef struct NpyIter_InternalOnly NpyIter;

构造和析构#

NpyIter* NpyIter_New(PyArrayObject* op, npy_uint32 flags, NPY_ORDER order, NPY_CASTING casting, PyArray_Descr* dtype, npy_intp a_ndim, npy_intp *axes, npy_intp buffersize)

为给定的 NumPy 数组对象 op 创建一个迭代器。

可以在 flags 中传递的标志是 NpyIter_MultiNew 中记录的全局标志和每个操作数标志的任意组合，除了 NPY_ITER_ALLOCATE。

可以将任何 NPY_ORDER 枚举值传递给 order。为了高效迭代，NPY_KEEPORDER 是最佳选择，其他顺序强制执行特定的迭代模式。

可以将任何 NPY_CASTING 枚举值传递给 casting。这些值包括 NPY_NO_CASTING、NPY_EQUIV_CASTING、NPY_SAFE_CASTING、NPY_SAME_KIND_CASTING 和 NPY_UNSAFE_CASTING。要允许发生类型转换，还必须启用复制或缓冲。

如果 dtype 不是 NULL，则它需要该数据类型。如果允许复制，它将在数据可转换时进行临时复制。如果启用了 UPDATEIFCOPY，它还将在迭代器销毁时使用另一个类型转换复制回数据。

如果 a_ndim 大于零，则还必须提供 axes。在这种情况下，axes 是 op 轴的一个大小为 a_ndim 的数组。axes 中的值 -1 表示 newaxis。在 axes 数组中，轴不能重复。

如果 buffersize 为零，则使用默认缓冲区大小，否则它指定要使用的缓冲区大小。建议使用 2 的幂作为缓冲区大小，例如 512 或 1024。

如果发生错误，则返回 NULL，否则返回分配的迭代器。

要创建类似于旧迭代器的迭代器，这应该可以工作。

iter = NpyIter_New(op, NPY_ITER_READWRITE,
                    NPY_CORDER, NPY_NO_CASTING, NULL, 0, NULL);

如果要使用对齐的 double 代码编辑数组，但顺序无关紧要，则可以使用以下方法。

dtype = PyArray_DescrFromType(NPY_DOUBLE);
iter = NpyIter_New(op, NPY_ITER_READWRITE |
                    NPY_ITER_BUFFERED |
                    NPY_ITER_NBO,
                    NPY_ITER_ALIGNED,
                    NPY_KEEPORDER,
                    NPY_SAME_KIND_CASTING,
                    dtype, 0, NULL);
Py_DECREF(dtype);

NpyIter* NpyIter_MultiNew(npy_intp niter, PyArrayObject** op, npy_uint32 flags, NPY_ORDER order, NPY_CASTING casting, npy_uint32 *op_flags, PyArray_Descr** op_dtypes, npy_intp oa_ndim, npy_intp **op_axes, npy_intp buffersize)

为广播在 op 中提供的 niter 个数组对象创建迭代器。

对于正常用法，将 oa_ndim 用于 0，并将 op_axes 用于 NULL。有关这些参数的说明，请参见下文，这些参数允许自定义手动广播以及重新排序和省略轴。

可以将任何 NPY_ORDER 枚举值传递给 order。为了高效迭代，NPY_KEEPORDER 是最佳选择，其他顺序强制执行特定的迭代模式。当使用 NPY_KEEPORDER 时，如果还希望确保沿某个轴不反转迭代，则应传递标志 NPY_ITER_DONT_NEGATE_STRIDES。

可以将任何 NPY_CASTING 枚举值传递给 casting。这些值包括 NPY_NO_CASTING、NPY_EQUIV_CASTING、NPY_SAFE_CASTING、NPY_SAME_KIND_CASTING 和 NPY_UNSAFE_CASTING。要允许发生类型转换，还必须启用复制或缓冲。

如果 op_dtypes 不是 NULL，则它为每个 op[i] 指定数据类型或 NULL。

参数 oa_ndim 在非零时，指定将使用自定义广播迭代的维度数。如果提供了它，则还必须提供 op_axes。这两个参数允许您详细控制操作数数组的轴如何匹配在一起并进行迭代。在 op_axes 中，您必须提供一个大小为 niter 的数组，其中包含类型为 npy_intp 的大小为 oa_ndim 的数组指针。如果 op_axes 中的条目为 NULL，则将应用正常的广播规则。在 op_axes[j][i] 中存储了 op[j] 的有效轴或 -1（表示 newaxis）。在每个 op_axes[j] 数组中，轴不能重复。以下示例是如何将正常的广播应用于 3 维数组、2 维数组、1 维数组和标量。

npy_intp oa_ndim = 3;               /* # iteration axes */
npy_intp op0_axes[] = {0, 1, 2};    /* 3-D operand */
npy_intp op1_axes[] = {-1, 0, 1};   /* 2-D operand */
npy_intp op2_axes[] = {-1, -1, 0};  /* 1-D operand */
npy_intp op3_axes[] = {-1, -1, -1}  /* 0-D (scalar) operand */
npy_intp *op_axes[] = {op0_axes, op1_axes, op2_axes, op3_axes};

如果 buffersize 为零，则使用默认缓冲区大小，否则它指定要使用的缓冲区大小。建议使用 2 的幂作为缓冲区大小，例如 512 或 1024。

如果发生错误，则返回 NULL，否则返回分配的迭代器。

可以在 flags 中传递的标志，应用于整个迭代器，如下所示：

NPY_ITER_C_INDEX、NPY_ITER_F_INDEX

使迭代器跟踪与 C 或 Fortran 顺序匹配的索引。这些选项是互斥的。

NPY_ITER_COORDS

使迭代器跟踪数组坐标。这可以防止迭代器合并轴以生成更大的内部循环。

NPY_ITER_NO_INNER_ITERATION

使迭代器跳过最内层循环的迭代，允许迭代器的用户处理它。

此标志与 NPY_ITER_C_INDEX、NPY_ITER_F_INDEX 和 NPY_ITER_COORDS 不兼容。

NPY_ITER_DONT_NEGATE_STRIDES

当为 order 参数指定 NPY_KEEPORDER 时，这仅会影响迭代器。默认情况下，使用 NPY_KEEPORDER 时，迭代器会反转步长为负的轴，以便以正向遍历内存。这会禁用此步骤。如果您想使用轴的基础内存顺序，但不想反转轴，请使用此标志。例如，这就是 numpy.ravel(a, order='K') 的行为。

NPY_ITER_COMMON_DTYPE

使迭代器将所有操作数转换为公共数据类型，该类型根据 ufunc 类型提升规则计算得出。必须设置每个操作数的标志，以便允许适当的转换，并且必须启用复制或缓冲。

如果事先知道公共数据类型，请不要使用此标志。相反，为所有操作数设置请求的数据类型。

NPY_ITER_REFS_OK

指示可以接受并在迭代器中使用具有引用类型（对象数组或包含对象类型的结构化数组）的数组。如果启用了此标志，则调用方必须确保检查 NpyIter_IterationNeedsAPI(iter) 是否为真，在这种情况下，它可能不会在迭代期间释放 GIL。

NPY_ITER_ZEROSIZE_OK

指示应允许大小为零的数组。由于典型的迭代循环不能自然地处理大小为零的数组，因此您必须在进入迭代循环之前检查 IterSize 是否不为零。

NPY_ITER_REDUCE_OK

允许具有步长为零且大小大于一的维度的可写操作数。请注意，此类操作数必须是读/写操作数。

启用缓冲时，这也切换到一种特殊的缓冲模式，该模式会根据需要缩短循环长度，以避免践踏正在减少的值。

请注意，如果您想对自动分配的输出进行归约，则必须使用 NpyIter_GetOperandArray 获取其引用，然后在执行迭代循环之前将每个值设置为归约单元。在缓冲归约的情况下，这意味着您还必须指定标志 NPY_ITER_DELAY_BUFALLOC，然后在初始化分配的操作数以准备缓冲区后重置迭代器。

NPY_ITER_RANGED

启用对完整 iterindex 范围 [0, NpyIter_IterSize(iter)) 的子范围的迭代的支持。使用函数 NpyIter_ResetToIterIndexRange 指定迭代范围。

此标志仅可在启用 NPY_ITER_BUFFERED 时与 NPY_ITER_NO_INNER_ITERATION 一起使用。这是因为在没有缓冲的情况下，内部循环始终是最内层迭代维的大小，并且允许它被分割将需要特殊处理，有效地使其更像缓冲版本。

NPY_ITER_BUFFERED

使迭代器存储缓冲数据，并使用缓冲来满足数据类型、对齐和字节序要求。要缓冲操作数，请不要指定 NPY_ITER_COPY 或 NPY_ITER_UPDATEIFCOPY 标志，因为它们将覆盖缓冲。缓冲对于使用迭代器的 Python 代码特别有用，因为它允许一次处理更大的数据块以摊销 Python 解释器的开销。

如果与 NPY_ITER_NO_INNER_ITERATION 一起使用，则调用方的内部循环可能会获得比没有缓冲时更大的数据块，因为步长的布局方式。

请注意，如果操作数被赋予标志 NPY_ITER_COPY 或 NPY_ITER_UPDATEIFCOPY，则优先进行复制而不是缓冲。当数组被广播以便需要复制元素以获得恒定的步长时，仍然会发生缓冲。

在正常缓冲中，每个内部循环的大小等于缓冲区大小，或者如果指定了 NPY_ITER_GROWINNER，则可能更大。如果启用了 NPY_ITER_REDUCE_OK 并发生归约，则内部循环可能会变小，具体取决于归约的结构。

NPY_ITER_GROWINNER

启用缓冲时，这允许在不需要缓冲时内部循环的大小增长。如果您正在对所有数据进行直通传递，而不是任何具有用于每个内部循环的小型缓存友好临时值数组的操作，则此选项最适合使用。

NPY_ITER_DELAY_BUFALLOC

启用缓冲时，这会延迟缓冲区的分配，直到调用其中一个 NpyIter_Reset* 函数。此标志的存在是为了避免在为多线程迭代创建缓冲迭代器的多个副本时浪费地复制缓冲数据。

此标志的另一个用途是设置归约操作。在创建迭代器并且迭代器自动分配归约输出（确保使用 READWRITE 访问权限）之后，可以将其值初始化为归约单元。使用 NpyIter_GetOperandArray 获取对象。然后，调用 NpyIter_Reset 以分配缓冲区并将其填充初始值。

可在 op_flags[i] 中传递的标志，其中 0 <= i < niter

NPY_ITER_READWRITE、NPY_ITER_READONLY、NPY_ITER_WRITEONLY

指示迭代器的用户将如何读取或写入 op[i]。每个操作数必须精确指定这些标志之一。

NPY_ITER_COPY

如果 op[i] 不满足构造函数标志和参数指定的数据类型或对齐要求，则允许对其进行复制。

NPY_ITER_UPDATEIFCOPY

触发 NPY_ITER_COPY，并且当数组操作数被标记为写入并且被复制时，在迭代器被销毁时，会导致副本中的数据被复制回 op[i]。

如果操作数被标记为只写并且需要复制，则会创建一个未初始化的临时数组，然后在销毁时将其复制回 op[i]，而不是执行不必要的复制操作。

NPY_ITER_NBO、NPY_ITER_ALIGNED、NPY_ITER_CONTIG

使迭代器提供 op[i] 的数据，该数据采用本机字节序、根据 dtype 要求对齐、连续或任意组合。

默认情况下，迭代器会生成指向提供的数组的指针，这些指针可能是对齐的或未对齐的，并且具有任何字节序。如果未启用复制或缓冲并且操作数数据不满足约束条件，则会引发错误。

连续约束仅适用于内部循环，后续的内部循环可能具有任意指针更改。

如果请求的数据类型采用非本机字节序，则 NBO 标志会覆盖它，并且请求的数据类型会被转换为本机字节序。

NPY_ITER_ALLOCATE

这是用于输出数组的，并且需要设置标志 NPY_ITER_WRITEONLY。如果 op[i] 为 NULL，则创建一个具有最终广播维度的新数组，以及与迭代器的迭代顺序匹配的布局。

当 op[i] 为 NULL 时，请求的数据类型 op_dtypes[i] 也可能为 NULL，在这种情况下，它会根据标记为可读的数组的数据类型自动生成。生成 dtype 的规则与 UFuncs 相同。需要注意的是在所选 dtype 中处理字节序。如果只有一个输入，则按原样使用输入的 dtype。否则，如果组合了多个输入 dtype，则输出将采用本机字节序。

使用此标志分配后，调用方可以通过调用 NpyIter_GetOperandArray 并获取返回的 C 数组中的第 i 个对象来检索新数组。调用方必须对其调用 Py_INCREF 以声明对该数组的引用。

NPY_ITER_NO_SUBTYPE

用于 NPY_ITER_ALLOCATE，此标志禁用为输出分配数组子类型，强制其成为直通 ndarray。

TODO：也许最好引入一个函数 NpyIter_GetWrappedOutput 并删除此标志？

NPY_ITER_NO_BROADCAST

确保输入或输出与迭代维度完全匹配。

NPY_ITER_WRITEABLE_REFERENCES

默认情况下，如果迭代器具有可写操作数，其中数据类型涉及 Python 引用，则迭代器在创建时会失败。添加此标志表示使用迭代器的代码知道这种可能性并正确处理它。

NpyIter *NpyIter_Copy(NpyIter *iter)

创建给定迭代器的副本。此函数主要用于启用数据的多线程迭代。

TODO：将此移动到有关多线程迭代的部分。

多线程迭代的推荐方法是首先创建一个迭代器，其标志为 NPY_ITER_NO_INNER_ITERATION、NPY_ITER_RANGED、NPY_ITER_BUFFERED、NPY_ITER_DELAY_BUFALLOC，以及可能还有 NPY_ITER_GROWINNER。为每个线程创建此迭代器的副本（第一个迭代器减去一个）。然后，获取迭代索引范围 [0, NpyIter_GetIterSize(iter)) 并将其分成任务，例如使用 TBB parallel_for 循环。当线程获得要执行的任务时，它会通过调用 NpyIter_ResetToIterIndexRange 并遍历整个范围来使用其迭代器副本。

在多线程代码或不持有 Python GIL 的代码中使用迭代器时，必须注意仅调用在此上下文中安全的函数。 NpyIter_Copy 无法在没有 Python GIL 的情况下安全调用，因为它会增加 Python 引用。 Reset* 和一些其他函数可以通过传入 errmsg 参数作为非 NULL 来安全调用，以便函数通过它传递错误而不是设置 Python 异常。

int NpyIter_UpdateIter(NpyIter *iter, npy_intp i, npy_uint32 op_flags, NPY_CASTING casting, PyArray_Descr *dtype) 未实现

更新迭代器中的第 i 个操作数，使其可能具有新的数据类型或更严格的标志属性。此操作的用例是允许自动分配根据标准 NumPy 类型提升规则确定输出数据类型，然后在此函数中转换输入以及可能在处理期间自动转换输出到不同的数据类型。

此操作仅可在将 NPY_ITER_COORDS 作为标志传递给迭代器时执行。如果不需要坐标，请在不再需要调用 NpyIter_UpdateIter 时调用函数 NpyIter_RemoveCoords()。

如果第 i 个操作数已被复制，则会抛出错误。要避免这种情况，请省略所有标志，除了任何操作数的读/写指示符外，这些操作数稍后将对 NpyIter_UpdateIter 进行调用。

可在 op_flags 中传递的标志为 NPY_ITER_COPY、NPY_ITER_UPDATEIFCOPY、NPY_ITER_NBO、NPY_ITER_ALIGNED、NPY_ITER_CONTIG。

int NpyIter_RemoveAxis(NpyIter *iter, npy_intp axis)

从迭代中删除一个轴。这要求为迭代器创建设置了 NPY_ITER_COORDS，并且如果启用了缓冲或正在跟踪索引则不起作用。此函数还会将迭代器重置为其初始状态。

这对于设置累积循环很有用，例如。迭代器可以首先使用所有维度（包括累积轴）创建，以便正确创建输出。然后，可以删除累积轴，并以嵌套方式执行计算。

警告：此函数可能会更改迭代器的内部内存布局。迭代器的任何缓存函数或指针都必须重新获取！

返回 NPY_SUCCEED 或 NPY_FAIL。

int NpyIter_RemoveCoords(NpyIter *iter)

如果迭代器具有坐标，则此操作会剥离对坐标的支持，并执行如果不需要坐标则可能的进一步迭代优化。此函数还会将迭代器重置为其初始状态。

警告：此函数可能会更改迭代器的内部内存布局。迭代器的任何缓存函数或指针都必须重新获取！

调用此函数后，NpyIter_HasCoords(iter) 将返回 false。

返回 NPY_SUCCEED 或 NPY_FAIL。

int NpyIter_RemoveInnerLoop(NpyIter *iter)

如果使用了 UpdateIter/RemoveCoords，您可能希望指定标志 NPY_ITER_NO_INNER_ITERATION。此标志不允许与 NPY_ITER_COORDS 一起使用，因此提供此函数以便在调用 NpyIter_RemoveCoords 后启用此功能。此函数还会将迭代器重置为其初始状态。

警告：此函数会更改迭代器的内部逻辑。迭代器的任何缓存函数或指针都必须重新获取！

返回 NPY_SUCCEED 或 NPY_FAIL。

int NpyIter_Deallocate(NpyIter *iter)

释放迭代器对象。这还会释放任何进行的复制，并在必要时触发 UPDATEIFCOPY 行为。

返回 NPY_SUCCEED 或 NPY_FAIL。

int NpyIter_Reset(NpyIter *iter, char **errmsg)

将迭代器重置回其初始状态，回到迭代范围的开头。

返回 NPY_SUCCEED 或 NPY_FAIL。如果 errmsg 不为 NULL，则当返回 NPY_FAIL 时不会设置 Python 异常。相反，*errmsg 将被设置为错误消息。当 errmsg 不为 NULL 时，可以在不持有 Python GIL 的情况下安全地调用此函数。

int NpyIter_ResetToIterIndexRange(NpyIter *iter, npy_intp istart, npy_intp iend, char **errmsg)

重置迭代器并将其限制在 iterindex 范围 [istart, iend) 内。有关如何将其用于多线程迭代的说明，请参见 NpyIter_Copy。这要求将标志 NPY_ITER_RANGED 传递给迭代器构造函数。

如果要同时重置 iterindex 范围和基本指针，可以执行以下操作以避免额外的缓冲区复制（复制此代码时，请务必添加返回代码错误检查）。

/* Set to a trivial empty range */
NpyIter_ResetToIterIndexRange(iter, 0, 0);
/* Set the base pointers */
NpyIter_ResetBasePointers(iter, baseptrs);
/* Set to the desired range */
NpyIter_ResetToIterIndexRange(iter, istart, iend);

返回 NPY_SUCCEED 或 NPY_FAIL。如果 errmsg 不为 NULL，则当返回 NPY_FAIL 时不会设置 Python 异常。相反，*errmsg 将被设置为错误消息。当 errmsg 不为 NULL 时，可以在不持有 Python GIL 的情况下安全地调用此函数。

int NpyIter_ResetBasePointers(NpyIter *iter, char **baseptrs, char **errmsg)

将迭代器重置回其初始状态，但使用 baseptrs 中的值作为数据，而不是正在迭代的数组中的指针。此函数旨在与 op_axes 参数一起使用，由具有两个或多个迭代器的嵌套迭代代码使用。

返回 NPY_SUCCEED 或 NPY_FAIL。如果 errmsg 不为 NULL，则当返回 NPY_FAIL 时不会设置 Python 异常。相反，*errmsg 将被设置为错误消息。当 errmsg 不为 NULL 时，可以在不持有 Python GIL 的情况下安全地调用此函数。

待办：将以下内容移至有关嵌套迭代器的专用部分。

为嵌套迭代创建迭代器需要谨慎。所有迭代器操作数必须完全匹配，否则对 NpyIter_ResetBasePointers 的调用将无效。这意味着不应随意使用自动复制和输出分配。可以通过使用所有转换参数启用其中一个迭代器来创建它，然后使用 NpyIter_GetOperandArray 函数获取已分配的操作数并将其传递给其余迭代器的构造函数，从而仍然使用迭代器的自动数据转换和强制转换功能。

警告：在为嵌套迭代创建迭代器时，代码不得在不同的迭代器中多次使用维度。如果这样做，嵌套迭代将在迭代期间产生越界指针。

警告：在为嵌套迭代创建迭代器时，缓冲只能应用于最内层迭代器。如果将缓冲迭代器用作 baseptrs 的源，它将指向一个小缓冲区而不是数组，并且内部迭代将无效。

使用嵌套迭代器的模式如下

NpyIter *iter1, *iter1;
NpyIter_IterNext_Fn iternext1, iternext2;
char **dataptrs1;

/*
 * With the exact same operands, no copies allowed, and
 * no axis in op_axes used both in iter1 and iter2.
 * Buffering may be enabled for iter2, but not for iter1.
 */
iter1 = ...; iter2 = ...;

iternext1 = NpyIter_GetIterNext(iter1);
iternext2 = NpyIter_GetIterNext(iter2);
dataptrs1 = NpyIter_GetDataPtrArray(iter1);

do {
    NpyIter_ResetBasePointers(iter2, dataptrs1);
    do {
        /* Use the iter2 values */
    } while (iternext2(iter2));
} while (iternext1(iter1));

int NpyIter_GotoCoords(NpyIter *iter, npy_intp *coords)

调整迭代器以指向由 coords 指向的 ndim 个坐标。如果未跟踪坐标、坐标超出范围或禁用内部循环迭代，则返回错误。

返回 NPY_SUCCEED 或 NPY_FAIL。

int NpyIter_GotoIndex(NpyIter *iter, npy_intp index)

调整迭代器以指向指定的 index。如果使用标志 NPY_ITER_C_INDEX 构造迭代器，则 index 是 C 序索引，如果使用标志 NPY_ITER_F_INDEX 构造迭代器，则 index 是 Fortran 序索引。如果未跟踪任何索引、索引超出范围或禁用内部循环迭代，则返回错误。

返回 NPY_SUCCEED 或 NPY_FAIL。

npy_intp NpyIter_GetIterSize(NpyIter *iter)

返回正在迭代的元素数量。这是形状中所有维度的乘积。

npy_intp NpyIter_GetReduceBlockSizeFactor(NpyIter *iter) 未实现

这提供了一个因子，该因子必须除以用于范围迭代的块大小，才能安全地对约简进行多线程处理。如果迭代器没有约简，则返回 1。

当使用范围迭代对约简进行多线程处理时，有两种可能的约简方法

如果对一个小输出进行大规模约简，则为每个线程创建一个初始化为约简单元的临时数组，然后让每个线程将其约简到其临时数组中。完成后，将这些临时数组组合在一起。可以通过观察 NpyIter_GetReduceBlockSizeFactor 返回一个较大的值（例如 NpyIter_GetIterSize 的一半或三分之一）来检测这种情况。还应检查输出是否较小以确保。

如果对一个大输出进行许多小规模约简，并且约简维度是内部维度，则 NpyIter_GetReduceBlockSizeFactor 将返回一个较小的数字，并且只要为多线程选择的块大小为 NpyIter_GetReduceBlockSizeFactor(iter)*n（对于某些 n），则操作将是安全的。

不好的情况是约简维度是迭代器中最外层的循环。例如，如果有一个形状为 (3,1000,1000) 的 C 序数组，并且在维度 0 上进行约简，则 NpyIter_GetReduceBlockSizeFactor 将为 NPY_KEEPORDER 或 NPY_CORDER 迭代顺序返回一个等于 NpyIter_GetIterSize 的大小。虽然这对 CPU 缓存不利，但也许将来可以提供另一种顺序可能性，例如 NPY_REDUCEORDER，它将约简轴推到内部循环，但在其他方面与 NPY_KEEPORDER 相同。

npy_intp NpyIter_GetIterIndex(NpyIter *iter)

获取迭代器的 iterindex，这是一个与迭代器迭代顺序匹配的索引。

void NpyIter_GetIterIndexRange(NpyIter *iter, npy_intp *istart, npy_intp *iend)

获取正在迭代的 iterindex 子范围。如果未指定 NPY_ITER_RANGED，则始终返回范围 [0, NpyIter_IterSize(iter))。

int NpyIter_GotoIterIndex(NpyIter *iter, npy_intp iterindex)

调整迭代器以指向指定的 iterindex。IterIndex 是一个与迭代器迭代顺序匹配的索引。如果 iterindex 超出范围、启用了缓冲或禁用了内部循环迭代，则返回错误。

返回 NPY_SUCCEED 或 NPY_FAIL。

int NpyIter_HasInnerLoop(NpyIter *iter)

如果迭代器处理内部循环，则返回 1；如果调用者需要处理内部循环，则返回 0。这由构造函数标志 NPY_ITER_NO_INNER_ITERATION 控制。

int NpyIter_HasCoords(NpyIter *iter)

如果使用 NPY_ITER_COORDS 标志创建迭代器，则返回 1，否则返回 0。

int NpyIter_HasIndex(NpyIter *iter)

如果使用 NPY_ITER_C_INDEX 或 NPY_ITER_F_INDEX 标志创建迭代器，则返回 1，否则返回 0。

int NpyIter_IsBuffered(NpyIter *iter)

如果使用 NPY_ITER_BUFFERED 标志创建迭代器，则返回 1，否则返回 0。

int NpyIter_IsGrowInner(NpyIter *iter)

如果使用 NPY_ITER_GROWINNER 标志创建迭代器，则返回 1，否则返回 0。

npy_intp NpyIter_GetBufferSize(NpyIter *iter)

如果迭代器已缓冲，则返回正在使用的缓冲区的大小，否则返回 0。

npy_intp NpyIter_GetNDim(NpyIter *iter)

返回正在迭代的维数。如果在迭代器构造函数中未请求坐标，则此值可能小于原始对象中的维数。

npy_intp NpyIter_GetNIter(NpyIter *iter)

返回正在迭代的对象数量。

npy_intp *NpyIter_GetAxisStrideArray(NpyIter *iter, npy_intp axis)

获取指定轴的步长数组。要求迭代器跟踪坐标，并且未启用缓冲。

当您希望以某种方式匹配操作数轴，然后使用 NpyIter_RemoveAxis 手动处理它们的处理时，可以使用此方法。在移除轴之前调用此函数，可以获取手动处理的步长。

错误时返回 NULL。

int NpyIter_GetShape(NpyIter *iter, npy_intp *outshape)

在 outshape 中返回迭代器的广播形状。这只能在支持坐标的迭代器上调用。

返回 NPY_SUCCEED 或 NPY_FAIL。

PyArray_Descr **NpyIter_GetDescrArray(NpyIter *iter)

这会返回一个指向正在迭代的对象的 niter 个数据类型 Descr 的指针。结果指向 iter，因此调用者不会获得对 Descr 的任何引用。

可以在迭代循环之前缓存此指针，调用 iternext 不会更改它。

PyObject **NpyIter_GetOperandArray(NpyIter *iter)

这会返回一个指向正在迭代的 niter 个操作数 PyObject 的指针。结果指向 iter，因此调用者不会获得对 PyObject 的任何引用。

PyObject *NpyIter_GetIterView(NpyIter *iter, npy_intp i)

这会返回对一个新的 ndarray 视图的引用，该视图是数组 NpyIter_GetOperandArray() 中第 i 个对象的视图，其维度和步长与内部优化的迭代模式匹配。此视图的 C 序迭代等效于迭代器的迭代顺序。

例如，如果使用单个数组作为输入创建迭代器，并且可以重新排列其所有轴，然后将其折叠成单个步长迭代，则这将返回一个一维数组的视图。

void NpyIter_GetReadFlags(NpyIter *iter, char *outreadflags)

填充 niter 个标志。如果可以从 op[i] 读取，则将 outreadflags[i] 设置为 1，否则设置为 0。

void NpyIter_GetWriteFlags(NpyIter *iter, char *outwriteflags)

填充 niter 个标志。如果可以写入 op[i]，则将 outwriteflags[i] 设置为 1，否则设置为 0。

迭代函数#

NpyIter_IterNext_Fn NpyIter_GetIterNext(NpyIter *iter, char **errmsg)

返回用于迭代的函数指针。此函数可能会计算函数指针的专用版本，而不是将其存储在迭代器结构中。因此，为了获得良好的性能，需要将函数指针保存在变量中，而不是在每次循环迭代中检索它。

如果发生错误，则返回 NULL。如果 errmsg 不为 NULL，则当返回 NPY_FAIL 时不会设置 Python 异常。相反，*errmsg 将被设置为错误消息。当 errmsg 不为 NULL 时，可以在不持有 Python GIL 的情况下安全地调用此函数。

典型的循环结构如下

NpyIter_IterNext_Fn iternext = NpyIter_GetIterNext(iter, NULL);
char **dataptr = NpyIter_GetDataPtrArray(iter);

do {
    /* use the addresses dataptr[0], ... dataptr[niter-1] */
} while(iternext(iter));

当指定 NPY_ITER_NO_INNER_ITERATION 时，典型的内部循环结构如下

NpyIter_IterNext_Fn iternext = NpyIter_GetIterNext(iter, NULL);
char **dataptr = NpyIter_GetDataPtrArray(iter);
npy_intp *stride = NpyIter_GetInnerStrideArray(iter);
npy_intp *size_ptr = NpyIter_GetInnerLoopSizePtr(iter), size;
npy_intp iiter, niter = NpyIter_GetNIter(iter);

do {
    size = *size_ptr;
    while (size--) {
        /* use the addresses dataptr[0], ... dataptr[niter-1] */
        for (iiter = 0; iiter < niter; ++iiter) {
            dataptr[iiter] += stride[iiter];
        }
    }
} while (iternext());

请注意，我们正在使用迭代器内部的 dataptr 数组，而不是将值复制到本地临时变量。这是可能的，因为当调用 iternext() 时，这些指针将被覆盖为新值，而不是增量更新。

如果正在使用编译时固定的缓冲区（同时使用 NPY_ITER_BUFFERED 和 NPY_ITER_NO_INNER_ITERATION 标志），则内部大小也可以用作信号。当 iternext() 返回 false 时，大小保证变为零，从而启用以下循环结构。请注意，如果使用此结构，则不应将 NPY_ITER_GROWINNER 作为标志传递，因为它会在某些情况下导致更大的大小

/* The constructor should have buffersize passed as this value */
#define FIXED_BUFFER_SIZE 1024

NpyIter_IterNext_Fn iternext = NpyIter_GetIterNext(iter, NULL);
char **dataptr = NpyIter_GetDataPtrArray(iter);
npy_intp *stride = NpyIter_GetInnerStrideArray(iter);
npy_intp *size_ptr = NpyIter_GetInnerLoopSizePtr(iter), size;
npy_intp i, iiter, niter = NpyIter_GetNIter(iter);

/* One loop with a fixed inner size */
size = *size_ptr;
while (size == FIXED_BUFFER_SIZE) {
    /*
     * This loop could be manually unrolled by a factor
     * which divides into FIXED_BUFFER_SIZE
     */
    for (i = 0; i < FIXED_BUFFER_SIZE; ++i) {
        /* use the addresses dataptr[0], ... dataptr[niter-1] */
        for (iiter = 0; iiter < niter; ++iiter) {
            dataptr[iiter] += stride[iiter];
        }
    }
    iternext();
    size = *size_ptr;
}

/* Finish-up loop with variable inner size */
if (size > 0) do {
    size = *size_ptr;
    while (size--) {
        /* use the addresses dataptr[0], ... dataptr[niter-1] */
        for (iiter = 0; iiter < niter; ++iiter) {
            dataptr[iiter] += stride[iiter];
        }
    }
} while (iternext());

NpyIter_GetCoords_Fn NpyIter_GetGetCoords(NpyIter *iter, char **errmsg)

返回一个用于获取迭代器坐标的函数指针。如果迭代器不支持坐标，则返回NULL。建议在迭代循环之前将此函数指针缓存到局部变量中。

如果发生错误，则返回 NULL。如果 errmsg 不为 NULL，则当返回 NPY_FAIL 时不会设置 Python 异常。相反，*errmsg 将被设置为错误消息。当 errmsg 不为 NULL 时，可以在不持有 Python GIL 的情况下安全地调用此函数。

char **NpyIter_GetDataPtrArray(NpyIter *iter)

此函数返回指向niter个数据指针的指针。如果未指定NPY_ITER_NO_INNER_ITERATION，则每个数据指针都指向迭代器的当前数据项。如果没有指定内部迭代，则它指向内部循环的第一个数据项。

此指针可以在迭代循环之前缓存，调用iternext不会改变它。此函数可以在不持有Python GIL的情况下安全调用。

npy_intp *NpyIter_GetIndexPtr(NpyIter *iter)

此函数返回指向正在跟踪的索引的指针，如果未跟踪任何索引则返回NULL。只有在构造期间指定了标志NPY_ITER_C_INDEX或NPY_ITER_F_INDEX之一时才能使用它。

当使用标志NPY_ITER_NO_INNER_ITERATION时，代码需要知道执行内部循环的参数。这些函数提供了这些信息。

npy_intp *NpyIter_GetInnerStrideArray(NpyIter *iter)

返回一个指向niter个步长的数组的指针，每个迭代对象一个，供内部循环使用。

此指针可以在迭代循环之前缓存，调用iternext不会改变它。此函数可以在不持有Python GIL的情况下安全调用。

npy_intp* NpyIter_GetInnerLoopSizePtr(NpyIter *iter)

返回指向内部循环应执行的迭代次数的指针。

此地址可以在迭代循环之前缓存，调用iternext不会改变它。值本身可能会在迭代期间发生变化，特别是如果启用了缓冲。此函数可以在不持有Python GIL的情况下安全调用。

void NpyIter_GetInnerFixedStrideArray(NpyIter *iter, npy_intp *out_strides)

获取一个步长数组，这些步长是固定的，或者在整个迭代过程中不会改变。对于可能发生变化的步长，在步长中放置值NPY_MAX_INTP。

一旦迭代器准备就绪（如果使用了NPY_DELAY_BUFALLOC，则在重置后），调用此函数以获取可用于选择快速内部循环函数的步长。例如，如果步长为0，则表示内部循环始终可以将其值加载到一个变量中，然后在整个循环中使用该变量；或者如果步长等于itemsize，则可以使用该操作数的连续版本。

此函数可以在不持有Python GIL的情况下安全调用。

示例#

使用迭代器的复制函数。order参数用于控制分配结果的内存布局。

如果输入是引用类型，则此函数将失败。要解决此问题，代码必须更改为专门处理可写引用，并将NPY_ITER_WRITEABLE_REFERENCES添加到标志中。

/* NOTE: This code has not been compiled/tested */
PyObject *CopyArray(PyObject *arr, NPY_ORDER order)
{
    NpyIter *iter;
    NpyIter_IterNext_Fn iternext;
    PyObject *op[2], *ret;
    npy_uint32 flags;
    npy_uint32 op_flags[2];
    npy_intp itemsize, *innersizeptr, innerstride;
    char **dataptrarray;

    /*
     * No inner iteration - inner loop is handled by CopyArray code
     */
    flags = NPY_ITER_NO_INNER_ITERATION;
    /*
     * Tell the constructor to automatically allocate the output.
     * The data type of the output will match that of the input.
     */
    op[0] = arr;
    op[1] = NULL;
    op_flags[0] = NPY_ITER_READONLY;
    op_flags[1] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE;

    /* Construct the iterator */
    iter = NpyIter_MultiNew(2, op, flags, order, NPY_NO_CASTING,
                            op_flags, NULL, 0, NULL);
    if (iter == NULL) {
        return NULL;
    }

    /*
     * Make a copy of the iternext function pointer and
     * a few other variables the inner loop needs.
     */
    iternext = NpyIter_GetIterNext(iter);
    innerstride = NpyIter_GetInnerStrideArray(iter)[0];
    itemsize = NpyIter_GetDescrArray(iter)[0]->elsize;
    /*
     * The inner loop size and data pointers may change during the
     * loop, so just cache the addresses.
     */
    innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);
    dataptrarray = NpyIter_GetDataPtrArray(iter);

    /*
     * Note that because the iterator allocated the output,
     * it matches the iteration order and is packed tightly,
     * so we don't need to check it like the input.
     */
    if (innerstride == itemsize) {
        do {
            memcpy(dataptrarray[1], dataptrarray[0],
                                    itemsize * (*innersizeptr));
        } while (iternext(iter));
    } else {
        /* Should specialize this further based on item size... */
        npy_intp i;
        do {
            npy_intp size = *innersizeptr;
            char *src = dataaddr[0], *dst = dataaddr[1];
            for(i = 0; i < size; i++, src += innerstride, dst += itemsize) {
                memcpy(dst, src, itemsize);
            }
        } while (iternext(iter));
    }

    /* Get the result from the iterator object array */
    ret = NpyIter_GetOperandArray(iter)[1];
    Py_INCREF(ret);

    if (NpyIter_Deallocate(iter) != NPY_SUCCEED) {
        Py_DECREF(ret);
        return NULL;
    }

    return ret;
}

Python lambda UFunc示例#

为了展示新的迭代器如何允许在纯Python中定义高效的类似UFunc的函数，我们演示了函数luf，它使lambda表达式像UFunc一样工作。这与numexpr库非常相似，但只需要几行代码。

首先，这是luf函数的定义。

def luf(lamdaexpr, *args, **kwargs):
    """Lambda UFunc

        e.g.
        c = luf(lambda i,j:i+j, a, b, order='K',
                            casting='safe', buffersize=8192)

        c = np.empty(...)
        luf(lambda i,j:i+j, a, b, out=c, order='K',
                            casting='safe', buffersize=8192)
    """

    nargs = len(args)
    op = args + (kwargs.get('out',None),)
    it = np.nditer(op, ['buffered','no_inner_iteration'],
            [['readonly','nbo_aligned']]*nargs +
                            [['writeonly','allocate','no_broadcast']],
            order=kwargs.get('order','K'),
            casting=kwargs.get('casting','safe'),
            buffersize=kwargs.get('buffersize',0))
    while not it.finished:
        it[-1] = lamdaexpr(*it[:-1])
        it.iternext()

    return it.operands[-1]

然后，通过使用luf而不是直接的Python表达式，我们可以从更好的缓存行为中获得一些性能提升。

In [2]: a = np.random.random((50,50,50,10))
In [3]: b = np.random.random((50,50,1,10))
In [4]: c = np.random.random((50,50,50,1))

In [5]: timeit 3*a+b-(a/c)
1 loops, best of 3: 138 ms per loop

In [6]: timeit luf(lambda a,b,c:3*a+b-(a/c), a, b, c)
10 loops, best of 3: 60.9 ms per loop

In [7]: np.all(3*a+b-(a/c) == luf(lambda a,b,c:3*a+b-(a/c), a, b, c))
Out[7]: True

Python加法示例#

迭代器已经基本编写完成并公开给Python。为了了解它的行为，让我们看看我们可以用np.add ufunc做什么。即使不更改NumPy的核心，我们也可以使用迭代器来创建更快的加法函数。

Python公开提供了两种迭代接口，一种遵循Python迭代器协议，另一种模仿C风格的do-while模式。在大多数情况下，原生Python方法更好，但如果您需要迭代器的坐标或索引，请使用C风格模式。

以下是我们如何使用Python迭代器协议编写iter_add函数。

def iter_add_py(x, y, out=None):
    addop = np.add

    it = np.nditer([x,y,out], [],
                [['readonly'],['readonly'],['writeonly','allocate']])

    for (a, b, c) in it:
        addop(a, b, c)

    return it.operands[2]

以下是相同的函数，但遵循C风格模式。

def iter_add(x, y, out=None):
    addop = np.add

    it = np.nditer([x,y,out], [],
                [['readonly'],['readonly'],['writeonly','allocate']])

    while not it.finished:
        addop(it[0], it[1], it[2])
        it.iternext()

    return it.operands[2]

此函数的一些值得注意的要点

  • 将np.add缓存为局部变量以减少命名空间查找。

  • 输入是只读的，输出是只写的，如果为None，将自动分配。

  • 使用np.add的out参数以避免额外的复制。

让我们创建一些测试变量，并对该函数以及内置的np.add进行计时。

In [1]: a = np.arange(1000000,dtype='f4').reshape(100,100,100)
In [2]: b = np.arange(10000,dtype='f4').reshape(1,100,100)
In [3]: c = np.arange(10000,dtype='f4').reshape(100,100,1)

In [4]: timeit iter_add(a, b)
1 loops, best of 3: 7.03 s per loop

In [5]: timeit np.add(a, b)
100 loops, best of 3: 6.73 ms per loop

慢了1000倍，这显然不是很好。迭代器的一个旨在帮助加快内部循环的功能是标志no_inner_iteration。这与旧迭代器的PyArray_IterAllButAxis相同，但略微更智能。让我们修改iter_add以使用此功能。

def iter_add_noinner(x, y, out=None):
    addop = np.add

    it = np.nditer([x,y,out], ['no_inner_iteration'],
                [['readonly'],['readonly'],['writeonly','allocate']])

    for (a, b, c) in it:
        addop(a, b, c)

    return it.operands[2]

性能显著提升。

In[6]: timeit iter_add_noinner(a, b)
100 loops, best of 3: 7.1 ms per loop

性能基本上与内置函数一样好！事实证明，这是因为迭代器能够合并最后两个维度，从而导致每次执行100次10000个元素的加法。如果内部循环没有变得那么大，性能就不会提高那么多。让我们使用c而不是b来了解它是如何工作的。

In[7]: timeit iter_add_noinner(a, c)
10 loops, best of 3: 76.4 ms per loop

它仍然比7秒好得多，但仍然比内置函数差10倍以上。在这里，内部循环有100个元素，它迭代了10000次。如果我们在C中编码，我们的性能已经与内置性能一样好，但在Python中，开销太大。

这让我们想到了迭代器的另一个功能，它能够为我们提供迭代内存的视图。它提供的视图结构化，以便以C顺序处理它们（就像内置的NumPy代码一样），提供与迭代器本身相同的访问顺序。有效地，我们使用迭代器来解决良好的内存访问模式，然后使用其他NumPy机制来有效地执行它。让我们再次修改iter_add。

def iter_add_itview(x, y, out=None):
    it = np.nditer([x,y,out], [],
                [['readonly'],['readonly'],['writeonly','allocate']])

    (a, b, c) = it.itviews
    np.add(a, b, c)

    return it.operands[2]

现在性能与内置函数非常接近。

In [8]: timeit iter_add_itview(a, b)
100 loops, best of 3: 6.18 ms per loop

In [9]: timeit iter_add_itview(a, c)
100 loops, best of 3: 6.69 ms per loop

现在让我们回到与新迭代器的最初动机类似的情况。以下是Fortran内存顺序而不是C内存顺序的相同计算。

In [10]: a = np.arange(1000000,dtype='f4').reshape(100,100,100).T
In [12]: b = np.arange(10000,dtype='f4').reshape(100,100,1).T
In [11]: c = np.arange(10000,dtype='f4').reshape(1,100,100).T

In [39]: timeit np.add(a, b)
10 loops, best of 3: 34.3 ms per loop

In [41]: timeit np.add(a, c)
10 loops, best of 3: 31.6 ms per loop

In [44]: timeit iter_add_itview(a, b)
100 loops, best of 3: 6.58 ms per loop

In [43]: timeit iter_add_itview(a, c)
100 loops, best of 3: 6.33 ms per loop

如您所见，内置函数的性能显着下降，但我们新编写的加法函数基本保持了相同的性能。作为最后的测试，让我们尝试几个连锁的加法。

In [4]: timeit np.add(np.add(np.add(a,b), c), a)
1 loops, best of 3: 99.5 ms per loop

In [9]: timeit iter_add_itview(iter_add_itview(iter_add_itview(a,b), c), a)
10 loops, best of 3: 29.3 ms per loop

此外，只是为了检查它是否执行了相同操作，

In [22]: np.all(
   ....: iter_add_itview(iter_add_itview(iter_add_itview(a,b), c), a) ==
   ....: np.add(np.add(np.add(a,b), c), a)
   ....: )

Out[22]: True

图像合成示例重温#

为了激发兴趣，我们举了一个对两张图像执行“叠加”合成操作的例子。现在让我们看看如何使用新的迭代器编写该函数。

为了参考，这里有一个原始函数和一些随机图像数据。

In [5]: rand1 = np.random.random(1080*1920*4).astype(np.float32)
In [6]: rand2 = np.random.random(1080*1920*4).astype(np.float32)
In [7]: image1 = rand1.reshape(1080,1920,4).swapaxes(0,1)
In [8]: image2 = rand2.reshape(1080,1920,4).swapaxes(0,1)

In [3]: def composite_over(im1, im2):
  ....:     ret = (1-im1[:,:,-1])[:,:,np.newaxis]*im2
  ....:     ret += im1
  ....:     return ret

In [4]: timeit composite_over(image1,image2)
1 loops, best of 3: 1.39 s per loop

这是重写为使用新迭代器的相同函数。请注意，添加可选输出参数有多么容易。

In [5]: def composite_over_it(im1, im2, out=None, buffersize=4096):
  ....:     it = np.nditer([im1, im1[:,:,-1], im2, out],
  ....:                     ['buffered','no_inner_iteration'],
  ....:                     [['readonly']]*3+[['writeonly','allocate']],
  ....:                     op_axes=[None,[0,1,np.newaxis],None,None],
  ....:                     buffersize=buffersize)
  ....:     while not it.finished:
  ....:         np.multiply(1-it[1], it[2], it[3])
  ....:         it[3] += it[0]
  ....:         it.iternext()
  ....:     return it.operands[3]

In [6]: timeit composite_over_it(image1, image2)
1 loops, best of 3: 197 ms per loop

即使与使用直接NumPy和C顺序数组的最佳先前尝试相比，速度也大幅提升！通过调整缓冲区大小，我们可以看到速度是如何提高的，直到我们在内部循环中达到CPU缓存的限制。

In [7]: timeit composite_over_it(image1, image2, buffersize=2**7)
1 loops, best of 3: 1.23 s per loop

In [8]: timeit composite_over_it(image1, image2, buffersize=2**8)
1 loops, best of 3: 699 ms per loop

In [9]: timeit composite_over_it(image1, image2, buffersize=2**9)
1 loops, best of 3: 418 ms per loop

In [10]: timeit composite_over_it(image1, image2, buffersize=2**10)
1 loops, best of 3: 287 ms per loop

In [11]: timeit composite_over_it(image1, image2, buffersize=2**11)
1 loops, best of 3: 225 ms per loop

In [12]: timeit composite_over_it(image1, image2, buffersize=2**12)
1 loops, best of 3: 194 ms per loop

In [13]: timeit composite_over_it(image1, image2, buffersize=2**13)
1 loops, best of 3: 180 ms per loop

In [14]: timeit composite_over_it(image1, image2, buffersize=2**14)
1 loops, best of 3: 192 ms per loop

In [15]: timeit composite_over_it(image1, image2, buffersize=2**15)
1 loops, best of 3: 280 ms per loop

In [16]: timeit composite_over_it(image1, image2, buffersize=2**16)
1 loops, best of 3: 328 ms per loop

In [17]: timeit composite_over_it(image1, image2, buffersize=2**17)
1 loops, best of 3: 345 ms per loop

最后，为了再次检查它是否正常工作，我们可以比较这两个函数。

In [18]: np.all(composite_over(image1, image2) ==
    ...:        composite_over_it(image1, image2))
Out[18]: True

使用NumExpr进行图像合成#

作为迭代器的测试，numexpr已得到增强，允许使用迭代器而不是其内部广播代码。首先，让我们使用numexpr实现合成操作。

In [22]: def composite_over_ne(im1, im2, out=None):
   ....:     ima = im1[:,:,-1][:,:,np.newaxis]
   ....:     return ne.evaluate("im1+(1-ima)*im2")

In [23]: timeit composite_over_ne(image1,image2)
1 loops, best of 3: 1.25 s per loop

这比直接的NumPy操作更好，但不是很好。切换到numexpr的迭代器版本，我们获得了比使用迭代器的纯Python函数大幅提升。请注意，这在双核机器上。

In [29]: def composite_over_ne_it(im1, im2, out=None):
   ....:     ima = im1[:,:,-1][:,:,np.newaxis]
   ....:     return ne.evaluate_iter("im1+(1-ima)*im2")

In [30]: timeit composite_over_ne_it(image1,image2)
10 loops, best of 3: 67.2 ms per loop

In [31]: ne.set_num_threads(1)
In [32]: timeit composite_over_ne_it(image1,image2)
10 loops, best of 3: 91.1 ms per loop
本页内容