Optimizing Iterator/UFunc Performance

作者:Mark Wiebe <mwwiebe@gmail.com>
内容类型:text / x-rst
创建:2010年11月25日

Abstract

这NEP提出用一个新的迭代器替换NumPy迭代器和多迭代器,设计更灵活,允许更多的缓存友好的数据访问。新的迭代器也包含了大部分的核心ufunc功能,使得很容易获得当前的ufunc优势在不精确适合ufunc模型的上下文中。主要优点包括:

  • 自动重排序以找到高速缓存友好的访问模式
  • 标准和可定制的广播
  • 自动类型/字节顺序/对齐转换
  • 可选缓冲以最小化转换内存使用
  • 可选输出数组,当不提供时具有自动分配
  • 自动输出或通用类型选择

这个迭代器设计的很大一部分已经实现了有希望的结果。构建开销略大(a.flat:0.5us,nditer(a):1.4us,广播(a,b):1.4us,nditer([a,b]):2.2us)例如,已经可以在迭代器中改进纯Python代码中的内置NumPy机制的性能。一个示例重写np.add,使用一些Fortran连续数组获得四倍的改进,另一个将图像合成代码从1.4s改进到180ms。

实现尝试考虑在NumPy 2.0重构器中做出的设计决策,使其未来与libndarray的集成相对简单。

Motivation

NumPy默认为从UFuncs返回C-contiguous数组。当处理结构不同的数据时,这可能导致极差的内存访问模式。一个简单的时序示例说明了这一点,将Fortran连续数组添加到一起的性能达到了8倍以上。所有的时间是使用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中的图像缓冲区的简单代码。

Image Compositing Example

对于更现实的示例,考虑图像缓冲器。图像通常以Fortran连续的顺序存储,并且颜色通道可以被视为结构化的“RGB”类型或长度为三的额外维度。所得到的内存布局既不是C也不是Fortran连续的,但是由于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默认,性能大约是7倍。

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)的单个存储器块,而不是单独的图像和α。这在图像处理应用程序中更为典型,下面是使用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之后的示例,接近文档的底部。

Improving Cache-Coherency

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

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

为了理解如何实现这一点,我们首先考虑在形状已经被规范化用于广播之后的所有输入的步幅。通过确定一组步幅是否兼容和/或模糊,我们可以确定使一致性最大化的输出存储器布局。

在广播中,首先通过前置单个维度将输入形状变换为广播形状,然后创建广播步幅,其中任何奇异维度的步幅被设置为零。

步幅也可以是负的,并且在某些情况下,这可以被归一化以适合下面的讨论。如果特定轴的所有步幅为负或零,则在适当调整基本数据指针之后可以否定该尺寸的步幅。

这里有一个例子,说明如何使用C连续布局的三个输入导致广播步幅。为了简化,这些示例使用了itemize 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)

兼容步幅 - 如果存在轴的排列,使得步幅对于集合中的每个步幅都减小,则不包括0的条目,则一组步幅是兼容的。

上面的例子满足身份置换的定义。在动机图像示例中,如果我们分离颜色和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的内存布局来解决它,但仍然与输入步幅兼容。

Output Layout Selection Algorithm

我们想要产生的输出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

Coalescing Dimensions

在许多情况下,存储器布局允许使用一维环而不是跟踪迭代器内的多个坐标。当数据是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

Inner Loop Specialization

专业化纯粹由内部循环函数处理,所以这种优化是独立于其他。一些专门化已经完成了,就像reduce操作一样。这个想法在http://projects.scipy.org/numpy/wiki/ProjectIdeas中提到,“使用内部函数(SSE指令)来加速NumPy中的低级循环。

这里有一些用于双参数函数的可能性,包括add / subtract / multiply / divide的重要情况。

  • 第一个或第二个参数是单个值(即0步幅值),并且不对输出进行别名。arr = arr + 1; arr = 1 + arr
    • 可以加载常量一次,而不是每次从内存重新加载它
  • 步幅匹配数据类型的大小。C-或Fortran连续数据
    • 可以做一个简单的循环,而不使用步幅
  • 步长与数据类型的大小匹配,并且它们都是16字节对齐(或者不同于以相同偏移量对齐的16字节)
    • 可以使用SSE一次处理多个值
  • 第一输入和输出是相同的单个值(即,缩减操作)。
    • 这已经专门用于现有代码中的许多UFunc

