Commit 75d17ffe authored by Liam Healy's avatar Liam Healy

Use grid:gref instead of maref

parent bf2413a3
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2010-06-29 22:15:23EDT basis-splines.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT basis-splines.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -125,18 +125,18 @@
(let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float))
(yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
(sample rng :gaussian :sigma sigma))))
(setf (maref x i) xi
(maref y i) yi
(maref w i) (/ (expt sigma 2)))))
(setf (grid:gref x i) xi
(grid:gref y i) yi
(grid:gref w i) (/ (expt sigma 2)))))
;; Uniform breakpoints [0, 15]
(uniform-knots 0.0d0 15.0d0 bw)
;; Fit matrix
(dotimes (i ndata)
;; Compute B_j
(evaluate bw (maref x i) :b B)
(evaluate bw (grid:gref x i) :b B)
;; Fill in row i of X
(dotimes (j ncoeffs)
(setf (maref Xmatrix i j) (maref B j))))
(setf (grid:gref Xmatrix i j) (grid:gref B j))))
;; Do the fit
(linear-mfit Xmatrix y c w nil cov mw)
;; Return the smoothed curve
......
;; Copy grid:foreign-arrays to/from CL arrays
;; Liam Healy 2009-02-11 19:28:44EST copy-cl.lisp
;; Time-stamp: <2010-06-29 22:51:21EDT copy-cl.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT copy-cl.lisp>
;;
;; Copyright 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -42,7 +42,7 @@
(unless (equal (dimensions source) (array-dimensions destination))
(error 'grid:array-mismatch))
(loop for i below (dim0 source)
do (setf (aref destination i) (maref source i)))
do (setf (aref destination i) (grid:gref source i)))
destination)
(defmethod grid:copy-to-destination ((source matrix) (destination array))
......@@ -50,14 +50,14 @@
(error 'grid:array-mismatch))
(loop for i below (dim0 source) do
(loop for j below (dim1 source) do
(setf (aref destination i j) (maref source i j))))
(setf (aref destination i j) (grid:gref source i j))))
destination)
(defmethod grid:copy-to-destination ((source array) (destination mvector))
(unless (equal (array-dimensions source) (dimensions destination))
(error 'grid:array-mismatch))
(loop for i below (length source)
do (setf (maref destination i) (aref source i)))
do (setf (grid:gref destination i) (aref source i)))
destination)
(defmethod grid:copy-to-destination ((source array) (destination matrix))
......@@ -65,7 +65,7 @@
(error 'grid:array-mismatch))
(loop for i below (array-dimension source 0) do
(loop for j below (array-dimension source 1) do
(setf (maref destination i j) (aref source i j))))
(setf (grid:gref destination i j) (aref source i j))))
destination)
(defmethod grid:copy-to-destination ((source grid:foreign-array) (destclass (eql 'array)))
......
;; Permutations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2010-06-30 18:00:28EDT permutation.lisp>
;; Time-stamp: <2010-06-30 19:53:22EDT permutation.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -333,7 +333,7 @@
while (nth-value 1 (permutation-previous perm)))))
(save-test permutation
(let ((perm-1 (make-permutation 4 t))) ;maref
(let ((perm-1 (make-permutation 4 t)))
(grid:gref perm-1 2))
(let ((perm-1 (make-permutation 4 t))) ;grid:contents
(grid:contents perm-1))
......
......@@ -343,9 +343,9 @@ area of memory.</p>
</table>
</center>
<p>
Individual elements are obtained using <code>maref</code> (analogous
Individual elements are obtained using <code>grid:gref</code> (analogous
to Lisp's <code>aref</code>), and are set
with <code>setf maref</code>. A complete CL array may
with <code>setf grid:gref</code>. A complete CL array may
be extracted with the function <code>#'cl-array</code>.
</p>
<p>Copying grid:foreign-arrays is performed with the function <code>copy</code>.
......@@ -472,7 +472,7 @@ and arrays used internally or for function return.
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
<small>
Time-stamp: <2010-06-29 22:51:17EDT index.html>
Time-stamp: <2010-06-30 19:57:28EDT index.html>
</small>
<!-- hhmts end -->
</div>
......
;; Example FFT: transform a pulse (using the "clean" fft interface)
;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24
;; Time-stamp: <2010-06-29 22:15:25EDT example.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT example.lisp>
;;
;; Copyright 2009 Sumant Oemrawsingh, Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -42,10 +42,10 @@
(assert (and (integerp dimension) (> dimension 20)))
(let ((pulse (grid:make-foreign-array element-type :dimensions dimension))
(init-value (coerce 1 element-type)))
(setf (maref pulse 0) init-value)
(setf (grid:gref pulse 0) init-value)
(loop for i from 1 to 10
do (setf (maref pulse i) init-value
(maref pulse (- dimension i)) init-value))
do (setf (grid:gref pulse i) init-value
(grid:gref pulse (- dimension i)) init-value))
(forward-fourier-transform pulse)))
(save-test
......@@ -80,7 +80,7 @@
:dimensions (list (* stride dimension)))))
(loop for i from 0 below (* stride dimension) by stride
do
(setf (maref vec i)
(setf (grid:gref vec i)
(if (subtypep element-type 'complex)
(coerce (complex (urand) (urand)) element-type)
(complex (coerce (urand) element-type)))))
......@@ -93,8 +93,8 @@
(grid:component-float-type (element-type complex-vector))
:dimensions (dimensions complex-vector))))
(loop for i below (total-size complex-vector) do
(setf (maref real-vector i)
(realpart (maref complex-vector i))))
(setf (grid:gref real-vector i)
(realpart (grid:gref complex-vector i))))
real-vector))
(defun size-vector-real (vector &key (stride 1))
......
......@@ -76,12 +76,12 @@
(zero-pos (if (evenp n) (- split 1) split)))
;; Positive frequencies
(loop for i from zero-pos to (- n 1)
do (setf (maref shifted (* i stride))
(maref vector (* (- i zero-pos) stride))))
do (setf (grid:gref shifted (* i stride))
(grid:gref vector (* (- i zero-pos) stride))))
;; Negative frequencies
(loop for i from (+ split 1) to (- n 1)
do (setf (maref shifted (* (- i split 1) stride))
(maref vector (* i stride))))
do (setf (grid:gref shifted (* (- i split 1) stride))
(grid:gref vector (* i stride))))
shifted))
(export 'fft-inverse-shift)
......@@ -96,10 +96,10 @@
(zero-pos (if (evenp n) (- split 1) split)))
;; Positive frequencies
(loop for i from zero-pos to (- n 1)
do (setf (maref inv-shifted (* (- i zero-pos) stride))
(maref vector (* i stride))))
do (setf (grid:gref inv-shifted (* (- i zero-pos) stride))
(grid:gref vector (* i stride))))
;; Negative frequencies
(loop for i from (+ split 1) to (- n 1)
do (setf (maref inv-shifted (* i stride))
(maref vector (* (- i split 1) stride))))
do (setf (grid:gref inv-shifted (* i stride))
(grid:gref vector (* (- i split 1) stride))))
inv-shifted))
;; Discrete Hankel Transforms.
;; Liam Healy, Sat Dec 8 2007 - 16:50
;; Time-stamp: <2010-06-29 22:15:22EDT hankel.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT hankel.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -81,7 +81,7 @@
(let ((hank (make-hankel 128 0.0d0 100.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128
do (setf (maref in n)
do (setf (grid:gref in n)
(/ (1+ (expt (sample-x-hankel hank n) 2)))))
(copy (apply-hankel hank in) 'array))
;; Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
......@@ -89,13 +89,13 @@
(let ((hank (make-hankel 128 1.0d0 20.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128
do (setf (maref in n) (exp (- (sample-x-hankel hank n)))))
do (setf (grid:gref in n) (exp (- (sample-x-hankel hank n)))))
(copy (apply-hankel hank in) 'array))
;; Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
(let ((hank (make-hankel 128 1.0d0 1.0d0))
(in (grid:make-foreign-array 'double-float :dimensions 128)))
(loop for n from 0 below 128
do (setf (maref in n)
do (setf (grid:gref in n)
(let ((x (sample-x-hankel hank n)))
(* x (- 1 (expt x 2))))))
(copy (apply-hankel hank in) 'array)))
;; N-tuples
;; Liam Healy Sat Feb 3 2007 - 12:53
;; Time-stamp: <2009-12-27 09:47:40EST ntuple.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT ntuple.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -235,7 +235,7 @@
(close-ntuple ntuple)
;; Check the histogram to see if it matches ntuple-example-values
(dotimes (row 100 t)
(unless (= (maref histo row) (aref answer row))
(unless (= (grid:gref histo row) (aref answer row))
(return nil))))))
(save-test ntuple
......
;; Updating and accessing histogram elements.
;; Liam Healy, Mon Jan 1 2007 - 14:43
;; Time-stamp: <2009-12-27 09:48:05EST updating-accessing.lisp>
;; Time-stamp: <2010-06-30 21:02:26EDT updating-accessing.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -42,9 +42,9 @@
histograms for a small range of a larger dataset, ignoring the values
outside the range of interest.")
(defmfun maref ((histogram histogram) i &optional i2 element-type)
(defmfun grid:gref ((histogram histogram) &rest indices)
"gsl_histogram_get"
(((mpointer histogram) :pointer) (i sizet))
(((mpointer histogram) :pointer) ((first indices) sizet))
:definition :method
:c-return :double
:documentation ; FDL
......@@ -127,22 +127,22 @@
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(maref histo 1))
(grid:gref histo 1))
(let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(maref histo 2))
(grid:gref histo 2))
(let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(maref histo 6))
(grid:gref histo 6))
(let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
(increment histo 6.9d0 2.0d0)
(maref histo 16))
(grid:gref histo 16))
(let ((histo (make-histogram 10)))
(set-ranges-uniform histo 0.0d0 10.0d0)
(increment histo 2.7d0)
......
;; Generate a lambda that calls the user function; will be called by callback.
;; Liam Healy
;; Time-stamp: <2010-06-29 22:49:08EDT funcallable.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT funcallable.lisp>
;;
;; Copyright 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -71,7 +71,7 @@
array, or a scalar."
(if (parse-callback-argspec argspec 'dimensions)
(if (eql (parse-callback-argspec argspec 'array-type) :foreign-array)
`(maref
`(grid:gref
,foreign-variable-name
,@(let ((dims (value-from-dimensions argspec dimension-values)))
(if (= (length dims) 2) ; matrix
......
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2010-06-29 22:15:24EDT spline-example.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT spline-example.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -35,7 +35,7 @@
(loop for i from 0.0d0 below 10.0d0
collect (+ i (cos (expt i 2))))))
(spline (make-spline +cubic-spline-interpolation+ xarr yarr)))
(loop for xi from (maref xarr 0) below (maref xarr 9) by step
(loop for xi from (grid:gref xarr 0) below (grid:gref xarr 9) by step
collect (list xi (evaluate spline xi)))))
(defun evaluate-integral-example (&optional (intervals 4))
......
;; BLAS level 1, Vector operations
;; Liam Healy, Wed Apr 26 2006 - 15:23
;; Time-stamp: <2010-06-29 22:15:23EDT blas1.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT blas1.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -264,7 +264,7 @@
(sines (array-default 8 t))
(cosines (array-default 8 t)))
(loop for i below 8 do
(setf (maref sines i) (sin (maref angles i)))
(setf (maref cosines i) (cos (maref angles i))))
(setf (grid:gref sines i) (sin (grid:gref angles i)))
(setf (grid:gref cosines i) (cos (grid:gref angles i))))
(givens-rotation v1 v2 cosines sines)
(list (cl-array v1) (cl-array v2))))
;; Cholesky Decomposition
;; Liam Healy, Wed May 3 2006 - 16:38
;; Time-stamp: <2010-06-29 22:51:20EDT cholesky.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT cholesky.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -98,7 +98,7 @@
(let ((decomp (cholesky-decomposition (copy matrix))))
(dotimes (row (dim0 matrix) decomp)
(loop for col from (1+ row) below (dim1 matrix) do
(setf (maref decomp row col) 0.0d0)))
(setf (grid:gref decomp row col) 0.0d0)))
(matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans)))
(defun test-cholesky-invert-dim (matrix)
......
;; Tridiagonal and Bidiagonal matrices
;; Liam Healy, Thu May 4 2006 - 15:43
;; Time-stamp: <2010-06-29 22:15:24EDT diagonal.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT diagonal.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -255,11 +255,11 @@
(subdiag (mvec (1- n) :initial-element 0.5d0)))
(setf
;; i=0 matrix elements
(maref diag 0) 1d0
(maref superdiag 0) 0d0
(maref b 0) (coerce (expt (1- n) 2) 'double-float)
(grid:gref diag 0) 1d0
(grid:gref superdiag 0) 0d0
(grid:gref b 0) (coerce (expt (1- n) 2) 'double-float)
;; i=n-1 matrix elements
(maref diag (1- n)) 1d0
(maref subdiag (- n 2)) 0d0
(maref b (1- n)) 0d0)
(grid:gref diag (1- n)) 1d0
(grid:gref subdiag (- n 2)) 0d0
(grid:gref b (1- n)) 0d0)
(solve-tridiagonal diag superdiag subdiag b x))))
;; QR decomposition
;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
;; Time-stamp: <2010-06-29 22:51:19EDT qr.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT qr.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -258,7 +258,7 @@
dim1 nil))
(qr1
(create-matrix
(lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j))))
(lambda (i j) (+ (grid:gref matrix i j) (* (grid:gref u i) (grid:gref v j))))
dim0 dim1))
(qr2 (copy matrix))
(w (grid:make-foreign-array 'double-float :dimensions dim0 :initial-element 0)))
......
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2010-06-29 22:15:22EDT ode-system.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT ode-system.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -43,11 +43,11 @@
(,ctime (grid:make-foreign-array 'double-float :dimensions 1))
(,cstep (grid:make-foreign-array 'double-float :dimensions 1)))
(symbol-macrolet
((,time (maref ,ctime 0))
(,step-size (maref ,cstep 0))
((,time (grid:gref ,ctime 0))
(,step-size (grid:gref ,cstep 0))
,@(loop for symb in dependent
for i from 0
collect `(,symb (maref ,dep ,i))))
collect `(,symb (grid:gref ,dep ,i))))
(flet ((next-step ()
(apply-evolution
evolve ,ctime ,dep ,cstep control stepperobj ,max-time)))
......
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
;; Time-stamp: <2010-06-29 22:15:24EDT series-acceleration.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT series-acceleration.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -95,8 +95,8 @@
(let ((levin (make-levin maxterms))
(array (grid:make-foreign-array 'double-float :dimensions maxterms)))
(dotimes (n maxterms)
(setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
(incf sum (maref array n)))
(setf (grid:gref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
(incf sum (grid:gref array n)))
(multiple-value-bind (accelerated-sum error)
(accelerate array levin)
(when print-explanation
......
;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
;; Time-stamp: <2010-06-29 22:15:22EDT linear-least-squares.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT linear-least-squares.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -315,15 +315,15 @@
(loop for i from 0 below (dim0 x)
do
(format t "data: ~12,5f ~12,5f ~12,5f~&"
(maref x i)
(maref y i)
(/ (maref w i))))
(grid:gref x i)
(grid:gref y i)
(/ (grid:gref w i))))
(loop for i from -30 below 130 by 10 ; don't print everything
for
xf = (+ (maref x 0)
xf = (+ (grid:gref x 0)
(* (/ i 100)
(- (maref x (1- (dim0 x)))
(maref x 0))))
(- (grid:gref x (1- (dim0 x)))
(grid:gref x 0))))
do
(multiple-value-bind (yf yferr)
(linear-estimate xf c0 c1 cov00 cov01 cov11)
......@@ -353,26 +353,26 @@
(w (grid:make-foreign-array 'double-float :dimensions n)))
(loop for i from 0
for row in data do
(setf (maref X i 0) 1.0d0
(maref X i 1) (first row)
(maref X i 2) (expt (first row) 2)
(maref y i) (second row)
(maref w i) (/ (expt (third row) 2))))
(setf (grid:gref X i 0) 1.0d0
(grid:gref X i 1) (first row)
(grid:gref X i 2) (expt (first row) 2)
(grid:gref y i) (second row)
(grid:gref w i) (/ (expt (third row) 2))))
(multiple-value-bind (parameters cov chisq)
(linear-mfit X y 3 w)
(when print-details
(format t "Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2~&"
(maref parameters 0) (maref parameters 1) (maref parameters 2))
(grid:gref parameters 0) (grid:gref parameters 1) (grid:gref parameters 2))
(format t "Covariance matrix:~&")
(format t "~10,8f ~10,8f ~10,8f~&"
(maref cov 0 0) (maref cov 0 1) (maref cov 0 2))
(grid:gref cov 0 0) (grid:gref cov 0 1) (grid:gref cov 0 2))
(format t "~10,8f ~10,8f ~10,8f~&"
(maref cov 1 0) (maref cov 1 1) (maref cov 1 2))
(grid:gref cov 1 0) (grid:gref cov 1 1) (grid:gref cov 1 2))
(format t "~10,8f ~10,8f ~10,8f~&"
(maref cov 2 0) (maref cov 2 1) (maref cov 2 2))
(grid:gref cov 2 0) (grid:gref cov 2 1) (grid:gref cov 2 2))
(format t "Chisq = ~10,6f~&" chisq))
(values
(maref parameters 0) (maref parameters 1) (maref parameters 2)
(grid:gref parameters 0) (grid:gref parameters 1) (grid:gref parameters 2)
chisq))))
(save-test linear-least-squares
......
;; Multivariate minimization.
;; Liam Healy <Tue Jan 8 2008 - 21:28>
;; Time-stamp: <2010-06-29 22:49:05EDT minimization-multi.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT minimization-multi.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -384,13 +384,13 @@
(when print-steps
(let ((x (solution minimizer)))
(format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~40t~8,3f~&"
iter (maref x 0) (maref x 1)
iter (grid:gref x 0) (grid:gref x 1)
(function-value minimizer)
size)))
finally
(return
(let ((x (solution minimizer)))
(values (maref x 0) (maref x 1) (function-value minimizer))))))))
(values (grid:gref x 0) (grid:gref x 1) (function-value minimizer))))))))
;;; Example using derivatives, taking a vector argument.
;;; Note that these functions are written to read objects of
......@@ -400,8 +400,8 @@
(defun paraboloid-vector (gsl-vector)
"A paraboloid function of two arguments, given in GSL manual Sec. 35.4.
This version takes a vector-double-float argument."
(let ((x (maref gsl-vector 0))
(y (maref gsl-vector 1))
(let ((x (grid:gref gsl-vector 0))
(y (grid:gref gsl-vector 1))
(dp0 (aref *paraboloid-center* 0))
(dp1 (aref *paraboloid-center* 1)))
(+ (* 10 (expt (- x dp0) 2))
......@@ -410,13 +410,13 @@
(defun paraboloid-derivative
(arguments-gv-pointer derivative-gv-pointer)
(let ((x (maref arguments-gv-pointer 0))
(y (maref arguments-gv-pointer 1))
(let ((x (grid:gref arguments-gv-pointer 0))
(y (grid:gref arguments-gv-pointer 1))
(dp0 (aref *paraboloid-center* 0))
(dp1 (aref *paraboloid-center* 1)))
(setf (maref derivative-gv-pointer 0)
(setf (grid:gref derivative-gv-pointer 0)
(* 20 (- x dp0))
(maref derivative-gv-pointer 1)
(grid:gref derivative-gv-pointer 1)
(* 40 (- y dp1)))))
(defun paraboloid-and-derivative
......@@ -447,12 +447,12 @@
(when print-steps
(let ((x (solution minimizer)))
(format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~&"
iter (maref x 0) (maref x 1)
iter (grid:gref x 0) (grid:gref x 1)
(function-value minimizer))))
finally
(return
(let ((x (solution minimizer)))
(values (maref x 0) (maref x 1) (function-value minimizer)))))))
(values (grid:gref x 0) (grid:gref x 1) (function-value minimizer)))))))
(defun paraboloid-derivative-scalar (x y)
(let ((dp0 (aref *paraboloid-center* 0))
......@@ -486,12 +486,12 @@
(when print-steps
(let ((x (solution minimizer)))
(format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~&"
iter (maref x 0) (maref x 1)
iter (grid:gref x 0) (grid:gref x 1)
(function-value minimizer))))
finally
(return
(let ((x (solution minimizer)))
(values (maref x 0) (maref x 1) (function-value minimizer)))))))
(values (grid:gref x 0) (grid:gref x 1) (function-value minimizer)))))))
......
;; Nonlinear least squares fitting.
;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
;; Time-stamp: <2010-06-29 22:49:06EDT nonlinear-least-squares.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT nonlinear-least-squares.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -301,7 +301,7 @@
(let ((arr (grid:make-foreign-array 'double-float :dimensions number-of-observations))
(rng (make-random-number-generator +mt19937+ 0)))
(dotimes (i number-of-observations arr)
(setf (maref arr i)
(setf (grid:gref arr i)
(+ 1 (* 5 (exp (* -1/10 i)))
(sample rng :gaussian :sigma 0.1d0)))))
:sigma
......@@ -311,32 +311,32 @@
(defun exponential-residual (x f)
"Compute the negative of the residuals with the exponential model
for the nonlinear least squares example."
(let ((A (maref x 0))
(lambda (maref x 1))
(b (maref x 2)))
(let ((A (grid:gref x 0))
(lambda (grid:gref x 1))
(b (grid:gref x 2)))
(symbol-macrolet
((y (exponent-fit-data-y *nlls-example-data*))
(sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i (exponent-fit-data-n *nlls-example-data*))
(setf (maref f i)
(setf (grid:gref f i)
;; the difference model - observation = - residual
(/ (- (+ (* A (exp (* (- lambda) i))) b) (maref y i))
(maref sigma i)))))))
(/ (- (+ (* A (exp (* (- lambda) i))) b) (grid:gref y i))
(grid:gref sigma i)))))))
(defun exponential-residual-derivative (x jacobian)
"Compute the partial derivatives of the negative of the
residuals with the exponential model
for the nonlinear least squares example."
(let ((A (maref x 0))
(lambda (maref x 1)))
(let ((A (grid:gref x 0))
(lambda (grid:gref x 1)))
(symbol-macrolet
((sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i (exponent-fit-data-n *nlls-example-data*))
(let ((e (exp (* (- lambda) 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)))))))
(s (grid:gref sigma i)))
(setf (grid:gref jacobian i 0) (/ e s)
(grid:gref jacobian i 1) (* -1 i A (/ e s))
(grid:gref jacobian i 2) (/ s)))))))
(defun exponential-residual-fdf (x f jacobian)
"Compute the function and partial derivatives of the negative of the
......@@ -363,8 +363,8 @@
'(exponential-residual
exponential-residual-derivative exponential-residual-fdf)
init nil)))
(macrolet ((fitx (i) `(maref (solution fit) ,i))
(err (i) `(sqrt (maref covariance ,i ,i))))
(macrolet ((fitx (i) `(grid:gref (solution fit) ,i))
(err (i) `(sqrt (grid:gref covariance ,i ,i))))
(when print-steps
(format t "iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g~&"
0 (fitx 0) (fitx 1) (fitx 2)
......
;;; Multivariate roots.
;;; Liam Healy 2008-01-12 12:49:08
;;; Time-stamp: <2010-06-29 22:51:17EDT roots-multi.lisp>
;;; Time-stamp: <2010-06-30 19:57:28EDT roots-multi.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -413,15 +413,15 @@
(when print-steps
(format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
iter
(maref argval 0)
(maref argval 1)
(maref fnval 0)
(maref fnval 1)))
(grid:gref argval 0)
(grid:gref argval 1)
(grid:gref fnval 0)
(grid:gref fnval 1)))
finally (return
(values (maref argval 0)
(maref argval 1)
(maref fnval 0)
(maref fnval 1))))))
(values (grid:gref argval 0)
(grid:gref argval 1)
(grid:gref fnval 0)
(grid:gref fnval 1))))))
(defun rosenbrock-df (arg0 arg1)
"The partial derivatives of the Rosenbrock functions."
......@@ -447,10 +447,10 @@
(when print-steps
(format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
iter
(maref argval 0)
(maref argval 1)
(maref fnval 0)
(maref fnval 1)))))
(grid:gref argval 0)