srfi-179

TOC »

srfi-179

Nonempty Intervals and Generalized Arrays (Updated)

Chicken implementation and documentation notes

NOTE: These are just differences to the reference implementation, not deviations from the SRFI specification. If any behavior is different from the SRFI, please contact the egg maintainer.

NOTE: This document is mostly the original SRFI html wrangled into wiki syntax using pandoc and manual edits. If anything looks wrong, please contact the egg maintainer.

Complex number vector operations (c64, c128) are provided by srfi-160
specialized-array-default-safe? and specialized-aray-default-mutable? are implemented in terms of chicken's make-parameter.
Checkers for flonum storage classes (f32, f64, c64, c128) have been relaxed according to acceptable values to chicken's srfi-4 -set! procedures.

This SRFI specifies an array mechanism for Scheme. Arrays as defined here are quite general; at their most basic, an array is simply a mapping, or function, from multi-indices of exact integers $i_0,\ldots,i_{d-1}$ to Scheme values. The set of multi-indices $i_0,\ldots,i_{d-1}$ that are valid for a given array form the domain of the array. In this SRFI, each array's domain consists of the cross product of nonempty intervals of exact integers $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ of $\mathbb Z^d$, $d$-tuples of integers. Thus, we introduce a data type called $d$-intervals, or more briefly intervals, that encapsulates this notion. (We borrow this terminology from, e.g., Elias Zakon's Basic Concepts of Mathematics.) Specialized variants of arrays are specified to provide portable programs with efficient representations for common use cases.

Rationale

This SRFI was motivated by a number of somewhat independent notions, which we outline here and which are explained below.

Provide a general API (Application Program Interface) that specifies the minimal required properties of any given array, without requiring any specific implementation strategy from the programmer for that array.
Provide a single, efficient implementation for dense arrays (which we call specialized arrays).
Provide useful array transformations by exploiting the algebraic structure of affine one-to-one mappings on multi-indices.
Separate the routines that specify the work to be done (array-map, array-outer-product, etc.) from the routines that actually do the work (array-copy, array-assign!, array-fold, etc.). This approach avoids temporary intermediate arrays in computations.
Encourage bulk processing of arrays rather than word-by-word operations.

This SRFI differs from the finalized SRFI 122 in the following ways:

The procedures specialized-array-default-mutable?, interval-for-each, interval-cartesian-product, interval-rotate and array-elements-in-order?, array-outer-product, array-tile, array-rotate, array-reduce, array-assign!, array-ref, array-set!, and specialized-array-reshape have been added together with some examples.
Global variables f8-storage-class and f16-storage-class have been added.
Homogeneous storage classes must be implemented using homogeneous vectors, or be defined as #f.
The procedure make-interval now takes one or two arguments.
Specialized arrays can be mutable or immutable; the default, which can be changed, is mutable. Shared arrays inherit safety and mutability from source arrays.
The discussion of Haar transforms as examples of separable transforms has been corrected.
The documentation has a few more examples of image processing algorithms.
Some matrix examples have been added to this document.

Overview

Bawden-style arrays

In a 1993 post to the news group comp.lang.scheme, Alan Bawden gave a simple implementation of multi-dimensional arrays in R4RS scheme. The only constructor of new arrays required specifying an initial value, and he provided the three low-level primitives array-ref, array-set!, and array?, as well as make-array and make-shared-array. His arrays were defined on rectangular intervals in $\mathbb Z^d$ of the form $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$. I'll note that his function array-set! put the value to be entered into the array at the front of the variable-length list of indices that indicate where to place the new value. He offered an intriguing way to "share" arrays in the form of a routine make-shared-array that took a mapping from a new interval of indices into the domain of the array to be shared. His implementation incorporated what he called an indexer, which was a function from the interval $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$ to an interval $[0,N)$, where the body of the array consisted of a single Scheme vector of length $N$. Bawden called the mapping specified in make-shared-array linear, but I prefer the term affine, as I explain later.

Mathematically, Bawden's arrays can be described as follows. We'll use the vector notation $\vec i$ for a multi-index $i_0,\ldots,i_{d-1}$. (Multi-indices correspond to Scheme values.) Arrays will be denoted by capital letters $A,B,\ldots$, the domain of the array $A$ will be denoted by $D_A$, and the indexer of $A$, mapping $D_A$ to the interval $[0,N)$, will be denoted by $I_A$. Initially, Bawden constructs $I_A$ such that $I_A(\vec i)$ steps consecutively through the values $0,1,\ldots,N-1$ as $\vec i$ steps through the multi-indices $(l_0,\ldots,l_{d-2},l_{d-1})$, $(l_0,\ldots,l_{d-2},l_{d-1}+1)$, $\ldots$, $(l_0,\ldots,l_{d-2}+1,l_{d-1})$, etc., in lexicographical order, which means that if $\vec i$ and $\vec j$ are two multi-indices, then $\vec i<\vec j$ iff the first coordinate $k$ where $\vec i$ and $\vec j$ differ satisfies $\vec i_k<\vec j_k$.

In make-shared-array, Bawden allows you to specify a new $r$-dimensional interval $D_B$ as the domain of a new array $B$, and a mapping $T_{BA}:D_B\to D_A$ of the form $T_{BA}(\vec i)=M\vec i+\vec b$; here $M$ is a $d\times r$ matrix of integer values and $\vec b$ is a $d$-vector. So this mapping $T_{BA}$ is affine, in that $T_{BA}(\vec i)-T_{BA}(\vec j)=M(\vec i-\vec j)$ is linear (in a linear algebra sense) in $\vec i-\vec j$. The new indexer of $B$ satisfies $I_B(\vec i)=I_A(T_{BA}(\vec i))$.

A fact Bawden exploits in the code, but doesn't point out in the short post, is that $I_B$ is again an affine map, and indeed, the composition of any two affine maps is again affine.

Our extensions of Bawden-style arrays

We incorporate Bawden-style arrays into this SRFI, but extend them in one minor way that we find useful.

We introduce the notion of a storage class, an object that contains functions that manipulate, store, check, etc., different types of values. A generic-storage-class can manipulate any Scheme value, whereas, e.g., a u1-storage-class can store only the values 0 and 1 in each element of a body.

We also require that our affine maps be one-to-one, so that if $\vec i\neq\vec j$ then $T(\vec i)\neq T(\vec j)$. Without this property, modifying the $\vec i$th component of $A$ would cause the $\vec j$th component to change.

Common transformations on Bawden-style arrays

Requiring the transformations $T_{BA}:D_B\to D_A$ to be affine may seem esoteric and restricting, but in fact many common and useful array transformations can be expressed in this way. We give several examples below:

