Beautiful Code [237]
A straightforward approach is to start at the end of the counter, or coordinate, array and work backward. At each index position, the code checks to see whether the coordinate is currently smaller than ni-1. If it is, it just adds 1 to that coordinate position and adds strides[i] to the memory address of the current value pointer. Whenever this happens, the for loop can break early, and the counter increment is done.
If the ith coordinate is greater than or equal to ni-1, it needs to be reset to 0 and (ni-1)xstrides[i] must be subtracted from the memory address of the current value pointer (to move back to the beginning of that dimension). In this case, the previous index position is checked.
All the necessary information can be represented in a structure we'll call it. The contents are:
coords
The coordinate index, N
dims_m1
The index of the highest element ni-1 for each dimension
strides
The stride in each dimension
backstrides
The amount to move back in order to return from the end of each dimension to the beginning: (ni-1) x strides[i]
nd_m1
The number of dimensions
currptr
A pointer to memory for the current position in the array
The code for the counter-tracking can then be written in C as follows:
for (i=it->nd_m1; i>=0; i--) {
if (it->coords[i] < it->dims_m1[i]) {
it->coords[i]++;
it->dataptr += it->strides[i];
break;
}
else {
it->coords[i] = 0;
it->dataptr -= it->backstrides[i];
}
}
This implementation uses the break statement and a for loop. We could instead have used a while statement and a flag indicating whether to continue looping:
done = 0;
i = it->nd_m1;
while (!done || i>=0) {
if (it->coords[i] < it->dims_m1[i]) {
it->coords[i]++;
it->dataptr += it->strides[i];
done = 1;
}
else {
it->coords[i] = 0;
it->dataptr -= it->backstrides[i];
}
i--;
}
Part of the reason I chose the for loop implementation is that the while loop looks a lot like the for loop (initialize counter, check against a value, decrement the counter), anyway. I typically reserve while loops for situations where the iteration requires more than a single iteration index. A bigger reason for choosing the for loop version, however, is that this code snippet implementing the counter increment will be used as a macro inside of every iterator loop. I wanted to avoid defining the extra variable done.
19.4.5. Iterator Structure
We are now in a position to understand the entire structure of the NumPy iterator. It's represented as the following struct in C:
typedef struct {
PyObject_HEAD
int nd_m1;
npy_intp index, size;
npy_intp coords[NPY_MAXDIMS];
npy_intp dims_m1[NPY_MAXDIMS];
npy_intp strides[NPY_MAXDIMS];
npy_intp backstrides[NPY_MAXDIMS];
npy_intp factors[NPY_MAXDIMS];
PyArrayObject *ao;
char *dataptr;
npy_bool contiguous;
} PyArrayIterObject;
The arrays in this structure (coords, dims_m1, strides, backstrides, and factors) are fixed-size arrays with dimensions controlled by the NPY_MAXDIMS constant. This choice was made to simplify memory management. However, it does limit the number of dimensions that can be used. It could easily be handled differently by dynamically allocating the needed memory when the iterator is created; such a change would not alter the fundamental behavior.
The npy_intp variables are integers just large enough to hold a pointer for the platform. npy_bool is a flag that should be either TRUE or FALSE. The PyObject_HEAD part of the structure contains the required portion that all Python objects must contain.
All of the variables have been explained before, but for clarity they are:
nd_m1
One less than the number of dimensions of the array: N-1.
index
A running counter indicating which element the iterator is currently at in the array. This counter runs from 0 to size-1.
size
The total number of elements in the array: n1 x n2 x…x nN.
coords
An array of N integers providing the counter, or the N-dimensional