diff options
| -rw-r--r-- | doc/neps/new-iterator-ufunc.rst | 1462 |
1 files changed, 1462 insertions, 0 deletions
diff --git a/doc/neps/new-iterator-ufunc.rst b/doc/neps/new-iterator-ufunc.rst new file mode 100644 index 000000000..59de93acd --- /dev/null +++ b/doc/neps/new-iterator-ufunc.rst @@ -0,0 +1,1462 @@ +:Title: Optimizing Iterator/UFunc Performance +:Author: Mark Wiebe <mwwiebe@gmail.com> +:Content-Type: text/x-rst +:Created: 25-Nov-2010 + +******** +Abstract +******** + +This NEP proposes to replace the NumPy iterator and multi-iterator +with a single new iterator, designed to be more flexible and allow for +more cache-friendly data access. The new iterator also subsumes much +of the core ufunc functionality, making it easy to get the current +ufunc benefits in contexts which don't precisely fit the ufunc mold. +Key benefits include: + +* automatic reordering to find a cache-friendly access pattern +* standard and customizable broadcasting +* automatic type/byte-order/alignment conversions +* optional buffering to minimize conversion memory usage +* optional output arrays +* automatic output or common type selection + +A large fraction of this iterator design has already been implemented with +promising results. Construction overhead is slightly greater (a.flat: +0.5 us, newiter(a): 1.4 us and broadcast(a,b): 1.4 us, newiter([a,b]): +2.2 us), but, as shown in an example, it is already possible to improve +on the performance of the built-in NumPy mechanisms in pure Python code +together with the iterator. One example rewrites np.add, getting a +four times improvement with some Fortran-contiguous arrays, and +another improves image compositing code from 1.4s to 180ms. + +The implementation attempts to take into account +the design decisions made in the NumPy 2.0 refactor, to make its future +integration into libndarray relatively simple. + +********** +Motivation +********** + +NumPy defaults to returning C-contiguous arrays from UFuncs. This can +result in extremely poor memory access patterns when dealing with data +that is structured differently. A simple timing example illustrates +this with a more than eight times performance hit from adding +Fortran-contiguous arrays together. All timings are done using Numpy +2.0dev (Nov 22, 2010) on an Athlon 64 X2 4200+, with a 64-bit OS.:: + + 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 + +In this case, it is simple to recover the performance by switching to +a view of the memory, adding, then reshaping back. To further examine +the problem and see how it isn’t always as trivial to work around, +let’s consider simple code for working with image buffers in NumPy. + +Image Compositing Example +========================= + +For a more realistic example, consider an image buffer. Images are +generally stored in a Fortran-contiguous order, and the colour +channel can be treated as either a structured 'RGB' type or an extra +dimension of length three. The resulting memory layout is neither C- +nor Fortran-contiguous, but is easy to work with directly in NumPy, +because of the flexibility of the ndarray. This appears ideal, because +it makes the memory layout compatible with typical C or C++ image code, +while simultaneously giving natural access in Python. Getting the color +of pixel (x,y) is just ‘image[x,y]’. + +The performance of this layout in NumPy turns out to be very poor. +Here is code which creates two black images, and does an ‘over’ +compositing operation on them.:: + + 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 + +If we give up the convenient layout, and use the C-contiguous default, +the performance is about seven times better.:: + + 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 + +But this is not all, since it turns out that broadcasting the alpha +channel is exacting a performance price as well. If we use an alpha +channel with 3 values instead of one, we get:: + + 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 + +For a final comparison, let’s see how it performs when we use +one-dimensional arrays to ensure just a single loop does the +calculation.:: + + 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 + +To get a reference performance number, I implemented this simple operation +straightforwardly in C (careful to use the same compile options as NumPy). +If I emulated the memory allocation and layout of the Python code, the +performance was roughly 0.3 seconds, very much in line with NumPy’s +performance. Combining the operations into one pass reduced the time +to roughly 0.15 seconds. + +A slight variation of this example is to use a single memory block +with four channels (1920,1080,4) instead of separate image and alpha. +This is more typical in image processing applications, and here’s how +that looks with a C-contiguous layout.:: + + 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 + +To see the improvements that implementation of the new iterator as +proposed can produce, go to the example continued after the +proposed API, near the bottom of the document. + +************************* +Improving Cache-Coherency +************************* + +In order to get the best performance from UFunc calls, the pattern of +memory reads should be as regular as possible. Modern CPUs attempt to +predict the memory read/write pattern and fill the cache ahead of time. +The most predictable pattern is for all the inputs and outputs to be +sequentially processed in the same order. + +I propose that by default, the memory layout of the UFunc outputs be as +close to that of the inputs as possible. Whenever there is an ambiguity +or a mismatch, it defaults to a C-contiguous layout. + +To understand how to accomplish this, we first consider the strides of +all the inputs after the shapes have been normalized for broadcasting. +By determining whether a set of strides are compatible and/or ambiguous, +we can determine an output memory layout which maximizes coherency. + +In broadcasting, the input shapes are first transformed to broadcast +shapes by prepending singular dimensions, then the broadcast strides +are created, where any singular dimension’s stride is set to zero. + +Strides may be negative as well, and in certain cases this can be +normalized to fit the following discussion. If all the strides for a +particular axis are negative or zero, the strides for that dimension +can be negated after adjusting the base data pointers appropriately. + +Here's an example of how three inputs with C-contiguous layouts result in +broadcast strides. To simplify things, the examples use an itemsize of 1. + +================== ======== ======= ======= +Input shapes: (5,3,7) (5,3,1) (1,7) +Broadcast shapes: (5,3,7) (5,3,1) (1,1,7) +Broadcast strides: (21,7,1) (3,1,0) (0,0,1) +================== ======== ======= ======= + +*Compatible Strides* - A set of strides are compatible if there exists +a permutation of the axes such that the strides are decreasing for every +stride in the set, excluding entries that are zero. + +The example above satisfies the definition with the identity permutation. +In the motivation image example, the strides are slightly different if +we separate the colour and alpha information or not. The permutation +which demonstrates compatibility here is the transposition (0,1). + +============================= ===================== ===================== +Input/Broadcast shapes: Image (1920, 1080, 3) Alpha (1920, 1080, 1) +Broadcast strides (separate): (3,5760,1) (1,1920,0) +Broadcast strides (together): (4,7680,1) (4,7680,0) +============================= ===================== ===================== + +*Ambiguous Strides* - A set of compatible strides are ambiguous if +more than one permutation of the axes exists such that the strides are +decreasing for every stride in the set, excluding entries that are zero. + +This typically occurs when every axis has a 0-stride somewhere in the +set of strides. The simplest example is in two dimensions, as follows. + +================== ===== ===== +Broadcast shapes: (1,3) (5,1) +Broadcast strides: (0,1) (1,0) +================== ===== ===== + +There may, however, be unambiguous compatible strides without a single +input forcing the entire layout, as in this example: + +================== ======= ======= +Broadcast shapes: (1,3,4) (5,3,1) +Broadcast strides: (0,4,1) (3,1,0) +================== ======= ======= + +In the face of ambiguity, we have a choice to either completely throw away +the fact that the strides are compatible, or try to resolve the ambiguity +by adding an additional constraint. I think the appropriate choice +is to resolve it by picking the memory layout closest to C-contiguous, +but still compatible with the input strides. + +Output Layout Selection Algorithm +================================= + +The output ndarray memory layout we would like to produce is as follows: + +=============================== ============================================= +Consistent/Unambiguous strides: The single consistent layout +Consistent/Ambiguous strides: The consistent layout closest to C-contiguous +Inconsistent strides: C-contiguous +=============================== ============================================= + +Here is pseudo-code for an algorithm to compute the permutation for the +output layout.:: + + 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 +********************* + +In many cases, the memory layout allows for the use of a one-dimensional +loop instead of tracking multiple coordinates within the iterator. +The existing code already exploits this when the data is C-contiguous, +but since we're reordering the axes, we can apply this optimization +more generally. + +Once the iteration strides have been sorted to be monotonically +decreasing, any dimensions which could be coalesced are side by side. +If for all the operands, incrementing by strides[i+1] shape[i+1] times +is the same as incrementing by strides[i], or strides[i+1]*shape[i+1] == +strides[i], dimensions i and i+1 can be coalesced into a single dimension. + +Here is pseudo-code for coalescing.:: + + # 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 +************************* + +Specialization is handled purely by the inner loop function, so this +optimization is independent of the others. Some specialization is +already done, like for the reduce operation. The idea is mentioned in +http://projects.scipy.org/numpy/wiki/ProjectIdeas, “use intrinsics +(SSE-instructions) to speed up low-level loops in NumPy.” + +Here are some possibilities for two-argument functions, +covering the important cases of add/subtract/multiply/divide. + +* The first or second argument is a single value (i.e. a 0 stride + value) and does not alias the output. arr = arr + 1; arr = 1 + arr + + * Can load the constant once instead of reloading it from memory every time + +* The strides match the size of the data type. C- or + Fortran-contiguous data, for example + + * Can do a simple loop without using strides + +* The strides match the size of the data type, and they are + both 16-byte aligned (or differ from 16-byte aligned by the same offset) + + * Can use SSE to process multiple values at once + +* The first input and the output are the same single value + (i.e. a reduction operation). + + * This is already specialized for many UFuncs in the existing code + +The above cases are not generally mutually exclusive, for example a +constant argument may be combined with SSE when the strides match the +data type size, and reductions can be optimized with SSE as well. + +********************** +Implementation Details +********************** + +Except for inner loop specialization, the discussed +optimizations significantly affect ufunc_object.c and the +PyArrayIterObject/PyArrayMultiIterObject used to do the broadcasting. +In general, it should be possible to emulate the current behavior where it +is desired, but I believe the default should be to produce and manipulate +memory layouts which will give the best performance. + +To support the new cache-friendly behavior, we introduce a new +option ‘K’ (for “keep”) for any ``order=`` parameter. + +The proposed ‘order=’ flags become as follows: + +=== ===================================================================================== +‘C’ C-contiguous layout +‘F’ Fortran-contiguous layout +‘A’ ‘F’ if the input(s) have a Fortran-contiguous layout, ‘C’ otherwise (“Any Contiguous”) +‘K’ a layout equivalent to ‘C’ followed by some permutation of the axes, as close to the layout of the input(s) as possible (“Keep Layout”) +=== ===================================================================================== + +Perhaps a good strategy is to first implement the capabilities discussed +here without changing the defaults. Once they are implemented and +well-tested, the defaults can change from ``order='C'`` to ``order='K'`` +everywhere appropriate. UFuncs additionally should gain an ``order=`` +parameter to control the layout of their output(s). + +Iterator Rewrite +================ + +Based on an analysis of the code, it appears that refactoring the existing +iteration objects to implement these optimizations is prohibitively +difficult. Additionally, some usage of the iterator requires modifying +internal values or flags, so code using the iterator would have to +change anyway. Thus we propose creating a new iterator object which +subsumes the existing iterator functionality and expands it to account +for the optimizations. + +High level goals for the replacement iterator include: + +* Small memory usage and a low number of memory allocations. +* Simple cases (like flat arrays) should have very little overhead. +* Combine single and multiple iteration into one object. + +Capabilities that should be provided to user code: + +* Iterate in C, Fortran, or “Fastest” (default) order. +* Track a C-style or Fortran-style flat index if requested + (existing iterator always tracks a C-style index). This can be done + independently of the iteration order. +* Track the coordinates if requested (the existing iterator requires + manually changing an internal iterator flag to guarantee this). +* Skip iteration of the last internal dimension so that it can be + processed with an inner loop. +* Jump to a specific coordinate in the array. +* Iterate an arbitrary subset of axes (to support, for example, reduce + with multiple axes at once). +* Ability to automatically allocate output parameters if a NULL input + is provided, These outputs should have a memory layout matching + the iteration order, and are the mechanism for the ``order='K'`` + support. +* Automatic copying and/or buffering of inputs which do not satisfy + type/byte-order/alignment requirements. The caller's iteration inner + loop should be the same no matter what buffering or copying is done. + +Notes for implementation: + +* User code must never touch the inside of the iterator. This allows + for drastic changes of the internal memory layout in the future, if + higher-performance implementation strategies are found. +* Use a function pointer instead of a macro for iteration. + This way, specializations can be created for the common cases, + like when ndim is small, for different flag settings, and when the + number of arrays iterated is small. Also, an iteration pattern + can be prescribed that makes a copy of the function pointer first + to allow the compiler to keep the function pointer + in a register. +* Dynamically create the memory layout, to minimize the number of + cache lines taken up by the iterator (for LP64, + sizeof(PyArrayIterObject) is about 2.5KB, and a binary operation + like plus needs three of these for the Multi-Iterator). +* Isolate the C-API object from Python reference counting, so that + it can be used naturally from C. The Python object then becomes + a wrapper around the C iterator. This is analogous to the + PEP 3118 design separation of Py_buffer and memoryview. + +Proposed Iterator Memory Layout +=============================== + +The following struct describes the iterator memory. All items +are packed together, which means that different values of the flags, +ndim, and niter will produce slightly different layouts. :: + + 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]; + }; + +The array of axis_data structs is ordered to be in increasing rapidity +of increment updates. If the ``perm`` is the identity, this means it’s +reversed from the C-order. This is done so data items touched +most often are closest to the beginning of the struct, where the +common properties are, resulting in increased cache coherency. +It also simplifies the iternext call, while making getcoord and +related functions slightly more complicated. + +Proposed Iterator API +===================== + +The existing iterator API includes functions like PyArrayIter_Check, +PyArray_Iter* and PyArray_ITER_*. The multi-iterator array includes +PyArray_MultiIter*, PyArray_Broadcast, and PyArray_RemoveSmallest. The +new iterator design replaces all of this functionality with a single object +and associated API. One goal of the new API is that all uses of the +existing iterator should be replaceable with the new iterator without +significant effort. + +The C-API naming convention chosen is based on the one in the numpy-refactor +branch, where libndarray has the array named ``NpyArray`` and functions +named ``NpyArray_*``. The iterator is named ``NpyIter`` and functions are +named ``NpyIter_*``. + +The Python exposure has the iterator named ``np.newiter``. One possible +release strategy for this iterator would be to release a 1.X (1.6?) version +with the iterator added, but not used by the NumPy code. Then, 2.0 can +be release with it fully integrated. If this strategy is chosen, the +naming convention and API should be finalized as much as possible before +the 1.X release. I would suggest the name ``np.iter`` within Python, +as it is currently unused. + +In addition to the performance goals set out for the new iterator, +it appears the API can be refactored to better support some common +NumPy programming idioms. + +By moving some functionality currently in the UFunc code into the +iterator, it should make it easier for extension code which wants +to emulate UFunc behavior in cases which don't quite fit the +UFunc paradigm. In particular, emulating the UFunc buffering behavior +is not a trivial enterprise. + +Old -> New Iterator API Conversion +---------------------------------- + +For the regular iterator: + +=============================== ============================================= +``PyArray_IterNew`` ``NpyIter_New`` +``PyArray_IterAllButAxis`` ``NpyIter_New`` + ``axes`` parameter **or** + Iterator flag ``NPY_ITER_NO_INNER_ITERATION`` +``PyArray_BroadcastToShape`` **NOT SUPPORTED** (but could be, if needed) +``PyArrayIter_Check`` Will need to add this in Python exposure +``PyArray_ITER_RESET`` ``NpyIter_Reset`` +``PyArray_ITER_NEXT`` Function pointer from ``NpyIter_GetIterNext`` +``PyArray_ITER_DATA`` ``NpyIter_GetDataPtrArray`` +``PyArray_ITER_GOTO`` ``NpyIter_GotoCoords`` +``PyArray_ITER_GOTO1D`` ``NpyIter_GotoIndex`` +``PyArray_ITER_NOTDONE`` Return value of ``iternext`` function pointer +=============================== ============================================= + +For the multi-iterator: + +=============================== ============================================= +``PyArray_MultiIterNew`` ``NpyIter_MultiNew`` +``PyArray_MultiIter_RESET`` ``NpyIter_Reset`` +``PyArray_MultiIter_NEXT`` Function pointer from ``NpyIter_GetIterNext`` +``PyArray_MultiIter_DATA`` ``NpyIter_GetDataPtrArray`` +``PyArray_MultiIter_NEXTi`` **NOT SUPPORTED** (always lock-step iteration) +``PyArray_MultiIter_GOTO`` ``NpyIter_GotoCoords`` +``PyArray_MultiIter_GOTO1D`` ``NpyIter_GotoIndex`` +``PyArray_MultiIter_NOTDONE`` Return value of ``iternext`` function pointer +``PyArray_Broadcast`` Handled by ``NpyIter_MultiNew`` +``PyArray_RemoveSmallest`` Iterator flag ``NPY_ITER_NO_INNER_ITERATION`` +=============================== ============================================= + +For other API calls: + +=============================== ============================================= +``PyArray_ConvertToCommonType`` Iterator flag ``NPY_ITER_COMMON_DATA_TYPE`` +=============================== ============================================= + + +Iterator Pointer Type +--------------------- + +The iterator structure is internally generated, but a type is still needed +to provide warnings and/or errors when the wrong type is passed to +the API. We do this with a typedef of an incomplete struct + +``typedef struct NpyIter_InternalOnly NpyIter;`` + + +Construction and Destruction +---------------------------- + +``NpyIter* NpyIter_New(PyArrayObject* op, npy_uint32 flags, NPY_ORDER order, PyArray_Descr* dtype, npy_intp a_ndim, npy_intp *axes, npy_intp buffersize)`` + + Creates an iterator for the given numpy array object ``op``. + + Flags that may be passed in ``flags`` are any combination + of the global and per-operand flags documented in + ``NpyIter_MultiNew``, except for ``NPY_ITER_ALLOCATE``. + + Any of the ``NPY_ORDER`` enum values may be passed to ``order``. For + efficient iteration, ``NPY_KEEPORDER`` is the best option, and the other + orders enforce the particular iteration pattern. + + If ``dtype`` isn't ``NULL``, then it requires that data type. + If copying is allowed, it will make a temporary copy if the data + is castable. If ``UPDATEIFCOPY`` is enabled, it will also copy + the data back with another cast upon iterator destruction. + + If ``a_ndim`` is greater than zero, ``axes`` must also be provided. + In this case, ``axes`` is an ``a_ndim``-sized array of ``op``'s axes. + A value of -1 in ``axes`` means ``newaxis``. Within the ``axes`` + array, axes may not be repeated. + + If ``buffersize`` is zero, a default buffer size is used, + otherwise it specifies how big of a buffer to use. Buffers + which are powers of 2 such as 512 or 1024 are recommended. + + Returns NULL if there is an error, otherwise returns the allocated + iterator. + + To make an iterator similar to the old iterator, this should work.:: + + iter = NpyIter_New(op, NPY_ITER_READWRITE, + NPY_CORDER, NULL, 0, NULL); + + If you want to edit an array with aligned ``double`` code, + but the order doesn't matter, you would use this.:: + + dtype = PyArray_DescrFromType(NPY_DOUBLE); + iter = NpyIter_New(op, NPY_ITER_READWRITE | + NPY_ITER_BUFFERED | + NPY_ITER_SAME_KIND_CASTS | + NPY_ITER_NBO_ALIGNED, + NPY_KEEPORDER, + dtype, 0, NULL); + Py_DECREF(dtype); + +``NpyIter* NpyIter_MultiNew(npy_intp niter, PyArrayObject** op, npy_uint32 flags, NPY_ORDER order, npy_uint32 *op_flags, PyArray_Descr** op_dtypes, npy_intp oa_ndim, npy_intp **op_axes, npy_intp buffersize)`` + + Creates an iterator for broadcasting the ``niter`` array objects provided + in ``op``. + + For normal usage, use 0 for ``oa_ndim`` and NULL for ``op_axes``. + See below for a description of these parameters, which allow for + custom manual broadcasting as well as reordering and leaving out axes. + + Any of the ``NPY_ORDER`` enum values may be passed to ``order``. For + efficient iteration, ``NPY_KEEPORDER`` is the best option, and the other + orders enforce the particular iteration pattern. + + If ``op_dtypes`` isn't ``NULL``, it specifies a data type or ``NULL`` + for each ``op[i]``. + + If the parameter ``oa_ndim`` is not 0, you must also specify + ``op_axes``. These parameters let you control in detail how the + axes of the operand arrays get matched together and iterated. + In ``op_axes``, you must provide an array of ``niter`` pointers + to ``oa_ndim``-sized arrays of type ``npy_intp``. If an entry + in ``op_axes`` is NULL, normal broadcasting rules will apply. + In ``op_axes[j][i]`` is stored either a valid axis of ``op[j]``, or + -1 which means ``newaxis``. Within each ``op_axes[j]`` array, axes + may not be repeated. The following example is how normal broadcasting + applies to a 3-D array, a 2-D array, a 1-D array and a scalar.:: + + npy_intp oa_ndim = 3; 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}; + + If ``buffersize`` is zero, a default buffer size is used, + otherwise it specifies how big of a buffer to use. Buffers + which are powers of 2 such as 512 or 1024 are recommended. + + Returns NULL if there is an error, otherwise returns the allocated + iterator. + + Flags that may be passed in ``flags``, applying to the whole + iterator, are: + + ``NPY_ITER_C_ORDER_INDEX``, ``NPY_ITER_F_ORDER_INDEX`` + + Causes the iterator to track an index matching C or + Fortran order. These options are mutually exclusive. + + ``NPY_ITER_COORDS`` + + Causes the iterator to track array coordinates. + This prevents the iterator from coalescing axes to + produce bigger inner loops. + + ``NPY_ITER_NO_INNER_ITERATION`` + + Causes the iterator to skip iteration of the innermost + loop, allowing the user of the iterator to handle it. + + This flag is incompatible with ``NPY_ITER_C_ORDER_INDEX``, + ``NPY_ITER_F_ORDER_INDEX``, and ``NPY_ITER_COORDS``. + + ``NPY_ITER_COMMON_DATA_TYPE`` + + Causes the iterator to convert all the operands to a common + data type, calculated based on the ufunc type promotion rules. + The flags for each operand must be set so that the appropriate + casting is permitted, and copying or buffering must be enabled. + + If the common data type is known ahead of time, don't use this + flag. Instead, set the requested dtype for all the operands. + + ``NPY_ITER_OFFSETS`` + + Causes the iterator to produce offsets in the addresses + returned by ``NpyIter_GetDataPtrArray`` instead of + pointers. If used, cast the DataPtrArray value from + ``char **`` to ``npy_intp *``. + + This option is useful, together with the ``op_axes`` parameter, + for nested loops with two or more iterators. + The inner loop iterator should be made without + ``NPY_ITER_OFFSETS``, and any outer loop iterators should + be made with it. The sum of the inner loop pointers and + the outer loop offsets produces the innermost element addresses. + + This flag is incompatible with copying or buffering inputs. + + ``NPY_ITER_BUFFERED`` **PARTIALLY IMPLEMENTED** + + Causes the iterator to store buffering data, and use buffering + to satisfy data type, alignment, and byte-order requirements. + To buffer an operand, do not specify the ``NPY_ITER_COPY`` + or ``NPY_ITER_UPDATEIFCOPY`` flags, because they will + override buffering. Buffering is especially useful for Python + code using the iterator, allowing for larger chunks + of data at once to amortize the Python interpreter overhead. + + If used with ``NPY_ITER_NO_INNER_ITERATION``, the inner loop + for the caller may get larger chunks than would be possible + without buffering, because of how the strides are laid out. + + Note that if an operand is given the flag ``NPY_ITER_COPY`` + or ``NPY_ITER_UPDATEIFCOPY``, a copy will be made in preference + to buffering. Buffering will still occur when the array was + broadcast so elements need to be duplicated to get a constant + stride. + + When an operand is write buffered, it must either be an + aligned singleton, so buffering can be skipped for the operand, + or must match the dimensions of the iterator broadcast shape. + This is because the write back would otherwise overwrite + values multiple times to the same output, and accumulation + wouldn't work correctly. + + ``NPY_ITER_BUFFERED_GROWINNER`` + + Enables buffering as ``NPY_ITER_BUFFERED`` does, but when + buffering is not needed it grows the size of the inner loop + as much as possible. This option is best used if you're doing + a straight pass through all the data, rather than + anything with small cache-friendly arrays of + temporary values for each inner loop. + + Flags that may be passed in ``op_flags[i]``, where ``0 <= i < niter``: + + ``NPY_ITER_READWRITE``, ``NPY_ITER_READONLY``, ``NPY_ITER_WRITEONLY`` + + Indicate how the user of the iterator will read or write + to ``op[i]``. Exactly one of these flags must be specified + per operand. + + ``NPY_ITER_COPY`` + + Allow a copy of ``op[i]`` to be made if it does not + meet the data type or alignment requirements as specified + by the constructor flags and parameters. + + ``NPY_ITER_UPDATEIFCOPY`` + + Triggers ``NPY_ITER_COPY``, and when an array operand + is flagged for writing and is copied, causes the data + in a copy to be copied back to ``op[i]`` when the iterator + is destroyed. + + If the operand is flagged as write-only and a copy is needed, + an uninitialized temporary array will be created and then copied + to back to ``op[i]`` on destruction, instead of doing + the unecessary copy operation. + + ``NPY_ITER_NBO_ALIGNED`` + + Causes the iterator to provide data for ``op[i]`` + that is in native byte order and aligned according to + the dtype requirements. + + By default, the iterator produces pointers into the + arrays provided, which may be aligned or unaligned, and + with any byte order. If copying or buffering is not + enabled and an alignment fix or byte swapping is required, + an error will be raised. + + If the requested data type is non-native byte order, this + flag overrides it and the requested data type is converted to be + native byte order. + + ``NPY_ITER_SAFE_CASTS``, ``NPY_ITER_SAME_KIND_CASTS`` + ``NPY_ITER_UNSAFE_CASTS`` + + ``SAFE_CASTS`` means casting which the function + ``PyArray_CanCastSafely`` says is permitted. + ``SAME_KIND_CASTS`` adds to that any casting between types + of the same scalar kind, like ``double`` to ``float``. + ``UNSAFE_CASTS`` allows any casting for which a cast function + exists. + + To allow the casts to occur, copying or buffering must + also be enabled. + + ``NPY_ITER_ALLOCATE`` + + This is for output arrays, and requires that the flag + ``NPY_ITER_WRITEONLY`` be set. If ``op[i]`` is NULL, + creates a new array with the final broadcast dimensions, + and a layout matching the iteration order of the iterator. + + When ``op[i]`` is NULL, the requested data type + ``op_dtypes[i]`` may be NULL as well, in which case it is + automatically generated from the dtypes of the arrays which + are flagged as readable. The rules for generating the dtype + are the same is for UFuncs. Of special note is handling + of byte order in the selected dtype. If there is exactly + one input, the input's dtype is used as is. Otherwise, + if more than one input dtypes are combined together, the + output will be in native byte order. + + After being allocated with this flag, the caller may retrieve + the new array by calling ``NpyIter_GetObjectArray`` and + getting the i-th object in the returned C array. The caller + must call Py_INCREF on it to claim a reference to the array. + + ``NPY_ITER_NO_SUBTYPE`` + + For use with ``NPY_ITER_ALLOCATE``, this flag disables + allocating an array subtype for the output, forcing + it to be a straight ndarray. + + TODO: Maybe it would be better to introduce a function + ``NpyIter_GetWrappedOutput`` and remove this flag? + + ``NPY_ITER_WRITEABLE_REFERENCES`` + + By default, the iterator fails on creation if the iterator + has a writeable operand where the data type involves Python + references. Adding this flag indicates that the code using + the iterator is aware of this possibility and handles it + correctly. + + +``int NpyIter_UpdateIter(NpyIter *iter, npy_intp i, npy_uint32 op_flags, PyArray_Descr *dtype)`` **UNIMPLEMENTED** + + Updates the i-th operand within the iterator to possibly have a new + data type or more restrictive flag attributes. A use-case for + this is to allow the automatic allocation to determine an + output data type based on the standard NumPy type promotion rules, + then use this function to convert the inputs and possibly the + automatic output to a different data type during processing. + + This operation can only be done if ``NPY_ITER_COORDS`` was passed + as a flag to the iterator. If coordinates are not needed, + call the function ``NpyIter_RemoveCoords()`` once no more calls to + ``NpyIter_UpdateIter`` are needed. + + If the i-th operand has already been copied, an error is thrown. To + avoid this, leave all the flags out except the read/write indicators + for any operand that later has ``NpyIter_UpdateIter`` called on it. + + The flags that may be passed in ``op_flags`` are + ``NPY_ITER_COPY``, ``NPY_ITER_UPDATEIFCOPY``, + ``NPY_ITER_NBO_ALIGNED``, ``NPY_ITER_SAFE_CASTS``, + ``NPY_ITER_SAME_KIND_CASTS``, and ``NPY_ITER_UNSAFE_CASTS``. + +``int NpyIter_RemoveCoords(NpyIter *iter)`` + + If the iterator has coordinates, this strips support for them, and + does further iterator optimizations that are possible if coordinates + are not needed. This function also resets the iterator to its initial + state. + + **WARNING**: This function may change the internal memory layout of + the iterator. Any cached functions or pointers from the iterator + must be retrieved again! + + After calling this function, ``NpyIter_HasCoords(iter)`` will + return false. + + Returns ``NPY_SUCCEED`` or ``NPY_FAIL``. + +``int NpyIter_RemoveInnerLoop(NpyIter *iter)`` + + If UpdateIter/RemoveCoords was used, you may want to specify the + flag ``NPY_ITER_NO_INNER_ITERATION``. This flag is not permitted + together with ``NPY_ITER_COORDS``, so this function is provided + to enable the feature after ``NpyIter_RemoveCoords`` is called. + This function also resets the iterator to its initial state. + + **WARNING**: This function changes the internal logic of the iterator. + Any cached functions or pointers from the iterator must be retrieved + again! + + Returns ``NPY_SUCCEED`` or ``NPY_FAIL``. + +``int NpyIter_Deallocate(NpyIter *iter)`` + + Deallocates the iterator object. This additionally frees any + copies made, triggering UPDATEIFCOPY behavior where necessary. + + Returns ``NPY_SUCCEED`` or ``NPY_FAIL``. + +``void NpyIter_Reset(NpyIter *iter)`` + + Resets the iterator back to its initial state, at the beginning + of the arrays. + +``int NpyIter_GotoCoords(NpyIter *iter, npy_intp *coords)`` + + Adjusts the iterator to point to the ``ndim`` coordinates + pointed to by ``coords``. If the iterator was constructed without + the ``NPY_ITER_COORDS`` flag, or the coordinates are out of + bounds, returns an error. + + Returns ``NPY_SUCCEED`` or ``NPY_FAIL``. + +``int NpyIter_GotoIndex(NpyIter *iter, npy_intp index)`` + + Adjusts the iterator to point to the ``index`` specified. + If the iterator was constructed with the flag + ``NPY_ITER_C_ORDER_INDEX``, ``index`` is the C-order index, + and if the iterator was constructed with the flag + ``NPY_ITER_F_ORDER_INDEX``, ``index`` is the Fortran-order + index. Returns an error if there is no index being tracked, + or the index is out of bounds. + + Returns ``NPY_SUCCEED`` or ``NPY_FAIL``. + +``int NpyIter_HasInnerLoop(NpyIter *iter`` + + Returns 1 if the iterator handles the inner loop, + or 0 if the caller needs to handle it. This is controlled + by the constructor flag ``NPY_ITER_NO_INNER_ITERATION``. + +``int NpyIter_HasCoords(NpyIter *iter)`` + + Returns 1 if the iterator was created with the + ``NPY_ITER_COORDS`` flag, 0 otherwise. + +``int NpyIter_HasIndex(NpyIter *iter)`` + + Returns 1 if the iterator was created with the + ``NPY_ITER_C_ORDER_INDEX`` or ``NPY_ITER_F_ORDER_INDEX`` + flag, 0 otherwise. + +``int NpyIter_HasOffsets(NpyIter *iter)`` + + Returns 1 if the iterator was created with the + ``NPY_ITER_OFFSETS`` flag, 0 otherwise. + +``npy_intp NpyIter_GetNDim(NpyIter *iter)`` + + Returns the number of dimensions being iterated. If coordinates + were not requested in the iterator constructor, this value + may be smaller than the number of dimensions in the original + objects. + +``npy_intp NpyIter_GetNIter(NpyIter *iter)`` + + Returns the number of objects being iterated. + +``npy_intp NpyIter_GetIterSize(NpyIter *iter)`` + + Returns the number of times the iterator will iterate + starting from a new or reset state. + +``int NpyIter_GetShape(NpyIter *iter, npy_intp *outshape)`` + + Returns the broadcast shape of the iterator in ``outshape``. + This can only be called on an iterator which supports coordinates. + + Returns ``NPY_SUCCEED`` or ``NPY_FAIL``. + +``PyArray_Descr **NpyIter_GetDescrArray(NpyIter *iter)`` + + This gives back a pointer to the ``niter`` data type Descrs for + the objects being iterated. The result points into ``iter``, + so the caller does not gain any references to the Descrs. + + This pointer may be cached before the iteration loop, calling + ``iternext`` will not change it. + +``PyObject **NpyIter_GetObjectArray(NpyIter *iter)`` + + This gives back a pointer to the ``niter`` operand PyObjects + that are being iterated. The result points into ``iter``, + so the caller does not gain any references to the PyObjects. + +``PyObject *NpyIter_GetIterView(NpyIter *iter, npy_intp i)`` + + This gives back a reference to a new ndarray view, which is a view + into the i-th object in the array ``NpyIter_GetObjectArray()``, + whose dimensions and strides match the internal optimized + iteration pattern. A C-order iteration of this view is equivalent + to the iterator's iteration order. + + For example, if an iterator was created with a single array as its + input, and it was possible to rearrange all its axes and then + collapse it into a single strided iteration, this would return + a view that is a one-dimensional array. + +``void NpyIter_GetReadFlags(NpyIter *iter, char *outreadflags)`` + + Fills ``niter`` flags. Sets ``outreadflags[i]`` to 1 if + ``op[i]`` can be read from, and to 0 if not. + +``void NpyIter_GetWriteFlags(NpyIter *iter, char *outwriteflags)`` + + Fills ``niter`` flags. Sets ``outwriteflags[i]`` to 1 if + ``op[i]`` can be written to, and to 0 if not. + +Functions For Iteration +----------------------- + +``NpyIter_IterNext_Fn NpyIter_GetIterNext(NpyIter *iter)`` + + Returns a function pointer for iteration. A specialized version + of the function pointer may be calculated by this function + instead of being stored in the iterator structure. Thus, to + get good performance, it is required that the function pointer + be saved in a variable rather than retrieved for each loop iteration. + +``NpyIter_GetCoords_Fn NpyIter_GetGetCoords(NpyIter *iter)`` + + Returns a function pointer for getting the coordinates + of the iterator. Returns NULL if the iterator does not + support coordinates. It is recommended that this function + pointer be cached in a local variable before the iteration + loop. + +``char **NpyIter_GetDataPtrArray(NpyIter *iter)`` + + This gives back a pointer to the ``niter`` data pointers. If + ``NPY_ITER_NO_INNER_ITERATION`` was not specified, each data + pointer points to the current data item of the iterator. If + inner iteration was specified, it points to the first data + item of the inner loop. + + This pointer may be cached before the iteration loop, calling + ``iternext`` will not change it. + +``npy_intp *NpyIter_GetIndexPtr(NpyIter *iter)`` + + This gives back a pointer to the index being tracked, or NULL + if no index is being tracked. It is only useable if one of + the flags ``NPY_ITER_C_ORDER_INDEX`` or ``NPY_ITER_F_ORDER_INDEX`` + were specified during construction. + +When the flag ``NPY_ITER_NO_INNER_ITERATION`` is used, the code +needs to know the parameters for doing the inner loop. These +functions provide that information. + +``npy_intp *NpyIter_GetInnerStrideArray(NpyIter *iter)`` + + Returns a pointer to an array of the ``niter`` strides, + one for each iterated object, to be used by the inner loop. + + This pointer may be cached before the iteration loop, calling + ``iternext`` will not change it. + +``npy_intp* NpyIter_GetInnerLoopSizePtr(NpyIter *iter)`` + + Returns a pointer to the number of iterations the + inner loop should execute. + + This address may be cached before the iteration loop, calling + ``iternext`` will not change it. The value itself may change + during iteration, in particular if buffering is enabled. + +Examples +-------- + +A copy function using the iterator. The ``order`` parameter +is used to force the memory layout of the allocated +result. + +If the input is a reference type, this function will fail. +To fix this, the code must be changed to specially handle writeable +references, and add ``NPY_ITER_WRITEABLE_REFERENCES`` to the flags.:: + + /* 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, 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_GetObjectArray(iter)[1]; + Py_INCREF(ret); + + if (NpyIter_Deallocate(iter) != NPY_SUCCEED) { + Py_DECREF(ret); + return NULL; + } + + return ret; + } + +Python Example +-------------- + +The iterator has been mostly written and exposed to Python. To +see how it behaves, let's see what we can do with the np.add ufunc. +Even without changing the core of NumPy, we will be able to use +the iterator to make a faster add function. + +The Python exposure supplies two iteration interfaces, one which +follows the Python iterator protocol, and another which mirrors the +C-style do-while pattern. The native Python approach is better +in most cases, but if you need the iterator's coordinates or +index, use the C-style pattern. + +Here is how we might write an ``iter_add`` function, using the +Python iterator protocol.:: + + def iter_add_py(x, y, out=None): + addop = np.add + + it = np.newiter([x,y,out], [], + [['readonly'],['readonly'],['writeonly','allocate']]) + + for (a, b, c) in it: + addop(a, b, c) + + return it.operands[2] + +Here is the same function, but following the C-style pattern.:: + + def iter_add(x, y, out=None): + addop = np.add + + it = np.newiter([x,y,out], [], + [['readonly'],['readonly'],['writeonly','allocate']]) + + while not it.finished: + addop(it[0], it[1], it[2]) + it.iternext() + + return it.operands[2] + +Some noteworthy points about this function: + +* Cache np.add as a local variable to reduce namespace lookups +* Inputs are readonly, output is writeonly, and will be allocated + automatically if it is None. +* Uses np.add's out parameter to avoid an extra copy. + +Let's create some test variables, and time this function as well as the +built-in 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 + +At a thousand times slower, this is clearly not very good. One feature +of the iterator, designed to help speed up the inner loops, is the flag +``no_inner_iteration``. This is the same idea as the old iterator's +``PyArray_IterAllButAxis``, but slightly smarter. Let's modify +``iter_add`` to use this feature.:: + + def iter_add_noinner(x, y, out=None): + addop = np.add + + it = np.newiter([x,y,out], ['no_inner_iteration'], + [['readonly'],['readonly'],['writeonly','allocate']]) + + for (a, b, c) in it: + addop(a, b, c) + + return it.operands[2] + +The performance improves dramatically.:: + + In[6]: timeit iter_add_noinner(a, b) + 100 loops, best of 3: 7.1 ms per loop + +The performance is basically as good as the built-in function! It +turns out this is because the iterator was able to coalesce the last two +dimensions, resulting in 100 adds of 10000 elements each. If the +inner loop doesn't become as large, the performance doesn't improve +as dramatically. Let's use ``c`` instead of ``b`` to see how this works.:: + + In[7]: timeit iter_add_noinner(a, c) + 10 loops, best of 3: 76.4 ms per loop + +It's still a lot better than seven seconds, but still over ten times worse +than the built-in function. Here, the inner loop has 100 elements, +and its iterating 10000 times. If we were coding in C, our performance +would already be as good as the built-in performance, but in Python +there is too much overhead. + +This leads us to another feature of the iterator, its ability to give +us views of the iterated memory. The views it gives us are structured +so that processing them in C-order, like the built-in NumPy code does, +gives the same access order as the iterator itself. Effectively, we +are using the iterator to solve for a good memory access pattern, then +using other NumPy machinery to efficiently execute it. Let's +modify ``iter_add`` once again.:: + + def iter_add_itview(x, y, out=None): + it = np.newiter([x,y,out], [], + [['readonly'],['readonly'],['writeonly','allocate']]) + + (a, b, c) = it.itviews + np.add(a, b, c) + + return it.operands[2] + +Now the performance pretty closely matches the built-in function's.:: + + 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 + +Let us now step back to a case similar to the original motivation for the +new iterator. Here are the same calculations in Fortran memory order instead +Of C memory order.:: + + 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 + +As you can see, the performance of the built-in function dropped +significantly, but our newly-written add function maintained essentially +the same performance. As one final test, let's try several adds chained +together.:: + + 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 + +Also, just to check that it's doing the same thing,:: + + 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 +----------------------------------- + +For motivation, we had an example that did an 'over' composite operation +on two images. Now let's see how we can write the function with +the new iterator. + +Here is one of the original functions, for reference, and some +random image data.:: + + 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 + +Here's the same function, rewritten to use a new iterator. Note how +easy it was to add an optional output parameter.:: + + In [5]: def composite_over_it(im1, im2, out=None, buffersize=4096): + ....: it = np.newiter([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 + +A big speed improvement, over even the best previous attempt using +straight NumPy and a C-order array! By playing with the buffer size, we can +see how the speed improves until we hit the limits of the CPU cache +in the inner loop.:: + + 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 + +And finally, to double check that it's working, we can compare the two +functions.:: + + In [18]: np.all(composite_over(image1, image2) == + ...: composite_over_it(image1, image2)) + Out[18]: True + |