上述情况通常不是相互排斥的,例如,当步幅匹配数据类型大小时,常数自变量可以与SSE组合,并且还可以使用SSE优化缩减。

Implementation Details

除了内循环特化,讨论的优化显着影响ufunc_object.c和用于做广播的PyArrayIterObject / PyArrayMultiIterObject。一般来说,应该可以模拟当前需要的行为,但我相信默认应该是产生和操作内存布局,这将提供最好的性能。

为了支持新的缓存友好行为,我们为任何order=参数引入一个新的选项'K'(用于“保持”)。

建议的“order =”标志变成如下:

'C' C连续布局
'F' Fortran连续布局
'一个' 如果输入具有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;

Iterator Rewrite

基于对代码的分析,似乎重构现有的迭代对象以实现这些优化是非常困难的。另外,迭代器的一些用法需要修改内部值或标志,因此使用迭代器的代码将不得不改变。因此,我们建议创建一个新的迭代器对象,它包含现有的迭代器功能,并将其扩展以考虑优化。

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

  • 小的内存使用和少量的内存分配。
  • 简单的情况(像平数组)应该有非常少的开销。
  • 将单次和多次迭代合并为一个对象。

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

  • 以C,Fortran或“最快”(默认)顺序迭代。
  • 如果需要,跟踪C风格或Fortran风格的扁平索引(现有迭代器总是跟踪C风格索引)。这可以独立于迭代顺序来完成。
  • 如果请求,跟踪坐标(现有的迭代器需要手动更改内部迭代器标志以保证这一点)。
  • 跳过最后一个内部维度的迭代,以便可以使用内部循环进行处理。
  • 跳转到数组中的特定坐标。
  • 迭代轴的任意子集(以支持,例如,一次减少多个轴)。
  • 如果提供了NULL输入,则能够自动分配输出参数,这些输出应具有与迭代顺序匹配的内存布局,并且是order='K'支持的机制。
  • 自动复制和/或缓冲不满足类型/字节顺序/对齐要求的输入。调用者的迭代内循环应该是相同的,无论什么缓冲或复制完成。

实施注意事项:

  • 用户代码绝不能触及迭代器的内部。如果找到更高性能的实现策略,这允许将来内部存储器布局的急剧变化。
  • 使用函数指针代替宏进行迭代。这样,可以为常见情况创建专门化,例如,当ndim很小时,针对不同的标志设置,以及当迭代的数组数量很小时。此外,可以规定一个迭代模式,该模式首先使函数指针复制,以允许编译器将函数指针保存在寄存器中。
  • 动态创建内存布局,以最小化迭代器占用的高速缓存行数(对于LP64,sizeof(PyArrayIterObject)大约为2.5KB,而一个二进制操作,例如加需要三个迭代器)。
  • 将C-API对象与Python引用计数隔离,以便可以从C中自然使用。然后,Python对象成为C iterator的包装器。这类似于PEP 3118设计分离的Py_buffer和memoryview。

Proposed Iterator Memory Layout

以下结构描述了迭代器内存。所有项目都打包在一起,这意味着不同的标志值,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];
        {intp indexstride;} #if (flags&FLAGS_HASINDEX);
        /* The operand pointers and index values for this axis */
        char* ptr[niter];
        {intp index;} #if (flags&FLAGS_HASINDEX);
    }[ndim];
};

轴数据结构的数组被排序为增量更新的快速增长。如果perm是身份,这意味着它与C顺序相反。这样做使得触摸最频繁的数据项最接近结构的开始,其中公共属性是,导致增加的高速缓存一致性。它还简化了iternext调用,同时使getcoord和相关函数稍微复杂一点。

Proposed Iterator 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缓冲行为不是一个微不足道的企业。

Old -> New Iterator API Conversion

对于常规迭代器:

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

Iterator Pointer Type

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

