Commit 38ddfa1a authored by Liam Healy's avatar Liam Healy

GSL arrays created by GSL functions

Made the generic function #'contents-from-pointer that takes a pointer
to a GSL struct and produces a list of the contents suitable for the
:initial-contents argument to #'make-marray.  Added an argument
:from-pointer to #'make-marray that then creates a new marray with the
same contents.  Some functions in solve-minimize-fit return a pointer
to a GSL vector with double-floats.  Copying in this way is necessary
because GSL has malloced the data in a non-CL place.  Sometimes, this
results in the copying of data whose only further use will be to
extract a pointer to pass to another GSL function.  Added a method
to #'mpointer that permits directly passing a pointer, but this isn't
useful for arrays because there are many element types so the
functions have been defined as methods.
parent 74d8995a
;; Get/set array or elements: cl-array, maref
;; Liam Healy 2008-08-27 22:43:10EDT maref.lisp
;; Time-stamp: <2008-12-28 16:32:54EST maref.lisp>
;; Time-stamp: <2008-12-28 18:09:19EST maref.lisp>
;; $Id: $
(in-package :gsl)
......@@ -41,32 +41,6 @@
(declare (ignore array-rank element-type))
object))
;;; Some functions in solve-minimize-fit return a pointer to a GSL
;;; vector with double-floats. #'cl-array will turn 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.
(defmethod cl-array ((pointer #.+foreign-pointer-class+)
&optional (array-rank 1) (element-type 'double-float))
(if (= array-rank 2)
;; Matrix
(let* ((dim1 (cffi:foreign-slot-value pointer 'gsl-matrix-c 'size1))
(dim2 (cffi:foreign-slot-value pointer 'gsl-matrix-c 'size2))
(array (make-marray (list dim1 dim2) element-type)))
;; Copy over from the C side
(loop for i below dim1
do (loop for j below dim2 do
(setf (aref array i j) (maref pointer i j))))
array)
;; Vector
(let* ((size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))
(array (make-marray size element-type)))
;; Copy over from the C side
(loop for i below size
do (setf (aref array i) (maref pointer i)))
array)))
;;;;****************************************************************************
;;;; Get or set elements of the array: maref, (setf maref)
;;;;****************************************************************************
......
;; A "marray" is an array in both GSL and CL
;; Liam Healy 2008-04-06 21:23:41EDT
;; Time-stamp: <2008-12-28 16:43:44EST marray.lisp>
;; Time-stamp: <2008-12-28 18:54:20EST marray.lisp>
;; $Id$
(in-package :gsl)
......@@ -92,18 +92,34 @@
(export 'make-marray)
(defun make-marray
(element-type &rest keys &key dimensions initial-contents &allow-other-keys)
(element-type &rest keys &key dimensions initial-contents from-pointer
&allow-other-keys)
"Make a GSLL array with the given element type,
:dimensions, :initial-contents and/or :initial-element."
(apply #'make-instance
(data-class-name
(if
(or
(and dimensions (listp dimensions) (eql (length dimensions) 2))
(and initial-contents (listp (first initial-contents))))
'matrix 'vector)
element-type)
keys))
:dimensions, :initial-contents and/or :initial-element.
If a pointer to a GSL object is given in :from-pointer, create
an object with duplicate contents; if a matrix, :dimensions must be set to 2."
;; Some functions in solve-minimize-fit return a pointer to a GSL
;; vector of double-floats. With the :from-pointer argument, this
;; function turn 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.
(if from-pointer
(make-marray element-type
:initial-contents
(contents-from-pointer
from-pointer
(if (eql dimensions 2) 'gsl-matrix-c 'gsl-vector-c)
element-type))
(apply #'make-instance
(data-class-name
(if
(or
(and dimensions (listp dimensions) (eql (length dimensions) 2))
(and initial-contents (listp (first initial-contents))))
'matrix 'vector)
element-type)
keys)))
(defun hashm-numeric-code (n)
"Get the appropriate element type for the numeric code n"
......
;; Matrices
;; Liam Healy 2008-04-15 21:57:52EDT matrix.lisp
;; Time-stamp: <2008-12-28 16:55:04EST matrix.lisp>
;; Time-stamp: <2008-12-28 18:01:25EST matrix.lisp>
;; $Id$
(in-package :gsl)
......@@ -25,6 +25,16 @@
;;; Define all supported matrix subclasses
#.(data-defclass 'matrix 'matrix)
(defmethod contents-from-pointer
(pointer (struct-type (eql 'gsl-matrix-c))
&optional (element-type 'double-float))
(let ((dim1 (cffi:foreign-slot-value pointer struct-type 'size1))
(dim2 (cffi:foreign-slot-value pointer struct-type 'size2)))
;; Copy over from the C side
(loop for i below dim1
collect (loop for j below dim2
collect (maref pointer i j element-type)))))
;;;;****************************************************************************
;;;; Mathematical
;;;;****************************************************************************
......
;; Vectors
;; Liam Healy 2008-04-13 09:39:02EDT vector.lisp
;; Time-stamp: <2008-12-28 16:55:05EST vector.lisp>
;; Time-stamp: <2008-12-28 18:01:59EST vector.lisp>
;; $Id$
(in-package :gsl)
......@@ -24,6 +24,12 @@
;;; Define all supported mvector subclasses
#.(data-defclass 'vector 'mvector)
(defmethod contents-from-pointer
(pointer (struct-type (eql 'gsl-vector-c))
&optional (element-type 'double-float))
(loop for i below (cffi:foreign-slot-value pointer struct-type 'size)
collect (maref pointer i nil element-type)))
;;;;****************************************************************************
;;;; Function definitions
;;;;****************************************************************************
......
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2008-12-26 19:16:05EST mobject.lisp>
;; Time-stamp: <2008-12-28 18:25:33EST mobject.lisp>
;; $Id$
;;; GSL objects are represented in GSLL as and instance of a 'mobject.
......@@ -139,3 +139,17 @@
(defconstant +foreign-pointer-class+ (class-name (class-of (cffi:null-pointer)))
"The class in which foreign pointers fall. This will be assumed to be a
GSL vector or matrix.")
;;; Some functions in solve-minimize-fit return a pointer to a GSL
;;; vector with double-floats. This function will return a contents
;;; form suitable for make-marray. 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.
(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 mpointer ((object #.+foreign-pointer-class+))
object)
;; Multivariate minimization.
;; Liam Healy <Tue Jan 8 2008 - 21:28>
;; Time-stamp: <2008-12-26 18:33:28EST minimization-multi.lisp>
;; Time-stamp: <2008-12-28 18:14:04EST minimization-multi.lisp>
;; $Id$
(in-package :gsl)
......@@ -122,16 +122,16 @@
(defmfun mfminimizer-x (minimizer)
"gsl_multimin_fminimizer_x"
(((mpointer minimizer) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The current best estimate of the location of the minimum.")
(defmfun mfdfminimizer-x (minimizer)
"gsl_multimin_fdfminimizer_x"
(((mpointer minimizer) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The current best estimate of the location of the minimum.")
......@@ -159,8 +159,8 @@
(defmfun mfdfminimizer-gradient (minimizer)
"gsl_multimin_fdfminimizer_gradient"
(((mpointer minimizer) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The current best estimate of the gradient for the minimizer.")
......@@ -177,7 +177,7 @@
(defmfun min-test-gradient (gradient absolute-error)
"gsl_multimin_test_gradient"
((gradient :pointer) (absolute-error :double))
(((mpointer gradient) :pointer) (absolute-error :double))
:c-return :success-continue
:documentation ; FDL
"Test the norm of the gradient against the
......
;; Nonlinear least squares fitting.
;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
;; Time-stamp: <2008-12-26 18:29:44EST nonlinear-least-squares.lisp>
;; Time-stamp: <2008-12-28 18:49:26EST nonlinear-least-squares.lisp>
;; $Id$
(in-package :gsl)
......@@ -129,16 +129,16 @@
(defmfun current-parameters-ffit (solver)
"gsl_multifit_fsolver_position"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The current best-fit parameters.")
(defmfun current-parameters-fdffit (solver)
"gsl_multifit_fdfsolver_position"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The current best-fit parameters.")
......@@ -295,13 +295,13 @@
(make-exponent-fit-data
:n *number-of-observations*
:y
(let ((arr (make-marray *number-of-observations* 'double-float)))
(let ((arr (make-marray 'double-float :dimensions *number-of-observations*)))
(let ((rng (make-random-number-generator *mt19937* 0)))
(dotimes (i *number-of-observations* arr)
(setf (aref arr i)
(setf (maref arr i)
(+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0))))))
:sigma
(make-marray *number-of-observations* 'double-float :initial-element 0.1d0))))
(make-marray 'double-float :dimensions *number-of-observations* :initial-element 0.1d0))))
(defun exponential-residual (x f)
"Compute the negative of the residuals with the exponential model
......@@ -310,13 +310,13 @@
(lambda (maref x 1))
(b (maref x 2)))
(symbol-macrolet
((y (exponent-fit-data-y *nlls-example-data*))
(sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i *number-of-observations*)
(setf (maref f i)
;; the difference model - observation = - residual
(/ (- (+ (* A (exp (* (- lambda) i))) b) (aref y i))
(aref sigma i)))))))
((y (exponent-fit-data-y *nlls-example-data*))
(sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i *number-of-observations*)
(setf (maref f i)
;; the difference model - observation = - residual
(/ (- (+ (* A (exp (* (- lambda) i))) b) (maref y i))
(maref sigma i)))))))
(defun exponential-residual-derivative (x jacobian)
"Compute the partial derivatives of the negative of the
......@@ -328,7 +328,7 @@
((sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i *number-of-observations*)
(let ((e (exp (* (- lambda) i)))
(s (aref sigma i)))
(s (maref sigma i)))
(setf (maref jacobian i 0) (/ e s)
(maref jacobian i 1) (* -1 i A (/ e s))
(maref jacobian i 2) (/ s)))))))
......@@ -346,8 +346,7 @@
(defun norm-f (fit)
"Find the norm of the fit function f."
;; Fix this
(euclidean-norm (fdffit-slot fit 'f)))
(euclidean-norm (make-marray 'double-float :from-pointer (fdffit-slot fit 'f))))
(defun solve-nonlinear-least-squares-example ()
(let* ((init #m(1.0d0 0.0d0 0.0d0))
......@@ -355,12 +354,12 @@
(make-marray 'double-float
:dimensions
(list *number-of-parameters* *number-of-parameters*)))
(fit (make-nonlinear-fdffit
*levenberg-marquardt*
*number-of-observations*
*number-of-parameters*
exponential-residual
init)))
(fit (mpointer (make-nonlinear-fdffit
*levenberg-marquardt*
*number-of-observations*
*number-of-parameters*
exponential-residual
init))))
(macrolet ((fitx (i) `(maref (fdffit-slot fit 'x) ,i))
(err (i) `(sqrt (maref covariance ,i ,i))))
(format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g"
......@@ -391,3 +390,4 @@
;;; Run example:
;;; (nlls-setup)
;;; (solve-nonlinear-least-squares-example)
;;; (5.045357801443204d0 0.10404905892045835d0 1.0192487061031013d0)
;;; Multivariate roots.
;;; Liam Healy 2008-01-12 12:49:08
;;; Time-stamp: <2008-12-26 18:37:15EST roots-multi.lisp>
;;; Time-stamp: <2008-12-28 18:25:31EST roots-multi.lisp>
;;; $Id$
(in-package :gsl)
......@@ -118,48 +118,48 @@
(defmfun mfsolver-root (solver)
"gsl_multiroot_fsolver_root"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The current estimate of the root for the solver.")
(defmfun mfdfsolver-root (solver)
"gsl_multiroot_fdfsolver_root"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation
"The current estimate of the root for the solver.")
(defmfun mfsolver-f (solver)
"gsl_multiroot_fsolver_f"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The function value f(x) at the current estimate x of the root for the solver.")
(defmfun mfdfsolver-f (solver)
"gsl_multiroot_fdfsolver_f"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The function value f(x) at the current estimate x of the root for the solver.")
(defmfun mfsolver-dx (solver)
"gsl_multiroot_fsolver_dx"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The last step dx taken by the solver.")
(defmfun mfdfsolver-dx (solver)
"gsl_multiroot_fsolver_dx"
(((mpointer solver) :pointer))
:c-return :pointer
:return (:c-return)
:c-return (crtn :pointer)
:return ((make-marray 'double-float :from-pointer crtn))
:documentation ; FDL
"The last step dx taken by the solver.")
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment