From 42cf5a7b2739e00ed18f0d315b502aee6df4bffb Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Sat, 24 Jan 2009 12:03:22 -0500 Subject: [PATCH] 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. --- solve-minimize-fit/generic.lisp | 34 +++- solve-minimize-fit/linear-least-squares.lisp | 60 ++++--- .../nonlinear-least-squares.lisp | 170 ++++++++++-------- solve-minimize-fit/roots-multi.lisp | 29 +-- tests/linear-least-squares.lisp | 4 +- tests/nonlinear-least-squares.lisp | 5 +- 6 files changed, 167 insertions(+), 135 deletions(-) diff --git a/solve-minimize-fit/generic.lisp b/solve-minimize-fit/generic.lisp index b887ea0a..8c33c5a7 100644 --- a/solve-minimize-fit/generic.lisp +++ b/solve-minimize-fit/generic.lisp @@ -1,10 +1,14 @@ ;; Generic functions for optimization ;; 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: $ (in-package :gsl) +;;;;**************************************************************************** +;;;; Generic functions for solve-minimize-fit objects +;;;;**************************************************************************** + (defgeneric iterate (object) (:documentation "Take the next iteration step for this object.")) @@ -19,3 +23,31 @@ (defgeneric last-step (object) (:documentation ; FDL "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)) + diff --git a/solve-minimize-fit/linear-least-squares.lisp b/solve-minimize-fit/linear-least-squares.lisp index c9881689..095a04ee 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: <2009-01-19 23:23:11EST linear-least-squares.lisp> +;; Time-stamp: <2009-01-24 09:24:25EST linear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -137,8 +137,27 @@ "multi-dimensional root solver with function only" "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 - (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" (((mpointer model) :pointer) ((mpointer observations) :pointer) (tolerance :double) @@ -161,7 +180,11 @@ zero singular value (to machine precision) are discarded from the fit.") (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" (((mpointer model) :pointer) ((mpointer observations) :pointer) (tolerance :double) @@ -189,22 +212,11 @@ tolerance tolerance, and the effective rank is returned as the second value.") -(defun size-array (array-or-size) - (if (numberp array-or-size) - array-or-size - (dim0 array-or-size))) - (defmfun weighted-linear-mfit (model weight observations parameters-or-size &optional - (covariance - (make-marray - 'double-float - :dimensions - (let ((s (size-array parameters-or-size))) (list s s)))) - (workspace - (make-fit-workspace - (dim0 observations) (size-array parameters-or-size))) + (covariance (default-covariance parameters-or-size)) + (workspace (default-lls-workspace observations parameters-or-size)) &aux (parameters (vdf parameters-or-size))) "gsl_multifit_wlinear" @@ -231,7 +243,11 @@ discarded from the fit.") (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" (((mpointer model) :pointer) ((mpointer weight) :pointer) @@ -274,7 +290,7 @@ ;;;; 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." ;; Results not given in manual so not verified yet. (let ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0)) @@ -318,7 +334,7 @@ collect (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 coefficients of x^0, x^1, x^2 for the best fit, and the chi squared." @@ -335,7 +351,7 @@ (maref w i) (/ (expt (third row) 2)))) (multiple-value-bind (parameters cov chisq) (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~&" (maref parameters 0) (maref parameters 1) (maref parameters 2)) (format t "Covariance matrix:~&") @@ -351,5 +367,5 @@ chisq)))) (save-test linear-least-squares - (univariate-linear-least-squares-example nil) - (mv-linear-least-squares-example (mv-linear-least-squares-data) nil)) + (linear-least-squares-univariate-example nil) + (linear-least-squares-multivariate-example (mv-linear-least-squares-data) nil)) diff --git a/solve-minimize-fit/nonlinear-least-squares.lisp b/solve-minimize-fit/nonlinear-least-squares.lisp index fca32989..a952572f 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: <2009-01-19 22:10:45EST nonlinear-least-squares.lisp> +;; Time-stamp: <2009-01-24 11:46:27EST nonlinear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -66,14 +66,16 @@ ;;;;**************************************************************************** (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." + ;; 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) - (dimensions sizet) + (number-of-observations sizet) + (number-of-parameters sizet) (parameters :pointer)) (cffi:defcstruct gsl-fdffit-function - ;; See "The definition of a function and its derivatives for nonlinear least squares fitting in GSL." (function :pointer) @@ -85,24 +87,35 @@ (export '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. The CL functions name and derivative should be defined previously 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 - (defmcallback ,function :success-failure :pointer :pointer) - ,@(when df - `((defmcallback ,df :success-failure :pointer :pointer) - (defmcallback ,fdf :success-failure :pointer (:pointer :pointer)))) - ,(if df - `(defcbstruct (,function function ,df df ,fdf fdf) - gsl-fdffit-function - ((number-of-observations ,number-of-observations) - (number-of-parameters ,number-of-parameters))) - `(defcbstruct (,function function) - gsl-ffit-function - ((number-of-observations ,number-of-observations) - (number-of-parameters ,number-of-parameters)))))) + (defmcallback ,function :success-failure :pointer :pointer) + ,@(when df + `((defmcallback ,df :success-failure :pointer :pointer) + (defmcallback ,fdf :success-failure :pointer (:pointer :pointer)))) + ,(if df + `(defcbstruct (,function function ,df df ,fdf fdf) + gsl-fdffit-function + ((number-of-observations ,number-of-observations) + (number-of-parameters ,number-of-parameters))) + `(defcbstruct (,function function) + gsl-ffit-function + ((number-of-observations ,number-of-observations) + (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 @@ -297,26 +310,22 @@ ;;; See the GSL source tree, doc/examples/expfit.c for the functions ;;; 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) (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." - (setf - *nlls-example-data* - (make-exponent-fit-data - :n *number-of-observations* - :y - (let ((arr (make-marray 'double-float :dimensions *number-of-observations*)) - (rng (make-random-number-generator *mt19937* 0))) - (dotimes (i *number-of-observations* arr) - (setf (maref arr i) - (+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0))))) - :sigma - (make-marray 'double-float :dimensions *number-of-observations* :initial-element 0.1d0)))) + (make-exponent-fit-data + :n number-of-observations + :y + (let ((arr (make-marray 'double-float :dimensions number-of-observations)) + (rng (make-random-number-generator *mt19937* 0))) + (dotimes (i number-of-observations arr) + (setf (maref arr i) + (+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0))))) + :sigma + (make-marray + 'double-float :dimensions number-of-observations :initial-element 0.1d0))) (defun exponential-residual (x f) "Compute the negative of the residuals with the exponential model @@ -327,7 +336,7 @@ (symbol-macrolet ((y (exponent-fit-data-y *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) ;; the difference model - observation = - residual (/ (- (+ (* A (exp (* (- lambda) i))) b) (maref y i)) @@ -341,7 +350,7 @@ (lambda (maref x 1))) (symbol-macrolet ((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))) (s (maref sigma i))) (setf (maref jacobian i 0) (/ e s) @@ -355,53 +364,58 @@ (exponential-residual x f) (exponential-residual-derivative x jacobian)) -(def-fitting-functions - exponential-residual *number-of-observations* *number-of-parameters* - exponential-residual-derivative exponential-residual-fdf) +(def-fitting-functions exponential-residual 3 + exponential-residual-derivative exponential-residual-fdf) (defun norm-f (fit) "Find the norm of the fit function f." (euclidean-norm (function-value fit))) -(defun solve-nonlinear-least-squares-example (&optional (print-steps t)) - (let* ((init #m(1.0d0 0.0d0 0.0d0)) - (covariance - (make-marray 'double-float - :dimensions - (list *number-of-parameters* *number-of-parameters*))) - (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)))) - (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) - (norm-f fit))) - (loop for iter from 0 below 25 - until - (and (plusp iter) - (fit-test-delta (last-step fit) (mpointer (solution fit)) 1.0d-4 1.0d-4)) - do - (iterate fit) - (ls-covariance (jacobian fit) 0.0d0 covariance) - (when print-steps - (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))))) +(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)) + (number-of-parameters (number-of-parameters exponential-residual)) + (covariance + (make-marray 'double-float + :dimensions + (list number-of-parameters number-of-parameters))) + (fit (make-nonlinear-fdffit + method + number-of-observations + number-of-parameters + exponential-residual + init))) + (macrolet ((fitx (i) `(maref (solution fit) ,i)) + (err (i) `(sqrt (maref 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) + (norm-f fit))) + (loop for iter from 0 below 25 + until + (and (plusp iter) + (fit-test-delta (last-step fit) (mpointer (solution fit)) 1.0d-4 1.0d-4)) + do + (iterate fit) + (ls-covariance (jacobian fit) 0.0d0 covariance) (when print-steps - (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)))))))) + (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))))) + (when print-steps + (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))))))))) (save-test nonlinear-least-squares - (progn (nlls-setup) (solve-nonlinear-least-squares-example nil))) + (nonlinear-least-squares-example 40 *levenberg-marquardt* nil)) diff --git a/solve-minimize-fit/roots-multi.lisp b/solve-minimize-fit/roots-multi.lisp index 8abe22e0..59efda2c 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: <2009-01-22 22:19:34EST roots-multi.lisp> +;;; Time-stamp: <2009-01-24 09:12:18EST roots-multi.lisp> ;;; $Id$ (in-package :gsl) @@ -13,33 +13,6 @@ ;;; marrays instead. This would allow directly manipulation of ;;; 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 ;;;;**************************************************************************** diff --git a/tests/linear-least-squares.lisp b/tests/linear-least-squares.lisp index e6b5f226..9d636a78 100644 --- a/tests/linear-least-squares.lisp +++ b/tests/linear-least-squares.lisp @@ -8,11 +8,11 @@ 39601.99999999999d0 -19.9d0 0.009999999999999998d0 0.8d0) (MULTIPLE-VALUE-LIST - (UNIVARIATE-LINEAR-LEAST-SQUARES-EXAMPLE NIL))) + (LINEAR-LEAST-SQUARES-UNIVARIATE-EXAMPLE NIL))) (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.1824632487186013d0 0.1845715900137661d0 1.3031038153096723d0 16.006137036426168d0) (MULTIPLE-VALUE-LIST - (MV-LINEAR-LEAST-SQUARES-EXAMPLE + (LINEAR-LEAST-SQUARES-MULTIVARIATE-EXAMPLE (MV-LINEAR-LEAST-SQUARES-DATA) NIL)))) diff --git a/tests/nonlinear-least-squares.lisp b/tests/nonlinear-least-squares.lisp index 0dcac436..67903142 100644 --- a/tests/nonlinear-least-squares.lisp +++ b/tests/nonlinear-least-squares.lisp @@ -8,7 +8,4 @@ (LIST 5.045357801443204d0 0.10404905892045835d0 1.0192487061031013d0)) (MULTIPLE-VALUE-LIST - (PROGN - (NLLS-SETUP) - (SOLVE-NONLINEAR-LEAST-SQUARES-EXAMPLE NIL))))) - + (NONLINEAR-LEAST-SQUARES-EXAMPLE 40 *levenberg-marquardt* NIL)))) -- GitLab