typedef struct NpyIter_InternalOnly NpyIter;

Construction and Destruction

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_CASTINGNPY_EQUIV_CASTINGNPY_SAFE_CASTINGNPY_SAME_KIND_CASTINGNPY_UNSAFE_CASTING要允许强制转换,还必须启用复制或缓冲。

如果dtype不是NULL,那么它需要数据类型。如果允许复制,如果数据是可转换的,它将会进行临时复制。如果启用了UPDATEIFCOPY,它也将在迭代器销毁时将数据复制回另一个cast。

如果a_ndim大于零,还必须提供axes在这种情况下,axesop轴的a_ndim尺寸数组。axes中的值-1表示newaxisaxes数组中,轴不能重复。

如果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 标志, NPY_ORDER 顺序, NPY_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_CASTINGNPY_EQUIV_CASTINGNPY_SAFE_CASTINGNPY_SAME_KIND_CASTINGNPY_UNSAFE_CASTING要允许强制转换,还必须启用复制或缓冲。

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

参数oa_ndim在非零时指定将使用自定义广播进行迭代的维数。如果提供,还必须提供op_axes这两个参数让你详细控制操作数数组的轴如何匹配和迭代。op_axes中,您必须为npy_intp类型的oa_ndim -sized数组提供niter的数组。如果op_axes中的条目为NULL,则应用正常的广播规则。op_axes[j][i]中存储op[j]的有效轴,或者表示newaxis的-1。在每个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_INDEXNPY_ITER_F_INDEX

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

NPY_ITER_COORDS

导致迭代器跟踪数组坐标。这防止迭代器合并轴以产生更大的内循环。

NPY_ITER_NO_INNER_ITERATION

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

此标记与NPY_ITER_C_INDEXNPY_ITER_F_INDEXNPY_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_COPYNPY_ITER_UPDATEIFCOPY标志,因为它们将覆盖缓冲。缓冲对于使用迭代器的Python代码特别有用,允许同时使用更大的数据块来摊销Python解释器的开销。

如果与NPY_ITER_NO_INNER_ITERATION一起使用,调用者的内部循环可能获得比没有缓冲的情况下更大的块,因为如何布置步幅。

注意,如果操作数被赋予标志NPY_ITER_COPYNPY_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 &lt; niter

NPY_ITER_READWRITENPY_ITER_READONLYNPY_ITER_WRITEONLY

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

NPY_ITER_COPY

如果不满足构造函数标志和参数指定的数据类型或对齐要求,则允许创建op[i]的副本。

NPY_ITER_UPDATEIFCOPY

触发器NPY_ITER_COPY,当数组操作数被标记为写入并被复制时,当迭代器被销毁时,使副本中的数据被复制回op[i]

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

NPY_ITER_NBONPY_ITER_ALIGNEDNPY_ITER_CONTIG

使迭代器为本地字节顺序提供op[i]的数据,根据dtype要求,连续或任何组合进行对齐。

默认情况下,迭代器产生指向提供的数组的指针,可以是对齐的或不对齐的,以及任何字节顺序。如果未启用复制或缓冲,并且操作数数据不满足约束,则会出现错误。

连续约束仅应用于内循环,连续内循环可以具有任意的指针改变。

如果请求的数据类型是非本地字节顺序,NBO标志覆盖它,并且请求的数据类型被转换为本地字节顺序。

NPY_ITER_ALLOCATE

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

op[i]为NULL时,请求的数据类型op_dtypes[i]也可以为NULL,在这种情况下,它从数组的dtypes自动生成被标记为可读。生成dtype的规则与UFuncs相同。特别注意的是处理所选dtype中的字节顺序。如果只有一个输入,则使用输入的dtype。否则,如果多个输入dtypes组合在一起,输出将按原生字节顺序。

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

NPY_ITER_NO_SUBTYPE

对于与NPY_ITER_ALLOCATE一起使用,此标志禁用为输出分配数组子类型,强制其为直线型。

TODO:也许最好引入一个函数NpyIter_GetWrappedOutput并删除这个标志?

NPY_ITER_NO_BROADCAST

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

NPY_ITER_WRITEABLE_REFERENCES

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

