231: Intervals and Generalized Arrays231: Intervals and Generalized Arraysby Bradley J. Lucier
SRFI 231 is currently in final status. Here is an explanation of each status that a SRFI can hold. To provide input on this SRFI, please send email to srfi-231@nospamsrfi.schemers.org. To subscribe to the list, follow these instructions. You can access previous messages via the mailing list archive.
The NumPy array library is a popular library for scientific computing that extends the Python programming language. I'd like to compare some of the properties of NumPy arrays and SRFI 231 arrays. I should state from the outset that (1) I am not in any way an expert in Python or NumPy, and (2) this commentary is not meant in any way as a criticism of NumPy, which has had tremendous success. We refer in the following to the NumPy API Reference section on Array manipulation routines.
I believe that NumPy arrays are conceptually similar to SRFI 231's "specialized arrays", i.e., with elements stored in a one-dimensional "vector" of some element type accessed by addresses given by $a_0\times i_0+\cdots a_{n-1}\times i_{n-1}+b$ when an array index is $i_0,\ldots,i_{n-1}$ and $a_0,\ldots,a_{n-1},b$ are integer constants. SRFI 231 also has "generalized arrays", with a general mapping between multi-indices $i_0,\ldots,i_{n-1}$ and values, without a backing vector.
NumPy array indices begin at zero; SRFI 231 arrays have arbitrary (exact integer) lower and upper bounds, though one can specify only the upper bounds of the domain of an array and the lower bounds are assumed to be zero. Because array axes can have nonzero lower bounds, SRFI 231 contains an array-translate routine that translates the domain of valid indices. SRFI 231 may benefit from a routine to automate translate an array's domain to set all axis lower bounds to zero.
NumPy uses powerful index notation for specifying subarrays, notation along the lines of start:end:skip to specify the beginning index along an axis of a subarray, the end index, and the number of elements to skip between successive indices. In SRFI 231 array-extract specifies subarrays, and array-sample samples elements that are uniformly spaced in an array; the latter routine requires that an array's lower bounds all be zero. Combining array-extract and array-sample provides roughly the same functionality as NumPy's index notation.
NumPy sometimes has a series of routines that do similar things, but along different axes.
For example, there is this section in the API reference:
| Joining arrays | |
|---|---|
| concatenate([axis, out, dtype, casting]) | Join a sequence of arrays along an existing axis |
| stack(arrays[, axis, out, dtype, casting]) | Join a sequence of arrays along a new axis. |
| block(arrays) | Assemble an nd-array from nested lists of blocks. |
| vstack(tup, *[, dtype, casting]) | Stack arrays in sequence vertically (row wise). |
| hstack(tup, *[, dtype, casting]) | Stack arrays in sequence horizontally (column wise). |
| dstack(tup, *[, dtype, casting]) | Stack arrays in sequence depth wise (along the third axis). |
| column_stack(tup) | Stack 1-D arrays as columns into a 2-D array. |
| row_stack(tup, *[, dtype, casting]) | Stack arrays in sequence vertically (row wise) |
The fundamental NumPy routines here are concatenate, where the dimension of the result is the same as the dimension of the argument arrays, and stack, where the dimension of the result is one more than the dimension of the argument arrays. SRFI 231 has the associated routines array-append and array-stack.
hstack, vstack, and dstack either stack or concatenate lists of arrays, depending on dimension. (I am not kidding.)
SRFI 231 also has array-decurry, which takes as an argument an array of arrays and allows simultaneous "stacking" of arrays, adding multiple new axes. In this case, all the new axes are placed before the axes of the argument arrays; these new axes can be redistributed with SRFI 231's array-permute. I don't see a corresponding routine in NumPy, i.e., a multi-axis "stack" routine.
NumPy's block takes a nested list of arrays and assembles them into one large array of the same dimension, but for some reason doesn't have the same optional arguments as the other routines. I don't know whether there are other routines that can operate on nested lists in NumPy. SRFI 231's array-block, inspired by numpy.block, takes an array of arrays as an argument.
For "splitting" arrays, NumPy has the following routines:
| Splitting arrays | |
|---|---|
| split(ary, indices_or_sections[, axis]) | Split an array into multiple sub-arrays as views into ary. |
| array_split(ary, indices_or_sections[, axis]) | Split an array into multiple sub-arrays. |
| dsplit(ary, indices_or_sections) | Split array into multiple sub-arrays along the 3rd axis (depth). |
| hsplit(ary, indices_or_sections) | Split an array into multiple sub-arrays horizontally (column-wise). |
| vsplit(ary, indices_or_sections) | Split an array into multiple sub-arrays vertically (row-wise). |
As far as I can tell, all these routines split an array into parts along one axis direction, with split and array_split being the most general.
By contrast, SRFI 231's array-tile can cut an array along various axes simultaneously, returning an array of subarrays; one can reconstruct the original array (translated to zero lower bounds) with array-block.
The NumPy routines *split are undone by NumPy's concatenate. I know of no NumPy routine that will decompose an array into sections appropriate to be reconstructed with NumPy's block.
NumPy has powerful syntactic notation for considering various views and slices of arrays, but I know of no NumPy routine that decomposes an array in a way that is reconstituted with the NumPy *stack routine. In contrast, one can combine SRFI 231's array-permute and array-curry to decompose an array by dimension, which is reconstituted with SRFI 231's array-stack (if along only one dimension) or array-decurry (if along several dimensions at once).
NumPy has routines that it categorizes as "Rearranging elements":
| Rearranging elements | |
|---|---|
| flip(m[, axis]) | Reverse the order of elements in an array along the given axis. |
| fliplr(m) | Reverse the order of elements along axis 1 (left/right). |
| flipud(m) | Reverse the order of elements along axis 0 (up/down). |
| reshape(a, newshape[, order]) | Gives a new shape to an array without changing its data. |
| roll(a, shift[, axis]) | Roll array elements along a given axis. |
| rot90(m[, k, axes]) | Rotate an array by 90 degrees in the plane specified by axes. |
NumPy's various "flip" routines are achieved by array-reverse of SRFI 231; NumPy's "flip" routines operate on one axis at a time, SRFI 231's array-reverse operate on all axes at once. NumPy's reshape motivated, and is basically the same as, SRFI 231's specialized-array-reshape. SRFI 231 does not have builtin functionality mimicking NumPy's roll. NumPy's rot90 routine is a special case of combining SRFI 231's array-reverse and array-permute. Indeed, by combining one call to array-reverse with one call to array-permute, SRFI 231 can generate all "symmetries" of a multi-dimensional hypercube.
NumPy has what it calls "Transpose-like operations:
| Transpose-like operations | |
|---|---|
| moveaxis(a, source, destination) | Move axes of an array to new positions. |
| rollaxis(a, axis[, start]) | Roll the specified axis backwards, until it lies in a given position. |
| swapaxes(a, axis1, axis2) | Interchange two axes of an array. |
| ndarray.T | View of the transposed array. |
| transpose(a[, axes]) | Returns an array with axes transposed. |
All of NumPy's various transpose-like operations, the most general of which is moveaxis, are effected in SRFI 231 with array-permute. SRFI 231 has a number of auxiliary routines that specify common useful permutations (rearrangements) of axes, or indices, as they're called here: index-first, index-last, index-rotate, and index-swap, the results of which are used as arguments to array-permute.
One philosophy of SRFI 231 is to separate the routines that specify actions and the routines that carry out actions.
For example, one may find the maximum of a nonempty array A of numbers with (array-reduce max A). One has to examine every element of A to determine the maximum.
Say now you want to find the maximum of the squares of the elements of A. It's clear that you'll need to evaluate the squares of all the elements, but you don't need to store the squares all at one time in their own array with a backing store. So when one says (array-reduce max (array-map square A)), the array-map routine returns a "generalized array", which, in effect, extracts an element of A and squares it for each multi-index in the domain of A.
So SRFI 231 routines like array-copy, array-assign!, array-fold-left, array-reduce, array-every, array-any, etc., actually do something with all elements of an array, while, e.g., array-map, array-outer-product, and array-inner-product return generalized arrays that specify how to compute array elements later, as needed, for each multi-index, perhaps when passed to an "action" routine like array-reduce.
specialized-array-reshapeSRFI 231's specialized-array-reshape is basically the same as NumPy's reshape routine. I'd like to give a simpler example here of what it does.
Assume that you have a specialized array that you'd like to reshape into a one-dimensional specialized array.
In SRFI 231, indexer functions for one-dimensional specialized arrays, which map the single argument $i$ to an index into a vector where the array element is stored, must be affine, i.e., of the form $a\times i+b$. This implies, for example, that the array elements are stored in uniformly spaced boxes within the vector, with distance $|a|$ between boxes.
And if a multi-dimensional specialized array's elements are not stored with this pattern inside the backing vector, it can't be reshaped to a one-dimensional vector.
As an example, consider the array (define A (list->array (make-interval '#(4 4)) (iota 16))). Its elements are stored consecutively in a vector of length 16, and we see:
(array->list* A) => ((0 1 2 3) (4 5 6 7) (8 9 10 11) (12 13 14 15))
We then have
(array->list* (array-scale A '#(1 2))) => ((0 2) (4 6) (8 10) (12 14))
and we see that the locations of thse elements in the backing store are
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15) ^ ^ ^ ^ ^ ^ ^ ^
i.e., they are uniformly spaced in the backing vector, and this array can be reshaped to be one dimensional. On the other hand, we have
(array->list* (array-sample A '#(2 1))) => ((0 1 2 3) (8 9 10 11))
and the locations of these elements in the backing store are
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15) (^ ^ ^ ^ ^ ^ ^ ^
and these are not uniformly spaced in the backing store, and this array cannot be reshaped into a one-dimensional array.
There are similar requirements for arrays of other dimensions.