Commit 42cf5a7b authored by Liam Healy's avatar Liam Healy

Relocate definition, optional arguments, nonlinear least squares streamline

Moved gsl-mfunction, def-mfunction, and gsl-mfunction-fdf from
roots-multi.lisp to generic.lisp because they are common to several
files in the solve-minimize-fit module.  Made arguments default
instead of required in #'linear-mfit, linear-mfit-svd,
weighted-linear-mfit, weighted-linear-mfit-svd.  Most of these do not
have examples, so untested.  This required new functions
#'default-covariance, #'default-lls-workspace.  Renamed examples to
#'linear-least-squares-univariate-example and
#'linear-least-squares-multivariate-example, and changed to
'print-details.  Revised definitions in nonlinear least squares to
streamline data flow; there is no setup function or global
pre-declaration of the number of parameters or observations required
anymore.  The example function now called
nonlinear-least-squares-example will take the number of observations
and the fitter method as optional arguments.
parent 78b88074
;; Generic functions for optimization ;; Generic functions for optimization
;; Liam Healy 2009-01-03 12:59:07EST generic.lisp ;; Liam Healy 2009-01-03 12:59:07EST generic.lisp
;; Time-stamp: <2009-01-03 13:14:47EST generic.lisp> ;; Time-stamp: <2009-01-24 09:11:55EST generic.lisp>
;; $Id: $ ;; $Id: $
(in-package :gsl) (in-package :gsl)
;;;;****************************************************************************
;;;; Generic functions for solve-minimize-fit objects
;;;;****************************************************************************
(defgeneric iterate (object) (defgeneric iterate (object)
(:documentation "Take the next iteration step for this object.")) (:documentation "Take the next iteration step for this object."))
...@@ -19,3 +23,31 @@ ...@@ -19,3 +23,31 @@
(defgeneric last-step (object) (defgeneric last-step (object)
(:documentation ; FDL (:documentation ; FDL
"The last step dx taken by the solver.")) "The last step dx taken by the solver."))
;;;;****************************************************************************
;;;; Structures and macros for function definition
;;;;****************************************************************************
(cffi:defcstruct gsl-mfunction
;; See /usr/include/gsl/gsl_multiroots.h
"The definition of a function for multiroot finding in GSL."
(function :pointer)
(dimensions sizet)
(parameters :pointer))
(export 'def-mfunction)
(defmacro def-mfunction (name dimensions)
"Define a function for multivariate root solving."
`(def-single-function ,name :success-failure :double gsl-mfunction
,dimensions))
(cffi:defcstruct gsl-mfunction-fdf
;; See /usr/include/gsl/gsl_multiroots.h
"The definition of a function and its derivatives for multiroot
finding in GSL."
(function :pointer)
(df :pointer)
(fdf :pointer)
(dimensions sizet)
(parameters :pointer))
;; Linear least squares, or linear regression ;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp> ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
;; Time-stamp: <2009-01-19 23:23:11EST linear-least-squares.lisp> ;; Time-stamp: <2009-01-24 09:24:25EST linear-least-squares.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -137,8 +137,27 @@ ...@@ -137,8 +137,27 @@
"multi-dimensional root solver with function only" "multi-dimensional root solver with function only"
"Make a workspace for a multidimensional linear least-squares fit.") "Make a workspace for a multidimensional linear least-squares fit.")
(defun size-array (array-or-size)
(if (numberp array-or-size)
array-or-size
(dim0 array-or-size)))
(defun default-covariance (parameters-or-size)
(make-marray
'double-float
:dimensions
(let ((s (size-array parameters-or-size))) (list s s))))
(defun default-lls-workspace (observations parameters-or-size)
(make-fit-workspace
(dim0 observations) (size-array parameters-or-size)))
(defmfun linear-mfit (defmfun linear-mfit
(model observations parameters covariance tolerance workspace) (model observations parameters-or-size tolerance
&optional
(covariance (default-covariance parameters-or-size))
(workspace (default-lls-workspace observations parameters-or-size))
&aux (parameters (vdf parameters-or-size)))
"gsl_multifit_linear" "gsl_multifit_linear"
(((mpointer model) :pointer) ((mpointer observations) :pointer) (((mpointer model) :pointer) ((mpointer observations) :pointer)
(tolerance :double) (tolerance :double)
...@@ -161,7 +180,11 @@ ...@@ -161,7 +180,11 @@
zero singular value (to machine precision) are discarded from the fit.") zero singular value (to machine precision) are discarded from the fit.")
(defmfun linear-mfit-svd (defmfun linear-mfit-svd
(model observations parameters covariance tolerance workspace) (model observations parameters-or-size tolerance
&optional
(covariance (default-covariance parameters-or-size))
(workspace (default-lls-workspace observations parameters-or-size))
&aux (parameters (vdf parameters-or-size)))
"gsl_multifit_linear_svd" "gsl_multifit_linear_svd"
(((mpointer model) :pointer) ((mpointer observations) :pointer) (((mpointer model) :pointer) ((mpointer observations) :pointer)
(tolerance :double) (tolerance :double)
...@@ -189,22 +212,11 @@ ...@@ -189,22 +212,11 @@
tolerance tolerance, and the effective rank is returned as the tolerance tolerance, and the effective rank is returned as the
second value.") second value.")
(defun size-array (array-or-size)
(if (numberp array-or-size)
array-or-size
(dim0 array-or-size)))
(defmfun weighted-linear-mfit (defmfun weighted-linear-mfit
(model weight observations parameters-or-size (model weight observations parameters-or-size
&optional &optional
(covariance (covariance (default-covariance parameters-or-size))
(make-marray (workspace (default-lls-workspace observations parameters-or-size))
'double-float
:dimensions
(let ((s (size-array parameters-or-size))) (list s s))))
(workspace
(make-fit-workspace
(dim0 observations) (size-array parameters-or-size)))
&aux &aux
(parameters (vdf parameters-or-size))) (parameters (vdf parameters-or-size)))
"gsl_multifit_wlinear" "gsl_multifit_wlinear"
...@@ -231,7 +243,11 @@ ...@@ -231,7 +243,11 @@
discarded from the fit.") discarded from the fit.")
(defmfun weighted-linear-mfit-svd (defmfun weighted-linear-mfit-svd
(model weight observations parameters covariance tolerance workspace) (model weight observations parameters-or-size tolerance
&optional
(covariance (default-covariance parameters-or-size))
(workspace (default-lls-workspace observations parameters-or-size))
&aux (parameters (vdf parameters-or-size)))
"gsl_multifit_wlinear_svd" "gsl_multifit_wlinear_svd"
(((mpointer model) :pointer) (((mpointer model) :pointer)
((mpointer weight) :pointer) ((mpointer weight) :pointer)
...@@ -274,7 +290,7 @@ ...@@ -274,7 +290,7 @@
;;;; Examples ;;;; Examples
;;;;**************************************************************************** ;;;;****************************************************************************
(defun univariate-linear-least-squares-example (&optional (print-steps t)) (defun linear-least-squares-univariate-example (&optional (print-steps t))
"First example in Section 36.5 of the GSL manual." "First example in Section 36.5 of the GSL manual."
;; Results not given in manual so not verified yet. ;; Results not given in manual so not verified yet.
(let ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0)) (let ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0))
...@@ -318,7 +334,7 @@ ...@@ -318,7 +334,7 @@
collect collect
(list xd (+ y0 (gaussian rng sigma)) sigma)))) (list xd (+ y0 (gaussian rng sigma)) sigma))))
(defun mv-linear-least-squares-example (data &optional (print-steps t)) (defun linear-least-squares-multivariate-example (data &optional (print-details t))
"Second example in Section 36.5 of the GSL manual. Returns the "Second example in Section 36.5 of the GSL manual. Returns the
coefficients of x^0, x^1, x^2 for the best fit, and the chi coefficients of x^0, x^1, x^2 for the best fit, and the chi
squared." squared."
...@@ -335,7 +351,7 @@ ...@@ -335,7 +351,7 @@
(maref w i) (/ (expt (third row) 2)))) (maref w i) (/ (expt (third row) 2))))
(multiple-value-bind (parameters cov chisq) (multiple-value-bind (parameters cov chisq)
(weighted-linear-mfit X w y 3) (weighted-linear-mfit X w y 3)
(when print-steps (when print-details
(format t "Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2~&" (format t "Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2~&"
(maref parameters 0) (maref parameters 1) (maref parameters 2)) (maref parameters 0) (maref parameters 1) (maref parameters 2))
(format t "Covariance matrix:~&") (format t "Covariance matrix:~&")
...@@ -351,5 +367,5 @@ ...@@ -351,5 +367,5 @@
chisq)))) chisq))))
(save-test linear-least-squares (save-test linear-least-squares
(univariate-linear-least-squares-example nil) (linear-least-squares-univariate-example nil)
(mv-linear-least-squares-example (mv-linear-least-squares-data) nil)) (linear-least-squares-multivariate-example (mv-linear-least-squares-data) nil))
;; Nonlinear least squares fitting. ;; Nonlinear least squares fitting.
;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp ;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
;; Time-stamp: <2009-01-19 22:10:45EST nonlinear-least-squares.lisp> ;; Time-stamp: <2009-01-24 11:46:27EST nonlinear-least-squares.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -66,14 +66,16 @@ ...@@ -66,14 +66,16 @@
;;;;**************************************************************************** ;;;;****************************************************************************
(cffi:defcstruct gsl-ffit-function (cffi:defcstruct gsl-ffit-function
;; See /usr/include/gsl/gsl_multifit_nlin.h
"The definition of a function for nonlinear least squares fitting in GSL." "The definition of a function for nonlinear least squares fitting in GSL."
;; See gsl_multifit_function in /usr/include/gsl/gsl_multifit_nlin.h
;; Note that this is moot because GSL currently does not have any
;; derivative-free fitting functions.
(function :pointer) (function :pointer)
(dimensions sizet) (number-of-observations sizet)
(number-of-parameters sizet)
(parameters :pointer)) (parameters :pointer))
(cffi:defcstruct gsl-fdffit-function (cffi:defcstruct gsl-fdffit-function
;; See
"The definition of a function and its derivatives for nonlinear "The definition of a function and its derivatives for nonlinear
least squares fitting in GSL." least squares fitting in GSL."
(function :pointer) (function :pointer)
...@@ -85,10 +87,13 @@ ...@@ -85,10 +87,13 @@
(export 'def-fitting-functions) (export 'def-fitting-functions)
(defmacro def-fitting-functions (defmacro def-fitting-functions
(function number-of-observations number-of-parameters &optional df fdf) (function number-of-parameters &optional df fdf (number-of-observations 0))
"Setup functions for nonlinear least squares fitting. "Setup functions for nonlinear least squares fitting.
The CL functions name and derivative should be defined previously The CL functions name and derivative should be defined previously
with defuns." with defuns."
;; The number of observations doesn't have anything to do with the
;; fitting functions, only with the data. Therefore it defaults to
;; 0 and should be specified separately with the data.
`(progn `(progn
(defmcallback ,function :success-failure :pointer :pointer) (defmcallback ,function :success-failure :pointer :pointer)
,@(when df ,@(when df
...@@ -104,6 +109,14 @@ ...@@ -104,6 +109,14 @@
((number-of-observations ,number-of-observations) ((number-of-observations ,number-of-observations)
(number-of-parameters ,number-of-parameters)))))) (number-of-parameters ,number-of-parameters))))))
(export '(number-of-parameters number-of-observations))
(defun number-of-parameters (pointer)
(cffi:foreign-slot-value pointer 'gsl-fdffit-function 'number-of-parameters))
(defmacro number-of-observations (pointer)
`(cffi:foreign-slot-value ,pointer 'gsl-fdffit-function 'number-of-observations))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Iteration ;;;; Iteration
;;;;**************************************************************************** ;;;;****************************************************************************
...@@ -297,26 +310,22 @@ ...@@ -297,26 +310,22 @@
;;; See the GSL source tree, doc/examples/expfit.c for the functions ;;; See the GSL source tree, doc/examples/expfit.c for the functions
;;; and doc/examples/nlfit.c for the solver. ;;; and doc/examples/nlfit.c for the solver.
(defparameter *number-of-observations* 40)
(defparameter *number-of-parameters* 3)
(defstruct exponent-fit-data n y sigma) (defstruct exponent-fit-data n y sigma)
(defvar *nlls-example-data*) (defvar *nlls-example-data*)
(defun nlls-setup () (defun generate-nlls-data (&optional (number-of-observations 40))
"Create the data used in the nonlinear least squares fit example." "Create the data used in the nonlinear least squares fit example."
(setf
*nlls-example-data*
(make-exponent-fit-data (make-exponent-fit-data
:n *number-of-observations* :n number-of-observations
:y :y
(let ((arr (make-marray 'double-float :dimensions *number-of-observations*)) (let ((arr (make-marray 'double-float :dimensions number-of-observations))
(rng (make-random-number-generator *mt19937* 0))) (rng (make-random-number-generator *mt19937* 0)))
(dotimes (i *number-of-observations* arr) (dotimes (i number-of-observations arr)
(setf (maref arr i) (setf (maref arr i)
(+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0))))) (+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0)))))
:sigma :sigma
(make-marray 'double-float :dimensions *number-of-observations* :initial-element 0.1d0)))) (make-marray
'double-float :dimensions number-of-observations :initial-element 0.1d0)))
(defun exponential-residual (x f) (defun exponential-residual (x f)
"Compute the negative of the residuals with the exponential model "Compute the negative of the residuals with the exponential model
...@@ -327,7 +336,7 @@ ...@@ -327,7 +336,7 @@
(symbol-macrolet (symbol-macrolet
((y (exponent-fit-data-y *nlls-example-data*)) ((y (exponent-fit-data-y *nlls-example-data*))
(sigma (exponent-fit-data-sigma *nlls-example-data*))) (sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i *number-of-observations*) (dotimes (i (exponent-fit-data-n *nlls-example-data*))
(setf (maref f i) (setf (maref f i)
;; the difference model - observation = - residual ;; the difference model - observation = - residual
(/ (- (+ (* A (exp (* (- lambda) i))) b) (maref y i)) (/ (- (+ (* A (exp (* (- lambda) i))) b) (maref y i))
...@@ -341,7 +350,7 @@ ...@@ -341,7 +350,7 @@
(lambda (maref x 1))) (lambda (maref x 1)))
(symbol-macrolet (symbol-macrolet
((sigma (exponent-fit-data-sigma *nlls-example-data*))) ((sigma (exponent-fit-data-sigma *nlls-example-data*)))
(dotimes (i *number-of-observations*) (dotimes (i (exponent-fit-data-n *nlls-example-data*))
(let ((e (exp (* (- lambda) i))) (let ((e (exp (* (- lambda) i)))
(s (maref sigma i))) (s (maref sigma i)))
(setf (maref jacobian i 0) (/ e s) (setf (maref jacobian i 0) (/ e s)
...@@ -355,30 +364,35 @@ ...@@ -355,30 +364,35 @@
(exponential-residual x f) (exponential-residual x f)
(exponential-residual-derivative x jacobian)) (exponential-residual-derivative x jacobian))
(def-fitting-functions (def-fitting-functions exponential-residual 3
exponential-residual *number-of-observations* *number-of-parameters*
exponential-residual-derivative exponential-residual-fdf) exponential-residual-derivative exponential-residual-fdf)
(defun norm-f (fit) (defun norm-f (fit)
"Find the norm of the fit function f." "Find the norm of the fit function f."
(euclidean-norm (function-value fit))) (euclidean-norm (function-value fit)))
(defun solve-nonlinear-least-squares-example (&optional (print-steps t)) (defun nonlinear-least-squares-example
(&optional (number-of-observations 40)
(method *levenberg-marquardt*)
(print-steps t))
(let ((*nlls-example-data* (generate-nlls-data number-of-observations)))
(setf (number-of-observations exponential-residual) number-of-observations)
(let* ((init #m(1.0d0 0.0d0 0.0d0)) (let* ((init #m(1.0d0 0.0d0 0.0d0))
(number-of-parameters (number-of-parameters exponential-residual))
(covariance (covariance
(make-marray 'double-float (make-marray 'double-float
:dimensions :dimensions
(list *number-of-parameters* *number-of-parameters*))) (list number-of-parameters number-of-parameters)))
(fit (make-nonlinear-fdffit (fit (make-nonlinear-fdffit
*levenberg-marquardt* method
*number-of-observations* number-of-observations
*number-of-parameters* number-of-parameters
exponential-residual exponential-residual
init))) init)))
(macrolet ((fitx (i) `(maref (solution fit) ,i)) (macrolet ((fitx (i) `(maref (solution fit) ,i))
(err (i) `(sqrt (maref covariance ,i ,i)))) (err (i) `(sqrt (maref covariance ,i ,i))))
(when print-steps (when print-steps
(format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" (format t "iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g~&"
0 (fitx 0) (fitx 1) (fitx 2) 0 (fitx 0) (fitx 1) (fitx 2)
(norm-f fit))) (norm-f fit)))
(loop for iter from 0 below 25 (loop for iter from 0 below 25
...@@ -389,19 +403,19 @@ ...@@ -389,19 +403,19 @@
(iterate fit) (iterate fit)
(ls-covariance (jacobian fit) 0.0d0 covariance) (ls-covariance (jacobian fit) 0.0d0 covariance)
(when print-steps (when print-steps
(format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" (format t "iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g~&"
(1+ iter) (fitx 0) (fitx 1) (fitx 2) (1+ iter) (fitx 0) (fitx 1) (fitx 2)
(norm-f fit))) (norm-f fit)))
finally finally
(let* ((chi (norm-f fit)) (let* ((chi (norm-f fit))
(dof (- *number-of-observations* *number-of-parameters*)) (dof (- number-of-observations number-of-parameters))
(c (max 1.0d0 (/ chi (sqrt dof))))) (c (max 1.0d0 (/ chi (sqrt dof)))))
(when print-steps (when print-steps
(format t "~&chisq/dof = ~g" (/ (expt chi 2) dof)) (format t "chisq/dof = ~g~&" (/ (expt chi 2) dof))
(format t "~&A = ~,5f +/- ~,5f" (fitx 0) (* c (err 0))) (format t "A = ~,5f +/- ~,5f~&" (fitx 0) (* c (err 0)))
(format t "~&lambda = ~,5f +/- ~,5f" (fitx 1) (* c (err 1))) (format t "lambda = ~,5f +/- ~,5f~&" (fitx 1) (* c (err 1)))
(format t "~&b = ~,5f +/- ~,5f" (fitx 2) (* c (err 2)))) (format t "b = ~,5f +/- ~,5f~&" (fitx 2) (* c (err 2))))
(return (list (fitx 0) (fitx 1) (fitx 2)))))))) (return (list (fitx 0) (fitx 1) (fitx 2)))))))))
(save-test nonlinear-least-squares (save-test nonlinear-least-squares
(progn (nlls-setup) (solve-nonlinear-least-squares-example nil))) (nonlinear-least-squares-example 40 *levenberg-marquardt* nil))
;;; Multivariate roots. ;;; Multivariate roots.
;;; Liam Healy 2008-01-12 12:49:08 ;;; Liam Healy 2008-01-12 12:49:08
;;; Time-stamp: <2009-01-22 22:19:34EST roots-multi.lisp> ;;; Time-stamp: <2009-01-24 09:12:18EST roots-multi.lisp>
;;; $Id$ ;;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -13,33 +13,6 @@ ...@@ -13,33 +13,6 @@
;;; marrays instead. This would allow directly manipulation of ;;; marrays instead. This would allow directly manipulation of
;;; marrays by the user function. Notes Mon Jan 19 2009. ;;; marrays by the user function. Notes Mon Jan 19 2009.
;;;;****************************************************************************
;;;; Function definition
;;;;****************************************************************************
(cffi:defcstruct gsl-mfunction
;; See /usr/include/gsl/gsl_multiroots.h
"The definition of a function for multiroot finding in GSL."
(function :pointer)
(dimensions sizet)
(parameters :pointer))
(export 'def-mfunction)
(defmacro def-mfunction (name dimensions)
"Define a function for multivariate root solving."
`(def-single-function ,name :success-failure :double gsl-mfunction
,dimensions))
(cffi:defcstruct gsl-mfunction-fdf
;; See /usr/include/gsl/gsl_multiroots.h
"The definition of a function and its derivatives for multiroot
finding in GSL."
(function :pointer)
(df :pointer)
(fdf :pointer)
(dimensions sizet)
(parameters :pointer))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Initialization ;;;; Initialization
;;;;**************************************************************************** ;;;;****************************************************************************
......
...@@ -8,11 +8,11 @@ ...@@ -8,11 +8,11 @@
39601.99999999999d0 -19.9d0 39601.99999999999d0 -19.9d0
0.009999999999999998d0 0.8d0) 0.009999999999999998d0 0.8d0)
(MULTIPLE-VALUE-LIST (MULTIPLE-VALUE-LIST
(UNIVARIATE-LINEAR-LEAST-SQUARES-EXAMPLE NIL))) (LINEAR-LEAST-SQUARES-UNIVARIATE-EXAMPLE NIL)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 1.1824632487186013d0 0.1845715900137661d0 (LIST 1.1824632487186013d0 0.1845715900137661d0
1.3031038153096723d0 16.006137036426168d0) 1.3031038153096723d0 16.006137036426168d0)
(MULTIPLE-VALUE-LIST (MULTIPLE-VALUE-LIST
(MV-LINEAR-LEAST-SQUARES-EXAMPLE (LINEAR-LEAST-SQUARES-MULTIVARIATE-EXAMPLE
(MV-LINEAR-LEAST-SQUARES-DATA) NIL)))) (MV-LINEAR-LEAST-SQUARES-DATA) NIL))))
...@@ -8,7 +8,4 @@ ...@@ -8,7 +8,4 @@
(LIST 5.045357801443204d0 0.10404905892045835d0 (LIST 5.045357801443204d0 0.10404905892045835d0
1.0192487061031013d0)) 1.0192487061031013d0))
(MULTIPLE-VALUE-LIST (MULTIPLE-VALUE-LIST
(PROGN (NONLINEAR-LEAST-SQUARES-EXAMPLE 40 *levenberg-marquardt* NIL))))
(NLLS-SETUP)
(SOLVE-NONLINEAR-LEAST-SQUARES-EXAMPLE NIL)))))
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