Find file
Fetching contributors…
Cannot retrieve contributors at this time
642 lines (567 sloc) 22.7 KB
;;; array
;;; 1997 - 2001 Jussi Piitulainen
;;; --- Intro ---
;;; This interface to arrays is based on Alan Bawden's array.scm of
;;; 1993 (earlier version in the Internet Repository and another
;;; version in SLIB). This is a complete rewrite, to be consistent
;;; with the rest of Scheme and to make arrays independent of lists.
;;; Some modifications are due to discussion in srfi-25 mailing list.
;;; (array? obj)
;;; (make-array shape [obj]) changed arguments
;;; (shape bound ...) new
;;; (array shape obj ...) new
;;; (array-rank array) changed name back
;;; (array-start array dimension) new
;;; (array-end array dimension) new
;;; (array-ref array k ...)
;;; (array-ref array index) new variant
;;; (array-set! array k ... obj) changed argument order
;;; (array-set! array index obj) new variant
;;; (share-array array shape proc) changed arguments
;;; All other variables in this file have names in "array:".
;;; Should there be a way to make arrays with initial values mapped
;;; from indices? Sure. The current "initial object" is lame.
;;; Removed (array-shape array) from here. There is a new version
;;; in arlib though.
;;; --- Representation type dependencies ---
;;; The mapping from array indices to the index to the underlying vector
;;; is whatever array:optimize returns. The file "opt" provides three
;;; representations:
;;; mbda) mapping is a procedure that allows an optional argument
;;; tter) mapping is two procedures that takes exactly the indices
;;; ctor) mapping is a vector of a constant term and coefficients
;;; Choose one in "opt" to make the optimizer. Then choose the matching
;;; implementation of array-ref and array-set!.
;;; These should be made macros to inline them. Or have a good compiler
;;; and plant the package as a module.
;;; 1. Pick an optimizer.
;;; 2. Pick matching index representation.
;;; 3. Pick a record implementation; as-procedure is generic; syntax inlines.
;;; 3. This file is otherwise portable.
;;; --- Portable R5RS (R4RS and multiple values) ---
;;; (array? obj)
;;; returns #t if `obj' is an array and #t or #f otherwise.
(define (array? obj)
(array:array? obj))
;;; (make-array shape)
;;; (make-array shape obj)
;;; makes array of `shape' with each cell containing `obj' initially.
(define (make-array shape . rest)
(or (array:good-shape? shape)
(error "make-array: shape is not a shape"))
(apply array:make-array shape rest))
(define (array:make-array shape . rest)
(let ((size (array:size shape)))
(if (pair? rest)
(apply (lambda (o) (make-vector size o)) rest)
(make-vector size))
(if (= size 0)
(vector-ref (array:shape shape) 1))
(array:make-index shape)
(vector-ref (array:shape shape) 1)))
(array:shape->vector shape))))
;;; (shape bound ...)
;;; makes a shape. Bounds must be an even number of exact, pairwise
;;; non-decreasing integers. Note that any such array can be a shape.
(define (shape . bounds)
(let ((v (list->vector bounds)))
(or (even? (vector-length v))
(error (string-append "shape: uneven number of bounds: "
(array:list->string bounds))))
(let ((shp (array:make
(if (pair? bounds)
(vector 0 (quotient (vector-length v) 2)
0 2))))
(or (array:good-shape? shp)
(error (string-append "shape: bounds are not pairwise "
"non-decreasing exact integers: "
(array:list->string bounds))))
;;; (array shape obj ...)
;;; is analogous to `vector'.
(define (array shape . elts)
(or (array:good-shape? shape)
(error (string-append "array: shape " (array:thing->string shape)
" is not a shape")))
(let ((size (array:size shape)))
(let ((vector (list->vector elts)))
(or (= (vector-length vector) size)
(error (string-append "array: an array of shape "
(array:vector shape))
" has "
(number->string size)
" elements but got "
(number->string (vector-length vector))
" values: "
(array:list->string elts))))
(if (= size 0)
(vector-ref (array:shape shape) 1))
(array:make-index shape)
(vector-ref (array:shape shape) 1)))
(array:shape->vector shape)))))
;;; (array-rank array)
;;; returns the number of dimensions of `array'.
(define (array-rank array)
(quotient (vector-length (array:shape array)) 2))
;;; (array-start array k)
;;; returns the lower bound index of array along dimension k. This is
;;; the least valid index along that dimension if the dimension is not
;;; empty.
(define (array-start array d)
(vector-ref (array:shape array) (+ d d)))
;;; (array-end array k)
;;; returns the upper bound index of array along dimension k. This is
;;; not a valid index. If the dimension is empty, this is the same as
;;; the lower bound along it.
(define (array-end array d)
(vector-ref (array:shape array) (+ d d 1)))
;;; (share-array array shape proc)
;;; makes an array that shares elements of `array' at shape `shape'.
;;; The arguments to `proc' are indices of the result. The values of
;;; `proc' are indices of `array'.
;;; Todo: in the error message, should recognise the mapping and show it.
(define (share-array array subshape f)
(or (array:good-shape? subshape)
(error (string-append "share-array: shape "
(array:thing->string subshape)
" is not a shape")))
(let ((subsize (array:size subshape)))
(or (array:good-share? subshape subsize f (array:shape array))
(error (string-append "share-array: subshape "
(array:vector subshape))
" does not map into supershape "
(array:shape array))
" under mapping "
(vector-ref (array:shape subshape) 1)))))
(let ((g (array:index array)))
(array:vector array)
(if (= subsize 0)
(vector-ref (array:shape subshape) 1))
(lambda ks
(lambda () (apply f ks))
(lambda ks (array:vector-index g ks))))
(vector-ref (array:shape subshape) 1)))
(array:shape->vector subshape)))))
;;; --- Hrmph ---
;;; (array:share/index! ...)
;;; reuses a user supplied index object when recognising the
;;; mapping. The mind balks at the very nasty side effect that
;;; exposes the implementation. So this is not in the spec.
;;; But letting index objects in at all creates a pressure
;;; to go the whole hog. Arf.
;;; Use array:optimize-empty for an empty array to get a
;;; clearly invalid vector index.
;;; Surely it's perverse to use an actor for index here? But
;;; the possibility is provided for completeness.
(define (array:share/index! array subshape proc index)
(array:vector array)
(if (= (array:size subshape) 0)
(quotient (vector-length (array:shape array)) 2))
((if (vector? index)
(lambda (subindex)
(let ((superindex (proc subindex)))
(if (vector? superindex)
(quotient (vector-length (array:shape array)) 2)
(array:index array)
(quotient (vector-length (array:shape array)) 2)
(array:index array)
(array:vector superindex)
(array:index superindex)))))
(array:shape->vector subshape)))
(define (array:optimize/vector f v)
(let ((r (vector-length v)))
(do ((k 0 (+ k 1)))
((= k r))
(vector-set! v k 0))
(let ((n0 (f v))
(cs (make-vector (+ r 1)))
(apply (array:applier-to-vector (+ r 1))))
(vector-set! cs 0 n0)
(let wok ((k 0))
(if (< k r)
(let ((k1 (+ k 1)))
(vector-set! v k 1)
(let ((nk (- (f v) n0)))
(vector-set! v k 0)
(vector-set! cs k1 nk)
(wok k1)))))
(apply (array:maker r) cs))))
(define (array:optimize/actor f a)
(let ((r (array-end a 0))
(v (array:vector a))
(i (array:index a)))
(do ((k 0 (+ k 1)))
((= k r))
(vector-set! v (array:actor-index i k) 0))
(let ((n0 (f a))
(cs (make-vector (+ r 1)))
(apply (array:applier-to-vector (+ r 1))))
(vector-set! cs 0 n0)
(let wok ((k 0))
(if (< k r)
(let ((k1 (+ k 1))
(t (array:actor-index i k)))
(vector-set! v t 1)
(let ((nk (- (f a) n0)))
(vector-set! v t 0)
(vector-set! cs k1 nk)
(wok k1)))))
(apply (array:maker r) cs))))
;;; --- Internals ---
(define (array:shape->vector shape)
(let ((idx (array:index shape))
(shv (array:vector shape))
(rnk (vector-ref (array:shape shape) 1)))
(let ((vec (make-vector (* rnk 2))))
(do ((k 0 (+ k 1)))
((= k rnk)
(vector-set! vec (+ k k)
(vector-ref shv (array:shape-vector-index idx k 0)))
(vector-set! vec (+ k k 1)
(vector-ref shv (array:shape-vector-index idx k 1)))))))
;;; (array:size shape)
;;; returns the number of elements in arrays of shape `shape'.
(define (array:size shape)
(let ((idx (array:index shape))
(shv (array:vector shape))
(rnk (vector-ref (array:shape shape) 1)))
(do ((k 0 (+ k 1))
(s 1 (* s
(- (vector-ref shv (array:shape-vector-index idx k 1))
(vector-ref shv (array:shape-vector-index idx k 0))))))
((= k rnk) s))))
;;; (array:make-index shape)
;;; returns an index function for arrays of shape `shape'. This is a
;;; runtime composition of several variable arity procedures, to be
;;; passed to array:optimize for recognition as an affine function of
;;; as many variables as there are dimensions in arrays of this shape.
(define (array:make-index shape)
(let ((idx (array:index shape))
(shv (array:vector shape))
(rnk (vector-ref (array:shape shape) 1)))
(do ((f (lambda () 0)
(lambda (k . ks)
(+ (* s (- k (vector-ref
(array:shape-vector-index idx (- j 1) 0))))
(apply f ks))))
(s 1 (* s (- (vector-ref
(array:shape-vector-index idx (- j 1) 1))
(array:shape-vector-index idx (- j 1) 0)))))
(j rnk (- j 1)))
((= j 0)
;;; --- Error checking ---
;;; (array:good-shape? shape)
;;; returns true if `shape' is an array of the right shape and its
;;; elements are exact integers that pairwise bound intervals `[lo..hi)´.
(define (array:good-shape? shape)
(and (array:array? shape)
(let ((u (array:shape shape))
(v (array:vector shape))
(x (array:index shape)))
(and (= (vector-length u) 4)
(= (vector-ref u 0) 0)
(= (vector-ref u 2) 0)
(= (vector-ref u 3) 2))
(let ((p (vector-ref u 1)))
(do ((k 0 (+ k 1))
(true #t (let ((lo (vector-ref
(array:shape-vector-index x k 0)))
(hi (vector-ref
(array:shape-vector-index x k 1))))
(and true
(integer? lo)
(exact? lo)
(integer? hi)
(exact? hi)
(<= lo hi)))))
((= k p) true))))))
;;; (array:good-share? subv subsize mapping superv)
;;; returns true if the extreme indices in the subshape vector map
;;; into the bounds in the supershape vector.
;;; If some interval in `subv' is empty, then `subv' is empty and its
;;; image under `f' is empty and it is trivially alright. One must
;;; not call `f', though.
(define (array:good-share? subshape subsize f super)
(or (zero? subsize)
((sub (array:vector subshape))
(dex (array:index subshape))
(ck (lambda (k ks)
(if (zero? k)
(lambda () (apply f ks))
(lambda qs (array:good-indices? qs super)))
(and (ck (- k 1)
(cons (vector-ref
(- k 1)
(ck (- k 1)
(cons (- (vector-ref
(- k 1)
(let ((rnk (vector-ref (array:shape subshape) 1)))
(or (array:unchecked-share-depth? rnk)
(ck rnk '()))))))
;;; Check good-share on 10 dimensions at most. The trouble is,
;;; the cost of this check is exponential in the number of dimensions.
(define (array:unchecked-share-depth? rank)
(if (> rank 10)
(display `(warning: unchecked depth in share:
,rank subdimensions))
;;; (array:check-indices caller indices shape-vector)
;;; (array:check-indices.o caller indices shape-vector)
;;; (array:check-index-vector caller index-vector shape-vector)
;;; return if the index is in bounds, else signal error.
;;; Shape-vector is the internal representation, with
;;; b and e for dimension k at 2k and 2k + 1.
(define (array:check-indices who ks shv)
(or (array:good-indices? ks shv)
(error (array:not-in who ks shv))))
(define (array:check-indices.o who ks shv)
(or (array:good-indices.o? ks shv)
(error (array:not-in who (reverse (cdr (reverse ks))) shv))))
(define (array:check-index-vector who ks shv)
(or (array:good-index-vector? ks shv)
(error (array:not-in who (vector->list ks) shv))))
(define (array:check-index-actor who ks shv)
(let ((shape (array:shape ks)))
(or (and (= (vector-length shape) 2)
(= (vector-ref shape 0) 0))
(error "not an actor"))
(or (array:good-index-actor?
(vector-ref shape 1)
(array:vector ks)
(array:index ks)
(array:not-in who (do ((k (vector-ref shape 1) (- k 1))
(m '() (cons (vector-ref
(array:vector ks)
(array:index ks)
(- k 1)))
((= k 0) m))
(define (array:good-indices? ks shv)
(let ((d2 (vector-length shv)))
(do ((kp ks (if (pair? kp)
(cdr kp)))
(k 0 (+ k 2))
(true #t (and true (pair? kp)
(array:good-index? (car kp) shv k))))
((= k d2)
(and true (null? kp))))))
(define (array:good-indices.o? ks.o shv)
(let ((d2 (vector-length shv)))
(do ((kp ks.o (if (pair? kp)
(cdr kp)))
(k 0 (+ k 2))
(true #t (and true (pair? kp)
(array:good-index? (car kp) shv k))))
((= k d2)
(and true (pair? kp) (null? (cdr kp)))))))
(define (array:good-index-vector? ks shv)
(let ((r2 (vector-length shv)))
(and (= (* 2 (vector-length ks)) r2)
(do ((j 0 (+ j 1))
(k 0 (+ k 2))
(true #t (and true
(array:good-index? (vector-ref ks j) shv k))))
((= k r2) true)))))
(define (array:good-index-actor? r v i shv)
(and (= (* 2 r) (vector-length shv))
(do ((j 0 (+ j 1))
(k 0 (+ k 2))
(true #t (and true
(array:good-index? (vector-ref
(array:actor-index i j))
((= j r) true))))
;;; (array:good-index? index shape-vector 2d)
;;; returns true if index is within bounds for dimension 2d/2.
(define (array:good-index? w shv k)
(and (integer? w)
(exact? w)
(<= (vector-ref shv k) w)
(< w (vector-ref shv (+ k 1)))))
(define (array:not-in who ks shv)
(let ((index (array:list->string ks))
(bounds (array:shape-vector->string shv)))
(error (string-append who
": index " index
" not in bounds " bounds))))
(define (array:list->string ks)
(do ((index "" (string-append index (array:thing->string (car ks)) " "))
(ks ks (cdr ks)))
((null? ks) index)))
(define (array:shape-vector->string shv)
(do ((bounds "" (string-append bounds
(number->string (vector-ref shv t))
(number->string (vector-ref shv (+ t 1)))
" "))
(t 0 (+ t 2)))
((= t (vector-length shv)) bounds)))
(define (array:thing->string thing)
((number? thing) (number->string thing))
((symbol? thing) (string-append "#<symbol>" (symbol->string thing)))
((char? thing) "#<char>")
((string? thing) "#<string>")
((list? thing) (string-append "#" (number->string (length thing))
((pair? thing) "#<pair>")
((array? thing) "#<array>")
((vector? thing) (string-append "#" (number->string
(vector-length thing))
((procedure? thing) "#<procedure>")
(case thing
((()) "()")
((#t) "#t")
((#f) "#f")
;;; And to grok an affine map, vector->vector type. Column k of arr
;;; will contain coefficients n0 ... nm of 1 k1 ... km for kth value.
;;; These are for the error message when share fails.
(define (array:index-ref ind k)
(if (vector? ind)
(vector-ref ind k)
(array:vector ind)
(array:actor-index (array:index ind) k))))
(define (array:index-set! ind k o)
(if (vector? ind)
(vector-set! ind k o)
(array:vector ind)
(array:actor-index (array:index ind) k)
(define (array:index-length ind)
(if (vector? ind)
(vector-length ind)
(vector-ref (array:shape ind) 1)))
(define (array:map->string proc r)
(let* ((m (array:grok/arguments proc r))
(s (vector-ref (array:shape m) 3)))
(do ((i "" (string-append i c "k" (number->string k)))
(c "" ", ")
(k 1 (+ k 1)))
((< r k)
(do ((o "" (string-append o c (array:map-column->string m r k)))
(c "" ", ")
(k 0 (+ k 1)))
((= k s)
(string-append i " => " o)))))))
(define (array:map-column->string m r k)
(let ((v (array:vector m))
(i (array:index m)))
(let ((n0 (vector-ref v (array:vector-index i (list 0 k)))))
(let wok ((j 1)
(e (if (= n0 0) "" (number->string n0))))
(if (<= j r)
(let ((nj (vector-ref v (array:vector-index i (list j k)))))
(if (= nj 0)
(wok (+ j 1) e)
(let* ((nj (if (= nj 1) ""
(if (= nj -1) "-"
(string-append (number->string nj)
" "))))
(njkj (string-append nj "k" (number->string j))))
(if (string=? e "")
(wok (+ j 1) njkj)
(wok (+ j 1) (string-append e " + " njkj))))))
(if (string=? e "") "0" e))))))
(define (array:grok/arguments proc r)
(lambda (vec)
(lambda ()
(array:apply-to-vector r proc vec))
(make-vector r)))
(define (array:grok/index! proc in)
(let ((m (array:index-length in)))
(do ((k 0 (+ k 1)))
((= k m))
(array:index-set! in k 0))
(let* ((n0 (proc in))
(n (array:index-length n0)))
(let ((arr (make-array (shape 0 (+ m 1) 0 n)))) ; (*)
(do ((k 0 (+ k 1)))
((= k n))
(array-set! arr 0 k (array:index-ref n0 k))) ; (**)
(do ((j 0 (+ j 1)))
((= j m))
(array:index-set! in j 1)
(let ((nj (proc in)))
(array:index-set! in j 0)
(do ((k 0 (+ k 1)))
((= k n))
(array-set! arr (+ j 1) k (- (array:index-ref nj k) ; (**)
(array:index-ref n0 k))))))
;; (*) Should not use `make-array' and `shape' here
;; (**) Should not use `array-set!' here
;; Should use something internal to the library instead: either lower
;; level code (preferable but complex) or alternative names to these same.