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.4s 缩短至 180ms。

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

动机

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

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 通道而不是 1 个，我们会得到：

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 连续布局的输入如何导致广播步长。为了简化，示例使用 itemsize 为 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 连续时已经利用了这一点，但由于我们正在对轴进行重新排序，我们可以更通用地应用此优化。

一旦迭代步长被排序为单调递减，任何可以合并的维度都会并排。如果对于所有操作数，通过 strides[i+1] * shape[i+1] 增加与通过 strides[i] 增加相同，或者 strides[i+1] * shape[i+1] == strides[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

内部循环专门化

专门化完全由内部循环函数处理，因此此优化与其他优化无关。某些专门化已经完成，例如 reduce 操作。这个想法在 https://projects.scipy.org.cn/numpy/wiki/ProjectIdeas 中提到，“使用内部函数（SSE 指令）来加速 NumPy 中的底层循环。”

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

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

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

• 步长与数据类型的大小匹配。例如 C 或 Fortran 连续数据

  • 可以在不使用步长的情况下进行简单的循环

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

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

• 第一个输入和输出是相同的单个值（即一个归约操作）。

  • 这在现有代码中已经针对许多 UFunc 进行了专门化

上述情况通常不是互斥的，例如，当步长与数据类型大小匹配时，常量参数可以与 SSE 结合，归约也可以用 SSE 优化。

实现细节

除了内部循环专门化外，所讨论的优化显著影响了 ufunc_object.c 和用于进行广播的 PyArrayIterObject/PyArrayMultiIterObject。通常情况下，应该可以在需要的地方模拟当前行为，但我认为默认应倾向于生成和操作能提供最佳性能的内存布局。

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

拟议的 ‘order=’ 标志如下

‘C’	C 连续布局
‘F’	Fortran 连续布局
‘A’	如果输入具有 Fortran 连续布局则为 ‘F’，否则为 ‘C’（“任意连续”）
‘K’	等同于 ‘C’ 后跟某种轴排列的布局，尽可能接近输入(s)的布局（“保持布局”）

或者作为枚举

/* 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'。UFunc 此外应获得一个 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 版本可以将其完全集成。如果选择此策略，命名约定和 API 应在 1.X 发布前尽可能最终确定。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 传递错误类型时提供警告和/或错误。我们通过 incomplete struct 的 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-D 数组、2-D 数组、1-D 数组和标量。

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 类型提升规则计算得出。必须设置每个操作数的标志以便允许相应的转换，并且必须启用复制或缓冲。

如果事先知道通用数据类型，则不要使用此标志。相反，为所有操作数设置所需的 dtype。

NPY_ITER_REFS_OK

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

NPY_ITER_ZEROSIZE_OK

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

NPY_ITER_REDUCE_OK

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

当启用缓冲时，这也切换到一种特殊的缓冲模式，该模式根据需要减小循环长度，以免影响正在被归约的值。

请注意，如果您想对自动分配的输出进行归约，则必须使用 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 生成。生成 dtype 的规则与 UFunc 相同。特别值得注意的是所选 dtype 中字节序的处理。如果只有一个输入，则按原样使用输入的 dtype。否则，如果组合了多个输入 dtype，输出将以本地字节序形式呈现。

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

NPY_ITER_NO_SUBTYPE

与 NPY_ITER_ALLOCATE 一起使用，此标志禁用为输出分配数组子类型，强制其成为直接的 ndarray。

待办事项：引入函数 NpyIter_GetWrappedOutput 并删除此标志是否会更好？

NPY_ITER_NO_BROADCAST

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

NPY_ITER_WRITEABLE_REFERENCES

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

NpyIter *NpyIter_Copy(NpyIter *iter)

制作给定迭代器的副本。提供此函数主要是为了实现数据的多线程迭代。

待办事项：将其移至关于多线程迭代的章节。

多线程迭代的推荐方法是首先创建一个带有标志 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 引用。可以通过传递非 NULL 的 errmsg 参数来安全调用 Reset* 和一些其他函数，这样函数会将错误通过该参数传回，而不是设置 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 将返回一个小数，只要您选择用于多线程的块大小对于某个 n 是 NpyIter_GetReduceBlockSizeFactor(iter)*n，该操作就是安全的。

糟糕的情况是归约维度是迭代器中最外层的循环。例如，如果您有一个形状为 (3,1000,1000) 的 C 顺序数组，并且您在维度 0 上进行归约，对于 NPY_KEEPORDER 或 NPY_CORDER 迭代顺序，NpyIter_GetReduceBlockSizeFactor 将返回等于 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 数据类型 Descrs 的指针。结果指向 iter，因此调用者不会获得对 Descrs 的任何引用。

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

PyObject **NpyIter_GetOperandArray(NpyIter *iter)

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

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)

获取一组在整个迭代过程中固定不变的步长（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 风格模式。

这是我们如何编写 iter_add 函数的示例，使用了 Python 迭代器协议。

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

慢了一千倍，显然这并不理想。迭代器的一个特性是标志 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

它仍然比七秒要好得多，但仍然比内置函数慢十倍以上。在这里，内部循环有 100 个元素，并且迭代了 10000 次。如果我们用 C 语言编写，性能就已经可以达到内置函数的水平了，但在 Python 中，开销太大了。

这引出了迭代器的另一个特性：它能够提供已迭代内存的视图。它提供的视图结构使得像内置 NumPy 代码那样以 C 顺序处理它们时，可以获得与迭代器本身相同的访问顺序。实际上，我们正在利用迭代器来确定良好的内存访问模式，然后使用其他 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

重新审视图像合成示例

出于动机，我们有一个对两张图像进行“over”合成操作的示例。现在让我们看看如何用新的迭代器编写该函数。

这是其中一个原始函数供参考，以及一些随机图像数据。

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