;;;;"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))))))