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-contiguous 数组。这在处理结构不同的数据时可能导致极差的内存访问模式。一个简单的计时示例说明了这一点，将 Fortran-contiguous 数组相加会导致超过八倍的性能损失。所有计时均使用 NumPy 2.0dev（2010 年 11 月 22 日）在 Athlon 64 X2 4200+、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-contiguous 顺序存储，颜色通道可以被视为结构化的“RGB”类型或长度为 3 的额外维度。由此产生的内存布局既不是 C-contiguous 也不是 Fortran-contiguous，但由于 ndarray 的灵活性，可以直接在 NumPy 中使用。这似乎是理想的，因为它使内存布局与典型的 C 或 C++ 图像代码兼容，同时在 Python 中提供了自然的访问方式。获取像素 (x,y) 的颜色只需“image[x,y]”。

这种布局在 NumPy 中的性能非常差。以下代码创建了两个黑色图像，并对其进行“over”合成操作。

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-contiguous 默认值，性能大约会提高七倍。

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 个值而不是 1 个值的 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-contiguous 布局下的情况。

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-contiguous 布局。

要理解如何实现这一点，我们首先考虑所有输入在形状已标准化以进行广播后的步幅。通过确定一组步幅是否兼容和/或模糊，我们可以确定一个最大化一致性的输出内存布局。

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

步幅也可以是负数，在某些情况下可以将其规范化以适应以下讨论。如果某个特定轴的所有步幅都是负数或零，则在适当调整基本数据指针后，该维度的步幅可以取反。

以下是三个具有 C-contiguous 布局的输入如何导致广播步幅的示例。为简化起见，示例中使用 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-contiguous 但仍与输入步幅兼容的内存布局来解决它。

输出布局选择算法#

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

一致/明确的步幅

单一的一致布局

一致/模糊的步幅

最接近 C-contiguous 的一致布局

不一致的步幅

C-contiguous

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

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-contiguous 时已经利用了这一点，但由于我们正在重新排序轴，因此可以更普遍地应用此优化。

一旦迭代步幅已排序为单调递减，则可以合并的任何维度都将并排放置。如果对于所有操作数，按 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-contiguous 数据。

    • 可以执行一个简单的循环，而无需使用步幅。

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

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

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

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

上述情况通常并非互斥，例如，当步幅与数据类型大小匹配时，常数参数可以与 SSE 结合使用，并且归约也可以通过 SSE 进行优化。

实现细节#

除了内循环专门化，所讨论的优化显著影响了 ufunc_object.c 和用于广播的 PyArrayIterObject/PyArrayMultiIterObject。总的来说，应该可以模拟当前的期望行为，但我相信默认值应该是生成和操作能够提供最佳性能的内存布局。

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

提议的“order=”标志如下：

‘C’

C-contiguous 布局

‘F’

Fortran-contiguous 布局

‘A’

如果输入具有 Fortran-contiguous 布局，则为 'F'，否则为 'C'（“Any Contiguous”）。

‘K’

与 'C' 等效的布局，后跟轴的某种排列，尽可能接近输入布局（“Keep Layout”）。

或作为枚举

/* 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。此迭代器的一种可能发布策略是发布一个添加了迭代器但 NumPy 代码未使用的 1.X (1.6?) 版本。然后，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 是一个 a_ndim 大小、包含 op 的轴的数组。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 类型提升规则计算。必须设置每个操作数的标志，以便允许适当的类型转换，并且必须启用复制或缓冲。

如果通用数据类型已知，请勿使用此标志。而是为所有操作数设置请求的数据类型。

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_NO_INNER_ITERATION 一起使用，且必须启用 NPY_ITER_BUFFERED。这是因为在没有缓冲的情况下，内循环的大小始终是最内层迭代维度的大小，允许其被切分将需要特殊处理，从而使其更像缓冲版本。

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] 不满足构造函数标志和参数指定的数据类型或对齐要求，则允许创建 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 的规则与 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* 和其他一些函数可以通过传递非 NULL 的 errmsg 参数来安全调用，这样函数将通过它返回错误，而不是设置 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 的情况下安全调用该函数。

TODO：将以下内容移至嵌套迭代器的特殊部分。

创建嵌套迭代器需要一些注意。所有迭代器操作数必须完全匹配，否则对 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 将返回一个等于 NpyIter_GetIterSize 的大小，用于 NPY_KEEPORDER 或 NPY_CORDER 迭代顺序。虽然这不利于 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)

如果迭代器已缓冲，则返回 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 个数据类型描述符的指针。结果指向 iter 内部，因此调用者不会获得对描述符的任何引用。

此指针可以在迭代循环之前缓存，调用 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

慢了一千倍，这显然不太好。迭代器的一个特性，旨在帮助加速内循环，是标志 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 中有太多的开销。

这引出了迭代器的另一个特性，即它能够为我们提供迭代内存的视图。它给出的视图的结构是这样的：以 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

图像合成示例回顾#

作为动机，我们有一个对两张图像执行“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