Commit f6c5e260 authored by Liam Healy's avatar Liam Healy

Generic functions for solve-minimize-fit

Introduced in generic.lisp the generic functions iterate, solution,
function-value, last-step.  These replace regular functions with
solver-specific names.
parent d0b63022
......@@ -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 ===
......
;; 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))))
;; 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."))
;; 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))))))))
;; 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)))
......
;; 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))
......
;;; 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)
......
;; 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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment