;; A "marray" is an array in both GSL and CL ;; Liam Healy 2008-04-06 21:23:41EDT ;; Time-stamp: <2010-01-20 15:31:57EST marray.lisp> ;; ;; Copyright 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License ;; ;; This program 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. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public License ;; along with this program. If not, see . (in-package :gsl) ;;;;**************************************************************************** ;;;; The class marray and its construction ;;;;**************************************************************************** (export 'marray) (defclass marray (mobject c-array:foreign-array) ((block-pointer :initarg :block-pointer :reader block-pointer) (total-size :reader size)) (:documentation "A superclass for arrays represented in GSL and CL.")) ;;; We don't allocate or free the C array data, because that is ;;; handled by c-array:foreign-array. We can use the GSL functions ;;; *_alloc_from_block because they will allocate only the structure, ;;; but we must use CFFI to allocate the block structure; otherwise ;;; GSL would try to allocate the C array. We do not use any of the ;;; GSL free functions because they would free the C array data. (defmethod initialize-instance :after ((object marray) &rest initargs) (declare (ignore initargs)) ; required by &rest arg for c-array:foreign-array? ;; Need to check that all allocations succeeded. ;; Don't do anything if mpointer has been assigned (shouldn't happen) ;; or this is a permutation or combination (they have their own methods). (unless (and (slot-boundp object 'mpointer) (mpointer object)) (let ((blockptr (cffi:foreign-alloc 'gsl-block-c))) (setf (cffi:foreign-slot-value blockptr 'gsl-block-c 'size) (size object) (slot-value object 'block-pointer) blockptr) (let ((array-struct (alloc-from-block object blockptr))) (setf (slot-value object 'mpointer) array-struct) #-native (set-struct-array-pointer object) (tg:finalize object (lambda () (cffi:foreign-free blockptr) (cffi:foreign-free array-struct))))))) (defun set-struct-array-pointer (object) "Set the pointer in the 'data slot of the foreign structs to be the current c-pointer value. In non-native implementations this need be called only once when the marray is made. In native implementations it is called whenever mpointer is requested because of the possibility that a GC moved the pointer." (setf (cffi:foreign-slot-value (block-pointer object) 'gsl-block-c 'data) (c-pointer object) ;; alloc-from-block automatically copies over the data pointer ;; from the block to the vector/matrix; we must do that manually here (cffi:foreign-slot-value (slot-value object 'mpointer) (if (typep object 'matrix) 'gsl-matrix-c 'gsl-vector-c) 'data) (c-pointer object))) #+native (defmethod mpointer ((object marray)) "Compute the c-pointer of the array and place it in the GSL struct because the stored version is untrustworthy unless it was computed within the same native-pointer-protect form as this mpointer extraction." (set-struct-array-pointer object) (call-next-method)) (defmethod make-load-form ((object marray) &optional env) (declare (ignore env)) `(make-marray ',(element-type object) :dimensions ',(dimensions object) :initial-contents ',(loop for elt across (subseq (c-array:original-array object) (c-array:offset object) (+ (c-array:offset object) (size object))) collect elt))) ;;;;**************************************************************************** ;;;; Definition of specific data classes ;;;;**************************************************************************** (defparameter *class-element-type* nil "The mapping between the class name and the CL element type.") (defun data-class-name (category element-type) "The class name from the type of element." ;; e.g. (data-class-name 'vector '(unsigned-byte 8)) ;; -> VECTOR-UNSIGNED-BYTE-8 (if (member category '(vector matrix)) (intern (format nil "~a-~a" category (c-array:cl-single element-type :gsl)) :gsl) category)) (defun data-defclass (category superclass) "Define all the subclasses based on the known element types." (cons 'progn (mapcan (lambda (element-type-cl) (let* ((class-name (data-class-name category element-type-cl))) ;; Define the class `((defclass ,class-name (,superclass) ((element-type :initform ',element-type-cl :allocation :class))) ;; Push mapping onto *class-element-type* (pushnew ',(cons class-name element-type-cl) *class-element-type* :test #'equal) (export ',class-name)))) c-array:*array-element-types*))) ;;;;**************************************************************************** ;;;; Make data from either the dimensions provided or from the initial values ;;;;**************************************************************************** (export 'make-marray) (defun make-marray (class-or-element-type &rest keys &key dimensions initial-contents cl-array &allow-other-keys) "Make a GSLL array with the given element type, :dimensions, and :initial-contents, :initial-element or :cl-array. If the :cl-array is supplied, it should be a CL array generated with #'make-ffa." (apply #'make-instance (if (subtypep class-or-element-type 'marray) class-or-element-type (data-class-name (if (or (and dimensions (listp dimensions) (eql (length dimensions) 2)) (and initial-contents (listp (first initial-contents))) (and cl-array (eql (length (array-dimensions cl-array)) 2))) 'matrix 'vector) class-or-element-type)) keys)) (eval-when (:compile-toplevel :load-toplevel :execute) (pushnew 'marray grid:*grid-data-superclasses*)) (defmethod grid:make-grid-data ((type (eql 'marray)) dimensions rest-spec &key initial-element) (make-marray (grid:spec-scalar-p rest-spec) :dimensions dimensions :initial-element initial-element)) ;;; The reader macro #m will read a list of arguments, evaluating the ;;; contents, and construct a marray from it. If the symbol ^ occurs ;;; in the list, it indicates the object being made is a matrix, and ;;; this symbol indicates the end of each row. For compatibility ;;; reasons, a list of numeric values may be given for each row also. ;;; Examples: ;;; #m(1 2 (exp 1)) ;;; # ;;; #m(1 2 3 ^ pi (cos (/ pi 4)) -12) ;;; # ;;; #2m(2 1 3 -1) ;;; # (defun hashm-numeric-code (n) "Get the appropriate element type for the numeric code n" (case n ((nil 1) 'double-float) (2 '(complex double-float)) (3 'single-float) (4 '(complex single-float)) (7 '(signed-byte 8)) (8 '(unsigned-byte 8)) (15 '(signed-byte 16)) (16 '(unsigned-byte 16)) (31 '(signed-byte 32)) (32 '(unsigned-byte 32)) (63 '(signed-byte 64)) (64 '(unsigned-byte 64)))) (defun quote-numeric-list (s) (if (and (listp s) (numberp (first s))) `',s s)) (defvar *row-separator* '^) (export *row-separator*) (defun conslist (l) (cons 'list l)) (set-dispatch-macro-character #\# #\m (lambda (stream subchar arg) (declare (ignore subchar)) (read-char stream) (let ((list (read-delimited-list #\) stream))) `(make-marray ',(hashm-numeric-code arg) :initial-contents ,(if (find *row-separator* list) (conslist (mapcar 'conslist (cl-utilities:split-sequence *row-separator* list))) `(list ,@(mapcar 'quote-numeric-list list))))))) ;;;;**************************************************************************** ;;;; Copy to and from bare mpointers ;;;;**************************************************************************** (defgeneric contents-from-pointer (pointer struct-type &optional element-type) (:documentation "Create a contents list from the GSL object of type struct-type referenced by pointer.")) (defmethod c-array:copy-making-destination ((pointer #.+foreign-pointer-class+)) (foreign-pointer-method pointer ;; Default assumption when destination isn't given in #'copy is ;; that this should make a vector-double-float. (c-array:copy-to-destination pointer 'vector-double-float))) (defmethod c-array:copy-to-destination ((pointer #.+foreign-pointer-class+) (class-name symbol)) (foreign-pointer-method pointer (make-marray class-name :initial-contents (contents-from-pointer pointer (if (subtypep class-name 'matrix) 'gsl-matrix-c 'gsl-vector-c) (c-array:lookup-type class-name *class-element-type*))))) ;; Some functions in solve-minimize-fit return a pointer to a GSL ;; vector of double-floats. This function turns that into a ;; foreign-friendly array. There is no choice but to copy over the ;; data even on native implementations; because GSL is doing the ;; mallocing, the data are not CL-accessible. ;;;;**************************************************************************** ;;;; Convenient defaulting of defmfun arguments ;;;;**************************************************************************** (defun make-marray-or-default (default dimensions &optional dont-make-if-nil (element-type 'double-float) initial-element) "If default is T or dont-make-if-nil, make a marray and returned with the dimensions and element-type. If default is a marray, it is returned after being synced for non-native." (if (member default (list t dont-make-if-nil)) (if initial-element (make-marray element-type :dimensions dimensions :initial-element initial-element) (make-marray element-type :dimensions dimensions)) (progn #-native (when (typep default 'marray) (c-array:copy-cl-to-c default)) default)))