[bpt/guile.git] / module / slib / array.scm

;;;;"array.scm" Arrays for Scheme
; Copyright (C) 1993 Alan Bawden
;
; Permission to copy this software, to redistribute it, and to use it
; for any purpose is granted, subject to the following restrictions and
; understandings.
;
; 1.  Any copy made of this software must include this copyright notice
; in full.
;
; 2.  Users of this software agree to make their best efforts (a) to
; return to me any improvements or extensions that they make, so that
; these may be included in future releases; and (b) to inform me of
; noteworthy uses of this software.
;
; 3.  I have made no warrantee or representation that the operation of
; this software will be error-free, and I am under no obligation to
; provide any services, by way of maintenance, update, or otherwise.
;
; 4.  In conjunction with products arising from the use of this material,
; there shall be no use of my name in any advertising, promotional, or
; sales literature without prior written consent in each case.
;
; Alan Bawden
; MIT Room NE43-510
; 545 Tech. Sq.
; Cambridge, MA 02139
; Alan@LCS.MIT.EDU

(require 'record)

;(declare (usual-integrations))

(define array:rtd
  (make-record-type "Array"
    '(indexer		; Must be a -linear- function!
      shape		; Inclusive bounds: ((lower upper) ...)
      vector		; The actual contents
      )))

