diff --git a/documentation/status.text b/documentation/status.text index 3e158ac1bab424600383dd0a8ea7369d6d5cdcbc..49617c73d981554db3e582cd15ce43ceb353c12d 100644 --- a/documentation/status.text +++ b/documentation/status.text @@ -28,7 +28,9 @@ This gives an errorneous error value (4th returned value) when run from the shell; run inside SLIME gives the correct answer. === Other known bugs === - +Fail in CCL but not in SBCL: +(mv-linear-least-squares-example (mv-linear-least-squares-data)) +(multimin-example-fletcher-reeves) === GSL questions === diff --git a/gsll.asd b/gsll.asd index 64163579c62bf05c6b8a2889afdff3205f990855..ce6f4437c1fb919b6f11e79d4328d2e401d7c32d 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2008-12-31 21:10:22EST gsll.asd> +;; Time-stamp: <2009-01-03 14:10:02EST gsll.asd> ;; $Id$ (asdf:defsystem "gsll" @@ -179,10 +179,11 @@ (:module solve-minimize-fit :depends-on (init data random) :components - ((:file "roots-one") - (:file "minimization-one") - (:file "roots-multi" :depends-on (roots-one)) - (:file "minimization-multi") + ((:file "generic") + (:file "roots-one" :depends-on (generic)) + (:file "minimization-one" :depends-on (generic)) + (:file "roots-multi" :depends-on (roots-one generic)) + (:file "minimization-multi" :depends-on (generic)) (:file "linear-least-squares") - (:file "nonlinear-least-squares"))) + (:file "nonlinear-least-squares" :depends-on (generic)))) (:file "basis-splines" :depends-on (init data)))) diff --git a/solve-minimize-fit/generic.lisp b/solve-minimize-fit/generic.lisp new file mode 100644 index 0000000000000000000000000000000000000000..b887ea0aab071245b3ead0e9762088d37c77ecc4 --- /dev/null +++ b/solve-minimize-fit/generic.lisp @@ -0,0 +1,21 @@ +;; Generic functions for optimization +;; Liam Healy 2009-01-03 12:59:07EST generic.lisp +;; Time-stamp: <2009-01-03 13:14:47EST generic.lisp> +;; $Id: $ + +(in-package :gsl) + +(defgeneric iterate (object) + (:documentation "Take the next iteration step for this object.")) + +(defgeneric solution (object) + (:documentation + "The current value of the independent variable(s) that solves this object.")) + +(defgeneric function-value (object) + (:documentation + "The current value of the function that solves this object.")) + +(defgeneric last-step (object) + (:documentation ; FDL + "The last step dx taken by the solver.")) diff --git a/solve-minimize-fit/minimization-multi.lisp b/solve-minimize-fit/minimization-multi.lisp index f206fddbeb1f8f7d57aca6859eba8ec8cb2674cb..2be2669cc107947b32dab9faeff247eefd392e8b 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-12-30 10:17:15EST minimization-multi.lisp> +;; Time-stamp: <2009-01-03 13:27:12EST minimization-multi.lisp> ;; $Id$ (in-package :gsl) @@ -103,48 +103,54 @@ ;;;; Iteration ;;;;**************************************************************************** -(defmfun iterate-mfminimizer (minimizer) +(defmfun iterate ((minimizer multi-dimensional-minimizer-f)) "gsl_multimin_fminimizer_iterate" (((mpointer minimizer) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the minimizer. If the iteration encounters an unexpected problem then an error code will be returned.") -(defmfun iterate-mfdfminimizer (minimizer) +(defmfun iterate ((minimizer multi-dimensional-minimizer-fdf)) "gsl_multimin_fdfminimizer_iterate" (((mpointer minimizer) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the minimizer. If the iteration encounters an unexpected problem then an error code will be returned.") -(defmfun mfminimizer-x (minimizer) +(defmfun solution ((minimizer multi-dimensional-minimizer-f)) "gsl_multimin_fminimizer_x" (((mpointer minimizer) :pointer)) + :definition :method :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) +(defmfun solution ((minimizer multi-dimensional-minimizer-fdf)) "gsl_multimin_fdfminimizer_x" (((mpointer minimizer) :pointer)) + :definition :method :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 mfminimizer-minimum (minimizer) +(defmfun function-value ((minimizer multi-dimensional-minimizer-f)) "gsl_multimin_fminimizer_minimum" (((mpointer minimizer) :pointer)) + :definition :method :c-return :double :documentation ; FDL "The current best estimate of the value of the minimum.") -(defmfun mfdfminimizer-minimum (minimizer) +(defmfun function-value ((minimizer multi-dimensional-minimizer-fdf)) "gsl_multimin_fdfminimizer_minimum" (((mpointer minimizer) :pointer)) + :definition :method :c-return :double :documentation ; FDL "The current best estimate of the value of the minimum.") @@ -316,17 +322,17 @@ for iter from 0 below 100 while status do - (iterate-mfdfminimizer minimizer) + (iterate minimizer) (setf status (not (min-test-gradient (mfdfminimizer-gradient minimizer) 1.0d-3))) - (let ((x (mfdfminimizer-x minimizer))) + (let ((x (solution minimizer))) (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f" iter (maref x 0) (maref x 1) - (mfdfminimizer-minimum minimizer))) + (function-value minimizer))) finally (return - (let ((x (mfdfminimizer-x minimizer))) + (let ((x (solution minimizer))) (values (maref x 0) (maref x 1))))))) ;;; Because def-minimization-functions bind a symbol @@ -347,16 +353,16 @@ (loop with status = T and size for iter from 0 below 100 while status - do (iterate-mfminimizer minimizer) + do (iterate minimizer) (setf size (mfminimizer-size minimizer) status (not (min-test-size size 1.0d-2))) - (let ((x (mfminimizer-x minimizer))) + (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) - (mfminimizer-minimum minimizer) + (function-value minimizer) size)) finally (return - (let ((x (mfminimizer-x minimizer))) + (let ((x (solution minimizer))) (values (maref x 0) (maref x 1)))))))) diff --git a/solve-minimize-fit/minimization-one.lisp b/solve-minimize-fit/minimization-one.lisp index b29ca9ffb833cad4edafdf629edcbd24fa60c35c..1a6d60283cfd66f6f5fb145cef1c60b8f277ec63 100644 --- a/solve-minimize-fit/minimization-one.lisp +++ b/solve-minimize-fit/minimization-one.lisp @@ -1,6 +1,6 @@ ;; Univariate minimization ;; Liam Healy Tue Jan 8 2008 - 21:02 -;; Time-stamp: <2008-12-26 18:35:53EST minimization-one.lisp> +;; Time-stamp: <2009-01-03 15:42:59EST minimization-one.lisp> ;; $Id$ (in-package :gsl) @@ -47,9 +47,10 @@ ;;;; Iteration ;;;;**************************************************************************** -(defmfun iterate-fminimizer (minimizer) +(defmfun iterate ((minimizer one-dimensional-minimizer)) "gsl_min_fminimizer_iterate" (((mpointer minimizer) :pointer)) + :definition :method :c-return :success-continue :documentation ; FDL "Perform a single iteration of the minimizer. The following @@ -59,9 +60,10 @@ :FAILURE, the algorithm could not improve the current best approximation or bounding interval.") -(defmfun fminimizer-x-minimum (minimizer) +(defmfun solution ((minimizer one-dimensional-minimizer)) "gsl_min_fminimizer_x_minimum" (((mpointer minimizer) :pointer)) + :definition :method :c-return :double :documentation ; FDL "The current estimate of the position of the minimum for the minimizer.") @@ -80,9 +82,10 @@ :documentation ; FDL "The current upper bound of the interval for the minimizer.") -(defmfun fminimizer-f-minimum (minimizer) +(defmfun function-value ((minimizer one-dimensional-minimizer)) "gsl_min_fminimizer_f_minimum" (((mpointer minimizer) :pointer)) + :definition :method :c-return :double :documentation ; FDL "The value of the function at the current estimate of the minimum for the @@ -183,10 +186,10 @@ *brent-fminimizer* minimization-one-fn 2.0d0 0.0d0 6.0d0))) (format t "~&iter ~6t [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower") (loop for iter from 0 - for min = (fminimizer-x-minimum minimizer) + for min = (solution minimizer) for lower = (fminimizer-x-lower minimizer) for upper = (fminimizer-x-upper minimizer) - do (iterate-fminimizer minimizer) + do (iterate minimizer) while (and (< iter max-iter) ;; abs and rel error swapped in example? (not (min-test-interval lower upper 0.001d0 0.0d0))) diff --git a/solve-minimize-fit/nonlinear-least-squares.lisp b/solve-minimize-fit/nonlinear-least-squares.lisp index db14c26cbf773c8b3161bf5923b1b7238ff4e9d7..f77a390cf070a43a9cbeb283019e5d028331cb23 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-12-28 18:49:26EST nonlinear-least-squares.lisp> +;; Time-stamp: <2009-01-03 16:05:48EST nonlinear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -61,10 +61,6 @@ (dx :pointer) (state :pointer)) -(export '(fdffit-slot)) -(defun fdffit-slot (solver slot) - (cffi:foreign-slot-value solver 'gsl-fdffit-solver slot)) - ;;;;**************************************************************************** ;;;; The function to be minimized ;;;;**************************************************************************** @@ -112,36 +108,55 @@ ;;;; Iteration ;;;;**************************************************************************** -(defmfun iterate-ffit (solver) +(defmfun iterate ((solver nonlinear-ffit)) "gsl_multifit_fsolver_iterate" (((mpointer solver) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the solver. The solver maintains a current estimate of the best-fit parameters at all times. ") -(defmfun iterate-fdffit (solver) +(defmfun iterate ((solver nonlinear-fdffit)) "gsl_multifit_fdfsolver_iterate" (((mpointer solver) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the solver. The solver maintains a current estimate of the best-fit parameters at all times. ") -(defmfun current-parameters-ffit (solver) +(defmfun solution ((solver nonlinear-ffit)) "gsl_multifit_fsolver_position" (((mpointer solver) :pointer)) + :definition :method :c-return (crtn :pointer) :return ((make-marray 'double-float :from-pointer crtn)) :documentation ; FDL "The current best-fit parameters.") -(defmfun current-parameters-fdffit (solver) +(defmfun solution ((solver nonlinear-fdffit)) "gsl_multifit_fdfsolver_position" (((mpointer solver) :pointer)) + :definition :method :c-return (crtn :pointer) :return ((make-marray 'double-float :from-pointer crtn)) :documentation ; FDL "The current best-fit parameters.") +;;; Why doesn't GSL have functions to extract these values? +(defmethod function-value ((solver nonlinear-fdffit)) + (make-marray + 'double-float + :from-pointer + (cffi:foreign-slot-value (mpointer solver) 'gsl-fdffit-solver 'f))) + +(defmethod last-step ((solver nonlinear-fdffit)) + ;; Raw pointer, because we presume we're passing it on to another GSL function. + (cffi:foreign-slot-value (mpointer solver) 'gsl-fdffit-solver 'dx)) + +(defun jacobian (solver) + ;; Raw pointer, because we presume we're passing it on to another GSL function. + (cffi:foreign-slot-value (mpointer solver) 'gsl-fdffit-solver 'jacobian)) + ;;;;**************************************************************************** ;;;; Search stopping ;;;;**************************************************************************** @@ -346,7 +361,7 @@ (defun norm-f (fit) "Find the norm of the fit function f." - (euclidean-norm (make-marray 'double-float :from-pointer (fdffit-slot fit 'f)))) + (euclidean-norm (function-value fit))) (defun solve-nonlinear-least-squares-example () (let* ((init #m(1.0d0 0.0d0 0.0d0)) @@ -354,13 +369,13 @@ (make-marray 'double-float :dimensions (list *number-of-parameters* *number-of-parameters*))) - (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)) + (fit (make-nonlinear-fdffit + *levenberg-marquardt* + *number-of-observations* + *number-of-parameters* + exponential-residual + init))) + (macrolet ((fitx (i) `(maref (solution fit) ,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) @@ -368,12 +383,10 @@ (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)) + (fit-test-delta (last-step fit) (mpointer (solution fit)) 1.0d-4 1.0d-4)) do - (iterate-fdffit fit) - (ls-covariance (fdffit-slot fit 'jacobian) 0.0d0 covariance) + (iterate fit) + (ls-covariance (jacobian fit) 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)) diff --git a/solve-minimize-fit/roots-multi.lisp b/solve-minimize-fit/roots-multi.lisp index 5aea26a6e5ce248c19256c5f23767c55214a236b..f62b92f434b9685c5b9318ff754f1266c2c1293f 100644 --- a/solve-minimize-fit/roots-multi.lisp +++ b/solve-minimize-fit/roots-multi.lisp @@ -1,6 +1,6 @@ ;;; Multivariate roots. ;;; Liam Healy 2008-01-12 12:49:08 -;;; Time-stamp: <2008-12-30 09:57:11EST roots-multi.lisp> +;;; Time-stamp: <2009-01-03 16:07:02EST roots-multi.lisp> ;;; $Id$ (in-package :gsl) @@ -93,9 +93,10 @@ ;;;; Iteration ;;;;**************************************************************************** -(defmfun iterate-mfsolver (solver) +(defmfun iterate ((solver multi-dimensional-root-solver-f)) "gsl_multiroot_fsolver_iterate" (((mpointer solver) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the solver. The following errors may be signalled: 'bad-function-supplied, the iteration encountered a @@ -104,9 +105,10 @@ function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.") -(defmfun iterate-mfdfsolver (solver) +(defmfun iterate ((solver multi-dimensional-root-solver-fdf)) "gsl_multiroot_fdfsolver_iterate" (((mpointer solver) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the solver. The following errors may be signalled: 'bad-function-supplied, the iteration encountered a @@ -115,49 +117,55 @@ function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.") -(defmfun mfsolver-root (solver) +(defmfun solution ((solver multi-dimensional-root-solver-f)) "gsl_multiroot_fsolver_root" (((mpointer solver) :pointer)) + :definition :method :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) +(defmfun solution ((solver multi-dimensional-root-solver-fdf)) "gsl_multiroot_fdfsolver_root" (((mpointer solver) :pointer)) + :definition :method :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) +(defmfun function-value ((solver multi-dimensional-root-solver-f)) "gsl_multiroot_fsolver_f" (((mpointer solver) :pointer)) + :definition :method :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) +(defmfun function-value ((solver multi-dimensional-root-solver-fdf)) "gsl_multiroot_fdfsolver_f" (((mpointer solver) :pointer)) + :definition :method :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) +(defmfun last-step ((solver multi-dimensional-root-solver-f)) "gsl_multiroot_fsolver_dx" (((mpointer solver) :pointer)) + :definition :method :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) +(defmfun last-step ((solver multi-dimensional-root-solver-fdf)) "gsl_multiroot_fsolver_dx" (((mpointer solver) :pointer)) + :definition :method :c-return (crtn :pointer) :return ((make-marray 'double-float :from-pointer crtn)) :documentation ; FDL @@ -168,7 +176,7 @@ ;;;;**************************************************************************** ;;; The only place we need to pick apart the gsl_multiroot_fsolver -;;; struct is here. We could use mfsolver-dx etc., but then we'd have +;;; struct is here. We could use last-step etc., but then we'd have ;;; to discriminate on mfsolver vs. mfdfsolver. (cffi:defcstruct gsl-multiroot-fsolver ;; See /usr/include/gsl/gsl_multiroots.h @@ -379,7 +387,7 @@ (defun roots-multi-example () "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual." (let ((max-iter 1000)) - (let* ((vect #m(a -10.0d0 -5.0d0)) + (let* ((vect #m(-10.0d0 -5.0d0)) (solver (make-multi-dimensional-root-solver-f *hybrid-scaled* rosenbrock vect))) (loop for iter from 0 @@ -388,9 +396,9 @@ (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))) + (iterate solver) + (setf fnval (cl-array (function-value solver)) + argval (cl-array (solution solver))) (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g" iter (aref argval 0) @@ -433,19 +441,19 @@ (maref fnval 0) (maref fnval 1)))) (let ((max-iter 1000)) - (let* ((vect #m(a -10.0d0 -5.0d0)) + (let* ((vect #m(-10.0d0 -5.0d0)) (solver (make-multi-dimensional-root-solver-fdf *gnewton-mfdfsolver* rosenbrock-f vect))) (loop for iter from 0 - with fnval = (mfdfsolver-f solver) - and argval = (mfdfsolver-root solver) + with fnval = (function-value solver) + and argval = (solution 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)) + (iterate solver) + (setf fnval (function-value solver) + argval (solution solver)) (print-state iter argval fnval) finally (return (values (maref argval 0) diff --git a/solve-minimize-fit/roots-one.lisp b/solve-minimize-fit/roots-one.lisp index 9a093eb1436a56490a3427e6215e2e7b07d7f65c..3316bfef883ae88575563c01749bec429e767edf 100644 --- a/solve-minimize-fit/roots-one.lisp +++ b/solve-minimize-fit/roots-one.lisp @@ -1,6 +1,6 @@ ;; One-dimensional root solver. ;; Liam Healy -;; Time-stamp: <2008-12-30 10:17:14EST roots-one.lisp> +;; Time-stamp: <2009-01-03 13:08:27EST roots-one.lisp> ;; $Id$ (in-package :gsl) @@ -81,9 +81,10 @@ ;;;;**************************************************************************** ;; It appears that this is always returning :SUCCESS (0). -(defmfun iterate-fsolver (solver) +(defmfun iterate ((solver one-dimensional-root-solver-f)) "gsl_root_fsolver_iterate" (((mpointer solver) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the solver. The following errors may be signalled: 'bad-function-supplied, the iteration encountered a @@ -92,9 +93,10 @@ function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.") -(defmfun iterate-fdfsolver (solver) +(defmfun iterate ((solver one-dimensional-root-solver-fdf)) "gsl_root_fdfsolver_iterate" (((mpointer solver) :pointer)) + :definition :method :documentation ; FDL "Perform a single iteration of the solver. The following errors may be signalled: 'bad-function-supplied, the iteration encountered a @@ -103,16 +105,18 @@ function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.") -(defmfun fsolver-root (solver) +(defmfun solution ((solver one-dimensional-root-solver-f)) "gsl_root_fsolver_root" (((mpointer solver) :pointer)) + :definition :method :c-return :double :documentation ; FDL "The current estimate of the root for the solver.") -(defmfun fdfsolver-root (solver) +(defmfun solution ((solver one-dimensional-root-solver-fdf)) "gsl_root_fdfsolver_root" (((mpointer solver) :pointer)) + :definition :method :c-return :double :documentation ; FDL "The current estimate of the root for the solver.") @@ -324,10 +328,10 @@ *brent-fsolver* quadratic 0.0d0 5.0d0))) (format t "~&iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)") (loop for iter from 0 - for root = (fsolver-root solver) + for root = (solution solver) for lower = (fsolver-lower solver) for upper = (fsolver-upper solver) - do (iterate-fsolver solver) + do (iterate solver) while (and (< iter max-iter) (not (root-test-interval lower upper 0.0d0 0.001d0))) do @@ -353,7 +357,7 @@ (format t "~&iter ~6t ~8troot ~22terr ~34terr(est)") (loop for iter from 0 for oldroot = initial then root - for root = (progn (iterate-fdfsolver solver) (fdfsolver-root solver)) + for root = (progn (iterate solver) (solution solver)) while (and (< iter max-iter) (not (root-test-delta root oldroot 0.0d0 1.0d-5))) do