diff --git a/gsll.asd b/gsll.asd index fe67b96e6a05d251e7b359a6eb82136338ae97f1..d3dcf4385f98387cf56442c27b24b316e43d85a9 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2008-08-31 10:05:43EDT gsll.asd> +;; Time-stamp: <2008-08-31 12:49:38EDT gsll.asd> ;; $Id$ (asdf:defsystem "gsll" @@ -181,9 +181,8 @@ :components ((:file "roots-one") (:file "minimization-one") - ;(:file "roots-multi" :depends-on (roots-one)) - ;(:file "minimization-multi") - ;(:file "linear-least-squares") - ;(:file "nonlinear-least-squares") - )) + (:file "roots-multi" :depends-on (roots-one)) + (:file "minimization-multi") + (:file "linear-least-squares") + (:file "nonlinear-least-squares"))) (:file "basis-splines" :depends-on (init data)))) diff --git a/solve-minimize-fit/linear-least-squares.lisp b/solve-minimize-fit/linear-least-squares.lisp index 4b29a6d6b72509e0d3546f87168f64da66ba5889..6c4fdc91c37d39d6498d9aa01febf22a92aca46e 100644 --- a/solve-minimize-fit/linear-least-squares.lisp +++ b/solve-minimize-fit/linear-least-squares.lisp @@ -1,6 +1,6 @@ ;; Linear least squares, or linear regression ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp> -;; Time-stamp: <2008-03-09 19:30:56EDT linear-least-squares.lisp> +;; Time-stamp: <2008-08-31 12:44:33EDT linear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -17,9 +17,9 @@ (defmfun linear-fit (x y &optional (x-stride 1) (y-stride 1)) "gsl_fit_linear" - (((gsl-array x) :pointer) (x-stride size) - ((gsl-array y) :pointer) (y-stride size) - ((dim0 x) size) (c0 :double) (c1 :double) + (((c-pointer x) :pointer) (x-stride sizet) + ((c-pointer y) :pointer) (y-stride sizet) + ((dim0 x) sizet) (c0 :double) (c1 :double) (cov00 :double) (cov01 :double) (cov11 :double) (sumsq :double)) :documentation ; FDL @@ -38,10 +38,10 @@ (defmfun weighted-linear-fit (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1)) "gsl_fit_wlinear" - (((gsl-array x) :pointer) (x-stride size) - ((gsl-array weight) :pointer) (weight-stride size) - ((gsl-array y) :pointer) (y-stride size) - ((dim0 x) size) (c0 :double) (c1 :double) + (((c-pointer x) :pointer) (x-stride sizet) + ((c-pointer weight) :pointer) (weight-stride sizet) + ((c-pointer y) :pointer) (y-stride sizet) + ((dim0 x) sizet) (c0 :double) (c1 :double) (cov00 :double) (cov01 :double) (cov11 :double) (chisq :double)) :documentation ; FDL @@ -77,9 +77,9 @@ (defmfun multiplier-fit (x y &optional (x-stride 1) (y-stride 1)) "gsl_fit_mul" - (((gsl-array x) :pointer) (x-stride size) - ((gsl-array y) :pointer) (y-stride size) - ((dim0 x) size) (c1 :double) (cov11 :double) + (((c-pointer x) :pointer) (x-stride sizet) + ((c-pointer y) :pointer) (y-stride sizet) + ((dim0 x) sizet) (c1 :double) (cov11 :double) (sumsq :double)) :documentation ; FDL "The best-fit linear regression coefficient c1 of the model Y = c_1 @@ -93,10 +93,10 @@ (defmfun weighted-multiplier-fit (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1)) "gsl_fit_wmul" - (((gsl-array x) :pointer) (x-stride size) - ((gsl-array weight) :pointer) (weight-stride size) - ((gsl-array y) :pointer) (y-stride size) - ((dim0 x) size) (c1 :double) (cov11 :double) + (((c-pointer x) :pointer) (x-stride sizet) + ((c-pointer weight) :pointer) (weight-stride sizet) + ((c-pointer y) :pointer) (y-stride sizet) + ((dim0 x) sizet) (c1 :double) (cov11 :double) (chisq :double)) :documentation ; FDL "Compute the best-fit linear regression coefficient @@ -130,7 +130,7 @@ (defmfun allocate-fit-workspace (number-of-observations number-of-parameters) "gsl_multifit_linear_alloc" - ((number-of-observations size) (number-of-parameters size)) + ((number-of-observations sizet) (number-of-parameters sizet)) :c-return :pointer :index (letm fit-workspace) :documentation ; FDL @@ -147,9 +147,9 @@ (defmfun linear-mfit (model observations parameters covariance tolerance workspace) "gsl_multifit_linear" - (((pointer model) :pointer) ((pointer observations) :pointer) + (((mpointer model) :pointer) ((mpointer observations) :pointer) (tolerance :double) - ((pointer parameters) :pointer) (covariance :pointer) (chisq :double) + ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double) (workspace :pointer)) :documentation ; FDL "Compute the best-fit parameters c of the model @@ -168,10 +168,10 @@ (defmfun linear-mfit-svd (model observations parameters covariance tolerance workspace) "gsl_multifit_linear_svd" - (((pointer model) :pointer) ((pointer observations) :pointer) + (((mpointer model) :pointer) ((mpointer observations) :pointer) (tolerance :double) - (rank size) - ((pointer parameters) :pointer) (covariance :pointer) (chisq :double) + (rank sizet) + ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double) (workspace :pointer)) :return ((dcref chisq) (scref rank)) :documentation ; FDL @@ -195,11 +195,11 @@ (defmfun weighted-linear-mfit (model weight observations parameters covariance workspace) "gsl_multifit_wlinear" - (((pointer model) :pointer) - ((pointer weight) :pointer) - ((pointer observations) :pointer) - ((pointer parameters) :pointer) - ((pointer covariance) :pointer) (chisq :double) + (((mpointer model) :pointer) + ((mpointer weight) :pointer) + ((mpointer observations) :pointer) + ((mpointer parameters) :pointer) + ((mpointer covariance) :pointer) (chisq :double) (workspace :pointer)) :documentation ; FDL "Compute the best-fit parameters c of the weighted @@ -217,12 +217,12 @@ (defmfun weighted-linear-mfit-svd (model weight observations parameters covariance tolerance workspace) "gsl_multifit_wlinear_svd" - (((pointer model) :pointer) - ((pointer weight) :pointer) - ((pointer observations) :pointer) + (((mpointer model) :pointer) + ((mpointer weight) :pointer) + ((mpointer observations) :pointer) (tolerance :double) - (rank size) - ((pointer parameters) :pointer) (covariance :pointer) (chisq :double) + (rank sizet) + ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double) (workspace :pointer)) :return ((dcref chisq) (scref rank)) :documentation ; FDL @@ -243,8 +243,8 @@ (defmfun multi-linear-estimate (x coefficients covariance) "gsl_multifit_linear_est" - (((pointer x) :pointer) ((pointer coefficients) :pointer) - ((pointer covariance) :pointer) (y :double) (y-error :double)) + (((mpointer x) :pointer) ((mpointer coefficients) :pointer) + ((mpointer covariance) :pointer) (y :double) (y-error :double)) :documentation ; FDL "Use the best-fit multilinear regression coefficients and their covariance matrix to compute the fitted function value @@ -258,9 +258,9 @@ (defun univariate-linear-least-squares-example () "First example in Section 36.5 of the GSL manual." ;; Results not given in manual so not verified yet. - (letm ((x (vector-double-float #(1970.0d0 1980.0d0 1990.0d0 2000.0d0))) - (y (vector-double-float #(12.0d0 11.0d0 14.0d0 13.0d0))) - (w (vector-double-float #(0.1d0 0.2d0 0.3d0 0.4d0)))) + (letm ((x (vector-double-float (a 1970.0d0 1980.0d0 1990.0d0 2000.0d0))) + (y (vector-double-float (a 12.0d0 11.0d0 14.0d0 13.0d0))) + (w (vector-double-float (a 0.1d0 0.2d0 0.3d0 0.4d0)))) (multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq) (weighted-linear-fit x w y) (format t "~&Best fit: Y = ~8,5f + ~8,5f X" c0 c1) @@ -301,8 +301,8 @@ (defun mv-linear-least-squares-example (data) "Second example in Section 36.5 of the GSL manual." (letm ((n (length data)) chisq - (x (matrix-double-float n 3)) - (cov (matrix-double-float 3 3)) + (x (matrix-double-float (list n 3))) + (cov (matrix-double-float '(3 3))) (y (vector-double-float n)) (w (vector-double-float n)) (c (vector-double-float 3))) diff --git a/solve-minimize-fit/minimization-multi.lisp b/solve-minimize-fit/minimization-multi.lisp index 7130d4824b96a125f5cef01a7624a15e90685382..29f52c50dfd7ae0c996cad42eb73b63de5d2b971 100644 --- a/solve-minimize-fit/minimization-multi.lisp +++ b/solve-minimize-fit/minimization-multi.lisp @@ -1,6 +1,6 @@ ;; Multivariate minimization. ;; Liam Healy <Tue Jan 8 2008 - 21:28> -;; Time-stamp: <2008-03-23 17:50:49EDT minimization-multi.lisp> +;; Time-stamp: <2008-08-31 12:25:59EDT minimization-multi.lisp> ;; $Id$ (in-package :gsl) @@ -56,7 +56,7 @@ (defmfun allocate-mfminimizer (type dimension) "gsl_multimin_fminimizer_alloc" - ((type :pointer) (dimension size)) + ((type :pointer) (dimension sizet)) :c-return :pointer :export nil :index (letm mfminimizer) @@ -66,7 +66,7 @@ (defmfun allocate-mfdfminimizer (type dimension) "gsl_multimin_fdfminimizer_alloc" - ((type :pointer) (dimension size)) + ((type :pointer) (dimension sizet)) :c-return :pointer :export nil :index (letm mfdfminimizer) @@ -77,7 +77,7 @@ (defmfun set-mfminimizer (minimizer function initial step-size) "gsl_multimin_fminimizer_set" ((minimizer :pointer) (function :pointer) - ((pointer initial) :pointer) ((pointer step-size) :pointer)) + ((mpointer initial) :pointer) ((mpointer step-size) :pointer)) :export nil :index (letm mfminimizer) :documentation ; FDL @@ -91,7 +91,7 @@ (minimizer function-derivative initial step-size tolerance) "gsl_multimin_fdfminimizer_set" ((minimizer :pointer) (function-derivative :pointer) - ((pointer initial) :pointer) (step-size :double) + ((mpointer initial) :pointer) (step-size :double) (tolerance :double)) :export nil :index (letm mfdfminimizer) @@ -160,16 +160,16 @@ (defmfun mfminimizer-x (minimizer) "gsl_multimin_fminimizer_x" ((minimizer :pointer)) - :c-return (canswer :pointer) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The current best estimate of the location of the minimum.") (defmfun mfdfminimizer-x (minimizer) "gsl_multimin_fdfminimizer_x" ((minimizer :pointer)) - :c-return (canswer :pointer) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The current best estimate of the location of the minimum.") @@ -197,8 +197,8 @@ (defmfun mfdfminimizer-gradient (minimizer) "gsl_multimin_fdfminimizer_gradient" ((minimizer :pointer)) - :c-return (canswer :pointer) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The current best estimate of the gradient for the minimizer.") @@ -215,7 +215,7 @@ (defmfun min-test-gradient (gradient absolute-error) "gsl_multimin_test_gradient" - (((pointer gradient) :pointer) (absolute-error :double)) + ((gradient :pointer) (absolute-error :double)) :c-return :success-continue :documentation ; FDL "Test the norm of the gradient against the @@ -316,8 +316,8 @@ (defun parabaloid (gsl-vector-pointer) "A parabaloid function of two arguments, given in GSL manual Sec. 35.4." - (let ((x (vref gsl-vector-pointer 0)) - (y (vref gsl-vector-pointer 1)) + (let ((x (maref gsl-vector-pointer 0)) + (y (maref gsl-vector-pointer 1)) (dp0 (aref *parabaloid-center* 0)) (dp1 (aref *parabaloid-center* 1))) (+ (* 10 (expt (- x dp0) 2)) @@ -326,13 +326,13 @@ (defun parabaloid-derivative (arguments-gv-pointer derivative-gv-pointer) - (let ((x (vref arguments-gv-pointer 0)) - (y (vref arguments-gv-pointer 1)) + (let ((x (maref arguments-gv-pointer 0)) + (y (maref arguments-gv-pointer 1)) (dp0 (aref *parabaloid-center* 0)) (dp1 (aref *parabaloid-center* 1))) - (setf (vref derivative-gv-pointer 0) + (setf (maref derivative-gv-pointer 0) (* 20 (- x dp0)) - (vref derivative-gv-pointer 1) + (maref derivative-gv-pointer 1) (* 40 (- y dp1))))) (defun parabaloid-and-derivative @@ -345,7 +345,7 @@ parabaloid 2 parabaloid-derivative parabaloid-and-derivative) (defun multimin-example-fletcher-reeves () - (letm ((initial (vector-double-float #(5.0d0 7.0d0))) + (letm ((initial (vector-double-float (a 5.0d0 7.0d0))) (minimizer (mfdfminimizer *conjugate-fletcher-reeves* 2 parabaloid initial 0.01d0 1.0d-4))) @@ -375,7 +375,7 @@ (def-minimization-functions parabaloid-f 2) (defun multimin-example-nelder-mead () - (letm ((initial (vector-double-float #(5.0d0 7.0d0))) + (letm ((initial (vector-double-float (a 5.0d0 7.0d0))) (step-size (vector-double-float 2))) (set-all step-size 1.0d0) (letm ((minimizer diff --git a/solve-minimize-fit/nonlinear-least-squares.lisp b/solve-minimize-fit/nonlinear-least-squares.lisp index 03ed2e094d5d4f188261563df99f03c8206199b6..366a3722a3939482180ac2427329a9f02ebe3f5d 100644 --- a/solve-minimize-fit/nonlinear-least-squares.lisp +++ b/solve-minimize-fit/nonlinear-least-squares.lisp @@ -1,6 +1,6 @@ ;; Nonlinear least squares fitting. ;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp -;; Time-stamp: <2008-03-09 19:30:55EDT nonlinear-least-squares.lisp> +;; Time-stamp: <2008-08-31 15:13:29EDT nonlinear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -18,7 +18,7 @@ (defmfun allocate-ffit (solver-type number-of-observations number-of-parameters) "gsl_multifit_fsolver_alloc" - ((solver-type :pointer) (number-of-observations size) (number-of-parameters size)) + ((solver-type :pointer) (number-of-observations sizet) (number-of-parameters sizet)) :c-return :pointer :export nil :index (letm nonlinear-ffit) @@ -59,7 +59,7 @@ (defmfun allocate-fdffit (solver-type number-of-observations number-of-parameters) "gsl_multifit_fdfsolver_alloc" - ((solver-type :pointer) (number-of-observations size) (number-of-parameters size)) + ((solver-type :pointer) (number-of-observations sizet) (number-of-parameters sizet)) :c-return :pointer :export nil :index (letm nonlinear-fdffit) @@ -69,7 +69,7 @@ (defmfun set-fdffit (solver function initial-guess) "gsl_multifit_fdfsolver_set" - ((solver :pointer) (function :pointer) ((pointer initial-guess) :pointer)) + ((solver :pointer) (function :pointer) ((mpointer initial-guess) :pointer)) :documentation ; FDL "Initialize or reinitialize an existing solver to use the function and the initial guess.") @@ -114,7 +114,7 @@ ;; See /usr/include/gsl/gsl_multifit_nlin.h "The definition of a function for nonlinear least squares fitting in GSL." (function :pointer) - (dimensions size) + (dimensions sizet) (parameters :pointer)) (cffi:defcstruct gsl-fdffit-function @@ -124,8 +124,8 @@ (function :pointer) (df :pointer) (fdf :pointer) - (number-of-observations size) - (number-of-parameters size) + (number-of-observations sizet) + (number-of-parameters sizet) (parameters :pointer)) (export 'def-fitting-functions) @@ -170,16 +170,16 @@ (defmfun current-parameters-ffit (solver) "gsl_multifit_fsolver_position" ((solver :pointer)) - :c-return (canswer :pointer) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The current best-fit parameters.") (defmfun current-parameters-fdffit (solver) "gsl_multifit_fdfsolver_position" ((solver :pointer)) - :c-return (canswer :pointer) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The current best-fit parameters.") @@ -216,7 +216,7 @@ (defmfun fit-gradient (jacobian function-values gradient) "gsl_multifit_gradient" - ((jacobian :pointer) ((pointer function-values) :pointer) (gradient :pointer)) + ((jacobian :pointer) ((mpointer function-values) :pointer) (gradient :pointer)) :documentation ; FDL "Compute the gradient of \Phi(x) = (1/2) ||F(x)||^2 from the Jacobian matrix and the function values using @@ -279,7 +279,7 @@ (defmfun ls-covariance (jacobian relative-error covariance) "gsl_multifit_covar" - ((jacobian :pointer) (relative-error :double) ((pointer covariance) :pointer)) + ((jacobian :pointer) (relative-error :double) ((mpointer covariance) :pointer)) :return (covariance) :documentation ; FDL "Compute the covariance matrix of the best-fit parameters @@ -333,25 +333,25 @@ (make-exponent-fit-data :n *number-of-observations* :y - (let ((arr (make-array *number-of-observations* :element-type 'double-float))) + (let ((arr (make-array* *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* :element-type 'double-float :initial-element 0.1d0)))) + (make-array* *number-of-observations* 'double-float :initial-element 0.1d0)))) (defun exponential-residual (x f) "Compute the negative of the residuals with the exponential model for the nonlinear least squares example." - (let ((A (vref x 0)) - (lambda (vref x 1)) - (b (vref x 2))) + (let ((A (maref x 0)) + (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 (vref f i) + (setf (maref f i) ;; the difference model - observation = - residual (/ (- (+ (* A (exp (* (- lambda) i))) b) (aref y i)) (aref sigma i))))))) @@ -360,16 +360,16 @@ "Compute the partial derivatives of the negative of the residuals with the exponential model for the nonlinear least squares example." - (let ((A (vref x 0)) - (lambda (vref x 1))) + (let ((A (maref x 0)) + (lambda (maref x 1))) (symbol-macrolet ((sigma (exponent-fit-data-sigma *nlls-example-data*))) (dotimes (i *number-of-observations*) (let ((e (exp (* (- lambda) i))) (s (aref sigma i))) - (setf (mref jacobian i 0) (/ e s) - (mref jacobian i 1) (* -1 i A (/ e s)) - (mref jacobian i 2) (/ s))))))) + (setf (maref jacobian i 0) (/ e s) + (maref jacobian i 1) (* -1 i A (/ e s)) + (maref jacobian i 2) (/ s))))))) (defun exponential-residual-fdf (x f jacobian) "Compute the function and partial derivatives of the negative of the @@ -382,42 +382,47 @@ exponential-residual *number-of-observations* *number-of-parameters* exponential-residual-derivative exponential-residual-fdf) +(defun norm-f (fit) + "Find the norm of the fit function f." + (letm ((arr (vector-double-float (cl-array (fdffit-slot fit 'f))))) + (euclidean-norm arr))) + (defun solve-nonlinear-least-squares-example () - (letm ((init (vector-double-float #(1.0d0 0.0d0 0.0d0))) + (letm ((init (vector-double-float (a 1.0d0 0.0d0 0.0d0))) (covariance - (matrix-double-float *number-of-parameters* *number-of-parameters*)) + (matrix-double-float (list *number-of-parameters* *number-of-parameters*))) (fit (nonlinear-fdffit *levenberg-marquardt* *number-of-observations* *number-of-parameters* exponential-residual init))) - (macrolet ((fitx (i) `(vref (fdffit-slot fit 'x) ,i)) + (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" 0 (fitx 0) (fitx 1) (fitx 2) - (norm (make-data-from-pointer (fdffit-slot fit 'f)))) + (norm-f fit)) (loop for iter from 0 below 25 - until - (and (plusp iter) - (fit-test-delta - (fdffit-slot fit 'dx) (fdffit-slot fit 'x) - 1.0d-4 1.0d-4)) - do - (iterate-fdffit fit) - (ls-covariance (fdffit-slot fit 'jacobian) 0.0d0 covariance) - (format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" - (1+ iter) (fitx 0) (fitx 1) (fitx 2) - (norm (make-data-from-pointer (fdffit-slot fit 'f)))) - finally - (let* ((chi (norm (make-data-from-pointer (fdffit-slot fit 'f)))) - (dof (- *number-of-observations* *number-of-parameters*)) - (c (max 1.0d0 (/ chi (sqrt dof))))) - (format t "~&chisq/dof = ~g" (/ (expt chi 2) dof)) - (format t "~&A = ~,5f +/- ~,5f" (fitx 0) (* c (err 0))) - (format t "~&lambda = ~,5f +/- ~,5f" (fitx 1) (* c (err 1))) - (format t "~&b = ~,5f +/- ~,5f" (fitx 2) (* c (err 2))) - (return (list (fitx 0) (fitx 1) (fitx 2)))))))) + until + (and (plusp iter) + (fit-test-delta + (fdffit-slot fit 'dx) (fdffit-slot fit 'x) + 1.0d-4 1.0d-4)) + do + (iterate-fdffit fit) + (ls-covariance (fdffit-slot fit 'jacobian) 0.0d0 covariance) + (format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" + (1+ iter) (fitx 0) (fitx 1) (fitx 2) + (norm-f fit)) + finally + (let* ((chi (norm-f fit)) + (dof (- *number-of-observations* *number-of-parameters*)) + (c (max 1.0d0 (/ chi (sqrt dof))))) + (format t "~&chisq/dof = ~g" (/ (expt chi 2) dof)) + (format t "~&A = ~,5f +/- ~,5f" (fitx 0) (* c (err 0))) + (format t "~&lambda = ~,5f +/- ~,5f" (fitx 1) (* c (err 1))) + (format t "~&b = ~,5f +/- ~,5f" (fitx 2) (* c (err 2))) + (return (list (fitx 0) (fitx 1) (fitx 2)))))))) ;;; Run example: ;;; (nlls-setup) diff --git a/solve-minimize-fit/roots-multi.lisp b/solve-minimize-fit/roots-multi.lisp index 0c59a5b9f1ae136b6a719781b0b658c8ce1ac2a2..058772daa8c0f7b7f91f9886c2487e7876c900eb 100644 --- a/solve-minimize-fit/roots-multi.lisp +++ b/solve-minimize-fit/roots-multi.lisp @@ -1,16 +1,10 @@ ;;; Multivariate roots. ;;; Liam Healy 2008-01-12 12:49:08 -;;; Time-stamp: <2008-03-09 19:29:20EDT roots-multi.lisp> +;;; Time-stamp: <2008-08-23 23:04:50EDT roots-multi.lisp> ;;; $Id$ (in-package :gsl) -;;; I don't like using make-data-from-pointer, but it's the only way -;;; to have access to the GSL functions when given a pointer. -;;; Alternatively, I could provide the GSL pointer and then the only -;;; thing the user could use is vref, or, of course -;;; make-data-from-pointer. - ;;;;**************************************************************************** ;;;; Function definition ;;;;**************************************************************************** @@ -19,7 +13,7 @@ ;; See /usr/include/gsl/gsl_multiroots.h "The definition of a function for multiroot finding in GSL." (function :pointer) - (dimensions size) + (dimensions sizet) (parameters :pointer)) (export 'def-mfunction) @@ -36,7 +30,7 @@ (function :pointer) (df :pointer) (fdf :pointer) - (dimensions size) + (dimensions sizet) (parameters :pointer)) ;;;;**************************************************************************** @@ -59,7 +53,7 @@ (defmfun allocate-mfsolver (type dimension) "gsl_multiroot_fsolver_alloc" - ((type :pointer) (dimension size)) + ((type :pointer) (dimension sizet)) :c-return :pointer :export nil :index (letm mfsolver) @@ -69,7 +63,7 @@ (defmfun allocate-mfdfsolver (type dimension) "gsl_multiroot_fdfsolver_alloc" - ((type :pointer) (dimension size)) + ((type :pointer) (dimension sizet)) :c-return :pointer :export nil :index (letm mfdfsolver) @@ -79,7 +73,7 @@ (defmfun set-mfsolver (solver function initial) "gsl_multiroot_fsolver_set" - ((solver :pointer) (function :pointer) ((pointer initial) :pointer)) + ((solver :pointer) (function :pointer) ((mpointer initial) :pointer)) :export nil :index (letm mfsolver) :documentation ; FDL @@ -89,7 +83,7 @@ (defmfun set-mfdfsolver (solver function-derivative initial) "gsl_multiroot_fdfsolver_set" ((solver :pointer) (function-derivative :pointer) - ((pointer initial) :pointer)) + ((mpointer initial) :pointer)) :export nil :index (letm mfdfsolver) :documentation ; FDL @@ -157,48 +151,48 @@ (defmfun mfsolver-root (solver) "gsl_multiroot_fsolver_root" ((solver :pointer)) - :c-return (canswer :pointer) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The current estimate of the root for the solver.") (defmfun mfdfsolver-root (solver) "gsl_multiroot_fdfsolver_root" ((solver :pointer)) - :c-return (canswer gsl-vector-c) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation "The current estimate of the root for the solver.") (defmfun mfsolver-f (solver) "gsl_multiroot_fsolver_f" ((solver :pointer)) - :c-return (canswer gsl-vector-c) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :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" ((solver :pointer)) - :c-return (canswer gsl-vector-c) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :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" ((solver :pointer)) - :c-return (canswer gsl-vector-c) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The last step dx taken by the solver.") (defmfun mfdfsolver-dx (solver) "gsl_multiroot_fsolver_dx" ((solver :pointer)) - :c-return (canswer gsl-vector-c) - :return ((make-data-from-pointer canswer)) + :c-return :pointer + :return (:c-return) :documentation ; FDL "The last step dx taken by the solver.") @@ -419,44 +413,45 @@ ;;; The recommended alternative (defun rosenbrock (argument return) "Rosenbrock test function." - (setf (vref return 0) - (* *rosenbrock-a* (- 1 (vref argument 0))) - (vref return 1) - (* *rosenbrock-b* (- (vref argument 1) (expt (vref argument 0) 2))))) + (setf (maref return 0) + (* *rosenbrock-a* (- 1 (maref argument 0))) + (maref return 1) + (* *rosenbrock-b* (- (maref argument 1) (expt (maref argument 0) 2))))) (def-mfunction rosenbrock 2) (defun roots-multi-example () "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual." (let ((max-iter 1000)) - (letm ((vect (vector-double-float #(-10.0d0 -5.0d0)))) - (letm ((solver (mfsolver *hybrid-scaled* rosenbrock vect))) - (let ((fnval (mfsolver-f solver)) - (argval (mfsolver-root solver))) - (loop for iter from 0 - while (and (< iter max-iter) - (not (multiroot-test-residual solver 1.0d-7))) - do - (iterate-mfsolver solver) - (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)) - finally (return - (values (maref argval 0) - (maref argval 1) - (maref fnval 0) - (maref fnval 1))))))))) - + (letm ((vect (vector-double-float (a -10.0d0 -5.0d0))) + (solver (mfsolver *hybrid-scaled* rosenbrock vect))) + (loop for iter from 0 + with fnval and argval + while (and (< iter max-iter) + (not (multiroot-test-residual solver 1.0d-7))) + + do + (iterate-mfsolver solver) + (setf fnval (cl-array (mfsolver-f solver)) + argval (cl-array (mfsolver-root solver))) + (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g" + iter + (aref argval 0) + (aref argval 1) + (aref fnval 0) + (aref fnval 1)) + finally (return + (values (aref argval 0) + (aref argval 1) + (aref fnval 0) + (aref fnval 1))))))) (defun rosenbrock-df (argument jacobian) "The partial derivatives of the Rosenbrock functions." - (setf (mref jacobian 0 0) (- *rosenbrock-a*) - (mref jacobian 0 1) 0.0d0 - (mref jacobian 1 0) (* -2 *rosenbrock-b* (vref argument 0)) - (mref jacobian 1 1) *rosenbrock-b*)) + (setf (maref jacobian 0 0) (- *rosenbrock-a*) + (maref jacobian 0 1) 0.0d0 + (maref jacobian 1 0) (* -2 *rosenbrock-b* (maref argument 0)) + (maref jacobian 1 1) *rosenbrock-b*)) (defun rosenbrock-fdf (argument value jacobian) (rosenbrock argument value) @@ -481,20 +476,22 @@ (maref fnval 0) (maref fnval 1)))) (let ((max-iter 1000)) - (letm ((vect (vector-double-float #(-10.0d0 -5.0d0)))) + (letm ((vect (vector-double-float (a -10.0d0 -5.0d0)))) (letm ((solver (mfdfsolver *gnewton-mfdfsolver* rosenbrock-f vect))) - (let ((fnval (mfdfsolver-f solver)) - (argval (mfdfsolver-root solver))) - (loop for iter from 0 - while (and (< iter max-iter) - (not (multiroot-test-residual solver 1.0d-7))) - initially (print-state iter argval fnval) - do - (iterate-mfdfsolver solver) - (print-state iter argval fnval) - finally (return - (values (maref argval 0) - (maref argval 1) - (maref fnval 0) - (maref fnval 1)))))))))) + (loop for iter from 0 + with fnval = (mfdfsolver-f solver) + and argval = (mfdfsolver-root solver) + while (and (< iter max-iter) + (not (multiroot-test-residual solver 1.0d-7))) + initially (print-state iter argval fnval) + do + (iterate-mfdfsolver solver) + (setf fnval (mfdfsolver-f solver) + argval (mfdfsolver-root solver)) + (print-state iter argval fnval) + finally (return + (values (maref argval 0) + (maref argval 1) + (maref fnval 0) + (maref fnval 1)))))))))