Beautiful Code [240]
Alternatively, NumPy provides a multi-iterator object that can simplify dealing with several iterators at once. This multi-iterator object also handles the broadcasting functionality of NumPy automatically. Broadcasting is the name given to the feature in NumPy that allows arrays with different shapes to be used together in operations that are supposed to work element-by-element. For example, broadcasting allows a (4,1)-shaped array to be added to a (3)-shaped array resulting in a (4,3)-shaped array. Broadcasting also allows simultaneous iteration over a (4,1)-shaped array, a (3)-shaped array, and a (5,1,1)-shaped array to produce a broadcasted iteration covering the elements of a (5,4,3)-shaped array.
The rules of broadcasting are:
Arrays with fewer dimensions are treated as occupying the last dimensions of an array that has the full number of dimensions, so that all arrays have the same number of dimensions. The new, initial dimensions are filled in with 1s.
The length of each dimension in the final broadcast shape is the greatest length of that dimension in any of the arrays.
For each dimension, all inputs must either have the same number of elements as the broadcast result or a 1 as the number of elements.
Arrays with a single element in a particular dimension act as if that element were virtually copied to all positions during the iteration. In effect, the element is "broadcast" to the additional positions.
The key to the implementation of broadcasting consists of surprisingly simple modifications to the array iterators. With these alterations, standard iterator loops can be used to implement the resulting calculations in a straightforward way. The modifications needed are changes to the shape of the iterators (not the underlying array) and changes to the strides and backstrides. The shape stored in the iterator is changed to match the broadcast shape. The strides and backstrides for broadcast dimensions are changed to 0. With a stride of 0, the standard iterator does not actually move the data pointer to the element in memory as the index in that dimension proceeds. This creates the desired effect of broadcasting without actually copying the memory.
The following code illustrates usage of the multi-iterator object:
Code View: Scroll / Show All
PyObject *multi;
PyObject *in1, *in2;
double *i1p, *i2p, *op;
/* get in1 and in2 (assumed to be arrays of NPY_DOUBLE) */
/* first argument is the number of input arrays; the
next (variable number of) arguments are the
array objects */
multi = PyArray_MultiNew(2, in1, in2);
/* construct output array */
out = PyArray_SimpleNew(PyArray_MultiIter_NDIM(multi),
PyArray_MultiIter_DIMS(multi),
NPY_DOUBLE);
op = PyArray_DATA(out);
while(PyArray_MultiIter_NOTDONE(multi)) {
/* get (pointers to) the current value in each array */
i1p = PyArray_MultiIter_DATA(multi, 0);
i2p = PyArray_MultiIter_DATA(multi, 1);
/* perform the operation for this element */
*op = *ip1 + *ip2
op += 1; /* Advance output array pointer */
/* Advance all the input iterators */
PyArray_MultiIter_NEXT(multi);
}
The code is very similar to a standard iterator loop, except the multi-iterator handles adjustments of the input iterators to accomplish broadcasting, as well as incrementing all the other input iterators. This code handles the broadcasting automatically as part of the iterator processing, so that the addition of a (3,1)-shaped array to a (4)-shaped one will produce a (3,4)-shaped output array.
19.6.3. Anecdotes
The NumPy iterator is used throughout the NumPy code base to simplify the construction of N-dimensional loops. Having the iterator available allowed me to code algorithms for more general (noncontiguous) arrays. Normally, the difficulty of thinking about how to handle the noncontiguous