Restricting the domain of an array: If the domain of $B$, $D_B$, is a subset of the domain of $A$, then $T_{BA}(\vec i)=\vec i$ is a one-to-one affine mapping. We define array-extract to define this common operation; it's like looking at a rectangular sub-part of a spreadsheet. We use it to extract the common part of overlapping domains of three arrays in an image processing example below.
Tiling an array: For various reasons (parallel processing, optimizing cache localization, GPU programming, etc.) one may wish to process a large array as a number of subarrays of the same dimensions, which we call tiling the array. The routine array-tile returns a new array, each entry of which is a subarray extracted (in the sense of array-extract) from the input array.
Translating the domain of an array: If $\vec d$ is a vector of integers, then $T_{BA}(\vec i)=\vec i-\vec d$ is a one-to-one affine map of $D_B=\{\vec i+\vec d\mid \vec i\in D_A\}$ onto $D_A$. We call $D_B$ the translate of $D_A$, and we define array-translate to provide this operation.
Permuting the coordinates of an array: If $\pi$ permutes the coordinates of a multi-index $\vec i$, and $\pi^{-1}$ is the inverse of $\pi$, then $T_{BA}(\vec i)=\pi (\vec i)$ is a one-to-one affine map from $D_B=\{\pi^{-1}(\vec i)\mid \vec i\in D_A\}$ onto $D_A$. We provide array-permute for this operation. (The only nonidentity permutation of a two-dimensional spreadsheet turns rows into columns and vice versa.) We also provide array-rotate for the special permutations that rotate the axes. For example, in three dimensions we have the following three rotations: $i\ j\ k\to j\ k\ i$; $i\ j\ k\to k\ i\ j$; and the trivial (identity) rotation $i\ j\ k\to i\ j\ k$. The three-dimensional permutations that are not rotations are $i\ j\ k\to i\ k\ j$; $i\ j\ k\to j\ i\ k$; and $i\ j\ k\to k\ j\ i$.
Currying an array: Let's denote the cross product of two intervals $\text{Int}_1$ and $\text{Int}_2$ by $\text{Int}_1\times\text{Int}_2$; if $\vec j=(j_0,\ldots,j_{r-1})\in \text{Int}_1$ and $\vec i=(i_0,\ldots,i_{s-1})\in \text{Int}_2$, then $\vec j\times\vec i$, which we define to be $(j_0,\ldots,j_{r-1},i_0,\ldots,i_{s-1})$, is in $\text{Int}_1\times\text{Int}_2$. If $D_A=\text{Int}_1\times\text{Int}_2$ and $\vec j\in\text{Int}_1$, then $T_{BA}(\vec i)=\vec j\times\vec i$ is a one-to-one affine mapping from $D_B=\text{Int}_2$ into $D_A$. For each vector $\vec j$ we can compute a new array in this way; we provide array-curry for this operation, which returns an array whose domain is $\text{Int}_1$ and whose elements are themselves arrays, each of which is defined on $\text{Int}_2$. Currying a two-dimensional array would be like organizing a spreadsheet into a one-dimensional array of rows of the spreadsheet.
Traversing some indices in a multi-index in reverse order: Consider an array $A$ with domain $D_A=[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$. Fix $D_B=D_A$ and assume we're given a vector of booleans $F$ ($F$ for "flip?"). Then define $T_{BA}:D_B\to D_A$ by $i_j\to i_j$ if $F_j$ is #f and $i_j\to u_j+l_j-1-i_j$ if $F_j$ is #t. In other words, we reverse the ordering of the $j$th coordinate of $\vec i$ if and only if $F_j$ is true. $T_{BA}$ is an affine mapping from $D_B\to D_A$, which defines a new array $B$, and we can provide array-reverse for this operation. Applying array-reverse to a two-dimensional spreadsheet might reverse the order of the rows or columns (or both).
Uniformly sampling an array: Assume that $A$ is an array with domain $[0,u_1)\times\cdots\times[0,u_{d-1})$ (i.e., an interval all of whose lower bounds are zero). We'll also assume the existence of vector $S$ of scale factors, which are positive exact integers. Let $D_B$ be a new interval with $j$th lower bound equal to zero and $j$th upper bound equal to $\operatorname{ceiling}(u_j/S_j)$ and let $T_{BA}(\vec i)_j=i_j\times S_j$, i.e., the $j$th coordinate is scaled by $S_j$. ($D_B$ contains precisely those multi-indices that $T_{BA}$ maps into $D_A$.) Then $T_{BA}$ is an affine one-to-one mapping, and we provide interval-scale and array-sample for these operations.

We make several remarks. First, all these operations could have been computed by specifying the particular mapping $T_{BA}$ explicitly, so that these routines are simply "convenience" procedures. Second, because the composition of any number of affine mappings is again affine, accessing or changing the elements of a restricted, translated, curried, permuted array is no slower than accessing or changing the elements of the original array itself. Finally, we note that by combining array currying and permuting, say, one can come up with simple expressions of powerful algorithms, such as extending one-dimensional transforms to multi-dimensional separable transforms, or quickly generating two-dimensional slices of three-dimensional image data. Examples are given below.

Generalized arrays

Bawden-style arrays are clearly useful as a programming construct, but they do not fulfill all our needs in this area. An array, as commonly understood, provides a mapping from multi-indices $(i_0,\ldots,i_{d-1})$ of exact integers in a nonempty, rectangular, $d$-dimensional interval $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ (the domain of the array) to Scheme objects. Thus, two things are necessary to specify an array: an interval and a mapping that has that interval as its domain.

Since these two things are often sufficient for certain algorithms, we introduce in this SRFI a minimal set of interfaces for dealing with such arrays.

Specifically, an array specifies a nonempty, multi-dimensional interval, called its domain, and a mapping from this domain to Scheme objects. This mapping is called the getter of the array, accessed with the procedure array-getter; the domain of the array (more precisely, the domain of the array's getter) is accessed with the procedure array-domain.

If this mapping can be changed, the array is said to be mutable and the mutation is effected by the array's setter, accessed by the procedure array-setter. We call an object of this type a mutable array. Note: If an array does not have a setter, then we call it immutable even though the array's getter might not be a "pure" function, i.e., the value it returns may not depend solely on the arguments passed to the getter.

In general, we leave the implementation of generalized arrays completely open. They may be defined simply by closures, or they may have hash tables or databases behind an implementation, one may read the values from a file, etc.

In this SRFI, Bawden-style arrays are called specialized. A specialized array is an example of a mutable array.

Sharing generalized arrays

Even if an array $A$ is not a specialized array, then it could be "shared" by specifying a new interval $D_B$ as the domain of a new array $B$ and an affine map $T_{BA}:D_B\to D_A$. Each call to $B$ would then be computed as $B(\vec i)=A(T_{BA}(\vec i))$.

One could again "share" $B$, given a new interval $D_C$ as the domain of a new array $C$ and an affine transform $T_{CB}:D_C\to D_B$, and then each access $C(\vec i)=A(T_{BA}(T_{CB}(\vec i)))$. The composition $T_{BA}\circ T_{CB}:D_C\to D_A$, being itself affine, could be precomputed and stored as $T_{CA}:D_C\to D_A$, and $C(\vec i)=A(T_{CA}(\vec i))$ can be computed with the overhead of computing a single affine transformation.

So, if we wanted, we could share generalized arrays with constant overhead by adding a single layer of (multi-valued) affine transformations on top of evaluating generalized arrays. Even though this could be done transparently to the user, we do not do that here; it would be a compatible extension of this SRFI to do so. We provide only the routine specialized-array-share, not a more general array-share.

Certain ways of sharing generalized arrays, however, are relatively easy to code and not that expensive. If we denote (array-getter A) by A-getter, then if B is the result of array-extract applied to A, then (array-getter B) is simply A-getter. Similarly, if A is a two-dimensional array, and B is derived from A by applying the permutation $\pi((i,j))=(j,i)$, then (array-getter B) is (lambda (i j) (A-getter j i)). Translation and currying also lead to transformed arrays whose getters are relatively efficiently derived from A-getter, at least for arrays of small dimension.

Thus, while we do not provide for sharing of generalized arrays for general one-to-one affine maps $T$, we do allow it for the specific functions array-extract, array-translate, array-permute, array-curry, array-reverse, array-tile, array-rotate and array-sample, and we provide relatively efficient implementations of these functions for arrays of dimension no greater than four.

Array-map does not produce a specialized array

Daniel Friedman and David Wise wrote a famous paper CONS should not Evaluate its Arguments. In the spirit of that paper, our procedure array-map does not immediately produce a specialized array, but a simple immutable array, whose elements are recomputed from the arguments of array-map each time they are accessed. This immutable array can be passed on to further applications of array-map for further processing without generating the storage bodies for intermediate arrays.

We provide the procedure array-copy to transform a generalized array (like that returned by array-map) to a specialized, Bawden-style array, for which accessing each element again takes $O(1)$ operations.

Notational convention

If A is an array, then we generally define A_ to be (array-getter A) and A! to be (array-setter A). The latter notation is motivated by the general Scheme convention that the names of functions that modify the contents of data structures end in !, while the notation for the getter of an array is motivated by the TeX notation for subscripts. See particularly the Haar transform example.

Notes

Relationship to nonstrict arrays in Racket. It appears that what we call simply arrays in this SRFI are called nonstrict arrays in the math/array library of Racket, which in turn was influenced by an array proposal for Haskell. Our "specialized" arrays are related to Racket's "strict" arrays.
Indexers. The argument new-domain->old-domain to specialized-array-share is, conceptually, a multi-valued array.
Source of function names. The function array-curry gets its name from the curry operator in programming—we are currying the getter of the array and keeping careful track of the domains. interval-projections can be thought of as currying the characteristic function of the interval, encapsulated here as interval-contains-multi-index?.
Choice of functions on intervals. The choice of functions for both arrays and intervals was motivated almost solely by what I needed for arrays.
No empty intervals. This SRFI considers arrays over only nonempty intervals of positive dimension. The author of this proposal acknowledges that other languages and array systems allow either zero-dimensional intervals or empty intervals of positive dimension, but prefers to leave such empty intervals as possibly compatible extensions to the current proposal.
Multi-valued arrays. While this SRFI restricts attention to single-valued arrays, wherein the getter of each array returns a single value, allowing multi-valued immutable arrays would be a compatible extension of this SRFI.
No low-level specialized array constructor. While the author of the SRFI uses mainly (make-array ...), array-map, and array-copy to construct arrays, and while there are several other ways to construct arrays, there is no really low-level interface given for constructing specialized arrays (where one specifies a body, an indexer, etc.). It was felt that certain difficulties, some surmountable (such as checking that a given body is compatible with a given storage class) and some not (such as checking that an indexer is indeed affine), made a low-level interface less useful. At the same time, the simple (make-array ...) mechanism is so general, allowing one to specify getters and setters as general functions, as to cover nearly all needs.

Specification

Miscellaneous Functions: translation?, permutation?.
Intervals: make-interval, interval?, interval-dimension, interval-lower-bound, interval-upper-bound, interval-lower-bounds->list, interval-upper-bounds->list, interval-lower-bounds->vector, interval-upper-bounds->vector, interval=, interval-volume, interval-subset?, interval-contains-multi-index?, interval-projections, interval-for-each, interval-dilate, interval-intersect, interval-translate, interval-permute, interval-rotate, interval-scale, interval-cartesian-product.
Storage Classes: make-storage-class, storage-class?, storage-class-getter, storage-class-setter, storage-class-checker, storage-class-maker, storage-class-copier, storage-class-length, storage-class-default, generic-storage-class, s8-storage-class, s16-storage-class, s32-storage-class, s64-storage-class, u1-storage-class, u8-storage-class, u16-storage-class, u32-storage-class, u64-storage-class, f8-storage-class, f16-storage-class, f32-storage-class, f64-storage-class, c64-storage-class, c128-storage-class.
Arrays: make-array, array?, array-domain, array-getter, array-dimension, mutable-array?, array-setter, specialized-array-default-safe?, specialized-array-default-mutable?, make-specialized-array, specialized-array?, array-storage-class, array-indexer, array-body, array-safe?, array-elements-in-order?, specialized-array-share, array-copy, array-curry, array-extract, array-tile, array-translate, array-permute, array-rotate, array-reverse, array-sample, array-outer-product, array-map, array-for-each, array-fold, array-fold-right, array-reduce, array-any, array-every, array->list, list->array, array-assign!, array-ref, array-set!, specialized-array-reshape.

Miscellaneous Functions

This document refers to translations and permutations. A translation is a vector of exact integers. A permutation of dimension $n$ is a vector whose entries are the exact integers $0,1,\ldots,n-1$, each occurring once, in any order.

Procedures

translation?

translation? objectprocedure: Returns #t if object is a translation, and #f otherwise.

permutation?

permutation? objectprocedure: Returns #t if object is a permutation, and #f otherwise.

Intervals

An interval represents the set of all multi-indices of exact integers $i_0,\ldots,i_{d-1}$ satisfying $l_0\leq i_0<u_0,\ldots,l_{d-1}\leq i_{d-1}<u_{d-1}$, where the lower bounds $l_0,\ldots,l_{d-1}$ and the upper bounds $u_0,\ldots,u_{d-1}$ are specified multi-indices of exact integers. The positive integer $d$ is the dimension of the interval. It is required that $l_0<u_0,\ldots,l_{d-1}<u_{d-1}$.

Intervals are a data type distinct from other Scheme data types.

Procedures

make-interval

make-interval arg1 #!optional arg2procedure

Create a new interval. arg1 and arg2 (if given) are nonempty vectors (of the same length) of exact integers.

If arg2 is not given, then the entries of arg1 must be positive, and they are taken as the upper-bounds of the interval, and lower-bounds is set to a vector of the same length with exact zero entries.

If arg2 is given, then arg1 is taken to be lower-bounds and arg2 is taken to be upper-bounds, which must satisfy

 (< (vector-ref lower-bounds i) (vector-ref upper-bounds i))

for $0\leq i<{}$(vector-length lower-bounds). It is an error if lower-bounds and upper-bounds do not satisfy these conditions.

interval?

interval? objprocedure: Returns #t if obj is an interval, and #f otherwise.

interval-dimension

interval-dimension intervalprocedure

If interval is an interval built with

(make-interval lower-bounds upper-bounds)

then interval-dimension returns (vector-length lower-bounds). It is an error to call interval-dimension if interval is not an interval.

Interval bounds

interval-lower-bound interval iprocedure

interval-upper-bound interval iprocedure

If interval is an interval built with

(make-interval lower-bounds upper-bounds)

and i is an exact integer that satisfies

$0 \leq i<$ (vector-length lower-bounds),

then interval-lower-bound returns (vector-ref lower-bounds i) and interval-upper-bound returns (vector-ref upper-bounds i). It is an error to call interval-lower-bound or interval-upper-bound if interval and i do not satisfy these conditions.

Interval bounds as lists

interval-lower-bounds->list intervalprocedure

interval-upper-bounds->list intervalprocedure

If interval is an interval built with

(make-interval lower-bounds upper-bounds)

then interval-lower-bounds->list returns (vector->list lower-bounds) and interval-upper-bounds->list returns (vector->list upper-bounds). It is an error to call interval-lower-bounds->list or interval-upper-bounds->list if interval does not satisfy these conditions.

Interval bounds as vectors

interval-lower-bounds->vector intervalprocedure

interval-upper-bounds->vector intervalprocedure

If interval is an interval built with

(make-interval lower-bounds upper-bounds)

then interval-lower-bounds->vector returns a copy of lower-bounds and interval-upper-bounds->vector returns a copy of upper-bounds. It is an error to call interval-lower-bounds->vector or interval-upper-bounds->vector if interval does not satisfy these conditions.

interval-volume

interval-volume intervalprocedure

If interval is an interval built with

(make-interval lower-bounds upper-bounds)

then, assuming the existence of vector-map, interval-volume returns

(apply * (vector->list (vector-map - upper-bounds lower-bounds)))

It is an error to call interval-volume if interval does not satisfy this condition.

interval=

interval= interval1 interval2procedure

If interval1 and interval2 are intervals built with

(make-interval lower-bounds1 upper-bounds1)

and

(make-interval lower-bounds2 upper-bounds2)

respectively, then interval= returns

(and (equal? lower-bounds1 lower-bounds2) (equal? upper-bounds1 upper-bounds2))

It is an error to call interval= if interval1 or interval2 do not satisfy this condition.

interval-subset?

interval-subset? interval1 interval2procedure

If interval1 and interval2 are intervals of the same dimension $d$, then interval-subset? returns #t if

(>= (interval-lower-bound interval1 j) (interval-lower-bound interval2 j))

and

(<= (interval-upper-bound interval1 j) (interval-upper-bound interval2 j))

for all $0\leq j<d$, otherwise it returns #f. It is an error if the arguments do not satisfy these conditions.

interval-contains-multi-index?

interval-contains-multi-index? interval index-0 index-1 ...procedure

If interval is an interval with dimension $d$ and index-0, index-1, ..., is a multi-index of length $d$, then interval-contains-multi-index? returns #t if

(interval-lower-bound interval j) $\leq$ index-j $<$ (interval-upper-bound interval j)

for $0\leq j < d$, and #f otherwise.

It is an error to call interval-contains-multi-index? if interval and index-0,..., do not satisfy this condition.

interval-projections

interval-projections interval right-dimensionprocedure

Conceptually, interval-projections takes a $d$-dimensional interval $[l_0,u_0)\times [l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ and splits it into two parts

$[l_0,u_0)\times\cdots\times[l_{d-\text{right-dimension}-1},u_{d-\text{right-dimension}-1})$

and

$[l_{d-\text{right-dimension}},u_{d-\text{right-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$

This function, the inverse of Cartesian products or cross products of intervals, is used to keep track of the domains of curried arrays.

More precisely, if interval is an interval and right-dimension is an exact integer that satisfies 0 < right-dimension < d then interval-projections returns two intervals:

(values
 (make-interval
  (vector (interval-lower-bound interval 0)
          ...
          (interval-lower-bound interval
                                (- d right-dimension 1)))
  (vector (interval-upper-bound interval 0)
          ...
          (interval-upper-bound interval
                                (- d right-dimension 1))))
 (make-interval
  (vector (interval-lower-bound interval
                                (- d right-dimension))
          ...
          (interval-lower-bound interval
                                (- d 1)))
  (vector (interval-upper-bound interval
                                (- d right-dimension))
          ...
          (interval-upper-bound interval
                                (- d 1)))))

It is an error to call interval-projections if its arguments do not satisfy these conditions.

interval-for-each

interval-for-each f intervalprocedure

This routine assumes that interval is an interval and f is a routine whose domain includes elements of interval. It is an error to call interval-for-each if interval and f do not satisfy these conditions.

interval-for-each calls f with each multi-index of interval as arguments, all in lexicographical order.

interval-dilate

interval-dilate interval lower-diffs upper-diffsprocedure

If interval is an interval with lower bounds $\ell_0,\dots,\ell_{d-1}$ and upper bounds $u_0,\dots,u_{d-1}$, and lower-diffs is a vector of exact integers $L_0,\dots,L_{d-1}$ and upper-diffs is a vector of exact integers $U_0,\dots,U_{d-1}$, then interval-dilate returns a new interval with lower bounds $\ell_0+L_0,\dots,\ell_{d-1}+L_{d-1}$ and upper bounds $u_0+U_0,\dots,u_{d-1}+U_{d-1}$, as long as this is a nonempty interval. It is an error if the arguments do not satisfy these conditions.

Examples:

(interval=
 (interval-dilate (make-interval '#(100 100))
                  '#(1 1) '#(1 1))
 (make-interval '#(1 1) '#(101 101))) => #t
(interval=
 (interval-dilate (make-interval '#(100 100))
                  '#(-1 -1) '#(1 1))
 (make-interval '#(-1 -1) '#(101 101))) => #t
(interval=
 (interval-dilate (make-interval '#(100 100))
                  '#(0 0) '#(-50 -50))
 (make-interval '#(50 50))) => #t
(interval-dilate
 (make-interval '#(100 100))
 '#(0 0) '#(-500 -50)) => error

interval-intersect

interval-intersect interval-1 interval-2 ...procedure

If all the arguments are intervals of the same dimension and they have a nonempty intersection, then interval-intersect returns that intersection; otherwise it returns #f.

It is an error if the arguments are not all intervals with the same dimension.

interval-translate

interval-translate interval translationprocedure

If interval is an interval with lower bounds $\ell_0,\dots,\ell_{d-1}$ and upper bounds $u_0,\dots,u_{d-1}$, and translation is a translation with entries $T_0,\dots,T_{d-1}$ , then interval-translate returns a new interval with lower bounds $\ell_0+T_0,\dots,\ell_{d-1}+T_{d-1}$ and upper bounds $u_0+T_0,\dots,u_{d-1}+T_{d-1}$. It is an error if the arguments do not satisfy these conditions.

One could define (interval-translate interval translation) by (interval-dilate interval translation translation).

interval-permute

interval-permute interval permutationprocedure

The argument interval must be an interval, and the argument permutation must be a valid permutation with the same dimension as interval. It is an error if the arguments do not satisfy these conditions.

Heuristically, this function returns a new interval whose axes have been permuted in a way consistent with permutation. But we have to say how the entries of permutation are associated with the new interval.

We have chosen the following convention: If the permutation is $(\pi_0,\ldots,\pi_{d-1})$, and the argument interval represents the cross product $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$, then the result represents the cross product $[l_{\pi_0},u_{\pi_0})\times[l_{\pi_1},u_{\pi_1})\times\cdots\times[l_{\pi_{d-1}},u_{\pi_{d-1}})$.

For example, if the argument interval represents $[0,4)\times[0,8)\times[0,21)\times [0,16)$ and the permutation is #(3 0 1 2), then the result of (interval-permute interval permutation) will be the representation of $[0,16)\times [0,4)\times[0,8)\times[0,21)$.

interval-rotate

interval-rotate interval dimprocedure

Informally, (interval-rotate interval dim) rotates the axes of interval dim places to the left.

More precisely, (interval-rotate interval dim) assumes that interval is an interval and dim is an exact integer between 0 (inclusive) and (interval-dimension interval) (exclusive). It computes the permutation (vector dim ... (- (interval-dimension interval) 1) 0 ... (- dim 1)) (unless dim is zero, in which case it constructs the identity permutation) and returns (interval-permute interval permutation). It is an error if the arguments do not satisfy these conditions.

interval-scale

interval-scale interval scalesprocedure

If interval is a $d$-dimensional interval $[0,u_1)\times\cdots\times[0,u_{d-1})$ with all lower bounds zero, and scales is a length-$d$ vector of positive exact integers, which we'll denote by $\vec s$, then interval-scale returns the interval $[0,\operatorname{ceiling}(u_1/s_1))\times\cdots\times[0,\operatorname{ceiling}(u_{d-1}/s_{d-1}))$.

It is an error if interval and scales do not satisfy this condition.

interval-cartesian-product

interval-cartesian-product interval #!rest intervalsprocedure

Implements the Cartesian product of the intervals in (cons interval intervals). Returns

(make-interval (list->vector (apply append (map array-lower-bounds->list (cons interval intervals))))
               (list->vector (apply append (map array-upper-bounds->list (cons interval intervals)))))

It is an error if any argument is not an interval.

Storage classes

Conceptually, a storage-class is a set of functions to manage the backing store of a specialized array. The functions allow one to make a backing store, to get values from the store and to set new values, to return the length of the store, and to specify a default value for initial elements of the backing store. Typically, a backing store is a (heterogeneous or homogeneous) vector. A storage-class has a type distinct from other Scheme types.

Procedures

make-storage-class

make-storage-class getter setter checker maker copier length defaultprocedure

Here we assume the following relationships between the arguments of make-storage-class. Assume that the "elements" of the backing store are of some "type", either heterogeneous (all Scheme types) or homogeneous (of some restricted type).

(maker n value) returns a linearly addressed object containing n elements of value value.
copier may be #f or a procedure; if a procedure then if to and from were created by maker, then (copier to at from start end) copies elements from from beginning at start (inclusive) and ending at end (exclusive) to to beginning at at. It is assumed that all the indices involved are within the domain of from and to, as needed. The order in which the elements are copied is unspecified.
If v is an object created by (maker n value) and 0 <== i < n, then (getter v i) returns the current value of the i'th element of v, and (checker (getter v i))> #t.
If v is an object created by (maker n value), 0 <== i < n, and (checker val)> #t, then (setter v i val) sets the value of the i'th element of v to val.
If v is an object created by (maker n value) then (length v) returns n.

If the arguments do not satisfy these conditions, then it is an error to call make-storage-class.

Note that we assume that getter and setter generally take O(1) time to execute.

storage-class?

storage-class? mprocedure: Returns #t if m is a storage class, and #f otherwise.

Storage class accessors

storage-class-getter mprocedure

storage-class-setter mprocedure

storage-class-checker mprocedure

storage-class-maker mprocedure

storage-class-copier mprocedure

storage-class-length mprocedure

storage-class-default mprocedure

If m is an object created by

(make-storage-class getter setter checker maker copier length default)

then storage-class-getter returns getter, storage-class-setter returns setter, storage-class-checker returns checker, storage-class-maker returns maker, storage-class-copier returns copier, storage-class-length returns length, and storage-class-default returns default. Otherwise, it is an error to call any of these routines.

Global Variables

generic-storage-classconstant

s8-storage-classconstant

s16-storage-classconstant

s32-storage-classconstant

s64-storage-classconstant

u1-storage-classconstant

u8-storage-classconstant

u16-storage-classconstant

u32-storage-classconstant

u64-storage-classconstant

f8-storage-classconstant

f16-storage-classconstant

f32-storage-classconstant

f64-storage-classconstant

c64-storage-classconstant

c128-storage-classconstant

generic-storage-class is defined as if by

(define generic-storage-class
  (make-storage-class vector-ref
                      vector-set!
                      (lambda (arg) #t)
                      make-vector
                      vector-copy!
                      vector-length
                      #f))

Implementations shall define sX-storage-class for X=8, 16, 32, and 64 (which have default values 0 and manipulate exact integer values between $-2^{x-1}$ and $2^{X-1}-1$ inclusive), uX-storage-class for X=1, 8, 16, 32, and 64 (which have default values 0 and manipulate exact integer values between 0 and $2^X-1$ inclusive), fX-storage-class for X= 8, 16, 32, and 64 (which have default value 0.0 and manipulate 8-, 16-, 32-, and 64-bit floating-point numbers), and cX-storage-class for X= 64 and 128 (which have default value 0.0+0.0i and manipulate complex numbers with, respectively, 32- and 64-bit floating-point numbers as real and imaginary parts).

Implementations with an appropriate homogeneous vector type should define the associated global variable using make-storage-class, otherwise they shall define the associated global variable to #f.

Arrays

Arrays are a data type distinct from other Scheme data types.

Procedures

make-array

make-array interval getter #!optional setterprocedure

Assume first that the optional argument setter is not given.

If interval is an interval and getter is a function from interval to Scheme objects, then make-array returns an array with domain interval and getter getter.

It is an error to call make-array if interval and getter do not satisfy these conditions.

If now setter is specified, assume that it is a procedure such that getter and setter satisfy: If

(i1,...,in) $\neq$ (j1,...,jn)

are elements of interval and

(getter j1 ... jn) => x

then "after"

 (setter v i1 ... in)

we have

(getter j1 ... jn) => x

and

(getter i1,...,in) => v

Then make-array builds a mutable array with domain interval, getter getter, and setter setter. It is an error to call make-array if its arguments do not satisfy these conditions.

Example:

  (define a (make-array (make-interval '#(1 1) '#(11 11))
                        (lambda (i j)
                          (if (= i j)
                              1
                              0))))

defines an array for which (array-getter a) returns 1 when i=j and 0 otherwise.

Example:

(define a   ;; a sparse array
  (let ((domain
         (make-interval '#(1000000 1000000)))
        (sparse-rows
         (make-vector 1000000 '())))
    (make-array
     domain
     (lambda (i j)
       (cond ((assv j (vector-ref sparse-rows i))
              => cdr)
             (else
              0.0)))
     (lambda (v i j)
       (cond ((assv j (vector-ref sparse-rows i))
              => (lambda (pair)
                   (set-cdr! pair v)))
             (else
              (vector-set!
               sparse-rows
               i
               (cons (cons j v)
                     (vector-ref sparse-rows i)))))))))

(define a_ (array-getter a))
(define a! (array-setter a))

(a_ 12345 6789)  => 0.
(a_ 0 0) => 0.
(a! 1.0 0 0) => undefined
(a_ 12345 6789)  => 0.
(a_ 0 0) => 1.

array?

array? objprocedure: Returns #t if obj is an array and #f otherwise.

array-domain

array-domain arrayprocedure

If array is an array built by

(make-array interval getter [setter])

(with or without the optional setter argument) then array-domain returns interval. It is an error to call array-domain if array is not an array.

array-getter

array-getter arrayprocedure

If array is an array built by

(make-array interval getter [setter])

(with or without the optional setter argument) array-getter returns getter. It is an error to call array-getter if array is not an array.

Example:

(define a (make-array (make-interval '#(1 1) '#(11 11))
                      (lambda (i j)
                        (if (= i j)
                            1
                            0))))
(define a_ (array-getter a))

(a_ 3 3) => 1
(a_ 2 3) => 0
(a_ 11 0) => is an error

array-dimension

array-dimension arrayprocedure: Shorthand for (interval-dimension (array-domain array)). It is an error to call this function if array is not an array.

mutable-array?

mutable-array? objprocedure: Returns #t if obj is a mutable array and #f otherwise.

array-setter

array-setter arrayprocedure

If array is an array built by

(make-array interval getter setter)

then array-setter returns setter. It is an error to call array-setter if array is not a mutable array.

specialized-array-default-safe?parameter

With no argument, returns #t if newly constructed specialized arrays check the arguments of setters and getters by default, and #f otherwise.

If bool is #t then the next call to specialized-array-default-safe? will return #t; if bool is #f then the next call to specialized-array-default-safe? will return #f; otherwise it is an error.

specialized-array-default-mutable? #!optional boolparameter

With no argument, returns #t if newly constructed specialized arrays are mutable by default, and #f otherwise.

If bool is #t then the next call to specialized-array-default-mutable? will return #t; if bool is #f then the next call to specialized-array-default-mutable? will return #f; otherwise it is an error.

make-specialized-array

(make-specialized-array interval [ storage-class generic-storage-class ] [ safe? (specialized-array-default-safe?) ])procedure

Constructs a mutable specialized array from its arguments.

interval must be given as a nonempty interval. If given, storage-class must be a storage class; if it is not given it defaults to generic-storage-class. If given, safe? must be a boolean; if it is not given it defaults to the current value of (specialized-array-default-safe?).

The body of the result is constructed as

  ((storage-class-maker storage-class)
   (interval-volume interval)
   (storage-class-default storage-class))

The indexer of the resulting array is constructed as the lexicographical mapping of interval onto the interval [0,(interval-volume interval)).

If safe is #t, then the arguments of the getter and setter (including the value to be stored) of the resulting array are always checked for correctness.

After correctness checking (if needed), (array-getter array) is defined simply as

  (lambda multi-index
    ((storage-class-getter storage-class)
     (array-body array)
     (apply (array-indexer array) multi-index)))

and (array-setter array) is defined as

  (lambda (val . multi-index)
    ((storage-class-setter storage-class)
     (array-body array)
     (apply (array-indexer array) multi-index)
     val))

It is an error if the arguments of make-specialized-array do not satisfy these conditions.

Examples. A simple array that can hold any type of element can be defined with (make-specialized-array (make-interval '#(3 3))). If you find that you're using a lot of unsafe arrays of unsigned 16-bit integers, one could define

  (define (make-u16-array interval)
    (make-specialized-array interval u16-storage-class #f))

and then simply call, e.g., (make-u16-array (make-interval '#(3 3))).

specialized-array?

specialized-array? objprocedure: Returns #t if obj is a specialized-array, and #f otherwise. A specialized-array is an array.

Specialized array accessors

array-storage-class arrayprocedure

array-indexer arrayprocedure

array-body arrayprocedure

array-safe? arrayprocedure

array-storage-class returns the storage-class of array. array-safe? is true if and only if the arguments of (array-getter array) and (array-setter array) (including the value to be stored in the array) are checked for correctness.

(array-body array) is a linearly indexed, vector-like object (e.g., a vector, string, u8vector, etc.) indexed from 0.

(array-indexer array) is assumed to be a one-to-one, but not necessarily onto, affine mapping from (array-domain array) into the indexing domain of (array-body array).

Please see make-specialized-array for how (array-body array), etc., are used.

It is an error to call any of these routines if array is not a specialized array.

array-elements-in-order?

array-elements-in-order? Aprocedure

Assumes that A is a specialized array, in which case it returns #t if the elements of A are in order and stored adjacently in (array-body A) and #f otherwise.

It is an error if A is not a specialized array.

specialized-array-share array new-domain new-domain->old-domainprocedure

Constructs a new specialized array that shares the body of the specialized array array. Returns an object that is behaviorally equivalent to a specialized array with the following fields:

domain:        new-domain
storage-class: (array-storage-class array)
body:          (array-body array)
indexer:       (lambda multi-index
                 (call-with-values
                     (lambda ()
                       (apply new-domain->old-domain
                              multi-index))
                   (array-indexer array)))

The resulting array inherits its safety and mutability from array.

Note: It is assumed that the affine structure of the composition of new-domain->old-domain and (array-indexer array) will be used to simplify:

  (lambda multi-index
    (call-with-values
        (lambda ()
          (apply new-domain->old-domain multi-index))
      (array-indexer array)))

It is an error if array is not a specialized array, or if new-domain is not an interval, or if new-domain->old-domain is not a one-to-one affine mapping from new-domain to (array-domain array).

Example: One can apply a "shearing" operation to an array as follows:

(define a
  (array-copy
   (make-array (make-interval '#(5 10))
               list)))
(define b
  (specialized-array-share
   a
   (make-interval '#(5 5))
   (lambda (i j)
     (values i (+ i j)))))
;; Print the "rows" of b
(array-for-each (lambda (row)
                  (pretty-print (array->list row)))
                (array-curry b 1))

;; which prints
;; ((0 0) (0 1) (0 2) (0 3) (0 4))
;; ((1 1) (1 2) (1 3) (1 4) (1 5))
;; ((2 2) (2 3) (2 4) (2 5) (2 6))
;; ((3 3) (3 4) (3 5) (3 6) (3 7))
;; ((4 4) (4 5) (4 6) (4 7) (4 8))

This "shearing" operation cannot be achieved by combining the procedures array-extract, array-translate, array-permute, array-translate, array-curry, array-reverse, and array-sample.

array-copy

(array-copy array [ result-storage-class generic-storage-class ] [ new-domain #f ] [ mutable? (specialized-array-default-mutable?) ] [ safe? (specialized-array-default-safe?) ])procedure

Assumes that array is an array, result-storage-class is a storage class that can manipulate all the elements of array, new-domain is either #f or an interval with the same volume as (array-domain array), and mutable? and safe? are booleans.

If new-domain is #f, then it is set to (array-domain array).

The specialized array returned by array-copy can be defined conceptually by:

(list->array (array->list array)
             new-domain
             result-storage-class
             mutable?
             safe?)

It is an error if the arguments do not satisfy these conditions.

Note: If new-domain is not the same as (array-domain array), one can think of the resulting array as a reshaped version of array.

array-curry

array-curry array inner-dimensionprocedure

If array is an array whose domain is an interval $[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$, and inner-dimension is an exact integer strictly between $0$ and $d$, then array-curry returns an immutable array with domain $[l_0,u_0)\times\cdots\times[l_{d-\text{inner-dimension}-1},u_{d-\text{inner-dimension}-1})$, each of whose entries is in itself an array with domain $[l_{d-\text{inner-dimension}},u_{d-\text{inner-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$.

For example, if A and B are defined by

(define interval (make-interval '#(10 10 10 10)))
(define A (make-array interval list))
(define B (array-curry A 1))

(define A_ (array-getter A))
(define B_ (array-getter B))

(A_ i j k l) => (list i j k l)

then B is an immutable array with domain (make-interval '#(10 10 10)), each of whose elements is itself an (immutable) array and

(equal?
 (A_ i j k l)
 ((array-getter (B_ i j k)) l)) => #t

for all multi-indices i j k l in interval.

The subarrays are immutable, mutable, or specialized according to whether the array argument is immutable, mutable, or specialized.

More precisely, if

0 < inner-dimension < (interval-dimension (array-domain array))

then array-curry returns a result as follows.

If the input array is specialized, then array-curry returns

(call-with-values
    (lambda () (interval-projections (array-domain array)
                                     inner-dimension))
  (lambda (outer-interval inner-interval)
    (make-array
     outer-interval
     (lambda outer-multi-index
       (specialized-array-share
        array
        inner-interval
        (lambda inner-multi-index
          (apply values
                 (append outer-multi-index
                         inner-multi-index))))))))

Otherwise, if the input array is mutable, then array-curry returns

(call-with-values
    (lambda () (interval-projections (array-domain array)
                                     inner-dimension))
  (lambda (outer-interval inner-interval)
    (make-array
     outer-interval
     (lambda outer-multi-index
       (make-array
        inner-interval
        (lambda inner-multi-index
          (apply (array-getter array)
                 (append outer-multi-index
                         inner-multi-index)))
        (lambda (v . inner-multi-index)
          (apply (array-setter array)
                 v
                 (append outer-multi-index
                         inner-multi-index))))))))

Otherwise, array-curry returns

(call-with-values
    (lambda () (interval-projections (array-domain array)
                                     inner-dimension))
  (lambda (outer-interval inner-interval)
    (make-array
     outer-interval
     (lambda outer-multi-index
       (make-array
        inner-interval
        (lambda inner-multi-index
          (apply (array-getter array)
                 (append outer-multi-index
                         inner-multi-index))))))))

It is an error to call array-curry if its arguments do not satisfy these conditions.

If array is a specialized array, the subarrays of the result inherit their safety and mutability from array.

Note: Let's denote by B the result of (array-curry A k). While the result of calling (array-getter B) is an immutable, mutable, or specialized array according to whether A itself is immutable, mutable, or specialized, B is always an immutable array, where (array-getter B), which returns an array, is computed anew for each call. If (array-getter B) will be called multiple times with the same arguments, it may be useful to store these results in a specialized array for fast repeated access.

Please see the note in the discussion of array-tile.

Example:

(define a (make-array (make-interval '#(10 10))
                      list))
(define a_ (array-getter a))
(a_ 3 4)  => (3 4)
(define curried-a (array-curry a 1))
(define curried-a_ (array-getter curried-a))
((array-getter (curried-a_ 3)) 4)
                    => (3 4)

array-extract

array-extract array new-domainprocedure

Returns a new array with the same getter (and setter, if appropriate) of the first argument, defined on the second argument.

Assumes that array is an array and new-domain is an interval that is a sub-interval of (array-domain array). If array is a specialized array, then returns

  (specialized-array-share array
                           new-domain
                           values)

Otherwise, if array is a mutable array, then array-extract returns

  (make-array new-domain
              (array-getter array)
              (array-setter array))

Finally, if array is an immutable array, then array-extract returns

  (make-array new-domain
              (array-getter array))

It is an error if the arguments of array-extract do not satisfy these conditions.

If array is a specialized array, the resulting array inherits its mutability and safety from array.

array-tile

array-tile A Sprocedure

Assume that A is an array and S is a vector of positive, exact integers. The routine array-tile returns a new immutable array $T$, each entry of which is a subarray of A whose domain has sidelengths given (mostly) by the entries of S. These subarrays completely "tile" A, in the sense that every entry in A is an entry of precisely one entry of the result $T$.

More formally, if S is the vector $(s_0,\ldots,s_{d-1})$, and the domain of A is the interval $[l_0,u_0)\times\cdots\times [l_{d-1},u_{d-1})$, then $T$ is an immutable array with all lower bounds zero and upper bounds given by $$ \operatorname{ceiling}((u_0-l_0)/s_0),\ldots,\operatorname{ceiling}((u_{d-1}-l_{d-1})/s_{d-1}). $$ The $i_0,\ldots,i_{d-1}$ entry of $T$ is (array-extract A D_i) with the interval D_i given by $$ [l_0+i_0*s_0,\min(l_0+(i_0+1)s_0,u_0))\times\cdots\times[l_{d-1}+i_{d-1}*s_{d-1},\min(l_{d-1}+(i_{d-1}+1)s_{d-1},u_{d-1})). $$ (The "minimum" operators are necessary if $u_j-l_j$ is not divisible by $s_j$.) Thus, each entry of $T$ will be a specialized, mutable, or immutable array, depending on the type of the input array A.

It is an error if the arguments of array-tile do not satisfy these conditions.

If A is a specialized array, the subarrays of the result inherit safety and mutability from A.

Note: The routines array-tile and array-curry both decompose an array into subarrays, but in different ways. For example, if A is defined as (make-array (make-interval '#(10 10)) list), then (array-tile A '#(1 10)) returns an array with domain (make-interval '#(10 1)), each element of which is an array with domain (make-interval '#(1 10)) (i.e., a two-dimensional array whose elements are two-dimensional arrays), while (array-curry A 1) returns an array with domain (make-interval '#(10)), each element of which has domain (make-interval '#(10)) (i.e., a one-dimensional array whose elements are one-dimensional arrays).

array-translate

array-translate array translationprocedure

Assumes that array is a valid array, translation is a valid translation, and that the dimensions of the array and the translation are the same. The resulting array will have domain (interval-translate (array-domain array) translation).

If array is a specialized array, returns a new specialized array

(specialized-array-share
 array
 (interval-translate (array-domain array)
                     translation)
 (lambda multi-index
   (apply values
          (map -
               multi-index
               (vector->list translation)))))

that shares the body of array, as well as inheriting its safety and mutability.

If array is not a specialized array but is a mutable array, returns a new mutable array

(make-array
 (interval-translate (array-domain array)
                     translation)
 (lambda multi-index
   (apply (array-getter array)
          (map -
               multi-index
               (vector->list translation))))
 (lambda (val . multi-index)
   (apply (array-setter array)
          val
          (map -
               multi-index
               (vector->list translation)))))

that employs the same getter and setter as the original array argument.

If array is not a mutable array, returns a new array

(make-array
 (interval-translate (array-domain array)
                     translation)
 (lambda multi-index
   (apply (array-getter array)
          (map - multi-index (vector->list translation)))))

that employs the same getter as the original array.

It is an error if the arguments do not satisfy these conditions.

array-permute

array-permute array permutationprocedure

Assumes that array is a valid array, permutation is a valid permutation, and that the dimensions of the array and the permutation are the same. The resulting array will have domain (interval-permute (array-domain array) permutation).

We begin with an example. Assume that the domain of array is represented by the interval $[0,4)\times[0,8)\times[0,21)\times [0,16)$, as in the example for interval-permute, and the permutation is #(3 0 1 2). Then the domain of the new array is the interval $[0,16)\times [0,4)\times[0,8)\times[0,21)$.

So the multi-index argument of the getter of the result of array-permute must lie in the new domain of the array, the interval $[0,16)\times [0,4)\times[0,8)\times[0,21)$. So if we define old-getter as (array-getter array), the definition of the new array must be in fact

(make-array (interval-permute (array-domain array)
                              '#(3 0 1 2))
            (lambda (l i j k)
              (old-getter i j k l)))

So you see that if the first argument if the new getter is in $[0,16)$, then indeed the fourth argument of old-getter is also in $[0,16)$, as it should be. This is a subtlety that I don't see how to overcome. It is the listing of the arguments of the new getter, the lambda, that must be permuted.

Mathematically, we can define $\pi^{-1}$, the inverse of a permutation $\pi$, such that $\pi^{-1}$ composed with $\pi$ gives the identity permutation. Then the getter of the new array is, in pseudo-code, (lambda multi-index (apply old-getter ($\pi^{-1}$ multi-index))). We have assumed that $\pi^{-1}$ takes a list as an argument and returns a list as a result.

Employing this same pseudo-code, if array is a specialized array and we denote the permutation by $\pi$, then array-permute returns the new specialized array

(specialized-array-share array
                         (interval-permute (array-domain array) π)
                         (lambda multi-index
                           (apply values (π-1 multi-index))))

The resulting array shares the body of array, as well as its safety and mutability.

Again employing this same pseudo-code, if array is not a specialized array, but is a mutable-array, then array-permute returns the new mutable

(make-array (interval-permute (array-domain array) π)
            (lambda multi-index
              (apply (array-getter array)
                     (π-1 multi-index)))
            (lambda (val . multi-index)
              (apply (array-setter array)
                     val
                     (π-1 multi-index))))

which employs the setter and the getter of the argument to array-permute.

Finally, if array is not a mutable array, then array-permute returns

(make-array (interval-permute (array-domain array) π)
            (lambda multi-index
              (apply (array-getter array)
                     (π-1 multi-index))))

It is an error to call array-permute if its arguments do not satisfy these conditions.

array-rotate

array-rotate array dimprocedure

Informally, (array-rotate array dim) rotates the axes of array dim places to the left.

More precisely, (array-rotate array dim) assumes that array is an array and dim is an exact integer between 0 (inclusive) and (array-dimension array) (exclusive). It computes the permutation (vector dim ... (- (array-dimension array) 1) 0 ... (- dim 1)) (unless dim is zero, in which case it constructs the identity permutation) and returns (array-permute array permutation). It is an error if the arguments do not satisfy these conditions.

array-reverse

array-reverse array #!optional flip?procedure

We assume that array is an array and flip?, if given, is a vector of booleans whose length is the same as the dimension of array. If flip? is not given, it is set to a vector with length the same as the dimension of array, all of whose elements are #t.

array-reverse returns a new array that is specialized, mutable, or immutable according to whether array is specialized, mutable, or immutable, respectively. Informally, if (vector-ref flip? k) is true, then the ordering of multi-indices in the k'th coordinate direction is reversed, and is left undisturbed otherwise.

More formally, we introduce the function

(define flip-multi-index
  (let* ((domain (array-domain array))
         (lowers (interval-lower-bounds->list domain))
         (uppers (interval-upper-bounds->list domain)))
    (lambda (multi-index)
      (map (lambda (i_k flip?_k l_k u_k)
             (if flip?_k
                 (- (+ l_k u_k -1) i_k)
                 i_k))
           multi-index
           (vector->list flip?)
           lowers
           uppers))))

Then if array is specialized, array-reverse returns

(specialized-array-share
 array
 domain
 (lambda multi-index
   (apply values
          (flip-multi-index multi-index))))

and the result inherits the safety and mutability of array.

Otherwise, if array is mutable, then array-reverse returns

(make-array
 domain
 (lambda multi-index
   (apply (array-getter array)
          (flip-multi-index multi-index)))
   (lambda (v . multi-index)
     (apply (array-setter array)
            v
            (flip-multi-index multi-index)))))

Finally, if array is immutable, then array-reverse returns

(make-array
 domain
 (lambda multi-index
   (apply (array-getter array)
          (flip-multi-index multi-index)))))

It is an error if array and flip? don't satisfy these requirements.

The following example using array-reverse was motivated by a blog post by Joe Marshall.

(define (palindrome? s)
  (let ((n (string-length s)))
    (or (< n 2)
        (let* ((a
                ;; an array accessing the characters of s
                (make-array (make-interval (vector n))
                            (lambda (i)
                              (string-ref s i))))
               (ra
                ;; the array accessed in reverse order
                (array-reverse a))
               (half-domain
                (make-interval (vector (quotient n 2)))))
          (array-every
           char=?
           ;; the first half of s
           (array-extract a half-domain)
           ;; the reversed second half of s
           (array-extract ra half-domain))))))

(palindrome? "") => #t
(palindrome? "a") => #t
(palindrome? "aa") => #t
(palindrome? "ab") => #f
(palindrome? "aba") => #t
(palindrome? "abc") => #f
(palindrome? "abba") => #t
(palindrome? "abca") => #f
(palindrome? "abbc") => #f

array-sample

array-sample array scalesprocedure

We assume that array is an array all of whose lower bounds are zero, and scales is a vector of positive exact integers whose length is the same as the dimension of array.

Informally, if we construct a new matrix $S$ with the entries of scales on the main diagonal, then the $\vec i$th element of (array-sample array scales) is the $S\vec i$th element of array.

More formally, if array is specialized, then array-sample returns

(specialized-array-share
 array
 (interval-scale (array-domain array)
                 scales)
 (lambda multi-index
   (apply values
          (map * multi-index (vector->list scales)))))

with the result inheriting the safety and mutability of array.

Otherwise, if array is mutable, then array-sample returns

(make-array
 (interval-scale (array-domain array)
                 scales)
 (lambda multi-index
   (apply (array-getter array)
          (map * multi-index (vector->list scales))))
 (lambda (v . multi-index)
   (apply (array-setter array)
          v
          (map * multi-index (vector->list scales)))))

Finally, if array is immutable, then array-sample returns

(make-array
 (interval-scale (array-domain array)
                 scales)
 (lambda multi-index
   (apply (array-getter array)
          (map * multi-index (vector->list scales)))))

It is an error if array and scales don't satisfy these requirements.

array-outer-product

array-outer-product op array1 array2procedure

Implements the outer product of array1 and array2 with the operator op, similar to the APL function with the same name.

Assume that array1 and array2 are arrays and that op is a function of two arguments. Assume that (list-tail l n) returns the list remaining after the first n items of the list l have been removed, and (list-take l n) returns a new list consisting of the first n items of the list l. Then array-outer-product returns the immutable array

(make-array (interval-cartesian-product (array-domain array1)
                                        (array-domain array2))
            (lambda args
              (op (apply (array-getter array1) (list-take args (array-dimension array1)))
                  (apply (array-getter array2) (list-tail args (array-dimension array1))))))

This operation can be considered a partial inverse to array-curry. It is an error if the arguments do not satisfy these assumptions.

Note: You can see from the above definition that if C is (array-outer-product op A B), then each call to (array-getter C) will call op as well as (array-getter A) and (array-getter B). This implies that if all elements of C are eventually accessed, then (array-getter A) will be called (array-volume B) times; similarly (array-getter B) will be called (array-volume A) times.

This implies that if (array-getter A) is expensive to compute (for example, if it's returning an array, as does array-curry) then the elements of A should be pre-computed if necessary and stored in a specialized array, typically using array-copy, before that specialized array is passed as an argument to array-outer-product. In the examples below, the code for Gaussian elimination applies array-outer-product to shared specialized arrays, which are of course themselves specialized arrays; the code for matrix multiplication and inner-product applies array-outer-product to curried arrays, so we apply array-copy to the arguments before passage to array-outer-product.

array-map

array-map f array #!rest arraysprocedure

If array, (car arrays), ... all have the same domain and f is a procedure, then array-map returns a new immutable array with the same domain and getter

(lambda multi-index
  (apply f
         (map (lambda (g)
                (apply g multi-index))
              (map array-getter
                   (cons array arrays)))))

It is assumed that f is appropriately defined to be evaluated in this context.

It is expected that array-map and array-for-each will specialize the construction of

(lambda multi-index
  (apply f
         (map (lambda (g)
                (apply g multi-index))
              (map array-getter
                   (cons array
                         arrays)))))

It is an error to call array-map if its arguments do not satisfy these conditions.

Note: The ease of constructing temporary arrays without allocating storage makes it trivial to imitate, e.g., Javascript's map with index. For example, we can add the index to each element of an array a by

(array-map +
           a
           (make-array (array-domain a)
                       (lambda (i) i)))

or even

(make-array (array-domain a)
            (let ((a_ (array-getter a)))
              (lambda (i)
                (+ (a_ i) i))))

array-for-each

array-for-each f array #!rest arraysprocedure

If array, (car arrays), ... all have the same domain and f is an appropriate procedure, then array-for-each calls

(interval-for-each
 (lambda multi-index
   (apply f
          (map (lambda (g)
                 (apply g multi-index))
               (map array-getter
                    (cons array
                          arrays)))))
 (array-domain array))

In particular, array-for-each always walks the indices of the arrays in lexicographical order.

It is expected that array-map and array-for-each will specialize the construction of

(lambda multi-index
  (apply f
         (map (lambda (g)
                (apply g multi-index))
              (map array-getter
                   (cons array
                         arrays)))))

It is an error to call array-for-each if its arguments do not satisfy these conditions.

array-fold

array-fold kons knil arrayprocedure

If we use the defining relations for fold over lists from SRFI 1:

(fold kons knil lis)
    = (fold kons (kons (car lis) knil) (cdr lis))
(fold kons knil '())
    = knil

then (array-fold kons knil array) returns

(fold kons knil (array->list array))

It is an error if array is not an array, or if kons is not a procedure.

array-fold-right

array-fold-right kons knil arrayprocedure

If we use the defining relations for fold-right over lists from SRFI 1:

(fold-right kons knil lis)
    = (kons (car lis) (fold-right kons knil (cdr lis)))
(fold-right kons knil '())
    = knil

then (array-fold-right kons knil array) returns

(fold-right kons knil (array->list array))

It is an error if array is not an array, or if kons is not a procedure.

array-reduce

array-reduce op Aprocedure

We assume that A is an array and op is a procedure of two arguments that is associative, i.e., (op (op x y) z) is the same as (op x (op y z)).

Then (array-reduce op A) returns

(let ((box '())
      (A_ (array-getter A)))
  (interval-for-each
   (lambda args
     (if (null? box)
         (set! box (list (apply A_ args)))
         (set-car! box (op (car box)
                           (apply A_ args)))))
   (array-domain A))
  (car box))

The implementation is allowed to use the associativity of op to reorder the computations in array-reduce. It is an error if the arguments do not satisfy these conditions.

As an example, we consider the finite sum: $$ S_m=\sum_{k=1}^m \frac 1{k^2}. $$ One can show that $$ \frac 1 {m+1}<\frac{\pi^2}6-S_m<\frac 1m. $$ We attempt to compute this in floating-point arithmetic in two ways. In the first, we apply array-reduce to an array containing the terms of the series, basically a serial computation. In the second, we divide the series into blocks of no more than 1,000 consecutive terms, apply array-reduce to get a new sequence of terms, and repeat the process. The second way is approximately what might happen with GPU computing.

We find with $m=1{,}000{,}000{,}000$:

(define A (make-array (make-interval '#(1) '#(1000000001))
                      (lambda (k)
                        (fl/ (flsquare (inexact k))))))
(define (block-sum A)
  (let ((N (interval-volume (array-domain A))))
    (cond ((<= N 1000)
           (array-reduce fl+ A))
          ((<= N (square 1000))
           (block-sum (array-map block-sum
                                 (array-tile A (vector (integer-sqrt N))))))
          (else
           (block-sum (array-map block-sum
                                 (array-tile A (vector (quotient N 1000)))))))))
(array-reduce fl+ A) => 1.644934057834575
(block-sum A)        => 1.6449340658482325

Since $\pi^2/6\approx{}$1.6449340668482264, we see using the first method that the difference $\pi^2/6-{}$1.644934057834575${}\approx{}$9.013651380840315e-9 and with the second we have $\pi^2/6-{}$1.6449340658482325${}\approx{}$9.99993865491433e-10. The true difference should be between $\frac 1{1{,}000{,}000{,}001}\approx{}$9.99999999e-10 and $\frac 1{1{,}000{,}000{,}000}={}$1e-9. The difference for the first method is about 10 times too big, and, in fact, will not change further because any further terms, when added to the partial sum, are too small to increase the sum after rounding-to-nearest in double-precision IEEE-754 floating-point arithmetic.

array-any

array-any pred array1 array2 ...procedure

Assumes that array1, array2, etc., are arrays, all with the same domain, which we'll call interval. Also assumes that pred is a procedure that takes as many arguments as there are arrays and returns a single value.

array-any first applies (array-getter array1), etc., to the first element of interval in lexicographical order, to which values it then applies pred.

If the result of pred is not #f, then that result is returned by array-any. If the result of pred is #f, then array-any continues with the second element of interval, etc., returning the first nonfalse value of pred.

If pred always returns #f, then array-any returns #f.

If it happens that pred is applied to the results of applying (array-getter array1), etc., to the last element of interval, then this last call to pred is in tail position.

The functions (array-getter array1), etc., are applied only to those values of interval necessary to determine the result of array-any.

It is an error if the arguments do not satisfy these assumptions.

array-every

array-every pred array1 array2 ...procedure

array-every first applies (array-getter array1), etc., to the first element of interval in lexicographical order, to which values it then applies pred.

If the result of pred is #f, then that result is returned by array-every. If the result of pred is nonfalse, then array-every continues with the second element of interval, etc., returning the first value of pred that is #f.

If pred always returns a nonfalse value, then the last nonfalse value returned by pred is also returned by array-every.

If it happens that pred is applied to the results of applying (array-getter array1), etc., to the last element of interval, then this last call to pred is in tail position.

The functions (array-getter array1), etc., are applied only to those values of interval necessary to determine the result of array-every.

It is an error if the arguments do not satisfy these assumptions.

array->list

array->list arrayprocedure

Stores the elements of array into a newly allocated list in lexicographical order. It is an error if array is not an array.

It is guaranteed that (array-getter array) is called precisely once for each multi-index in (array-domain array) in lexicographical order.

list->array

(list->array l domain [ result-storage-class generic-storage-class ] [ mutable? (specialized-array-default-mutable?) ] [ safe? (specialized-array-default-safe?) ])procedure

Assumes that l is an list, domain is an interval with volume the same as the length of l, result-storage-class is a storage class that can manipulate all the elements of l, and mutable? and safe? are booleans.

Returns a specialized array with domain domain whose elements are the elements of the list l stored in lexicographical order. The result is mutable or safe depending on the values of mutable? and safe?.

It is an error if the arguments do not satisfy these assumptions, or if any element of l cannot be stored in the body of result-storage-class, and this last error shall be detected and raised.

array-assign!

array-assign! destination sourceprocedure

Assumes that destination is a mutable array and source is an array, and that the elements of source can be stored into destination.

The array destination must be compatible with source, in the sense that either destination and source have the same domain, or destination is a specialized array whose elements are stored adjacently and in order in its body and whose domain has the same volume as the domain of source.

Evaluates (array-getter source) on the multi-indices in (array-domain source) in lexicographical order, and assigns each value to the multi-index in destination in the same lexicographical order.

It is an error if the arguments don't satisfy these assumptions.

If assigning any element of destination affects the value of any element of source, then the result is undefined.

Note: If the domains of destination and source are not the same, one can think of destination as a reshaped copy of source.

array-ref

array-ref A i0 #!rest i-tailprocedure

Assumes that A is an array, and every element of (cons i0 i-tail) is an exact integer.

Returns (apply (array-getter A) i0 i-tail).

It is an error if A is not an array, or if the number of arguments specified is not the correct number for (array-getter A).

array-set!

array-set! A v i0 #!rest i-tailprocedure

Assumes that A is a mutable array, that v is a value that can be stored within that array, and that every element of (cons i0 i-tail) is an exact integer.

Returns (apply (array-setter A) v i0 i-tail).

It is an error if A is not a mutable array, if v is not an appropriate value to be stored in that array, or if the number of arguments specified is not the correct number for (array-setter A).

Note: In the sample implementation, because array-ref and array-set! take a variable number of arguments and they must check that A is an array of the appropriate type, programs written in a style using these functions, rather than the style in which 1D-Haar-loop is coded below, can take up to three times as long runtime.

Note: In the sample implementation, checking whether the multi-indices are exact integers and within the domain of the array, and checking whether the value is appropriate for storage into the array, is delegated to the underlying definition of the array argument. If the first argument is a safe specialized array, then these items are checked; if it is an unsafe specialized array, they are not. If it is a generalized array, it is up to the programmer whether to define the getter and setter of the array to check the correctness of the arguments.

specialized-array-reshape

(specialized-array-reshape array new-domain [ copy-on-failure? #f ])procedure

Assumes that array is a specialized array, new-domain is an interval with the same volume as (array-domain array), and copy-on-failure?, if given, is a boolean.

If there is an affine map that takes the multi-indices in new-domain to the cells in (array-body array) storing the elements of array in lexicographical order, returns a new specialized array, with the same body and elements as array and domain new-domain. The result inherits its mutability and safety from array.

If there is not an affine map that takes the multi-indices in new-domain to the cells storing the elements of array in lexicographical order and copy-on-failure? is #t, then returns a specialized array copy of array with domain new-domain, storage class (array-storage-class array), mutability (mutable-array? array), and safety (array-safe? array).

It is an error if these conditions on the arguments are not met.

Note: The code in the sample implementation to determine whether there exists an affine map from new-domain to the multi-indices of the elements of array in lexicographical order is modeled on the corresponding code in the Python library NumPy.

Note: In the sample implementation, if an array cannot be reshaped and copy-on-failure? is #f, an error is raised in tail position. An implementation might want to replace this error call with a continuable exception to give the programmer more flexibility.

Examples: Reshaping an array is not a Bawden-type array transform. For example, we use array-display defined below to see:

;;; The entries of A are the multi-indices of the locations

(define A (array-copy (make-array (make-interval '#(3 4)) list)))

(array-display A)

;;; Displays

;;; (0 0)   (0 1)   (0 2)   (0 3)
;;; (1 0)   (1 1)   (1 2)   (1 3)
;;; (2 0)   (2 1)   (2 2)   (2 3)

(array-display (array-rotate A 1))

;;; Displays

;;; (0 0)   (1 0)   (2 0)
;;; (0 1)   (1 1)   (2 1)
;;; (0 2)   (1 2)   (2 2)
;;; (0 3)   (1 3)   (2 3)

(array-display (specialized-array-reshape A (make-interval '#(4 3))))

;;; Displays

;;; (0 0)   (0 1)   (0 2)
;;; (0 3)   (1 0)   (1 1)
;;; (1 2)   (1 3)   (2 0)
;;; (2 1)   (2 2)   (2 3)

(define B (array-sample A '#(2 1)))

(array-display B)

;;; Displays

;;; (0 0)   (0 1)   (0 2)   (0 3)
;;; (2 0)   (2 1)   (2 2)   (2 3)

(array-display (specialized-array-reshape B (make-interval '#(8)))) => fails

(array-display (specialized-array-reshape B (make-interval '#(8)) #t))

;;; Displays

;;; (0 0)   (0 1)   (0 2)   (0 3)   (2 0)   (2 1)   (2 2)   (2 3)

The following examples succeed:

(specialized-array-reshape
 (array-copy (make-array (make-interval '#(2 1 3 1)) list))
 (make-interval '#(6)))
(specialized-array-reshape
 (array-copy (make-array (make-interval '#(2 1 3 1)) list))
 (make-interval '#(3 2)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)))
 (make-interval '#(6)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)))
 (make-interval '#(3 2)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t))
 (make-interval '#(3 2)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t))
 (make-interval '#(3 1 2 1)))
(specialized-array-reshape
 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list)) '#(#f #f #f #t)) '#(1 1 2 1))
 (make-interval '#(4)))
(specialized-array-reshape
 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list)) '#(#t #f #t #t)) '#(1 1 2 1))
 (make-interval '#(4)))

The following examples raise an exception:

(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#t #f #f #f))
 (make-interval '#(6)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#t #f #f #f))
 (make-interval '#(3 2)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #t #f))
 (make-interval '#(6)))
(specialized-array-reshape
 (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #t #t))
 (make-interval '#(3 2)))
(specialized-array-reshape
 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 3 1)) list)) '#(#f #f #f #t)) '#(1 1 2 1))
 (make-interval '#(4)) )
(specialized-array-reshape
 (array-sample (array-reverse (array-copy (make-array (make-interval '#(2 1 4 1)) list)) '#(#f #f #t #t)) '#(1 1 2 1))
 (make-interval '#(4)))

In the next examples, we start with vector fields, $100\times 100$ arrays of 4-vectors. In one example, we reshape each large array to $100\times 100\times2\times2$ vector fields (so we consider each 4-vector as a $2\times 2$ matrix), and multiply the $2\times 2$ matrices together. In the second example, we reshape each 4-vector to a $2\times 2$ matrix individually, and compare the times.

(define interval-flat (make-interval '#(100 100 4)))

(define interval-2x2  (make-interval '#(100 100 2 2)))

(define A (array-copy (make-array interval-flat (lambda args (random-integer 5)))))

(define B (array-copy (make-array interval-flat (lambda args (random-integer 5)))))

(define C (array-copy (make-array interval-flat (lambda args 0))))

(define (2x2-matrix-multiply-into! A B C)
  (let ((C! (array-setter C))
        (A_ (array-getter A))
        (B_ (array-getter B)))
    (C! (+ (* (A_ 0 0) (B_ 0 0))
           (* (A_ 0 1) (B_ 1 0)))
        0 0)
    (C! (+ (* (A_ 0 0) (B_ 0 1))
           (* (A_ 0 1) (B_ 1 1)))
        0 1)
    (C! (+ (* (A_ 1 0) (B_ 0 0))
           (* (A_ 1 1) (B_ 1 0)))
        1 0)
    (C! (+ (* (A_ 1 0) (B_ 0 1))
           (* (A_ 1 1) (B_ 1 1)))
        1 1)))

;;; Reshape A, B, and C to change all the 4-vectors to 2x2 matrices

(time
 (array-for-each 2x2-matrix-multiply-into!
                 (array-curry (specialized-array-reshape A interval-2x2) 2)
                 (array-curry (specialized-array-reshape B interval-2x2) 2)
                 (array-curry (specialized-array-reshape C interval-2x2) 2)))
;;; Displays

;;;    0.015186 secs real time
;;;    0.015186 secs cpu time (0.015186 user, 0.000000 system)
;;;    4 collections accounting for 0.004735 secs real time (0.004732 user, 0.000000 system)
;;;    46089024 bytes allocated
;;;    no minor faults
;;;    no major faults

;;; Reshape each 4-vector to a 2x2 matrix individually

(time
 (array-for-each (lambda (A B C)
                   (2x2-matrix-multiply-into!
                    (specialized-array-reshape A (make-interval '#(2 2)))
                    (specialized-array-reshape B (make-interval '#(2 2)))
                    (specialized-array-reshape C (make-interval '#(2 2)))))
                 (array-curry A 1)
                 (array-curry B 1)
                 (array-curry C 1)))

;;; Displays

;;;    0.039193 secs real time
;;;    0.039193 secs cpu time (0.039191 user, 0.000002 system)
;;;    6 collections accounting for 0.006855 secs real time (0.006851 user, 0.000001 system)
;;;    71043024 bytes allocated
;;;    no minor faults
;;;    no major faults

Implementation

We provide an implementation in Gambit Scheme; the nonstandard techniques used in the implementation are: DSSSL-style optional and keyword arguments; a unique object to indicate absent arguments; define-structure; and define-macro.

There is a git repository of this document, a sample implementation, a test file, and other materials.

Relationship to other SRFIs

Final SRFIs 25, 47, 58, and 63 deal with "Multi-dimensional Array Primitives", "Array", "Array Notation", and "Homogeneous and Heterogeneous Arrays", respectively. Each of these previous SRFIs deal with what we call in this SRFI specialized arrays. Many of the functions in these previous SRFIs have corresponding forms in this SRFI. For example, from SRFI 63, we can translate:

(array? obj)

: (array? obj)

(array-rank A)

: (array-dimension A)

(make-array prototype k1 ...)

: (make-specialized-array (make-interval (vector k1 ...)) storage-class).

(make-shared-array A mapper k1 ...)

: (specialized-array-share A (make-interval (vector k1 ...)) mapper)

(array-in-bounds? A index1 ...)

: (interval-contains-multi-index? (array-domain A) index1 ...)

(array-ref A k1 ...)

: (let ((A_ (array-getter A))) ... (A_ k1 ...) ... ) or (array-ref A k1 ...)

(array-set! A obj k1 ...)

: (let ((A! (array-setter A))) ... (A! obj k1 ...) ...) or (array-set! A obj k1 ...)

At the same time, this SRFI has some special features:

Intervals, used as the domains of arrays in this SRFI, are useful objects in their own rights, with their own procedures. We make a sharp distinction between the domains of arrays and the arrays themselves.
Intervals can have nonzero lower bounds in each dimension.
Intervals cannot be empty.
Arrays must have a getter, but may have no setter.
There are many predefined array transformations: array-extract, array-tile, array-translate, array-permute, array-rotate, array-sample, array-reverse.

Other examples

Image processing applications provided significant motivation for this SRFI.

Manipulating images in PGM format. On a system with eight-bit chars, one can write routines to read and write greyscale images in the PGM format of the netpbm package as follows. The lexicographical order in array-copy guarantees the the correct order of execution of the input procedures:

(define make-pgm   cons)
(define pgm-greys  car)
(define pgm-pixels cdr)

(define (read-pgm file)

  (define (read-pgm-object port)
    (skip-white-space port)
    (let ((o (read port)))
      ;; to skip the newline or next whitespace
      (read-char port)
      (if (eof-object? o)
          (error "reached end of pgm file")
          o)))

  (define (skip-to-end-of-line port)
    (let loop ((ch (read-char port)))
      (if (not (eq? ch #\newline))
          (loop (read-char port)))))

  (define (white-space? ch)
    (case ch
      ((#\newline #\space #\tab) #t)
      (else #f)))

  (define (skip-white-space port)
    (let ((ch (peek-char port)))
      (cond ((white-space? ch)
             (read-char port)
             (skip-white-space port))
            ((eq? ch #\#)
             (skip-to-end-of-line port)
             (skip-white-space port))
            (else #f))))

  ;; The image file formats defined in netpbm
  ;; are problematical because they read the data
  ;; in the header as variable-length ISO-8859-1 text,
  ;; including arbitrary whitespace and comments,
  ;; and then they may read the rest of the file
  ;; as binary data.
  ;; So we give here a solution of how to deal
  ;; with these subtleties in Gambit Scheme.

  (call-with-input-file
      (list path:          file
            char-encoding: 'ISO-8859-1
            eol-encoding:  'lf)
    (lambda (port)

      ;; We're going to read text for a while,
      ;; then switch to binary.
      ;; So we need to turn off buffering until
      ;; we switch to binary.

      (port-settings-set! port '(buffering: #f))

      (let* ((header (read-pgm-object port))
             (columns (read-pgm-object port))
             (rows (read-pgm-object port))
             (greys (read-pgm-object port)))

        ;; Now we switch back to buffering
        ;; to speed things up.

        (port-settings-set! port '(buffering: #t))

        (make-pgm
         greys
         (array-copy
          (make-array
           (make-interval (vector rows columns))
           (cond ((or (eq? header 'p5)
                      (eq? header 'P5))
                  ;; pgm binary
                  (if (< greys 256)
                      ;; one byte/pixel
                      (lambda (i j)
                        (char->integer
                         (read-char port)))
                      ;; two bytes/pixel,
                      ;;little-endian
                      (lambda (i j)
                        (let* ((first-byte
                                (char->integer
                                 (read-char port)))
                               (second-byte
                                (char->integer
                                 (read-char port))))
                          (+ (* second-byte 256)
                             first-byte)))))
                 ;; pgm ascii
                 ((or (eq? header 'p2)
                      (eq? header 'P2))
                  (lambda (i j)
                      (read port)))
                   (else
                    (error "not a pgm file"))))
          (if (< greys 256)
              u8-storage-class
              u16-storage-class)))))))

(define (write-pgm pgm-data file #!optional force-ascii)
  (call-with-output-file
      (list path:          file
            char-encoding: 'ISO-8859-1
            eol-encoding:  'lf)
    (lambda (port)
      (let* ((greys
              (pgm-greys pgm-data))
             (pgm-array
              (pgm-pixels pgm-data))
             (domain
              (array-domain pgm-array))
             (rows
              (fx- (interval-upper-bound domain 0)
                   (interval-lower-bound domain 0)))
             (columns
              (fx- (interval-upper-bound domain 1)
                   (interval-lower-bound domain 1))))
        (if force-ascii
            (display "P2" port)
            (display "P5" port))
        (newline port)
        (display columns port) (display  port)
        (display rows port) (newline port)
        (display greys port) (newline port)
        (array-for-each
         (if force-ascii
             (let ((next-pixel-in-line 1))
               (lambda (p)
                 (write p port)
                 (if (fxzero? (fxand next-pixel-in-line 15))
                     (begin
                       (newline port)
                       (set! next-pixel-in-line 1))
                     (begin
                       (display  port)
                       (set! next-pixel-in-line
                             (fx+ 1 next-pixel-in-line))))))
             (if (fx< greys 256)
                 (lambda (p)
                   (write-u8 p port))
                 (lambda (p)
                   (write-u8 (fxand p 255) port)
                   (write-u8 (fxarithmetic-shift-right p 8)
                             port))))
         pgm-array)))))

One can write a a routine to convolve an image with a filter as follows:

(define (array-convolve source filter)
  (let* ((source-domain
          (array-domain source))
         (S_
          (array-getter source))
         (filter-domain
          (array-domain filter))
         (F_
          (array-getter filter))
         (result-domain
          (interval-dilate
           source-domain
           ;; the left bound of an interval is an equality,
           ;; the right bound is an inequality, hence the
           ;; the difference in the following two expressions
           (vector-map -
                       (interval-lower-bounds->vector filter-domain))
           (vector-map (lambda (x)
                         (- 1 x))
                       (interval-upper-bounds->vector filter-domain)))))
    (make-array result-domain
                (lambda (i j)
                  (array-fold
                   (lambda (p q)
                     (+ p q))
                   0
                   (make-array
                    filter-domain
                    (lambda (k l)
                      (* (S_ (+ i k)
                             (+ j l))
                         (F_ k l))))))
                )))

together with some filters

(define sharpen-filter
  (list->array
   '(0 -1  0
    -1  5 -1
     0 -1  0)
   (make-interval '#(-1 -1) '#(2 2))))

(define edge-filter
  (list->array
   '(0 -1  0
    -1  4 -1
     0 -1  0)
   (make-interval '#(-1 -1) '#(2 2))))

Our computations might results in pixel values outside the valid range, so we define

(define (round-and-clip pixel max-grey)
  (max 0 (min (exact (round pixel)) max-grey)))

We can then compute edges and sharpen a test image as follows:

(define test-pgm (read-pgm "girl.pgm"))

(let ((greys (pgm-greys test-pgm)))
  (write-pgm
   (make-pgm
    greys
    (array-map (lambda (p)
                 (round-and-clip p greys))
               (array-convolve
                (pgm-pixels test-pgm)
                sharpen-filter)))
   "sharper-test.pgm"))

(let* ((greys (pgm-greys test-pgm))
       (edge-array
        (array-copy
         (array-map
          abs
          (array-convolve
           (pgm-pixels test-pgm)
           edge-filter))))
       (max-pixel
        (array-fold max 0 edge-array))
       (normalizer
        (inexact (/ greys max-pixel))))
  (write-pgm
   (make-pgm
    greys
    (array-map (lambda (p)
                 (- greys
                    (round-and-clip (* p normalizer) greys)))
               edge-array))
   "edge-test.pgm"))

Viewing two-dimensional slices of three-dimensional data. One example might be viewing two-dimensional slices of three-dimensional data in different ways. If one has a $1024 \times 512\times 512$ 3D image of the body stored as a variable body, then one could get 1024 axial views, each $512\times512$, of this 3D body by (array-curry body 2); or 512 median views, each $1024\times512$, by (array-curry (array-permute body '#(1 0 2)) 2); or finally 512 frontal views, each again $1024\times512$ pixels, by (array-curry (array-permute body '#(2 0 1)) 2); see Anatomical plane. Note that the first permutation is not a rotation—you want to have the head up in both the median and frontal views.

Calculating second differences of images. For another example, if a real-valued function is defined on a two-dimensional interval $I$, its second difference in the direction $d$ at the point $x$ is defined as $\Delta^2_df(x)=f(x+2d)-2f(x+d)+f(x)$, and this function is defined only for those $x$ for which $x$, $x+d$, and $x+2d$ are all in $I$. See the beginning of the section on "Moduli of smoothness" in these notes on wavelets and approximation theory for more details.

Using this definition, the following code computes all second-order forward differences of an image in the directions $d,2 d,3 d,\ldots$, defined only on the domains where this makes sense:

(define (all-second-differences image direction)
  (let ((image-domain (array-domain image)))
    (let loop ((i 1)
               (result '()))
      (let ((negative-scaled-direction
             (vector-map (lambda (j)
                           (* -1 j i))
                         direction))
            (twice-negative-scaled-direction
             (vector-map (lambda (j)
                           (* -2 j i))
                         direction)))
        (cond ((interval-intersect
                image-domain
                (interval-translate
                 image-domain
                 negative-scaled-direction)
                (interval-translate
                 image-domain
                 twice-negative-scaled-direction))
               =>
               (lambda (subdomain)
                 (loop
                  (+ i 1)
                  (cons
                   (array-copy
                    (array-map
                     (lambda (f_i f_i+d f_i+2d)
                       (+ f_i+2d
                          (* -2. f_i+d)
                          f_i))
                     (array-extract
                      image
                      subdomain)
                     (array-extract
                      (array-translate
                       image
                       negative-scaled-direction)
                      subdomain)
                     (array-extract
                      (array-translate
                       image
                       twice-negative-scaled-direction)
                      subdomain)))
                   result))))
              (else
               (reverse result)))))))

We can define a small synthetic image of size 8x8 pixels and compute its second differences in various directions:

(define image
 (array-copy
  (make-array (make-interval '#(8 8))
              (lambda (i j)
                (exact->inexact (+ (* i i) (* j j)))))))

(define (expose difference-images)
  (pretty-print (map (lambda (difference-image)
                       (list (array-domain difference-image)
                             (array->list difference-image)))
                     difference-images)))

(begin
  (display
   "\nSecond-differences in the direction $k\times (1,0)$:\n")
  (expose (all-second-differences image '#(1 0)))
  (display
   "\nSecond-differences in the direction $k\times (1,1)$:\n")
  (expose (all-second-differences image '#(1 1)))
  (display
   "\nSecond-differences in the direction $k\times (1,-1)$:\n")
  (expose (all-second-differences image '#(1 -1))))

On Gambit 4.8.5, this yields (after some hand editing):

Second-differences in the direction $k\times (1,0)$:
((#<##interval #2 lower-bounds: #(0 0) upper-bounds: #(6 8)>
 (2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
  2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
  2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.))
 (#<##interval #3 lower-bounds: #(0 0) upper-bounds: #(4 8)>
  (8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8.
   8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8.))
 (#<##interval #4 lower-bounds: #(0 0) upper-bounds: #(2 8)>
  (18. 18. 18. 18. 18. 18. 18. 18. 18.
   18. 18. 18. 18. 18. 18. 18.)))

Second-differences in the direction $k\times (1,1)$:
((#<##interval #5 lower-bounds: #(0 0) upper-bounds: #(6 6)>
  (4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
   4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.))
 (#<##interval #6 lower-bounds: #(0 0) upper-bounds: #(4 4)>
  (16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16.
   16. 16.))
 (#<##interval #7 lower-bounds: #(0 0) upper-bounds: #(2 2)>
  (36. 36. 36. 36.)))

Second-differences in the direction $k\times (1,-1)$:
((#<##interval #8 lower-bounds: #(0 2) upper-bounds: #(6 8)>
  (4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
   4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.))
 (#<##interval #9 lower-bounds: #(0 4) upper-bounds: #(4 8)>
  (16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16.
   16. 16.))
 (#<##interval #10 lower-bounds: #(0 6) upper-bounds: #(2 8)>
  (36. 36. 36. 36.)))

You can see that with differences in the direction of only the first coordinate, the domains of the difference arrays get smaller in the first coordinate while staying the same in the second coordinate, and with differences in the diagonal directions, the domains of the difference arrays get smaller in both coordinates.

Separable operators. Many multi-dimensional transforms in signal processing are separable, in that the multi-dimensional transform can be computed by applying one-dimensional transforms in each of the coordinate directions. Examples of such transforms include the Fast Fourier Transform and the Fast Hyperbolic Wavelet Transform. Each one-dimensional subdomain of the complete domain is called a pencil, and the same one-dimensional transform is applied to all pencils in a given direction. Given the one-dimensional array transform, one can define the multidimensional transform as follows:

(define (make-separable-transform 1D-transform)
  (lambda (a)
    (let ((n (array-dimension a)))
      (do ((d 0 (fx+ d 1)))
          ((fx= d n))
        (array-for-each
         1D-transform
         (array-curry (array-rotate a d) 1))))))

Here we have cycled through all rotations, putting each axis in turn at the end, and then applied 1D-transform to each of the pencils along that axis.

Wavelet transforms in particular are calculated by recursively applying a transform to an array and then downsampling the array; the inverse transform recursively downsamples and then applies a transform. So we define the following primitives:

(define (recursively-apply-transform-and-downsample transform)
  (lambda (a)
    (let ((sample-vector (make-vector (array-dimension a) 2)))
      (define (helper a)
        (if (fx< 1 (interval-upper-bound (array-domain a) 0))
            (begin
              (transform a)
              (helper (array-sample a sample-vector)))))
      (helper a))))

(define (recursively-downsample-and-apply-transform transform)
  (lambda (a)
    (let ((sample-vector (make-vector (array-dimension a) 2)))
      (define (helper a)
        (if (fx< 1 (interval-upper-bound (array-domain a) 0))
            (begin
              (helper (array-sample a sample-vector))
              (transform a))))
      (helper a))))

By adding a single loop that calculates scaled sums and differences of adjacent elements in a one-dimensional array, we can define various Haar wavelet transforms as follows:

(define (1D-Haar-loop a)
  (let ((a_ (array-getter a))
        (a! (array-setter a))
        (n (interval-upper-bound (array-domain a) 0)))
    (do ((i 0 (fx+ i 2)))
        ((fx= i n))
      (let* ((a_i               (a_ i))
             (a_i+1             (a_ (fx+ i 1)))
             (scaled-sum        (fl/ (fl+ a_i a_i+1) (flsqrt 2.0)))
             (scaled-difference (fl/ (fl- a_i a_i+1) (flsqrt 2.0))))
        (a! scaled-sum i)
        (a! scaled-difference (fx+ i 1))))))

(define 1D-Haar-transform
  (recursively-apply-transform-and-downsample 1D-Haar-loop))

(define 1D-Haar-inverse-transform
  (recursively-downsample-and-apply-transform 1D-Haar-loop))

(define hyperbolic-Haar-transform
  (make-separable-transform 1D-Haar-transform))

(define hyperbolic-Haar-inverse-transform
  (make-separable-transform 1D-Haar-inverse-transform))

(define Haar-transform
  (recursively-apply-transform-and-downsample
   (make-separable-transform 1D-Haar-loop)))

(define Haar-inverse-transform
  (recursively-downsample-and-apply-transform
   (make-separable-transform 1D-Haar-loop)))

We then define an image that is a multiple of a single, two-dimensional hyperbolic Haar wavelet, compute its hyperbolic Haar transform (which should have only one nonzero coefficient), and then the inverse transform:

(let ((image
       (array-copy
        (make-array (make-interval '#(4 4))
                    (lambda (i j)
                      (case i
                        ((0) 1.)
                        ((1) -1.)
                        (else 0.)))))))
  (display "
Initial image:
")
  (pretty-print (list (array-domain image)
                      (array->list image)))
  (hyperbolic-Haar-transform image)
  (display "\nArray of hyperbolic Haar wavelet coefficients: \n")
  (pretty-print (list (array-domain image)
                      (array->list image)))
  (hyperbolic-Haar-inverse-transform image)
  (display "\nReconstructed image: \n")
  (pretty-print (list (array-domain image)
                      (array->list image))))

This yields:

Initial image:
(#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)>
 (1. 1. 1. 1. -1. -1. -1. -1. 0. 0. 0. 0. 0. 0. 0. 0.))

Array of hyperbolic Haar wavelet coefficients:
(#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)>
 (0. 0. 0. 0. 2.8284271247461894 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.))

Reconstructed image:
(#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)>
 (.9999999999999996
  .9999999999999996
  .9999999999999996
  .9999999999999996
  -.9999999999999996
  -.9999999999999996
  -.9999999999999996
  -.9999999999999996
  0.
  0.
  0.
  0.
  0.
  0.
  0.
  0.))

In perfect arithmetic, this hyperbolic Haar transform is orthonormal, in that the sum of the squares of the elements of the image is the same as the sum of the squares of the hyperbolic Haar coefficients of the image. We can see that this is approximately true here.

We can apply the (nonhyperbolic) Haar transform to the same image as follows:

 (let ((image
       (array-copy
        (make-array (make-interval '#(4 4))
                    (lambda (i j)
                      (case i
                        ((0) 1.)
                        ((1) -1.)
                        (else 0.)))))))
  (display "\nInitial image:\n")
  (pretty-print (list (array-domain image)
                      (array->list image)))
  (Haar-transform image)
  (display "\nArray of Haar wavelet coefficients: \n")
  (pretty-print (list (array-domain image)
                      (array->list image)))
  (Haar-inverse-transform image)
  (display "\nReconstructed image: \n")
  (pretty-print (list (array-domain image)
                      (array->list image))))

This yields:

Initial image:
(#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)>
 (1. 1. 1. 1. -1. -1. -1. -1. 0. 0. 0. 0. 0. 0. 0. 0.))

Array of Haar wavelet coefficients:
(#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)>
 (0. 0. 0. 0. 1.9999999999999998 0. 1.9999999999999998 0. 0. 0. 0. 0. 0. 0. 0. 0.))

Reconstructed image:
(#<##interval #12 lower-bounds: #(0 0) upper-bounds: #(4 4)>
 (.9999999999999997
  .9999999999999997
  .9999999999999997
  .9999999999999997
  -.9999999999999997
  -.9999999999999997
  -.9999999999999997
  -.9999999999999997
  0.
  0.
  0.
  0.
  0.
  0.
  0.
  0.))

You see in this example that this particular image has two, not one, nonzero coefficients in the two-dimensional Haar transform, which is again approximately orthonormal.

Matrix multiplication and Gaussian elimination. While we have avoided conflating matrices and arrays, we give here some examples of matrix operations defined using operations from this SRFI.

Given a nonsingular square matrix $A$ we can overwrite $A$ with lower-triangular matrix $L$ with ones on the diagonal and upper-triangular matrix $U$ so that $A=LU$ as follows. (We assume "pivoting" isn't needed.) For example, if $$A=\begin{pmatrix} a_{11}&a_{12}&a_{13}\\ a_{21}&a_{22}&a_{23}\\ a_{31}&a_{32}&a_{33}\end{pmatrix}=\begin{pmatrix} 1&0&0\\ \ell_{21}&1&0\\ \ell_{31}&\ell_{32}&1\end{pmatrix}\begin{pmatrix} u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\ 0&0&u_{33}\end{pmatrix}$$ then $A$ is overwritten with $$ \begin{pmatrix} u_{11}&u_{12}&u_{13}\\ \ell_{21}&u_{22}&u_{23}\\ \ell_{31}&\ell_{32}&u_{33}\end{pmatrix}. $$ The code uses array-assign!, specialized-array-share, array-extract, and array-outer-product.

(define (LU-decomposition A)
  ;; Assumes the domain of A is [0,n)\times [0,n)
  ;; and that Gaussian elimination can be applied
  ;; without pivoting.
  (let ((n
         (interval-upper-bound (array-domain A) 0))
        (A_
         (array-getter A)))
    (do ((i 0 (fx+ i 1)))
        ((= i (fx- n 1)) A)
      (let* ((pivot
              (A_ i i))
             (column/row-domain
              ;; both will be one-dimensional
              (make-interval (vector (+ i 1))
                             (vector n)))
             (column
              ;; the column below the (i,i) entry
              (specialized-array-share A
                                       column/row-domain
                                       (lambda (k)
                                         (values k i))))
             (row
              ;; the row to the right of the (i,i) entry
              (specialized-array-share A
                                       column/row-domain
                                       (lambda (k)
                                         (values i k))))

             ;; the subarray to the right and
             ;; below the (i,i) entry
             (subarray
              (array-extract
               A (make-interval
                  (vector (fx+ i 1) (fx+ i 1))
                  (vector n         n)))))
        ;; Compute multipliers.
        (array-assign!
         column
         (array-map (lambda (x)
                      (/ x pivot))
                    column))
        ;; Subtract the outer product of i'th
        ;; row and column from the subarray.
        (array-assign!
         subarray
         (array-map -
                    subarray
                    (array-outer-product * column row)))))))

We then define a $4\times 4$ Hilbert matrix:

(define A
  (array-copy
   (make-array (make-interval '#(4 4))
               (lambda (i j)
                 (/ (+ 1 i j))))))

(define (array-display A)

  (define (display-item x)
    (display x) (display "\t"))

  (newline)
  (case (array-dimension A)
    ((1) (array-for-each display-item A) (newline))
    ((2) (array-for-each (lambda (row)
                           (array-for-each display-item row)
                           (newline))
                         (array-curry A 1)))
    (else
     (error "array-display can't handle > 2 dimensions: " A))))

(display "\nHilbert matrix:\n\n")

(array-display A)

;;; which displays:
;;; 1       1/2     1/3     1/4
;;; 1/2     1/3     1/4     1/5
;;; 1/3     1/4     1/5     1/6
;;; 1/4     1/5     1/6     1/7

(LU-decomposition A)

(display "\nLU decomposition of Hilbert matrix:\n\n")

(array-display A)

;;; which displays:
;;; 1       1/2     1/3     1/4
;;; 1/2     1/12    1/12    3/40
;;; 1/3     1       1/180   1/120
;;; 1/4     9/10    3/2     1/2800

We can now define matrix multiplication as follows to check our result:

;;; Functions to extract the lower- and upper-triangular
;;; matrices of the LU decomposition of A.

(define (L a)
  (let ((a_ (array-getter a))
        (d  (array-domain a)))
    (make-array
     d
     (lambda (i j)
       (cond ((= i j) 1)        ;; diagonal
             ((> i j) (a_ i j)) ;; below diagonal
             (else 0))))))      ;; above diagonal

(define (U a)
  (let ((a_ (array-getter a))
        (d  (array-domain a)))
    (make-array
     d
     (lambda (i j)
       (cond ((<= i j) (a_ i j)) ;; diagonal and above
             (else 0))))))       ;; below diagonal

(display "\nLower triangular matrix of decomposition of Hilbert matrix:\n\n")
(array-display (L A))

;;; which displays:
;;; 1       0       0       0
;;; 1/2     1       0       0
;;; 1/3     1       1       0
;;; 1/4     9/10    3/2     1


(display "\nUpper triangular matrix of decomposition of Hilbert matrix:\n\n")
(array-display (U A))

;;; which displays:
;;; 1       1/2     1/3     1/4
;;; 0       1/12    1/12    3/40
;;; 0       0       1/180   1/120
;;; 0       0       0       1/2800

;;; We'll define a brief, not-very-efficient matrix multiply routine.

(define (array-dot-product a b)
  (array-fold + 0 (array-map * a b)))

(define (matrix-multiply a b)
  (let ((a-rows
         (array-copy (array-curry a 1)))
        (b-columns
         (array-copy (array-curry (array-rotate b 1) 1))))
    (array-outer-product array-dot-product a-rows b-columns)))

;;; We'll check that the product of the result of LU
;;; decomposition of A is again A.

(define product (matrix-multiply (L A) (U A)))

(display "\nProduct of lower and upper triangular matrices \n")
(display "of LU decomposition of Hilbert matrix:\n\n")
(array-display product)

;;; which displays:
;;; 1       1/2     1/3     1/4
;;; 1/2     1/3     1/4     1/5
;;; 1/3     1/4     1/5     1/6
;;; 1/4     1/5     1/6     1/7

Inner products. One can define an APL-style inner product as

(define (inner-product A f g B)
  (array-outer-product
   (lambda (a b)
     (array-reduce f (array-map g a b)))
   (array-copy (array-curry A 1))
   (array-copy (array-curry (array-rotate B 1) 1))))

This routine differs from that found in APL in several ways: The arguments A and B must each have two or more dimensions, and the result is always an array, never a scalar.

We take some examples from the APLX Language Reference:

;; Examples from
;; http://microapl.com/apl_help/ch_020_020_880.htm

(define TABLE1
  (list->array
   '(1 2
     5 4
     3 0)
   (make-interval '#(3 2))))

(define TABLE2
  (list->array
   '(6 2 3 4
     7 0 1 8)
   (make-interval '#(2 4))))

(array-display (inner-product TABLE1 + * TABLE2))

;;; Displays
;;; 20      2       5       20
;;; 58      10      19      52
;;; 18      6       9       12

(define X   ;; a "row vector"
  (list->array '(1 3 5 7) (make-interval '#(1 4))))

(define Y   ;; a "column vector"
  (list->array '(2 3 6 7) (make-interval '#(4 1))))

(array-display (inner-product X + (lambda (x y) (if (= x y) 1 0)) Y))

;;; Displays
;;; 2

Acknowledgments

The SRFI author thanks Edinah K Gnang, John Cowan, Sudarshan S Chawathe, Jamison Hope, and Per Bothner for their comments and suggestions, and Arthur A. Gleckler, SRFI Editor, for his guidance and patience.

References

"multi-dimensional arrays in R5RS?", by Alan Bawden.
SRFI 4: Homogeneous Numeric Vector Datatypes, by Marc Feeley.
SRFI 25: Multi-dimensional Array Primitives, by Jussi Piitulainen.
SRFI 47: Array, by Aubrey Jaffer.
SRFI 58: Array Notation, by Aubrey Jaffer.
SRFI 63: Homogeneous and Heterogeneous Arrays, by Aubrey Jaffer.
SRFI 164: Enhanced multi-dimensional Arrays, by Per Bothner.

© 2016, 2018, 2020 Bradley J Lucier. All Rights Reserved.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice (including the next
paragraph) shall be included in all copies or substantial portions of the
Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

chickadee » srfi-179

Identifier search

TOC »

Contents »