NpyIter * NpyIter_Copy(NpyIter * iter)

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

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

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

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

int NpyIter_UpdateIter(NpyIter * iter, npy_intp i, npy_uint32 op_flags, NPY_CASTING cast, PyArray_Descr * dtype) UNIMPLEMENTED

更新迭代器中的第i个操作数,以可能具有新的数据类型或更具限制性的标志属性。这样做的一个用例是允许自动分配基于标准NumPy类型促销规则确定输出数据类型,然后使用此函数在处理期间将输入和可能的自动输出转换为不同的数据类型。

仅当NPY_ITER_COORDS作为标志传递给迭代器时,才能执行此操作。如果不需要坐标,则一旦不再需要对NpyIter_UpdateIter的调用,则调用函数NpyIter_RemoveCoords()

如果第i个操作数已经被复制,则抛出错误。为了避免这种情况,保留所有标志,除了以后调用NpyIter_UpdateIter的任何操作数的读/写指示符。

The flags that may be passed in op_flags are NPY_ITER_COPY, NPY_ITER_UPDATEIFCOPY, NPY_ITER_NBO, NPY_ITER_ALIGNED, NPY_ITER_CONTIG.

int NpyIter_RemoveAxis(NpyIter * iter, npy_intp axis) t0>

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

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

警告:此函数可能会更改迭代器的内部存储器布局。任何缓存的函数或来自迭代器的指针必须再次检索!

返回NPY_SUCCEEDNPY_FAIL

int NpyIter_RemoveCoords(NpyIter * iter)

如果迭代器有坐标,这条带支持它们,并且做进一步的迭代器优化,如果不需要坐标的话。此函数还将迭代器重置为其初始状态。

警告:此函数可能会更改迭代器的内部存储器布局。必须再次检索来自迭代器的任何缓存函数或指针!

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

返回NPY_SUCCEEDNPY_FAIL

int NpyIter_RemoveInnerLoop(NpyIter * iter)

如果使用UpdateIter / RemoveCoords,您可能需要指定标志NPY_ITER_NO_INNER_ITERATION此标志不能与NPY_ITER_COORDS一起使用,因此提供此功能以在调用NpyIter_RemoveCoords之后启用此功能。此函数还将迭代器重置为其初始状态。

警告:此函数更改迭代器的内部逻辑。必须再次检索来自迭代器的任何缓存函数或指针!

返回NPY_SUCCEEDNPY_FAIL

int NpyIter_Deallocate(NpyIter * iter)

取消分配迭代器对象。这另外释放所有副本,触发UPDATEIFCOPY行为在必要时。

返回NPY_SUCCEEDNPY_FAIL

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

在迭代范围的开始处将迭代器重置为其初始状态。

返回NPY_SUCCEEDNPY_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_SUCCEEDNPY_FAIL如果errmsg为非NULL,则在返回NPY_FAIL时,不会设置Python异常。相反,* errmsg设置为错误消息。当errmsg是非NULL时,可以安全地调用该函数,而不用保持Python GIL。

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

将迭代器重置为其初始状态,但使用baseptrs中的值替代正在迭代的数组的指针。此函数旨在与op_axes参数一起通过嵌套迭代代码与两个或多个迭代器一起使用。

返回NPY_SUCCEEDNPY_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) / t0>

调整迭代器以指向coords指向的ndim坐标。如果坐标未被跟踪,坐标超出边界或内循环迭代被禁用,则返回错误。

返回NPY_SUCCEEDNPY_FAIL

int NpyIter_GotoIndex(NpyIter * iter, npy_intp index) t0>

调整迭代器以指向指定的index如果迭代器被构造为具有标志NPY_ITER_C_INDEX,则index是C阶索引,并且如果迭代器被构造为具有标志NPY_ITER_F_INDEX index是Fortran顺序索引。如果没有正在跟踪的索引,索引超出边界或禁用内循环迭代,则返回错误。

返回NPY_SUCCEEDNPY_FAIL

npy_intp NpyIter_GetIterSize(NpyIter * iter)

返回要迭代的元素数。这是所有尺寸的形状的产物。

