;;;; Copyright (c) 2017, Jeremy Steward
;;;; All rights reserved.
;;;;
;;;; Redistribution and use in source and binary forms, with or without
;;;; modification, are permitted provided that the following conditions are met:
;;;;
;;;; 1. Redistributions of source code must retain the above copyright notice,
;;;; this list of conditions and the following disclaimer.
;;;;
;;;; 2. Redistributions in binary form must reproduce the above copyright notice,
;;;; this list of conditions and the following disclaimer in the documentation
;;;; and/or other materials provided with the distribution.
;;;;
;;;; 3. Neither the name of the copyright holder nor the names of its
;;;; contributors may be used to endorse or promote products derived from this
;;;; software without specific prior written permission.
;;;;
;;;; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
;;;; AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
;;;; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
;;;; ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
;;;; LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
;;;; CONSEQUENTIAL DAMAGES  INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
;;;; SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
;;;; INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
;;;; CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
;;;; ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
;;;; POSSIBILITY OF SUCH DAMAGE.

;; NOTE: Because of the way I include the files in the top-level module, some
;; definitions for helper procedures won't be used here. Check
;; `generalized-arrays/generalized-arrays-base.scm` for these procedures /
;; record-types, as they should be located there.

;;; The whole array

;; Returns a newly allocated array with the same bounds as the arrays; it is an
;; error if they do not all have the same bounds. For each valid index value,
;; proc is invoked, passing each corresponding element of the arrays. Whatever
;; proc returns becomes the value of the storage element corresponding to that
;; index in the result array. The order of invocations of proc is not
;; specified.
(define (array-map storage-class proc array #!rest arrays)
  (check-array 'array-map array)
  (for-each (lambda (a) (check-array 'array-map a))
            arrays)
  (let ((all-arrays (cons array arrays))
        (new-array (make-array storage-class (%array-shape array))))
    (let ((data (%array-storage-object new-array))
          (storage-set! (storage-class-mutator
                          (%array-storage-class new-array))))
      (%array-for-each-index
        (lambda (index)
          (let ((i (%array-index->storage-index array index)))
           (storage-set!
             data
             i
             (apply proc (map (lambda (a)
                                (let ((storage-ref (storage-class-accessor
                                                     (%array-storage-class a)))
                                      (storage-object (%array-storage-object a))
                                      (storage-index (%array-index->storage-index a index)))
                                  (storage-ref storage-object storage-index)))
                              all-arrays)))))
        (%array-lower-bound new-array)
        (%array-upper-bound new-array)))
    (unless (%array-mutable? array)
      (set! (%array-mutable? new-array) #f))
    new-array))


;; It is an error if the arrays do not all have the same bounds. For each valid
;; index value, proc is invoked, passing each corresponding element of the
;; arrays.  Whatever proc returns becomes the value of the storage element
;; corresponding to that index in the first specified array. The order of
;; invocations of proc is not specified. Returns an undefined value.
(define (array-map! proc array #!rest arrays)
  (check-array 'array-map! array)
  (for-each (lambda (x) (check-array 'array-map! x))
            arrays)
  (unless (array-mutable? array)
    (error "(array-map!): immutable arrays cannot be mapped using array-map!"))
  (let ((all-arrays (cons array arrays))
        (data (%array-storage-object array))
        (storage-set! (storage-class-mutator
                        (%array-storage-class array))))
    (%array-for-each-index
      (lambda (index)
        (let ((i (%array-index->storage-index array index)))
         (storage-set!
           data
           i
           (apply proc
                  (map (lambda (a)
                         (let ((storage-ref (storage-class-accessor
                                              (%array-storage-class a)))
                               (storage-object (%array-storage-object a))
                               (storage-index (%array-index->storage-index a index)))
                           (storage-ref storage-object storage-index)))
                       all-arrays)))))
      (%array-lower-bound array)
      (%array-upper-bound array))
    array))


(define (%array-fold proc seed array arrays)
  (let ((acc seed)
        (all-arrays (cons array arrays)))
    (%array-for-each-index
      (lambda (index)
        (let ((args (map (lambda (a)
                           (let ((storage-ref (storage-class-accessor
                                                (%array-storage-class a)))
                                 (storage-object (%array-storage-object a))
                                 (storage-index (%array-index->storage-index a index)))
                             (storage-ref storage-object storage-index)))
                         all-arrays)))
          (set! acc (apply proc acc args))))
      (%array-lower-bound array)
      (%array-upper-bound array))
    acc))

;; Returns a newly allocated array with the same bounds as the arrays; it is an
;; error if they do not all have the same bounds. For each valid index value,
;; proc is invoked in lexicographic order, passing each corresponding element
;; of the arrays and a seed, whose initial value is seed. Proc returns two
;; values, the value of the storage element corresponding to that index in the
;; result array, and the new value of the seed.
(define (array-fold proc seed array #!rest arrays)
  (check-array 'array-map! array)
  (for-each (lambda (x) (check-array 'array-map! x))
            arrays)
  (%array-fold proc seed array arrays))


(define (remove-axis-from-index index axis)
  (match (list index axis)
    ((#(i) 0) (vector 0))
    ((#(i j) 0) (vector j))
    ((#(i j) 1) (vector i))
    ((#(i j k) 0) (vector j k))
    ((#(i j k) 1) (vector i k))
    ((#(i j k) 2) (vector i j))
    ((#(i j k v) 0) (vector j k v))
    ((#(i j k v) 1) (vector i k v))
    ((#(i j k v) 2) (vector i j v))
    ((#(i j k v) 3) (vector i j k))
    (_
      (let ((new-index (make-vector (sub1 (vector-length index)))))
       (vector-fold (lambda (i x)
                      (cond
                        ((eq? i axis) (add1 i))
                        ((> i axis)
                         (vector-set! new-index (sub1 i) (vector-ref index i))
                         (add1 i))
                        (else
                          (vector-set! new-index i (vector-ref index i))
                          (add1 i))))
                    0
                    index)
       new-index))))

(define (%array-reduce proc seed array axis)
  (let* ((shape (%array-shape array))
         (new-array (make-array (%array-storage-class array)
                                (if (= (%array-rank array) 1)
                                  (vector 1)
                                  (remove-axis-from-index shape axis))
                                seed)))
    (%array-for-each-index
      (lambda (index)
        (let ((new-index (remove-axis-from-index index axis)))
         (%array-set! new-array
                      new-index
                      (proc (%array-ref new-array new-index)
                            (%array-ref array index)))))
      (%array-lower-bound array)
      (%array-upper-bound array))
    new-array))

;; Returns a newly allocated array whose rank is one less than the rank of
;; array, by combining all the groups of elements of length n (where the
;; default is the extent of axis) along axis using proc, a two-argument
;; procedure. The order and number of invocations of proc is unspecified. If
;; there is only one such element, it is unchanged. (APL reduce.)
(define (array-reduce proc seed array axis)
  (check-array 'array-reduce array)
  (check-natural-fixnum 'array-reduce axis 'axis)
  (unless (and (fx>= axis 0)
               (fx< axis (%array-rank array)))
    (error "(array-reduce): axis must be between 0 and the array rank." axis))
  (%array-reduce proc seed array axis))


(define (%array-cumulate proc seed array axis)
  (let* ((shape (%array-shape array))
         (new-array (make-array (%array-storage-class array) shape seed)))
    (%array-for-each-index
      (lambda (index)
        (let ((last-index (vector-copy index)))
          (unless (zero? (vector-ref index axis))
           (vector-set! last-index
                       axis
                       (sub1 (vector-ref last-index axis))))
          (%array-set! new-array
                       index
                       (proc (%array-ref new-array last-index)
                             (%array-ref array index)))))
      (%array-lower-bound array)
      (%array-upper-bound array))
    new-array))

;; Returns a newly allocated array whose bounds are the same as the bounds of
;; array. Each element along axis is constructed by reducing (as if by
;; array-reduce) successive prefixes of the elements of array along that axis.
;; (APL scan.)
(define (array-cumulate proc seed array axis)
  (check-array 'array-cumulate array)
  (check-natural-fixnum 'array-cumulate axis 'axis)
  (unless (and (fx>= axis 0)
               (fx< axis (%array-rank array)))
    (error "(array-cumulate): axis must be between 0 and the array rank." axis))
  (%array-cumulate proc seed array axis))

;; Returns an exact integer containing the number of elements in array that
;; satisfy pred.
(define (array-count pred array)
  (array-fold (lambda (len x)
                (if (pred x)
                  (add1 len)
                  len))
              0
              array))

;; Returns the index of the first element of array in lexicographic order that
;; satisfies pred.
(define (array-index pred array)
  (check-array 'array-index array)
  (call/cc
    (lambda (k)
      (%array-for-each-index
        (lambda (i)
          (let ((x (%array-ref array i)))
           (when (pred x) (k i))))
        (%array-lower-bound array)
        (%array-upper-bound array)))))

;; Returns the first element of array in lexicographic order that satisfies
;; pred.
(define (array-find pred array)
  (array-ref array (array-index pred array)))


(define (compress-index index booleans axis)
  (let* ((new-index (vector-copy index))
         (axis-value (vector-ref index axis))
         (nitems-skipped
           (vector-count not (subvector booleans 0 axis-value))))
    (vector-set! new-index axis (fx- axis-value nitems-skipped))
    new-index))

(define (%array-compress array booleans axis)
  (let ((new-axis-length (vector-count identity booleans))
        (new-shape (vector-copy (%array-shape array))))
    (vector-set! new-shape axis new-axis-length)
    (let ((new-array (make-array (%array-storage-class array) new-shape)))
     (%array-for-each-index
       (lambda (index)
         (when (vector-ref booleans (vector-ref index axis))
           (let ((compressed-index (compress-index index booleans axis)))
            (array-set! new-array compressed-index (array-ref array index)))))
       (%array-lower-bound array)
       (%array-upper-bound array))
     new-array)))

;; Returns an array with the same bounds as array except possibly along the
;; axis dimension. The array is sliced along axis and the elements of booleans
;; (a vector of boolean values) are used to decide which slices to incorporate
;; into the result: if the corresponding boolean is #t, the slice is
;; incorporated, otherwise not. (APL compress.)
(define (array-compress array booleans axis)
  (check-array 'array-compress array)
  (check-vector 'array-compress booleans)
  (unless (and (fx>= axis 0)
               (fx< axis (%array-rank array)))
    (error "(array-compress): axis must be between 0 and the array rank." axis))
  (unless (fx= (vector-length booleans)
               (vector-ref (%array-shape array) axis))
    (error "(array-compress): booleans must be same size as axis length"
           (vector-length booleans)
           axis))
  (%array-compress array booleans axis))


(define (expand-index index booleans axis)
  (let* ((new-index (vector-copy index))
         (axis-value (vector-ref index axis))
         (nitems-inserted
           (vector-count identity
                         (subvector booleans 0 axis-value))))
    (vector-set! new-index axis (fx+ axis-value nitems-inserted))
    new-index))

(define (%array-expand array booleans nil axis)
  (let ((new-axis-length (fx+ (vector-ref (%array-shape array) axis)
                              (vector-count identity booleans)))
        (new-shape (vector-copy (%array-shape array))))
    (vector-set! new-shape axis new-axis-length)
    (let ((new-array (make-array (%array-storage-class array) new-shape)))
     (%array-for-each-index
       (lambda (index)
         (let ((expanded-index (expand-index index booleans axis)))
          (array-set! new-array expanded-index (array-ref array index))
          (when (vector-ref booleans (vector-ref index axis))
            (vector-set! expanded-index axis
                         (add1 (vector-ref expanded-index axis)))
            (let ((nil-index (remove-axis-from-index index axis)))
             (array-set! new-array expanded-index (array-ref nil nil-index))))))
       (%array-lower-bound array)
       (%array-upper-bound array))
     new-array)))

;; Returns an array with the same bounds as array except possibly along the
;; axis dimension. Array is sliced along axis and the elements of booleans (a
;; vector of boolean values) are used to decide where, if anywhere, nil is to
;; be interpolated: if the corresponding boolean is #t, nil is interpolated,
;; otherwise the next slice is incorporated. It is an error if the size of
;; booleans is not to the extent of the axis dimension in the result. It is
;; also an error if nil does not have the same bounds as a slice. (APL expand.)
(define (array-expand array booleans nil axis)
  (check-array 'array-expand array)
  (check-array 'array-expand nil)
  (unless (vector= =
                   (%array-shape nil)
                   (remove-axis-from-index (%array-shape array) axis))
    (error "(array-expand): nil must be same shape as array along axis."
           (%array-shape nil)
           (remove-axis-from-index (%array-shape array) axis)))
  (check-vector 'array-expand booleans)
  (unless (and (fx>= axis 0)
               (fx< axis (%array-rank array)))
    (error "(array-expand): axis must be between 0 and the array rank." axis))
  (unless (fx= (vector-length booleans)
               (vector-ref (%array-shape array) axis))
    (error "(array-extend): booleans must be same size as axis length"
           (vector-length booleans)
           axis))
  (%array-expand array booleans nil axis))


(define (%array-rearrange array order axis)
  (let ((new-array (make-array (%array-storage-class array)
                               (%array-shape array)))
        (axis-length (vector-ref (%array-shape array) axis))
        (lower-bound (vector-copy (%array-lower-bound array)))
        (upper-bound (vector-copy (%array-upper-bound array)))
        (new-lower-bound (vector-copy (%array-lower-bound array)))
        (new-upper-bound (vector-copy (%array-upper-bound array))))
    (vector-for-each
      (lambda (i)
        (vector-set! lower-bound axis i)
        (vector-set! upper-bound axis (add1 i))
        (vector-set! new-upper-bound
                     axis
                     (add1 (vector-ref new-lower-bound axis)))
        (array-copy! (array-slice array lower-bound upper-bound)
                     (array-slice new-array new-lower-bound new-upper-bound))
        (vector-set! new-lower-bound
                     axis
                     (add1 (vector-ref new-lower-bound axis))))
      order)
    new-array))

;; Returns an array with the same bounds as array. Array is sliced along the
;; axis dimension, and the slices are reassembled in the order given by vector.
;; The slice whose number appears in the first element of vector appears first
;; in the result, and so on. It is an error if vector is not a vector of exact
;; integers.  (Generalized version of APL rotate.)
(define (array-rearrange array order axis)
  (check-array 'array-rearrange array)
  (check-vector 'array-rearrange order)
  (unless (and (fx>= axis 0)
               (fx< axis (%array-rank array)))
    (error "(array-rearrange): axis must be between 0 and the array rank." axis))
  (unless (and (fx= (vector-ref (%array-shape array) axis)
                    (vector-length order))
               (vector-every (lambda (x)
                               (fx< x (vector-ref (%array-shape array) axis)))
                             order))
    (error "(array-rearrange): All values in the order vector must be less than the axis length."))
  (%array-rearrange array order axis))


;;; Transformations

(define (%array-transform proc shape array)
  (let* ((new-array (make-array (%array-storage-class array)
                                shape))
         (new-data (%array-storage-object new-array))
         (data (%array-storage-object array)))
    (%array-for-each-index
      (lambda (index)
        (let ((from-i (%array-index->storage-index array (proc index)))
              (to-i (%array-index->storage-index new-array index))
              (storage-set! (storage-class-mutator
                              (%array-storage-class new-array)))
              (storage-ref (storage-class-accessor
                             (%array-storage-class array))))
          (storage-set! new-data to-i (storage-ref data from-i))))
      (%array-lower-bound new-array)
      (%array-upper-bound new-array))
    new-array))

;; The procedure proc specifies an affine function that returns an index of
;; array when given an index of the returned array. The array does not retain a
;; dependency on proc. (SRFI 25 share-array.)
(define (array-transform proc shape array)
  (check-array 'array-transform array)
  (%array-transform proc shape array))

;; Returns an array with the specified bounds whose elements in lexicographic
;; order are the same as the elements of array in lexicographic order. It is an
;; error if there are too many or too few elements. (APL reshape.)
(define (array-reshape array shape)
  (check-array 'array-reshape array)
  (unless (= (vector-fold fx* 1 (%array-shape array))
             (vector-fold fx* 1 shape))
    (error "(array-reshape): new shape must have equivalent number of elements as old shape."
           (vector-fold fx* 1 shape)
           (vector-fold fx* 1 (%array-shape array))))
  (let ((reshaped-array (array-copy array (%array-mutable? array))))
   (set! (%array-shape reshaped-array) shape)
   (set! (%array-stride reshaped-array)
         (shape->stride shape (row-major? array)))
   reshaped-array))

;; Returns a one-dimensional array which contains the diagonal elements of
;; array (that is, the elements whose indices are all the same integer).
(define (array-diagonal array)
  (check-array 'array-diagonal array)
  (if (= 1 (%array-rank array))
    array
    (%make-array (%array-storage-class array)
                 (vector (apply min (vector->list (%array-shape array))))
                 (vector (vector-fold fx+ 0 (%array-stride array)))
                 (%array-offset array)
                 (%array-mutable? array)
                 (%array-storage-object array))))

;; Returns an array with the same bounds as array, but whose elements on the
;; specified axis are reversed. (APL reverse.)
(define (array-reverse array axis)
  (check-array 'array-reverse array)
  (unless (and (fx>= axis 0)
               (fx< axis (%array-rank array)))
    (error "(array-rearrange): axis must be between 0 and the array rank." axis))
  (let ((limit (sub1 (vector-ref (%array-shape array) axis))))
   (array-transform (lambda (index)
                      (let ((j (vector-copy index)))
                       (vector-set! j
                                    axis
                                    (fx- limit (vector-ref j axis)))
                       j))
                    (array-shape array)
                    array)))

;; Returns an array whose axes appear in the reverse order of the axes of
;; array.  This implies that the upper and lower bound are the reverse of the
;; bounds of array. (APL monadic transpose.)
(define (array-transpose array)
  (check-array 'array-transpose array)
  (let ((shape (vector-copy (%array-shape array)))
        (stride (vector-copy (%array-stride array))))
   (vector-reverse! shape)
   (vector-reverse! stride)
   (%make-array (%array-storage-class array)
                shape
                stride
                (%array-offset array)
                (%array-mutable? array)
                (%array-storage-object array))))


(define (%array-rearrange-axes array order)
  (let ((shape (%array-shape array))
        (stride (%array-stride array)))
   (%make-array (%array-storage-class array)
               (vector-map (cute vector-ref shape <>) order)
               (vector-map (cute vector-ref stride <>) order)
               (%array-offset array)
               (%array-mutable? array)
               (%array-storage-object array))))

;; Returns an array whose axes are an arbitrary permutation of the axes of
;; array.  Vec specifies how to do the permutation: the axis whose number
;; appears in the first element of vector appears as the first axis of the
;; result, and so on. (APL dyadic transpose with integer-valued vector.)
(define (array-rearrange-axes array order)
  (check-array 'array-rearrange-axes array)
  (unless (= (vector-length order)
             (%array-rank array))
    (error "(array-rearrange-axes): order must be same length as rank of array."
           (vector-length order)
           (%array-rank array)))
  (unless (vector= fx=
                   (sort order fx<)
                   (vector-unfold values (%array-rank array)))
    (error "(array-rearrange-axes): all axes must be present in order vector."
           order))
  (%array-rearrange-axes array order))


(define (%array-slice array start end step)
  (%make-array (%array-storage-class array)
               (vector-map (compose inexact->exact
                                    ceiling
                                    /)
                           (vector-map fx- end start)
                           step)
               (vector-map fx* (%array-stride array) step)
               (%array-index->storage-index array start)
               (%array-mutable? array)
               (%array-storage-object array)))

;; Returns a smaller array with the same rank as array whose elements begin at
;; the "lower left" corner specified by start and end at the "upper right"
;; corner specified byend. Unlike array-copy, the result shares its storage
;; object with array. (APL take and drop.)
(define (array-slice array start end #!optional (step *UNSPECIFIED*))
  (check-array 'array-slice array)
  (check-array-index 'array-slice start)
  (check-array-index 'array-slice end)
  (unless (and (vector= fx>= start (%array-lower-bound array))
               (vector= fx< start (%array-upper-bound array))
               (vector= fx< start end))
    (error
      "(array-slice): start index must be greater-than-or-equal to the lower bound and less than the end / upper-bound."
      start
      end
      (%array-lower-bound array)
      (%array-upper-bound array)))
  (unless (vector= fx<= end (%array-upper-bound array))
    (error "(array-slice): end index cannot be greater than the upper-bound of array."
           end
           (%array-shape array)))
  (let ((step (if (eq? step *UNSPECIFIED*)
                (vector-map (lambda _ 1) start)
                step)))
    (unless (vector= >= step (vector-map (lambda _ 0) start))
      (error "(array-slice): step size must be greater than or equal to zero." step))
    (%array-slice array start end step)))


(define (squeeze-vector vec axes)
  (let ((limit (vector-length vec)))
   (let loop ((i 0)
              (acc '()))
     (cond
       ((fx>= i limit)
        (reverse-list->vector acc))
       ((vector-index (cute = <> i) axes)
        (loop (add1 i) acc))
       (else
         (loop (add1 i) (cons (vector-ref vec i) acc)))))))

(define (%array-squeeze array axes)
  (%make-array (%array-storage-class array)
               (squeeze-vector (%array-shape array) axes)
               (squeeze-vector (%array-stride array) axes)
               (%array-offset array)
               (%array-mutable? array)
               (%array-storage-object array)))

;; Returns an array with the ranks specified by the elements of vector removed
;; from array. It is an error if the extents of the specified ranks are not
;; equal to 1.  (NumPy squeeze)
(define (array-squeeze array #!optional (axes *UNSPECIFIED*))
  (check-array 'array-squeeze array)
  (if (eq? axes *UNSPECIFIED*)
    (let* ((shape (%array-shape array))
           (axes (list->vector
                  (filter (lambda (x) (= 1 (vector-ref shape x)))
                          (iota (vector-length shape))))))
      (%array-squeeze array axes))
    (begin
      (vector-for-each
        (lambda (x)
          (unless (= 1 (vector-ref (%array-shape array) x))
            (error "(array-squeeze): Can only squeeze axes with length of 1.")))
        axes)
      (%array-squeeze array axes))))


(define (unsqueeze-vector vec rank)
  (let* ((new-length (add1 (vector-length vec)))
         (result (make-vector new-length)))
    (let loop ((i 0))
     (cond
       ((fx> i new-length) result)
       ((fx= i rank)
        (vector-set! result i 1)
        (loop (add1 i)))
       (else
         (vector-set! result
                      i
                      (vector-ref vec
                                  (if (fx> rank i)
                                    i
                                    (sub1 i))))
         (loop (add1 i)))))))

(define (%array-unsqueeze array rank)
  (let ((shape (unsqueeze-vector (%array-shape array) rank)))
   (%make-array (%array-storage-class array)
                shape
                (shape->stride shape (row-major? array))
                (%array-offset array)
                (%array-mutable? array)
                (%array-storage-object array))))

;; Returns an array whose rank is one greater than the rank of array. This is
;; accomplished by inserting a new axis numbered 'axis' whose extent is (0:1).
;; (NumPy expand.)
(define (array-unsqueeze array rank)
  (check-array 'array-unsqueeze array)
  (unless (and (fx>= rank 0)
               (fx<= rank (%array-rank array)))
    (error "(array-unsqueeze): rank must be a value in [0, n] for an n-ranked array."
           rank))
  (%array-unsqueeze array rank))


;;; Copying and conversion

;; These procedures return arrays which do not share their storage objects with
;; any existing arrays.

;; When these procedures return arrays, both the array and the underlying
;; storage object are newly allocated

(define (%array-copy array mutable?)
  (let* ((new-array (make-array (%array-storage-class array)
                                (%array-shape array)))
         (new-data (%array-storage-object new-array))
         (data (%array-storage-object array)))
    (%array-for-each-index
      (lambda (index)
        (let ((from-i (%array-index->storage-index array index))
              (to-i (%array-index->storage-index new-array index))
              (storage-set! (storage-class-mutator
                              (%array-storage-class new-array)))
              (storage-ref (storage-class-accessor
                             (%array-storage-class array))))
          (storage-set! new-data to-i (storage-ref data from-i))))
      (%array-lower-bound array)
      (%array-upper-bound array))
    (set! (%array-mutable? new-array) mutable?)
    new-array))

;; Returns an array containing the elements of array.
;; The stride and offset are implementation-defined. The resulting array is
;; mutable if mutable? is true.  The returned array has the same storage class
;; as array.
(define (array-copy array #!optional (mutable? #t))
  (check-array 'array-copy array)
  (check-boolean 'array-copy mutable?)
  (%array-copy array mutable?))


(define (%array-copy! from to)
  (let ((from-data (%array-storage-object from))
        (from-offset (%array-offset from))
        (to-data (%array-storage-object to))
        (to-offset (%array-offset to)))
   (%array-for-each-index
    (lambda (index)
      (let ((from-i (%array-index->storage-index from index))
            (to-i (%array-index->storage-index to index))
            (storage-set! (storage-class-mutator
                            (%array-storage-class to)))
            (storage-ref (storage-class-accessor
                           (%array-storage-class from))))
       (storage-set! to-data
                     to-i
                     (storage-ref from-data from-i))))
    (%array-lower-bound from)
    (%array-upper-bound from)))
  to)

;; Copies the elements of array from the "to" argument to the "from" argument
;; starting at element "at".  It is an error if there are not enough elements
;; in "to" to make this possible. The storage objects of from and to need not
;; belong to the same storage classes.
(define (array-copy! from to #!optional (at *UNSPECIFIED*))
  (check-array 'array-copy! from)
  (check-array 'array-copy! to)
  (unless (array-mutable? to)
    (error "(array-copy!): Specified 'to' array is immutable."))
  (let* ((at (if (eq? at *UNSPECIFIED*) (%array-lower-bound to) at))
         (size-available (vector-map fx+ (%array-shape to) at)))
    (unless (vector= >= size-available (%array-shape from))
      (error "There is insufficient space available to copy into array."
             size-available
             (%array-shape from)))
    (%array-copy! from (array-slice to at (%array-shape to)))))


(define (%array-append storage-class axis arrays)
  (let* ((new-axis-length (fold (lambda (x k)
                                  (fx+ k (vector-ref (%array-shape x) axis)))
                                0
                                arrays))
         (new-shape (let ((shape (vector-copy (%array-shape (car arrays)))))
                      (vector-set! shape axis new-axis-length)
                      shape))
         (new-array (make-array storage-class new-shape)))
    (fold (lambda (array length)
            (let ((offset (%array-lower-bound array)))
              (vector-set! offset axis length)
              (array-copy!
               array
               (array-slice new-array
                            offset
                            (vector-map fx+ offset (%array-shape array)))))
            (fx+ length (vector-ref (%array-shape array) axis)))
          0
          arrays)
    new-array))

;; Returns a newly allocated array consisting of the arrays concatenated along
;; axis. It is an error unless the storage classes of the arrays are the same,
;; and the result has the same storage class. It is also an error unless the
;; bounds of all the arrays are the same, with the possible exception of axis.
;; The axisth element of the lower bound of the result is 0; the corresponding
;; element of the upper bound is the sum of the extents of the arrays.
(define (array-append axis array #!rest arrays)
  (check-array 'array-append array)
  (for-each
   (lambda (a)
     (check-array 'array-append a)
     (unless (vector= =
                      (remove-axis-from-index (%array-shape array) axis)
                      (remove-axis-from-index (%array-shape a) axis))
       (error "(array-append): shapes must be congruent for all arrays."
              (%array-shape a)
              (%array-shape array))))
   arrays)
  (unless (or (fx>= axis 0)
              (fx< axis (%array-rank array)))
    (error "(array-append): axis must be between 0 and array rank."))
  (%array-append (%array-storage-class array) axis (cons array arrays)))


;; Append repeat copies of array along axis axis, as if by array-append.
;; (Variant of NumPy repeat.)
(define (array-repeat array axis repeat)
  (apply array-append axis (make-list repeat array)))

;; Returns an array with the same bounds and elements as array, but whose
;; storage class is specified by storage-class.
(define (array-reclassify array storage-class)
  (check-array 'array-reclassify array)
  (let ((new-array (make-array storage-class (%array-shape array))))
    (array-copy array new-array)
    array))

;; Returns a newly allocated vector/list whose elements are also newly
;; constructed vectors/lists, and so on as far down as necessary to cover every
;; axis of the array. Bottom-level Scheme vectors/lists contain the elements of
;; array. Thus, if array has rank 1, the result is a vector/list; if the array
;; has rank 2, the result is a vector/list whose elements are vectors/lists,
;; and so on. As a special case, if array has rank 0, the sole element is
;; returned.
(define (array->nested-vector array)
  (check-array 'array->nested-vector array)
  (let loop ((subarray array))
    (let ((vec (make-vector (vector-ref (%array-shape subarray) 0))))
     (if (= (array-rank subarray) 1)
       (for-each (lambda (i)
                   (vector-set! vec i (%array-ref subarray (vector i))))
                 (iota (vector-length vec)))
       (let ((lower-bound (vector-copy (%array-lower-bound subarray)))
             (upper-bound (vector-copy (%array-shape subarray))))
        (for-each
         (lambda (i)
           (vector-set! lower-bound 0 i)
           (vector-set! upper-bound 0 (add1 i))
           (vector-set! vec
                        i
                        (loop
                         (%array-squeeze (array-slice subarray
                                                      lower-bound
                                                      upper-bound)
                                         (vector 0)))))
         (iota (vector-length vec)))))
     vec)))

(define (array->nested-list array)
  (check-array 'array->nested-list array)
  (let loop ((subarray array))
    (let ((dim-length (vector-ref (%array-shape subarray) 0)))
      (if (= (%array-rank subarray) 1)
        (map (compose (cute %array-ref subarray <>) vector)
             (iota dim-length))
        (let ((lower-bound (vector-copy (%array-lower-bound subarray)))
              (upper-bound (vector-copy (%array-shape subarray))))
          (map (lambda (i)
                 (vector-set! lower-bound 0 i)
                 (vector-set! upper-bound 0 (add1 i))
                 (loop
                  (array-squeeze (array-slice subarray
                                              lower-bound
                                              upper-bound)
                                 (vector 0))))
           (iota dim-length)))))))


;; Returns a newly allocated array of the storage-class with rank rank whose
;; elements are initialized from the vector nested-vector or list nested-list.
;; It is an error if this argument is not rectangular. As a special case, if
;; rank is 0, the sole element is nested-vector or nested-list, which in that
;; case need not be a Scheme vector/list.
(define (nested-vector->array nested-vector storage-class rank)
  (check-vector 'nested-vector->array nested-vector)
  (check-positive-fixnum 'nested-vector->array rank)
  (let loop ((vec nested-vector)
             (shape '())
             (i 0))
    (cond
      ((fx>= i rank)
       (array-tabulate storage-class
                       (lambda (index)
                         (vector-fold vector-ref nested-vector index))
                       (list->vector (reverse shape))
                       #t))
      (else
        (loop (vector-ref vec 0)
              (cons (vector-length vec) shape)
              (add1 i))))))


(define-inline (flatten-by-rank rank lst)
  (let loop ((i 1)
             (xs lst))
    (if (< i rank)
      (loop (add1 i) (concatenate xs))
      xs)))

(define (nested-list->array nested-list storage-class rank)
  (check-list 'nested-list->array nested-list)
  (check-positive-fixnum 'nested-list->array rank)
  (let loop ((lst nested-list)
             (shape '())
             (i 0))
    (cond
      ((fx>= i rank)
       (let* ((flat-list (flatten-by-rank rank nested-list))
              (vshape (list->vector (reverse shape)))
              (storage-object
               (make-storage-object storage-class
                                    (length flat-list))))
         (unless (= (storage-object-length storage-class
                                           storage-object)
                    (vector-fold fx* 1 vshape))
           (error "(nested-list->array): the lengths along each dimension in nested-list must be the same."))
         (for-each (lambda (elem i)
                     ((storage-class-mutator storage-class)
                      storage-object
                      i
                      elem))
                   flat-list
                   (iota (storage-object-length storage-class
                                                storage-object)))
         (%make-array storage-class
                      vshape
                      (shape->stride vshape)
                      0
                      #t
                      storage-object)))
      ((and (fx< (add1 i) rank)
            (not (pair? (car lst))))
       (error "Rank is inconsistent with array shape."
              rank
              (reverse (cons (length lst) shape))))
      (else
        (loop (car lst)
              (cons (length lst) shape)
              (add1 i))))))


;;; Inner and outer products

;; These procedures return arrays which do not share their storage objects with
;; any existing arrays.

;; Returns a newly allocated array using the storage-class whose bounds are
;; equal to the bounds of array1 without their last elements, concatenated with
;; the bounds of array2 without their first elements. It is an error if the
;; omitted upper bounds are not numerically equal; it is also an error if the
;; omitted lower bounds are not numerically equal. It is an error if both
;; arrays have rank 0.
;;
;; Each element of the result array results from applying proc2 to the
;; corresponding elements of the last vectors of array1 and the first vectors
;; of array2 and then reducing them with proc1 to a single value. The order and
;; number of invocations of the procedures is unspecified.
;;
;; In particular, if both arrays have rank 1, the last and first vectors are
;; the whole of the arrays, and the result has rank 0; if both arrays have rank
;; 2, the last vectors of array1 are the column-wise vectors, and the first
;; vectors of array2 are the row-wise vectors, and the result has rank 2. (APL
;; inner product.)
;;
;; Example: (array-inner-product vector-storage-class + * array1 array2), where
;; the arrays have rank 1, computes the usual dot product of two vectors.
(define (array-inner-product storage-class proc1 proc2 array1 array2)
  (void))

;; Returns a newly allocated array using the storage-class whose bounds are the
;; concatenation of the bounds of array1 and array2. Each element of the new
;; array is the result of applying proc to every element of array1 and every
;; element of array2. The order and number of invocations of proc is
;; unspecified. (APL outer product.)
(define (array-outer-product storage-class proc array1 array2)
  (void))