guile-ploy/ploy/as-array.scm

; (c) Daniel Llorens - 2012-2013
; Array memory order operations and array conversion.

; This library is free software; you can redistribute it and/or modify it under
; the terms of the GNU General Public License as published by the Free
; Software Foundation; either version 3 of the License, or (at your option) any
; later version.

; The order options are for interfacing with foreign libraries.
; In normal Guile array code one shouldn't be concerned about storage order.

(define-module (ploy as-array))
(import (ice-9 optargs) (srfi srfi-26) (srfi srfi-1) (srfi srfi-11)
        (srfi srfi-9) (srfi srfi-8) (ploy basic) (ploy assert) (ice-9 control))

(define (c-order? a)
  (reset
   (fold (lambda (i s c)
           (if (not (= (* c s) i))
             (shift k #f)
             (* c s)))
         1
         (reverse (shared-array-increments a))
         (cons 1 (reverse (array-dimensions a))))))

(define (fortran-order? a)
  (reset
   (fold (lambda (i s c)
           (if (not (= (* c s) i))
             (shift k #f)
             (* c s)))
         1
         (shared-array-increments a)
         (cons 1 (array-dimensions a)))))

(define (assignable? a)
  (reset
   (case (array-rank a)
     ((0) #t)
     ((1) (or (< (tally a) 2) (not (zero? (first (shared-array-increments a))))))
     (else
; ignore singleton axes, bail out on null axes.
      (let ((is (sort (filter (lambda (is)
                                (case (second is)
                                  ((0) (shift k #f))
                                  ((1) #f)
                                  (else #t)))
                              (zip (shared-array-increments a) (array-dimensions a)))
                      (lambda (a b) (< (magnitude (first a)) (magnitude (first b)))))))
        (let loop ((is is))
          (cond ((null? is) #t)
                ((zero? (first (car is))) #f)
                ((null? (cdr is)) #t)
                ((> (* (magnitude (first (car is))) (second (car is)))
                    (magnitude (first (cadr is))))
                 #f)
                (else (loop (cdr is))))))))))

(export c-order? fortran-order? assignable?)

; Convert possibly nested arraylike object to array.

; @todo (as-array (reshape 1 10) #:type 'f64) creates a full size 10 array. Is this what we want?
(define* (arraylike-dimensions A #:key rank)
  (let loop ((A A) (i 0) (dims '()) (nested-list? #t))
    (cond ((and rank (= i rank)) (values dims nested-list?))
          ((and rank (> i rank)) (throw 'as-array-dimensions-cannot-be-split-by-rank rank))
          ((array? A) (let ((dims (append dims (array-dimensions A)))
                            (root (shared-array-root A)))
                        (if (zero? (tally root))
                          (values dims #f)
                          (loop (array-ref root 0) (+ i (array-rank A)) dims #f))))
          ((list? A) (if (null? A)
                       (values (append dims '(0)) nested-list?)
                       (loop (car A) (+ i 1) (append dims (list (length A))) nested-list?)))
          ((and rank (< i rank)) (throw 'as-array-requested-rank-too-large rank i))
          (else (values dims nested-list?)))))

(define* (as-array A #:key type order rank unique? check?)
  (let*-values (((dims nested-list?) (arraylike-dimensions A #:rank rank))
                ((rank) (or rank (length dims))))
    (define (make-dest-array type)
      (case order
        ((c #f) (apply make-typed-array type *unspecified* dims))
        ((fortran) (apply transpose-array (apply make-typed-array type *unspecified* (reverse dims))
                          (iota rank (- rank 1) -1)))
        (else (throw 'arbitrary-order-not-implemented))))
    (cond
; never ever return #0().
     ((zero? rank) A)
; special case. @todo for order != c, list->array and then array-copy would still be faster.
     ((and nested-list? (or (eq? order 'c) (eq? order #f)))
      (list->typed-array (or type #t) rank A))
     ((and (array? A) (= (array-rank A) rank))
      (assert (not (zero? rank)) "BAD")
      (if (and (case order
                 ((#f) #t)
                 ((fortran) (fortran-order? A))
                 ((c) (c-order? A))
                 (else (throw 'arbitrary-order-not-implemented)))
               (or (not type) (eq? (array-type A) type))
               (not unique?))
        A
        (or check?
            (let ((B (make-dest-array (or type (array-type A)))))
              (array-copy! A B)
              B))))
; need to delve and convert.
     (else
      (and (not check?)
           (let ((B (make-dest-array (or type #t))))
             (let loopd ((A A) (B B) (rank rank))
               (if (list? A)
                 (let ((rank- (- rank 1)))
                   (if (zero? rank-)
                     (array-copy! (list->array 1 A) B)
                     (let loop ((A A) (i 0))
                       (unless (null? A)
                         (loopd (car A) (array-cell-ref B i) rank-)
                         (loop (cdr A) (+ i 1))))))
                 (let ((rank- (- rank (array-rank A))))
                   (if (zero? rank-)
                     (array-copy! A B)
                     (array-slice-for-each (array-rank A) (lambda (A B) (loopd (array-cell-ref A) (array-cell-ref B) rank-)) A B)))))
             B))))))

(export arraylike-dimensions as-array)