(define array:indexer (record-accessor array:rtd 'indexer))
(define array-shape (record-accessor array:rtd 'shape))
(define array:vector (record-accessor array:rtd 'vector))

(define array? (record-predicate array:rtd))

(define (array-rank obj)
  (if (array? obj) (length (array-shape obj)) 0))

(define (array-dimensions ra)
  (map (lambda (ind) (if (zero? (car ind)) (+ 1 (cadr ind)) ind))
       (array-shape ra)))

(define array:construct
  (record-constructor array:rtd '(shape vector indexer)))

(define (array:compute-shape specs)
  (map (lambda (spec)
	 (cond ((and (integer? spec)
		     (< 0 spec))
		(list 0 (- spec 1)))
	       ((and (pair? spec)
		     (pair? (cdr spec))
		     (null? (cddr spec))
		     (integer? (car spec))
		     (integer? (cadr spec))
		     (<= (car spec) (cadr spec)))
		spec)
	       (else (slib:error "array: Bad array dimension: " spec))))
       specs))

(define (make-array initial-value . specs)
  (let ((shape (array:compute-shape specs)))
    (let loop ((size 1)
	       (indexer (lambda () 0))
	       (l (reverse shape)))
      (if (null? l)
	  (array:construct shape
			   (make-vector size initial-value)
			   (array:optimize-linear-function indexer shape))
	  (loop (* size (+ 1 (- (cadar l) (caar l))))
		(lambda (first-index . rest-of-indices)
		  (+ (* size (- first-index (caar l)))
		     (apply indexer rest-of-indices)))
		(cdr l))))))

(define (make-shared-array array mapping . specs)
  (let ((new-shape (array:compute-shape specs))
	(old-indexer (array:indexer array)))
    (let check ((indices '())
		(bounds (reverse new-shape)))
      (cond ((null? bounds)
	     (array:check-bounds array (apply mapping indices)))
	    (else
	     (check (cons (caar bounds) indices) (cdr bounds))
	     (check (cons (cadar bounds) indices) (cdr bounds)))))
    (array:construct new-shape
		     (array:vector array)
		     (array:optimize-linear-function
		       (lambda indices
			 (apply old-indexer (apply mapping indices)))
		       new-shape))))

(define (array:in-bounds? array indices)
  (let loop ((indices indices)
	     (shape (array-shape array)))
    (if (null? indices)
	(null? shape)
	(let ((index (car indices)))
	  (and (not (null? shape))
	       (integer? index)
	       (<= (caar shape) index (cadar shape))
	       (loop (cdr indices) (cdr shape)))))))

(define (array:check-bounds array indices)
  (or (array:in-bounds? array indices)
      (slib:error "array: Bad indices for " array indices)))

(define (array-ref array . indices)
  (array:check-bounds array indices)
  (vector-ref (array:vector array)
	      (apply (array:indexer array) indices)))

(define (array-set! array new-value . indices)
  (array:check-bounds array indices)
  (vector-set! (array:vector array)
	       (apply (array:indexer array) indices)
	       new-value))

(define (array-in-bounds? array . indices)
  (array:in-bounds? array indices))

; Fast versions of ARRAY-REF and ARRAY-SET! that do no error checking,
; and don't cons intermediate lists of indices:

(define (array-1d-ref a i0)
  (vector-ref (array:vector a) ((array:indexer a) i0)))

(define (array-2d-ref a i0 i1)
  (vector-ref (array:vector a) ((array:indexer a) i0 i1)))

(define (array-3d-ref a i0 i1 i2)
  (vector-ref (array:vector a) ((array:indexer a) i0 i1 i2)))

(define (array-1d-set! a v i0)
  (vector-set! (array:vector a) ((array:indexer a) i0) v))

(define (array-2d-set! a v i0 i1)
  (vector-set! (array:vector a) ((array:indexer a) i0 i1) v))

(define (array-3d-set! a v i0 i1 i2)
  (vector-set! (array:vector a) ((array:indexer a) i0 i1 i2) v))

; STOP!  Do not read beyond this point on your first reading of
; this code -- you should simply assume that the rest of this file
; contains only the following single definition:
;
;   (define (array:optimize-linear-function f l) f)
;
; Of course everything would be pretty inefficient if this were really the
; case, but it isn't.  The following code takes advantage of the fact that
; you can learn everything there is to know from a linear function by
; simply probing around in its domain and observing its values -- then a
; more efficient equivalent can be constructed.

(define (array:optimize-linear-function f l)
  (let ((d (length l)))
    (cond
     ((= d 0)
      (array:0d-c (f)))
     ((= d 1)
      (let ((c (f 0)))
	(array:1d-c0 c (- (f 1) c))))
     ((= d 2)
      (let ((c (f 0 0)))
	(array:2d-c01 c (- (f 1 0) c) (- (f 0 1) c))))
     ((= d 3)
      (let ((c (f 0 0 0)))
	(array:3d-c012 c (- (f 1 0 0) c) (- (f 0 1 0) c) (- (f 0 0 1) c))))
     (else
      (let* ((v (map (lambda (x) 0) l))
	     (c (apply f v)))
	(let loop ((p v)
		   (old-val c)
		   (coefs '()))
	  (cond ((null? p)
		 (array:Nd-c* c (reverse coefs)))
		(else
		 (set-car! p 1)
		 (let ((new-val (apply f v)))
		   (loop (cdr p)
			 new-val
			 (cons (- new-val old-val) coefs)))))))))))

; 0D cases:

(define (array:0d-c c)
  (lambda () c))

; 1D cases:

(define (array:1d-c c)
  (lambda (i0) (+ c i0)))

(define (array:1d-0 n0)
  (cond ((= 1 n0) +)
	(else (lambda (i0) (* n0 i0)))))

(define (array:1d-c0 c n0)
  (cond ((= 0 c) (array:1d-0 n0))
	((= 1 n0) (array:1d-c c))
	(else (lambda (i0) (+ c (* n0 i0))))))

; 2D cases:

(define (array:2d-0 n0)
  (lambda (i0 i1) (+ (* n0 i0) i1)))

(define (array:2d-1 n1)
  (lambda (i0 i1) (+ i0 (* n1 i1))))

(define (array:2d-c0 c n0)
  (lambda (i0 i1) (+ c (* n0 i0) i1)))

(define (array:2d-c1 c n1)
  (lambda (i0 i1) (+ c i0 (* n1 i1))))

(define (array:2d-01 n0 n1)
  (cond ((= 1 n0) (array:2d-1 n1))
	((= 1 n1) (array:2d-0 n0))
	(else (lambda (i0 i1) (+ (* n0 i0) (* n1 i1))))))

(define (array:2d-c01 c n0 n1)
  (cond ((= 0 c) (array:2d-01 n0 n1))
	((= 1 n0) (array:2d-c1 c n1))
	((= 1 n1) (array:2d-c0 c n0))
	(else (lambda (i0 i1) (+ c (* n0 i0) (* n1 i1))))))

; 3D cases:

(define (array:3d-01 n0 n1)
  (lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) i2)))

(define (array:3d-02 n0 n2)
  (lambda (i0 i1 i2) (+ (* n0 i0) i1 (* n2 i2))))

(define (array:3d-12 n1 n2)
  (lambda (i0 i1 i2) (+ i0 (* n1 i1) (* n2 i2))))

(define (array:3d-c12 c n1 n2)
  (lambda (i0 i1 i2) (+ c i0 (* n1 i1) (* n2 i2))))

(define (array:3d-c02 c n0 n2)
  (lambda (i0 i1 i2) (+ c (* n0 i0) i1 (* n2 i2))))

(define (array:3d-c01 c n0 n1)
  (lambda (i0 i1 i2) (+ c (* n0 i0) (* n1 i1) i2)))

(define (array:3d-012 n0 n1 n2)
  (cond ((= 1 n0) (array:3d-12 n1 n2))
	((= 1 n1) (array:3d-02 n0 n2))
	((= 1 n2) (array:3d-01 n0 n1))
	(else (lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) (* n2 i2))))))

(define (array:3d-c012 c n0 n1 n2)
  (cond ((= 0 c) (array:3d-012 n0 n1 n2))
	((= 1 n0) (array:3d-c12 c n1 n2))
	((= 1 n1) (array:3d-c02 c n0 n2))
	((= 1 n2) (array:3d-c01 c n0 n1))
	(else (lambda (i0 i1 i2) (+ c (* n0 i0) (* n1 i1) (* n2 i2))))))

; ND cases:

(define (array:Nd-* coefs)
  (lambda indices (apply + (map * coefs indices))))

(define (array:Nd-c* c coefs)
  (cond ((= 0 c) (array:Nd-* coefs))
	(else (lambda indices (apply + c (map * coefs indices))))))
Commit	Line	Data
9ddacf86 KN	1	;;;;"array.scm" Arrays for Scheme
	2	; Copyright (C) 1993 Alan Bawden
	3	;
	4	; Permission to copy this software, to redistribute it, and to use it
	5	; for any purpose is granted, subject to the following restrictions and
	6	; understandings.
	7	;
	8	; 1. Any copy made of this software must include this copyright notice
	9	; in full.
	10	;
	11	; 2. Users of this software agree to make their best efforts (a) to
	12	; return to me any improvements or extensions that they make, so that
	13	; these may be included in future releases; and (b) to inform me of
	14	; noteworthy uses of this software.
	15	;
	16	; 3. I have made no warrantee or representation that the operation of
	17	; this software will be error-free, and I am under no obligation to
	18	; provide any services, by way of maintenance, update, or otherwise.
	19	;
	20	; 4. In conjunction with products arising from the use of this material,
	21	; there shall be no use of my name in any advertising, promotional, or
	22	; sales literature without prior written consent in each case.
	23	;
	24	; Alan Bawden
	25	; MIT Room NE43-510
	26	; 545 Tech. Sq.
	27	; Cambridge, MA 02139
	28	; Alan@LCS.MIT.EDU
	29
	30	(require 'record)
	31
	32	;(declare (usual-integrations))
	33
	34	(define array:rtd
	35	(make-record-type "Array"
	36	'(indexer ; Must be a -linear- function!
	37	shape ; Inclusive bounds: ((lower upper) ...)
	38	vector ; The actual contents
	39	)))
	40
	41	(define array:indexer (record-accessor array:rtd 'indexer))
	42	(define array-shape (record-accessor array:rtd 'shape))
	43	(define array:vector (record-accessor array:rtd 'vector))
	44
	45	(define array? (record-predicate array:rtd))
	46
	47	(define (array-rank obj)
	48	(if (array? obj) (length (array-shape obj)) 0))
	49
	50	(define (array-dimensions ra)
	51	(map (lambda (ind) (if (zero? (car ind)) (+ 1 (cadr ind)) ind))
	52	(array-shape ra)))
	53
	54	(define array:construct
	55	(record-constructor array:rtd '(shape vector indexer)))
	56
	57	(define (array:compute-shape specs)
	58	(map (lambda (spec)
	59	(cond ((and (integer? spec)
	60	(< 0 spec))
	61	(list 0 (- spec 1)))
	62	((and (pair? spec)
	63	(pair? (cdr spec))
	64	(null? (cddr spec))
65	(integer? (car spec))
66	(integer? (cadr spec))
67	(<= (car spec) (cadr spec)))
68	spec)
69	(else (slib:error "array: Bad array dimension: " spec))))
70	specs))
71
72	(define (make-array initial-value . specs)
73	(let ((shape (array:compute-shape specs)))
74	(let loop ((size 1)
75	(indexer (lambda () 0))
76	(l (reverse shape)))
77	(if (null? l)
78	(array:construct shape
79	(make-vector size initial-value)
80	(array:optimize-linear-function indexer shape))
81	(loop (* size (+ 1 (- (cadar l) (caar l))))
82	(lambda (first-index . rest-of-indices)
83	(+ (* size (- first-index (caar l)))
84	(apply indexer rest-of-indices)))
85	(cdr l))))))
86
87	(define (make-shared-array array mapping . specs)
88	(let ((new-shape (array:compute-shape specs))
89	(old-indexer (array:indexer array)))
90	(let check ((indices '())
91	(bounds (reverse new-shape)))
92	(cond ((null? bounds)
93	(array:check-bounds array (apply mapping indices)))
94	(else
95	(check (cons (caar bounds) indices) (cdr bounds))
96	(check (cons (cadar bounds) indices) (cdr bounds)))))
97	(array:construct new-shape
98	(array:vector array)
99	(array:optimize-linear-function
100	(lambda indices
101	(apply old-indexer (apply mapping indices)))
102	new-shape))))
103
104	(define (array:in-bounds? array indices)
105	(let loop ((indices indices)
106	(shape (array-shape array)))
107	(if (null? indices)
108	(null? shape)
109	(let ((index (car indices)))
110	(and (not (null? shape))
111	(integer? index)
112	(<= (caar shape) index (cadar shape))
113	(loop (cdr indices) (cdr shape)))))))
114
115	(define (array:check-bounds array indices)
116	(or (array:in-bounds? array indices)
117	(slib:error "array: Bad indices for " array indices)))
118
119	(define (array-ref array . indices)
120	(array:check-bounds array indices)
121	(vector-ref (array:vector array)
122	(apply (array:indexer array) indices)))
123
124	(define (array-set! array new-value . indices)
125	(array:check-bounds array indices)
126	(vector-set! (array:vector array)
127	(apply (array:indexer array) indices)
128	new-value))
129
130	(define (array-in-bounds? array . indices)
131	(array:in-bounds? array indices))
132
133	; Fast versions of ARRAY-REF and ARRAY-SET! that do no error checking,
134	; and don't cons intermediate lists of indices:
135
136	(define (array-1d-ref a i0)
137	(vector-ref (array:vector a) ((array:indexer a) i0)))
138
139	(define (array-2d-ref a i0 i1)
140	(vector-ref (array:vector a) ((array:indexer a) i0 i1)))
141
142	(define (array-3d-ref a i0 i1 i2)
143	(vector-ref (array:vector a) ((array:indexer a) i0 i1 i2)))
144
145	(define (array-1d-set! a v i0)
146	(vector-set! (array:vector a) ((array:indexer a) i0) v))
147
148	(define (array-2d-set! a v i0 i1)
149	(vector-set! (array:vector a) ((array:indexer a) i0 i1) v))
150
151	(define (array-3d-set! a v i0 i1 i2)
152	(vector-set! (array:vector a) ((array:indexer a) i0 i1 i2) v))
153
154	; STOP! Do not read beyond this point on your first reading of
155	; this code -- you should simply assume that the rest of this file
156	; contains only the following single definition:
157	;
158	; (define (array:optimize-linear-function f l) f)
159	;
160	; Of course everything would be pretty inefficient if this were really the
161	; case, but it isn't. The following code takes advantage of the fact that
162	; you can learn everything there is to know from a linear function by
163	; simply probing around in its domain and observing its values -- then a
164	; more efficient equivalent can be constructed.
165
166	(define (array:optimize-linear-function f l)
167	(let ((d (length l)))
168	(cond
169	((= d 0)
170	(array:0d-c (f)))
171	((= d 1)
172	(let ((c (f 0)))
173	(array:1d-c0 c (- (f 1) c))))
174	((= d 2)
175	(let ((c (f 0 0)))
176	(array:2d-c01 c (- (f 1 0) c) (- (f 0 1) c))))
177	((= d 3)
178	(let ((c (f 0 0 0)))
179	(array:3d-c012 c (- (f 1 0 0) c) (- (f 0 1 0) c) (- (f 0 0 1) c))))
180	(else
181	(let* ((v (map (lambda (x) 0) l))
182	(c (apply f v)))
183	(let loop ((p v)
184	(old-val c)
185	(coefs '()))
186	(cond ((null? p)
187	(array:Nd-c* c (reverse coefs)))
188	(else
189	(set-car! p 1)
190	(let ((new-val (apply f v)))
191	(loop (cdr p)
192	new-val
193	(cons (- new-val old-val) coefs)))))))))))
194
195	; 0D cases:
196
197	(define (array:0d-c c)
198	(lambda () c))
199
200	; 1D cases:
201
202	(define (array:1d-c c)
203	(lambda (i0) (+ c i0)))
204
205	(define (array:1d-0 n0)
206	(cond ((= 1 n0) +)
207	(else (lambda (i0) (* n0 i0)))))
208
209	(define (array:1d-c0 c n0)
210	(cond ((= 0 c) (array:1d-0 n0))
211	((= 1 n0) (array:1d-c c))
212	(else (lambda (i0) (+ c (* n0 i0))))))
213
214	; 2D cases:
215
216	(define (array:2d-0 n0)
217	(lambda (i0 i1) (+ (* n0 i0) i1)))
218
219	(define (array:2d-1 n1)
220	(lambda (i0 i1) (+ i0 (* n1 i1))))
221
222	(define (array:2d-c0 c n0)
223	(lambda (i0 i1) (+ c (* n0 i0) i1)))
224
225	(define (array:2d-c1 c n1)
226	(lambda (i0 i1) (+ c i0 (* n1 i1))))
227
228	(define (array:2d-01 n0 n1)
229	(cond ((= 1 n0) (array:2d-1 n1))
230	((= 1 n1) (array:2d-0 n0))
231	(else (lambda (i0 i1) (+ (* n0 i0) (* n1 i1))))))
232
233	(define (array:2d-c01 c n0 n1)
234	(cond ((= 0 c) (array:2d-01 n0 n1))
235	((= 1 n0) (array:2d-c1 c n1))
236	((= 1 n1) (array:2d-c0 c n0))
237	(else (lambda (i0 i1) (+ c (* n0 i0) (* n1 i1))))))
238
239	; 3D cases:
240
241	(define (array:3d-01 n0 n1)
242	(lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) i2)))
243
244	(define (array:3d-02 n0 n2)
245	(lambda (i0 i1 i2) (+ (* n0 i0) i1 (* n2 i2))))
246
247	(define (array:3d-12 n1 n2)
248	(lambda (i0 i1 i2) (+ i0 (* n1 i1) (* n2 i2))))
249
250	(define (array:3d-c12 c n1 n2)
251	(lambda (i0 i1 i2) (+ c i0 (* n1 i1) (* n2 i2))))
252
253	(define (array:3d-c02 c n0 n2)
254	(lambda (i0 i1 i2) (+ c (* n0 i0) i1 (* n2 i2))))
255
256	(define (array:3d-c01 c n0 n1)
257	(lambda (i0 i1 i2) (+ c (* n0 i0) (* n1 i1) i2)))
258
259	(define (array:3d-012 n0 n1 n2)
260	(cond ((= 1 n0) (array:3d-12 n1 n2))
261	((= 1 n1) (array:3d-02 n0 n2))
262	((= 1 n2) (array:3d-01 n0 n1))
263	(else (lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) (* n2 i2))))))
264
265	(define (array:3d-c012 c n0 n1 n2)
266	(cond ((= 0 c) (array:3d-012 n0 n1 n2))
267	((= 1 n0) (array:3d-c12 c n1 n2))
268	((= 1 n1) (array:3d-c02 c n0 n2))
269	((= 1 n2) (array:3d-c01 c n0 n1))
270	(else (lambda (i0 i1 i2) (+ c (* n0 i0) (* n1 i1) (* n2 i2))))))
271
272	; ND cases:
273
274	(define (array:Nd-* coefs)
275	(lambda indices (apply + (map * coefs indices))))
276
277	(define (array:Nd-c* c coefs)
278	(cond ((= 0 c) (array:Nd-* coefs))
279	(else (lambda indices (apply + c (map * coefs indices))))))