Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
226 lines (189 sloc) 9.17 KB
;;; -*- Mode:Lisp; Syntax:ANSI-Common-Lisp; Coding:utf-8 -*-
(in-package #:array-operations)
;;; coercing can be used with * forms
(defun coercing (element-type &optional (function #'identity))
"Return a function composed of a univariate function that coerces to ELEMENT-TYPE and function. When FUNCTION is not given, return a closure that coerces to ELEMENT-TYPE."
(compose (lambda (value) (coerce value element-type))
;;; creating arrays
(defun generate* (element-type function dimensions &optional arguments)
"Return an array with given DIMENSIONS and ELEMENT-TYPE, with elements
generated by calling FUNCTION with
- no arguments, when ARGUMENTS is nil
- the position (= row major index), when ARGUMENTS is :POSITION
- a list of subscripts, when ARGUMENTS is :SUBSCRIPTS
The traversal order is unspecified and may be nonlinear."
(let ((dimensions (ensure-dimensions dimensions)))
(aprog1 (make-array dimensions :element-type element-type)
(ecase arguments
(dotimes (position (array-total-size it))
(setf (row-major-aref it position)
(funcall function))))
(walk-subscripts (dimensions subscripts position)
(setf (row-major-aref it position) (funcall function position))))
(walk-subscripts-list (dimensions subscripts position)
(setf (row-major-aref it position)
(funcall function subscripts))))
(walk-subscripts-list (dimensions subscripts position)
(setf (row-major-aref it position)
(funcall function position subscripts))))))))
(defun generate (function dimensions &optional arguments)
(generate* t function dimensions arguments))
;;; permutations
;;; A permutation is a list of nonnegative, non-repeated integers, below some
;;; rank (of the array it is applied to).
(define-condition permutation-repeated-index (error)
((index :initarg :index)))
(define-condition permutation-invalid-index (error)
((index :initarg :index)))
(define-condition permutation-incompatible-rank (error)
(defun permutation-flags (permutation &optional (rank (length permutation)))
"Make a bit vector of flags with indexes from PERMUTATION, signalling errors
for invalid and repeated indices. NOT EXPORTED."
(aprog1 (make-array rank :element-type 'bit :initial-element 0)
(map nil (lambda (p)
(assert (and (integerp p) (< -1 p rank)) ()
'permutation-invalid-index :index p)
(assert (zerop (aref it p)) ()
'permutation-repeated-index :index p)
(setf (aref it p) 1))
(defun check-permutation (permutation
&optional (rank (length permutation) rank?))
"Check if PERMUTATION is a valid permutation (of the given RANK), and signal
an error if necessary."
(when rank?
(assert (= rank (length permutation)) ()
'permutation-incompatible-rank ))
(assert (every #'plusp (permutation-flags permutation)) ()
(defun complement-permutation (permutation rank)
"Return a list of increasing indices that complement PERMUTATION, ie form a
permutation when appended. Atoms are accepted and treated as lists of a
single element."
(loop for f across (permutation-flags (ensure-list permutation) rank)
for index from 0
when (zerop f)
collect index))
(defun complete-permutation (permutation rank)
"Return a completed version of permutation, appending it to its complement."
(let ((permutation (ensure-list permutation)))
(append permutation (complement-permutation permutation rank))))
(defun invert-permutation (permutation)
"Invert a permutation."
(check-permutation permutation)
(coerce (aprog1 (make-array (length permutation) :element-type 'fixnum)
(map nil (let ((index 0))
(lambda (p)
(setf (aref it p) index)
(incf index)))
(defun identity-permutation? (permutation
&optional (rank (length permutation)))
"Test if PERMUTATION is the identity permutation, ie a sequence of
consecutive integers starting at 0. Note that permutation is otherwise not
checked, ie it may not be a permutation."
(let ((index 0))
(lambda (p)
(prog1 (= index p)
(incf index)))
(= index rank))))
(defun permute (permutation array)
"Return ARRAY with the axes permuted by PERMUTATION, which is a sequence of
indexes. Specifically, an array A is transformed to B, where
B[b_1,...,b_n] = A[a_1,...,a_n] with b_i=a_{P[i]}
P is the permutation.
Array element type is preserved."
(let ((rank (array-rank array)))
(if (identity-permutation? permutation rank)
(let+ ((dimensions (array-dimensions array))
((&flet map-subscripts (subscripts-vector)
(map 'list (curry #'aref subscripts-vector) permutation))))
(check-permutation permutation rank)
(aprog1 (make-array (map-subscripts (coerce dimensions 'vector))
:element-type (array-element-type array))
(walk-subscripts (dimensions subscripts position)
(setf (apply #'aref it (map-subscripts subscripts))
(row-major-aref array position))))))))
;;; margin
(defun each* (element-type function array &rest other-arrays)
"Apply function to the array arguments elementwise, and return the result as
an array with the given ELEMENT-TYPE. Arguments are checked for dimension
(aprog1 (make-array (array-dimensions array) :element-type element-type)
(assert (apply #'same-dimensions? array other-arrays))
(apply #'map-into (flatten it) function
(flatten array) (mapcar #'flatten other-arrays))))
(defun each (function array &rest other-arrays)
"Like EACH*, with ELEMENT-TYPE T."
(apply #'each* t function array other-arrays))
(defun margin* (element-type function array inner
&optional (outer (complement-permutation inner
(array-rank array))))
"PERMUTE ARRAY with `(,@OUTER ,@INNER), split the inner subarrays, apply
FUNCTION to each, return the results in an array of dimensions OUTER, with the
(let ((outer (ensure-list outer)))
(each* element-type function
(split (permute (append outer (ensure-list inner)) array)
(length outer)))))
(defun margin (function array inner
&optional (outer (complement-permutation inner
(array-rank array))))
(margin* t function array inner outer))
;;; recycle
(defun recycle (object &key inner outer
(element-type (if (arrayp object)
(array-element-type object)
"Recycle elements of object, extending the dimensions by outer (repeating
OBJECT) and inner (repeating each element of OBJECT). When both INNER and
OUTER are nil, the OBJECT is returned as is. Non-array objects are intepreted
as rank 0 arrays, following the usual semantics."
(if (or inner outer)
(let ((inner (ensure-dimensions inner))
(outer (ensure-dimensions outer)))
(if (arrayp object)
(let ((dimensions (array-dimensions object)))
(aprog1 (make-array (append outer dimensions inner)
:element-type element-type)
(let* ((outer-size (product outer))
(size (product dimensions))
(inner-size (product inner))
(reshaped (reshape it (list outer-size size inner-size))))
(loop for outer-index below outer-size
do (loop for index below size
do (fill (sub reshaped outer-index index)
(row-major-aref object index)))))))
(make-array (append outer inner) :initial-element object
:element-type element-type)))
;;; outer produce
(defun outer* (element-type function &rest arrays)
"Generalized outer product of ARRAYS with FUNCTION. The resulting array has the concatenated dimensions of ARRAYS, and the given ELEMENT-TYPE."
(assert arrays)
(let* ((result (make-array (mapcan #'dims arrays) :element-type element-type))
(vectors (mapcar #'flatten arrays))
(flat-dimensions (mapcar #'length vectors))
(flat-result (reshape result flat-dimensions)))
(walk-subscripts (flat-dimensions subscripts position)
(setf (row-major-aref flat-result position)
(apply function (map 'list #'aref vectors subscripts))))
(defun outer (function &rest arrays)
"Like OUTER, with ELEMENT-TYPE t."
(apply #'outer* t function arrays))