npy_intp NpyIter_GetReduceBlockSizeFactor(NpyIter * iter) UNIMPLEMENTED

这提供了一个必须划分为用于安全多线程减少的范围迭代的块大小的因素。如果迭代器没有减少,它返回1。

当使用范围迭代到多线程减少时,有两种可能的方式来执行减少:

如果有一个大的减少到一个小的输出,使一个临时数组初始化为每个线程的缩减单位,然后让每个线程缩减为其临时。当完成时,将临时组合在一起。您可以通过观察NpyIter_GetReduceBlockSizeFactor返回一个大的值,例如NpyIter_GetIterSize的一半或三分之一来检测这种情况。你还应该检查输出是否只是为了确保。

如果对于大输出有很多小的减少,并且减少维度是内部维度,则NpyIter_GetReduceBlockSizeFactor将返回一个小数,并且只要多线程所选的块大小为NpyIter_GetReduceBlockSizeFactor(iter)*n,对于某些n,操作将是安全的。

坏的情况是当缩减维度是迭代器中最外面的循环时。例如,如果你有一个具有形状(3,1000,1000)的C阶数组,并且在维度0上减少,则NpyIter_GetReduceBlockSizeFactor将返回等于NpyIter_GetIterSize对于NPY_KEEPORDERNPY_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) t0>

调整迭代器以指向指定的iterindexIterIndex是与迭代器的迭代次序匹配的索引。如果iterindex超出边界,启用缓冲或禁用内循环迭代,则返回错误。

返回NPY_SUCCEEDNPY_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_INDEXNPY_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) / t0>

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

返回NPY_SUCCEEDNPY_FAIL

PyArray_Descr ** NpyIter_GetDescrArray(NpyIter * iter)

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

这个指针可以在迭代循环之前缓存,调用iternext将不会改变它。

PyObject ** NpyIter_GetOperandArray(NpyIter * iter)

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

PyObject * NpyIter_GetIterView(NpyIter * iter, npy_intp i) / t0>

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

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

void NpyIter_GetReadFlags(NpyIter * iter, char * outreadflags) / t0>

填充niter标志。如果op[i]可以读取,则将outreadflags[i]设置为1,如果不能读取,则设置为0。

void NpyIter_GetWriteFlags(NpyIter * iter, char * outwriteflags) / t0>

填充niter标志。如果op[i]可以写入,则将outwriteflags[i]设置为1,如果不能,则将其设置为0。

Functions For Iteration

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_BUFFEREDNPY_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_INDEXNPY_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) / t0>

获取固定的或在整个迭代过程中不会更改的步长数组。对于可能改变的步幅,值NPY_MAX_INTP被放置在步幅。

一旦迭代器准备好迭代(在使用NPY_DELAY_BUFALLOC的复位之后),调用它以获得可用于选择快速内循环函数的步长。例如,如果stride为0,那意味着内循环总是可以将其值加载到变量中一次,然后在整个循环中使用该变量,或者如果stride等于itemsize,则可以使用该操作数的连续版本。

这个函数可以安全地调用而无需持有Python GIL。

Examples

使用迭代器的复制函数。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 Example

为了显示新的迭代器如何允许在纯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 Addition Example

迭代器主要是编写并暴露给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]

关于这个功能的一些值得注意的点:

  • 输入为只读,输出为writeonly,如果为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个元素100。如果内部循环没有变得那么大,性能不会显着改善。让我们使用c而不是b来了解它是如何工作的。

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

它仍然比7秒钟好多了,但仍然比内置函数差十多倍。这里,内循环有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

正如你可以看到,内置函数的性能明显下降,但我们新编写的add函数保持基本相同的性能。作为最后一个测试,让我们尝试几个添加链接在一起。

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

Image Compositing Example Revisited

为了动机,我们有一个例子,对两个图像做了一个'over'复合操作。现在让我们看看如何使用新的迭代器编写函数。

这里是原始功能之一,用于参考和一些随机图像数据。

In [5]: rand1 = np.random.random_sample(1080*1920*4).astype(np.float32)
In [6]: rand2 = np.random.random_sample(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

Image Compositing With 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