Commit 7c66d57b authored by Liam Healy's avatar Liam Healy

New class foreign-array, new names: marray, make-marray, mobject.lisp

New class foreign-array represents foreign (C) arrays, without GSL or
mobject superclass.  New class name 'marray was 'gsl-data,
representing arrays in GSL.  These objects are made by #'make-marray,
formerely #'make-array*.  New file names: init/mobject.lisp was
init/gsl-objects.lisp and data/marray.lisp was data/data.lisp.
parent a9dda7c4
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2008-12-21 21:46:10EST basis-splines.lisp>
;; Time-stamp: <2008-12-26 10:25:46EST basis-splines.lisp>
;; $Id$
(in-package :gsl)
......@@ -112,13 +112,13 @@
(bw (make-basis-spline order nbreak))
(mw (fit-workspace ndata ncoeffs))
(rng (random-number-generator *mt19937* 0))
(B (make-array* 'double-float :dimensions ncoeffs))
(c (make-array* 'double-float :dimensions ncoeffs))
(cov (make-array* 'double-float :dimensions (list ncoeffs ncoeffs)))
(w (make-array* 'double-float :dimensions ndata))
(x (make-array* 'double-float :dimensions ndata))
(y (make-array* 'double-float :dimensions ndata))
(Xmatrix (make-array* 'double-float :dimensions (list ndata ncoeffs)))
(B (make-marray 'double-float :dimensions ncoeffs))
(c (make-marray 'double-float :dimensions ncoeffs))
(cov (make-marray 'double-float :dimensions (list ncoeffs ncoeffs)))
(w (make-marray 'double-float :dimensions ndata))
(x (make-marray 'double-float :dimensions ndata))
(y (make-marray 'double-float :dimensions ndata))
(Xmatrix (make-marray 'double-float :dimensions (list ndata ncoeffs)))
(sigma 0.1d0))
;; The data to be fitted.
(dotimes (i ndata)
......
;; Use the foreign-friendly arrays package.
;; Liam Healy 2008-03-22 15:40:08EDT
;; Time-stamp: <2008-11-30 18:44:50EST foreign-friendly.lisp>
;; Time-stamp: <2008-12-26 10:10:40EST foreign-friendly.lisp>
;; $Id$
;;; Foreign-friendly arrays (original implementation by Tamas Papp)
......@@ -116,10 +116,10 @@
,body))
#+sbcl
(defun c-pointer (gsl-data)
(defun c-pointer (marray)
"The pointer to the C array."
(cffi:inc-pointer
(sb-sys:vector-sap (original-array gsl-data)) (offset gsl-data)))
(sb-sys:vector-sap (original-array marray)) (offset marray)))
;;;;****************************************************************************
;;;; Pointer management
......
;; Get/set array or elements: cl-array, maref
;; Liam Healy 2008-08-27 22:43:10EDT maref.lisp
;; Time-stamp: <2008-12-07 17:08:25EST maref.lisp>
;; Time-stamp: <2008-12-26 10:25:50EST maref.lisp>
;; $Id: $
(in-package :gsl)
......@@ -18,7 +18,7 @@
;;; Both these functions take one of the following class arguments for
;;; the array:
;;; - a gsl-data
;;; - a marray
;;; - a Common Lisp array
;;; - a pointer to a GSL vector or matrix structure
......@@ -28,11 +28,11 @@
(defgeneric cl-array (object &optional array-rank element-type)
(:documentation
"The array as a CL native array. The object may be a gsl-data object,
"The array as a CL native array. The object may be a marray object,
a pointer to a GSL vector or matrix, or an ordinary CL array of one
or two dimensions. Optional arguments array-rank and element-type
are used only for pointers.")
(:method ((object gsl-data) &optional array-rank element-type)
(:method ((object marray) &optional array-rank element-type)
(declare (ignore array-rank element-type))
#-native (copy-c-to-cl object)
(slot-value object 'cl-array))
......@@ -53,7 +53,7 @@
;; Matrix
(let* ((dim1 (cffi:foreign-slot-value pointer 'gsl-matrix-c 'size1))
(dim2 (cffi:foreign-slot-value pointer 'gsl-matrix-c 'size2))
(array (make-array* (list dim1 dim2) element-type)))
(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
......@@ -61,7 +61,7 @@
array)
;; Vector
(let* ((size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))
(array (make-array* size element-type)))
(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)))
......@@ -73,9 +73,9 @@
(defgeneric maref (object index &optional index2 type)
(:documentation
"An element of the data. The object may be a gsl-data object, a pointer to
"An element of the data. The object may be a marray object, a pointer to
a GSL vector or matrix, or an ordinary CL array of one or two dimensions.")
(:method ((object gsl-data) index &optional index2 type)
(:method ((object marray) index &optional index2 type)
(declare (ignore type))
#-native (copy-c-to-cl object)
(if index2
......@@ -93,10 +93,10 @@
(defgeneric (setf maref) (value object index &optional index2 type)
(:documentation
"Set an element of the data. The object may be a gsl-data object,
"Set an element of the data. The object may be a marray object,
a pointer to a GSL vector or matrix, or an ordinary CL array
of one or two dimensions.")
(:method (value (object gsl-data) index &optional index2 type)
(:method (value (object marray) index &optional index2 type)
(declare (ignore type))
(if index2
(setf (aref (slot-value object 'cl-array) index index2) value)
......
;; "Data" is bulk arrayed data, like vectors, matrices, permutations,
;; combinations, or histograms.
;; A "marray" is an array in both GSL and CL
;; Liam Healy 2008-04-06 21:23:41EDT
;; Time-stamp: <2008-12-25 09:40:55EST data.lisp>
;; Time-stamp: <2008-12-26 10:25:50EST marray.lisp>
;; $Id$
(in-package :gsl)
;;;;****************************************************************************
;;;; The class gsl-data and element types that make subclasses
;;;; The class marray and element types that make subclasses
;;;;****************************************************************************
(defclass gsl-data (mobject)
(defclass foreign-array ()
((cl-array :documentation "The Lisp array.")
(block-pointer :initform nil :accessor block-pointer
:documentation "A pointer to the gsl-block-c.")
#-native
(c-pointer :accessor c-pointer :documentation "A pointer to the C array.")
(dimensions :reader dimensions)
......@@ -39,14 +36,19 @@
way. If a list of index sets, those indices disagree and the remainder
are correct."))
(:documentation
"A superclass for all data storage structures, such as vector, matrix,
etc."))
"A superclass for arrays represented in C and CL."))
(defclass marray (mobject foreign-array)
((block-pointer :initform nil :accessor block-pointer
:documentation "A pointer to the gsl-block-c."))
(:documentation
"A superclass for arrays represented in GSL and CL."))
(export '(dimensions total-size element-type))
;;; Allowable keys: :dimensions, :initial-contents, :initial-element.
(defmethod initialize-instance :after
((object gsl-data) &rest initargs &key &allow-other-keys)
((object marray) &rest initargs &key &allow-other-keys)
(with-slots (cl-array dimensions original-array offset total-size) object
(let ((ffa (apply #'make-ffa (element-type object) initargs)))
(setf cl-array ffa
......@@ -61,14 +63,14 @@
(cffi:foreign-type-size (cl-cffi (element-type object)))))))
(alloc-gsl-struct object))
(defmethod print-object ((object gsl-data) stream)
(defmethod print-object ((object marray) stream)
(print-unreadable-object (object stream :type t)
#-native (copy-c-to-cl object)
(princ (cl-array object) stream)))
(defmethod make-load-form ((object gsl-data) &optional env)
(defmethod make-load-form ((object marray) &optional env)
(declare (ignore env))
`(make-array*
`(make-marray
',(element-type object)
:initial-contents
',(loop for elt across
......@@ -125,8 +127,8 @@
;;;; Make data from either the dimensions provided or from the initial values
;;;;****************************************************************************
(export 'make-array*)
(defun make-array*
(export 'make-marray)
(defun make-marray
(element-type &rest keys &key dimensions initial-contents &allow-other-keys)
"Make a GSLL array with the given element type,
:dimensions, :initial-contents and/or :initial-element."
......@@ -162,7 +164,7 @@
(declare (ignore subchar))
(read-char stream)
(let ((list (read-delimited-list #\) stream)))
`(make-array* ',(hashm-numeric-code arg) :initial-contents ',list))))
`(make-marray ',(hashm-numeric-code arg) :initial-contents ',list))))
(defun component-type (eltype)
(cl-cffi (component-float-type eltype)))
......@@ -264,7 +266,7 @@
(defgeneric alloc-gsl-struct (object)
(:documentation "Allocate the GSL structure.")
(:method ((object gsl-data))
(:method ((object marray))
;; Need to check that all allocations succeeded.
(unless (block-pointer object)
(let ((blockptr (cffi:foreign-alloc 'gsl-block-c)))
......
;; Matrices
;; Liam Healy 2008-04-15 21:57:52EDT matrix.lisp
;; Time-stamp: <2008-11-15 13:35:20EST matrix.lisp>
;; Time-stamp: <2008-12-26 10:10:38EST matrix.lisp>
;; $Id$
(in-package :gsl)
......@@ -18,7 +18,7 @@
(block :pointer)
(owner :int))
(defclass matrix (gsl-data)
(defclass matrix (marray)
()
(:documentation "GSL matrices."))
......
;; Vectors
;; Liam Healy 2008-04-13 09:39:02EDT vector.lisp
;; Time-stamp: <2008-11-29 14:19:54EST vector.lisp>
;; Time-stamp: <2008-12-26 10:10:39EST vector.lisp>
;; $Id$
(in-package :gsl)
......@@ -17,7 +17,7 @@
(block :pointer)
(owner :int))
(defclass mvector (gsl-data)
(defclass mvector (marray)
()
(:documentation "GSL vectors."))
......
;; Eigenvectors and eigenvalues
;; Liam Healy, Sun May 21 2006 - 19:52
;; Time-stamp: <2008-12-25 11:28:44EST eigensystems.lisp>
;; Time-stamp: <2008-12-26 10:28:19EST eigensystems.lisp>
;; $Id$
(in-package :gsl)
......@@ -194,8 +194,8 @@
;;;;****************************************************************************
(defun eigenvalue-eigenvectors-example ()
(letm ((evecs (make-array* 'double-float :dimensions '(3 3)))
(evals (make-array* 'double-float :dimensions 3))
(letm ((evecs (make-marray 'double-float :dimensions '(3 3)))
(evals (make-marray 'double-float :dimensions 3))
(ws (eigen-symmv 3))
(mat #m((20.0d0 -10.0d0 0.0d0)
(-10.0d0 30.0d0 0.0d0)
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-11-30 23:17:52EST gsll.asd>
;; Time-stamp: <2008-12-26 10:54:32EST gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -15,7 +15,7 @@
:components
((:file "init")
(:file "conditions" :depends-on (init))
(:file "gsl-objects" :depends-on (init))
(:file "mobject" :depends-on (init))
(:file "element-types" :depends-on (init))
(:file "number-conversion" :depends-on (init))
(:file "interface"
......@@ -33,15 +33,15 @@
:depends-on (init)
:components
((:file "foreign-friendly")
(:file "data" :depends-on (foreign-friendly))
(:file "vector" :depends-on (data))
(:file "matrix" :depends-on (data vector))
(:file "maref" :depends-on (data vector matrix))
(:file "both" :depends-on (data vector matrix))
(:file "marray" :depends-on (foreign-friendly))
(:file "vector" :depends-on (marray))
(:file "matrix" :depends-on (marray vector))
(:file "maref" :depends-on (marray vector matrix))
(:file "both" :depends-on (marray vector matrix))
(:file "array-tests" :depends-on (both))
(:file "permutation" :depends-on (data))
(:file "permutation" :depends-on (marray))
(:file "combination" ; preload defgeneric
:depends-on (data matrix both permutation))))
:depends-on (marray matrix both permutation))))
(:file "polynomial" :depends-on (init data))
(:module special-functions
:depends-on (init)
......
;; Define examples.
;; Liam Healy 2008-09-07 21:00:48EDT generate-tests.lisp
;; Time-stamp: <2008-11-30 23:02:06EST generate-examples.lisp>
;; Time-stamp: <2008-12-26 10:25:47EST generate-examples.lisp>
;; $Id: $
;;; Define examples that can be displayed by users with the function
......@@ -78,7 +78,7 @@
(defun make-vector-from-pool (type length &optional (starting 0))
"Make a vector of the specified element type and length using the
pool data for the type and starting at the specified point in the pool."
(make-array*
(make-marray
length type
:initial-contents
(make-list-from-pool type length starting)))
......@@ -108,9 +108,9 @@
(let ((matrixp (listp spec)))
(if no-init
;; No initial values, just dimension
`(make-array* ',default-element-type :dimensions ',spec)
`(make-marray ',default-element-type :dimensions ',spec)
;; Initial contents
`(make-array*
`(make-marray
',default-element-type
:initial-contents
',(if matrixp
......
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2008-11-30 23:45:40EST spline-example.lisp>
;; Time-stamp: <2008-12-26 10:29:01EST spline-example.lisp>
;; $Id$
(in-package :gsl)
......@@ -20,8 +20,8 @@
(multiple-value-bind (xarr yarr)
(spline-example-arrays)
(letm ((acc (acceleration))
(cxarr (make-array* 'double-float :dimensions xarr))
(cyarr (make-array* 'double-float :dimensions yarr))
(cxarr (make-marray 'double-float :dimensions xarr))
(cyarr (make-marray 'double-float :dimensions yarr))
(spline (spline *cubic-spline-interpolation* cxarr cyarr)))
(loop for xi from (aref xarr 0) below (aref xarr 9) by 0.01d0
collect (list xi (evaluate-spline spline xi acc))))))
;; Exponential of a matrix
;; Liam Healy 2008-08-10 17:25:35EDT exponential.lisp
;; Time-stamp: <2008-12-07 18:31:37EST exponential.lisp>
;; Time-stamp: <2008-12-26 10:25:48EST exponential.lisp>
;; $Id: $
(in-package :gsl)
......@@ -25,7 +25,7 @@
;;; A matrix of the form [[0, 1], [-1, 0]]
;;; when exponentiated gives [[cos x, sin x], [-sin x, cos x]]
(letm ((mat #m((0.0d0 1.0d0) (-1.0d0 0.0d0)))
(exp (make-array* 'double-float :dimensions '(2 2))))
(exp (make-marray 'double-float :dimensions '(2 2))))
(cl-array (matrix-exponential mat exp)))
#2A((0.5403023058681384d0 0.841470984807895d0)
(-0.841470984807895d0 0.5403023058681384d0))
......
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2008-12-07 18:32:04EST lu.lisp>
;; Time-stamp: <2008-12-26 10:25:47EST lu.lisp>
;; $Id$
(in-package :gsl)
......@@ -136,11 +136,11 @@
(letm ((mmat mat)
(dim (array-dimension mat 0))
(per (make-permutation dim))
(inv (make-array* 'double-float :dimensions (list dim dim))))
(inv (make-marray 'double-float :dimensions (list dim dim))))
(LU-decomposition mmat per)
(lu-invert mmat per inv)
(cl-array inv)))
(save-test lu
(invert-matrix
(make-array* '(2 2) 'double-float :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0))))
(make-marray '(2 2) 'double-float :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0))))
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2008-12-25 11:42:39EST monte-carlo.lisp>
;; Time-stamp: <2008-12-26 10:28:19EST monte-carlo.lisp>
;; $Id$
(in-package :gsl)
......@@ -329,21 +329,21 @@
(defun random-walk-plain-example (&optional (nsamples 500000))
(letm ((ws (monte-carlo-plain 3))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-array* 'double-float :initial-contents (list pi pi pi)))
(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
(rng (random-number-generator *mt19937* 0)))
(monte-carlo-integrate-plain monte-carlo-g lower upper nsamples rng ws)))
(defun random-walk-miser-example (&optional (nsamples 500000))
(letm ((ws (monte-carlo-miser 3))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-array* 'double-float :initial-contents (list pi pi pi)))
(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
(rng (random-number-generator *mt19937* 0)))
(monte-carlo-integrate-miser monte-carlo-g lower upper nsamples rng ws)))
(defun random-walk-vegas-example (&optional (nsamples 500000))
(letm ((ws (monte-carlo-vegas 3))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-array* 'double-float :initial-contents (list pi pi pi)))
(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
(rng (random-number-generator *mt19937* 0)))
(monte-carlo-integrate-vegas monte-carlo-g lower upper nsamples rng ws)))
......
;; Polynomials
;; Liam Healy, Tue Mar 21 2006 - 18:33
;; Time-stamp: <2008-12-25 10:01:23EST polynomial.lisp>
;; Time-stamp: <2008-12-26 10:28:20EST polynomial.lisp>
;; $Id$
(in-package :gsl)
......@@ -155,7 +155,7 @@
real and imaginary parts."
(let ((len (total-size coefficients)))
;; Should this be making a complex array?
(letm ((answer (make-array* 'double-float :dimensions (* 2 (1- len))))
(letm ((answer (make-marray 'double-float :dimensions (* 2 (1- len))))
(ws (complex-workspace len)))
(values-list (polynomial-solve-ws coefficients ws answer)))))
......@@ -181,7 +181,7 @@
(save-test polynomial
(let ((xa #m(0.0d0 1.0d0 2.0d0 3.0d0))
(ya #m(2.5d0 7.2d0 32.7d0 91.0d0))
(dd (make-array* 'double-float :dimensions 4)))
(dd (make-marray 'double-float :dimensions 4)))
(divided-difference dd xa ya)
(list
(polynomial-eval-divided-difference dd xa 0.0d0)
......
;; Dirichlet distribution
;; Liam Healy, Sun Oct 29 2006
;; Time-stamp: <2008-11-30 23:17:06EST dirichlet.lisp>
;; Time-stamp: <2008-12-26 10:28:21EST dirichlet.lisp>
;; $Id$
(in-package :gsl)
......@@ -58,7 +58,7 @@
(save-test dirichlet
(letm ((rng (random-number-generator *mt19937* 0))
(alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
(theta (make-array* 'double-float :dimensions 4)))
(theta (make-marray 'double-float :dimensions 4)))
(dirichlet rng alpha theta)
(cl-array theta))
(letm ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
......
;; Multinomial distribution
;; Liam Healy, Sat Nov 25 2006 - 16:00
;; Time-stamp: <2008-11-30 23:19:48EST multinomial.lisp>
;; Time-stamp: <2008-12-26 10:28:20EST multinomial.lisp>
;; $Id$
(in-package :gsl)
......@@ -51,7 +51,7 @@
(save-test multinomial
(letm ((rng (random-number-generator *mt19937* 0))
(p #m(0.1d0 0.2d0 0.3d0 0.4d0))
(n (make-array* '(signed-byte 32) :dimensions 4)))
(n (make-marray '(signed-byte 32) :dimensions 4)))
(multinomial rng 8 p n)
(cl-array n))
(letm ((p #m(0.1d0 0.2d0 0.3d0 0.4d0))
......
;; Quasi-random sequences in arbitrary dimensions.
;; Liam Healy, Sun Jul 16 2006 - 15:54
;; Time-stamp: <2008-12-25 22:49:26EST quasi.lisp>
;; Time-stamp: <2008-12-26 10:28:22EST quasi.lisp>
;; $Id$
(in-package :gsl)
......@@ -82,7 +82,7 @@
(save-test quasi-random-number-generators
;; This example is given in the GSL documentation
(letm ((gen (quasi-random-number-generator 2 *sobol*))
(vec (make-array* 'double-float :dimensions 2)))
(vec (make-marray 'double-float :dimensions 2)))
(loop repeat 5
do (qrng-get gen vec)
append (coerce (cl-array vec) 'list))))
;; Shuffling and sampling
;; Liam Healy, Sat Dec 2 2006 - 18:40
;; Time-stamp: <2008-11-30 23:15:32EST shuffling-sampling.lisp>
;; Time-stamp: <2008-12-26 10:28:21EST shuffling-sampling.lisp>
;; $Id$
(in-package :gsl)
......@@ -63,11 +63,11 @@
(cl-array v1))
(letm ((rng (random-number-generator *mt19937* 0))
(v1 #31m(1 2 3 4 5 6 7 8))
(v2 (make-array* '(signed-byte 32) :dimensions 4)))
(v2 (make-marray '(signed-byte 32) :dimensions 4)))
(choose-random rng v2 v1)
(cl-array v2))
(letm ((rng (random-number-generator *mt19937* 0))
(v1 #31m(1 2 3 4 5 6 7 8))
(v2 (make-array* '(signed-byte 32) :dimensions 10)))
(v2 (make-marray '(signed-byte 32) :dimensions 10)))
(sample rng v2 v1)
(cl-array v2)))
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
;; Time-stamp: <2008-12-23 22:10:48EST series-acceleration.lisp>
;; Time-stamp: <2008-12-26 10:25:49EST series-acceleration.lisp>
;; $Id$
(in-package :gsl)
......@@ -127,7 +127,7 @@
(sum 0.0d0)
(zeta2 (/ (expt pi 2) 6)))
(letm ((levin (levin maxterms))
(array (make-array* 'double-float :dimensions maxterms)))
(array (make-marray 'double-float :dimensions maxterms)))
(dotimes (n maxterms)
(setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
(incf sum (maref array n)))
......
;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
;; Time-stamp: <2008-12-25 10:35:32EST linear-least-squares.lisp>
;; Time-stamp: <2008-12-26 10:28:25EST linear-least-squares.lisp>
;; $Id$
(in-package :gsl)
......@@ -322,11 +322,11 @@
(defun mv-linear-least-squares-example (data)
"Second example in Section 36.5 of the GSL manual."
(letm ((n (length data)) chisq
(x (make-array* 'double-float :dimensions (list n 3)))
(cov (make-array* 'double-float :dimensions '(3 3)))
(y (make-array* 'double-float :dimensions n))
(w (make-array* 'double-float :dimensions n))
(c (make-array* 'double-float :dimensions 3)))
(x (make-marray 'double-float :dimensions (list n 3)))
(cov (make-marray 'double-float :dimensions '(3 3)))
(y (make-marray 'double-float :dimensions n))
(w (make-marray 'double-float :dimensions n))
(c (make-marray 'double-float :dimensions 3)))
(loop for i from 0
for row in data do
(setf (maref X i 0) 1.0d0
......
;; Multivariate minimization.
;; Liam Healy <Tue Jan 8 2008 - 21:28>
;; Time-stamp: <2008-12-25 10:32:16EST minimization-multi.lisp>
;; Time-stamp: <2008-12-26 10:29:01EST minimization-multi.lisp>
;; $Id$
(in-package :gsl)
......@@ -411,7 +411,7 @@
(defun multimin-example-nelder-mead ()
(letm ((initial #m(5.0d0 7.0d0))
(step-size (make-array* 'double-float :dimensions 2)))
(step-size (make-marray 'double-float :dimensions 2)))
(set-all step-size 1.0d0)
(letm ((minimizer
(mfminimizer *simplex-nelder-mead* 2 parabaloid-f initial step-size)))
......
;; Nonlinear least squares fitting.
;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
;; Time-stamp: <2008-12-23 22:32:19EST nonlinear-least-squares.lisp>
;; Time-stamp: <2008-12-26 10:25:45EST nonlinear-least-squares.lisp>
;; $Id$
(in-package :gsl)
......@@ -355,13 +355,13 @@
(make-exponent-fit-data
:n *number-of-observations*
:y
(let ((arr (make-array* *number-of-observations* 'double-float)))
(let ((arr (make-marray *number-of-observations* 'double-float)))
(letm ((rng (random-number-generator *mt19937* 0)))
(dotimes (i *number-of-observations* arr)
(setf (aref arr i)
(+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0))))))
:sigma
(make-array* *number-of-observations* 'double-float :initial-element 0.1d0))))
(make-marray *number-of-observations* 'double-float :initial-element 0.1d0))))
(defun exponential-residual (x f)
"Compute the negative of the residuals with the exponential model
......@@ -412,7 +412,7 @@
(defun solve-nonlinear-least-squares-example ()
(letm ((init #m(1.0d0 0.0d0 0.0d0))
(covariance
(make-array* 'double-float
(make-marray 'double-float
:dimensions
(list *number-of-parameters* *number-of-parameters*)))
(fit (nonlinear-fdffit
......
;; Bessel functions
;; Liam Healy, Fri Mar 17 2006 - 18:42
;; Time-stamp: <2008-11-30 23:22:00EST bessel.lisp>
;; Time-stamp: <2008-12-26 10:28:24EST bessel.lisp>
;; $Id$
(in-package :gsl)
......@@ -460,61 +460,61 @@
(cylindrical-bessel-J0 4.0d0)
(cylindrical-bessel-J1 4.0d0)
(cylindrical-bessel-Jn 2 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(cylindrical-bessel-Jn-array 2.0d0 besarr 2)
(cl-array besarr))
(cylindrical-bessel-Y0 4.0d0)
(cylindrical-bessel-Y1 4.0d0)
(cylindrical-bessel-Yn 3 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(cylindrical-bessel-Yn-array 2.0d0 besarr 2)
(cl-array besarr))
(cylindrical-bessel-I0 4.0d0)
(cylindrical-bessel-I1 4.0d0)
(cylindrical-bessel-In 3 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(cylindrical-bessel-In-array 2.0d0 besarr 2)
(cl-array besarr))
(cylindrical-bessel-I0-scaled 4.0d0)
(cylindrical-bessel-I1-scaled 4.0d0)
(cylindrical-bessel-In-scaled 3 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(cylindrical-bessel-In-scaled-array 2.0d0 besarr 2)
(cl-array besarr))
(cylindrical-bessel-K0 4.0d0)
(cylindrical-bessel-K1 4.0d0)
(cylindrical-bessel-Kn 2 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(cylindrical-bessel-Kn-array 2.0d0 besarr 2)
(cl-array besarr))
(cylindrical-bessel-K0-scaled 4.0d0)
(cylindrical-bessel-K1-scaled 4.0d0)
(cylindrical-bessel-Kn-scaled 2 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(cylindrical-bessel-Kn-array 2.0d0 besarr 2)
(cl-array besarr))
(spherical-bessel-j0 4.0d0)
(spherical-bessel-j1 4.0d0)
(spherical-bessel-j2 4.0d0)
(spherical-bessel-jl 3 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(spherical-bessel-jl-array 4.0d0 besarr)
(cl-array besarr))
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(spherical-bessel-jl-steed-array 4.0d0 besarr)
(cl-array besarr))
(spherical-bessel-y0 4.0d0)
(spherical-bessel-y1 4.0d0)
(spherical-bessel-y2 4.0d0)
(spherical-bessel-yl 2 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(spherical-bessel-yl-array 4.0d0 besarr)
(cl-array besarr))
(spherical-bessel-i0-scaled 4.0d0)
(spherical-bessel-i1-scaled 4.0d0)
(spherical-bessel-i2-scaled 4.0d0)
(spherical-bessel-il-scaled 3 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(spherical-bessel-il-scaled-array 4.0d0 besarr)
(cl-array besarr))
(spherical-bessel-k0-scaled 4.0d0)
......@@ -522,7 +522,7 @@
(spherical-bessel-k2-scaled 4.0d0)
(spherical-bessel-kl-scaled 3 4.0d0)
(letm ((besarr (make-array* 'double-float :dimensions 4)))
(letm ((besarr (make-marray 'double-float :dimensions 4)))
(spherical-bessel-kl-scaled-array 4.0d0 besarr)
(cl-array besarr))
(bessel-jnu 3.0d0 4.0d0)
......
;; Coulumb functions
;; Liam Healy, Sat Mar 18 2006 - 23:23
;; Time-stamp: <2008-11-30 23:44:05EST coulomb.lisp>
;; Time-stamp: <2008-12-26 10:28:23EST coulomb.lisp>
;; $Id$
(in-package :gsl)
......@@ -128,16 +128,16 @@
(hydrogenicr-1 1.0d0 2.5d0)
(hydrogenicr 3 1 1.0d0 2.5d0)
(coulomb-wave-FG 0.0d0 1.0d0 2.0d0 0)
(letm ((arr (make-array* 'double-float :dimensions 3)))
(letm ((arr (make-marray 'double-float :dimensions 3)))
(coulomb-wave-F-array 0.0d0 1.0d0 2.0d0 arr)
(cl-array arr))
(coulomb-wave-fg 1.0d0 2.0d0 2.5d0 1)
(letm ((Farr (make-array* 'double-float :dimensions 3))
(Garr (make-array* 'double-float :dimensions 3)))
(letm ((Farr (make-marray 'double-float :dimensions 3))
(Garr (make-marray 'double-float :dimensions 3)))
(coulomb-wave-FG-array 1.5d0 1.0d0 1.0d0 Farr Garr)
(append (coerce (cl-array Farr) 'list) (coerce (cl-array Garr) 'list)))
(letm ((arr (make-array* 'double-float :dimensions 3)))