Beautiful Code [238]
dims_m1
An array of N integers providing one less than the number of elements along each dimension: n1-1, n2-1,…, nN-1.
strides
An array of N integers providing the number of bytes to skip when advancing to the next element in a particular dimension.
backstrides
An array of N integers providing the number of bytes to subtract when the internal index counter rolls from ni-1 to 0 in a particular dimension.
factors
An array of factors useful in rapidly calculating the mapping between the one-dimensional index and the N-dimensional coords array. This array is needed only if PyArray_ITER_GOTO1D is called.
ao
A pointer to the underlying array this iterator is built from.
datapr
A pointer to the (first byte of) the current value of the array.
contiguous
TRUE (1) if this iterator is for a contiguous array and FALSE (0) if otherwise. This is the same as (ao->flags & NPY_C_CONTIGUOUS). It's much simpler to find the next element in the array each time when the array is contiguous, so it is worth checking for.
Multidimensional Iterators in NumPy > Iterator Interface
19.5. Iterator Interface
The iterator is implemented in NumPy using a combination of macros and function calls. An iterator is created using the C-API function call it=PyArray_IterNew(ao). The check for iterator termination can be accomplished using the macro PyArray_ITER_NOTDONE(it). Finally, the next position in the iterator is accomplished using PyArray_ITER_NEXT(it), which is a macro to ensure that it occurs inline (avoiding the function call). Ideally, this macro would be an inline function because it is sufficiently complicated. However, because NumPy is written to ANSI C, which does not define inline functions, a macro is used. Finally, the pointer to the first byte of the current value can be obtained using PyArray_ITER_DATA(it), which avoids referencing the structure member dataptr directly (allowing for future name changes to the structure members).
An example of the iterator interface is the following code snippet, which computes the largest value in an N-dimensional array. We assume that the array is named ao, has elements of type double, and is correctly aligned:
#include double *currval, maxval=-DBL_MAX; PyObject *it; it = PyArray_IterNew(ao); while (PyArray_ITER_NOTDONE(it)) { currval = (double *)PyArray_ITER_DATA(it); if (*currval < maxval) maxval = *currval; PyArray_ITER_NEXT(it); } This code shows how relatively easy it is to construct a loop for a noncontiguous, N-dimensional array using the iterator structure. The simplicity of this code also illustrates the elegance of iterator abstraction. Notice how similar the code is to the simple iterator pseudocode shown at the beginning of the earlier section "Iterator Design." Consider also that this code works for arrays of arbitrary dimensions and arbitrary strides in each dimension, and you begin to appreciate the beauty of the multidimensional iterator. The iterator-based code is fast for both contiguous and noncontiguous arrays. However, the fastest contiguous-array loop is still something like: double *currval, maxval=-MAX_DOUBLE; int size; currval = (double *)PyArray_DATA(ao); size = PyArray_SIZE(ao); while (size--) { if (*currval > maxval) maxval = *currval; currval += 1; } The real benefit of the NumPy iterator is that it allows programmers to write contiguous-like code that is still fairly fast without worrying about whether their arrays are contiguous. It should be remembered that forcing a contiguous algorithm has performance costs as well because noncontiguous data must be copied to another array for processing. The speed difference between the NumPy iterator solution and the fastest contiguous-case solution could be largely eliminated if a remaining problem with the current NumPy iterator interface could be fixed. The problem is that the PyArray_ITER_NEXT macro checks each time through the loop whether the iterator can use the simplified contiguous approach. Ideally,