Commit 1082a0f5 authored by Liam Healy's avatar Liam Healy

Port remainder of solve-minimize-fit to ffa

Tests work and duplicate results from pre-ffa results.  This puts
completes the initial phase of the ffa port.
parent c6fa7363
;; Definition of GSLL system ;; Definition of GSLL system
;; Liam Healy ;; Liam Healy
;; Time-stamp: <2008-08-31 10:05:43EDT gsll.asd> ;; Time-stamp: <2008-08-31 12:49:38EDT gsll.asd>
;; $Id$ ;; $Id$
(asdf:defsystem "gsll" (asdf:defsystem "gsll"
...@@ -181,9 +181,8 @@ ...@@ -181,9 +181,8 @@
:components :components
((:file "roots-one") ((:file "roots-one")
(:file "minimization-one") (:file "minimization-one")
;(:file "roots-multi" :depends-on (roots-one)) (:file "roots-multi" :depends-on (roots-one))
;(:file "minimization-multi") (:file "minimization-multi")
;(:file "linear-least-squares") (:file "linear-least-squares")
;(:file "nonlinear-least-squares") (:file "nonlinear-least-squares")))
))
(:file "basis-splines" :depends-on (init data)))) (:file "basis-splines" :depends-on (init data))))
;; 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: <2008-03-09 19:30:56EDT linear-least-squares.lisp> ;; Time-stamp: <2008-08-31 12:44:33EDT linear-least-squares.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -17,9 +17,9 @@ ...@@ -17,9 +17,9 @@
(defmfun linear-fit (x y &optional (x-stride 1) (y-stride 1)) (defmfun linear-fit (x y &optional (x-stride 1) (y-stride 1))
"gsl_fit_linear" "gsl_fit_linear"
(((gsl-array x) :pointer) (x-stride size) (((c-pointer x) :pointer) (x-stride sizet)
((gsl-array y) :pointer) (y-stride size) ((c-pointer y) :pointer) (y-stride sizet)
((dim0 x) size) (c0 :double) (c1 :double) ((dim0 x) sizet) (c0 :double) (c1 :double)
(cov00 :double) (cov01 :double) (cov11 :double) (cov00 :double) (cov01 :double) (cov11 :double)
(sumsq :double)) (sumsq :double))
:documentation ; FDL :documentation ; FDL
...@@ -38,10 +38,10 @@ ...@@ -38,10 +38,10 @@
(defmfun weighted-linear-fit (defmfun weighted-linear-fit
(x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1)) (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
"gsl_fit_wlinear" "gsl_fit_wlinear"
(((gsl-array x) :pointer) (x-stride size) (((c-pointer x) :pointer) (x-stride sizet)
((gsl-array weight) :pointer) (weight-stride size) ((c-pointer weight) :pointer) (weight-stride sizet)
((gsl-array y) :pointer) (y-stride size) ((c-pointer y) :pointer) (y-stride sizet)
((dim0 x) size) (c0 :double) (c1 :double) ((dim0 x) sizet) (c0 :double) (c1 :double)
(cov00 :double) (cov01 :double) (cov11 :double) (cov00 :double) (cov01 :double) (cov11 :double)
(chisq :double)) (chisq :double))
:documentation ; FDL :documentation ; FDL
...@@ -77,9 +77,9 @@ ...@@ -77,9 +77,9 @@
(defmfun multiplier-fit (x y &optional (x-stride 1) (y-stride 1)) (defmfun multiplier-fit (x y &optional (x-stride 1) (y-stride 1))
"gsl_fit_mul" "gsl_fit_mul"
(((gsl-array x) :pointer) (x-stride size) (((c-pointer x) :pointer) (x-stride sizet)
((gsl-array y) :pointer) (y-stride size) ((c-pointer y) :pointer) (y-stride sizet)
((dim0 x) size) (c1 :double) (cov11 :double) ((dim0 x) sizet) (c1 :double) (cov11 :double)
(sumsq :double)) (sumsq :double))
:documentation ; FDL :documentation ; FDL
"The best-fit linear regression coefficient c1 of the model Y = c_1 "The best-fit linear regression coefficient c1 of the model Y = c_1
...@@ -93,10 +93,10 @@ ...@@ -93,10 +93,10 @@
(defmfun weighted-multiplier-fit (defmfun weighted-multiplier-fit
(x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1)) (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
"gsl_fit_wmul" "gsl_fit_wmul"
(((gsl-array x) :pointer) (x-stride size) (((c-pointer x) :pointer) (x-stride sizet)
((gsl-array weight) :pointer) (weight-stride size) ((c-pointer weight) :pointer) (weight-stride sizet)
((gsl-array y) :pointer) (y-stride size) ((c-pointer y) :pointer) (y-stride sizet)
((dim0 x) size) (c1 :double) (cov11 :double) ((dim0 x) sizet) (c1 :double) (cov11 :double)
(chisq :double)) (chisq :double))
:documentation ; FDL :documentation ; FDL
"Compute the best-fit linear regression coefficient "Compute the best-fit linear regression coefficient
...@@ -130,7 +130,7 @@ ...@@ -130,7 +130,7 @@
(defmfun allocate-fit-workspace (number-of-observations number-of-parameters) (defmfun allocate-fit-workspace (number-of-observations number-of-parameters)
"gsl_multifit_linear_alloc" "gsl_multifit_linear_alloc"
((number-of-observations size) (number-of-parameters size)) ((number-of-observations sizet) (number-of-parameters sizet))
:c-return :pointer :c-return :pointer
:index (letm fit-workspace) :index (letm fit-workspace)
:documentation ; FDL :documentation ; FDL
...@@ -147,9 +147,9 @@ ...@@ -147,9 +147,9 @@
(defmfun linear-mfit (defmfun linear-mfit
(model observations parameters covariance tolerance workspace) (model observations parameters covariance tolerance workspace)
"gsl_multifit_linear" "gsl_multifit_linear"
(((pointer model) :pointer) ((pointer observations) :pointer) (((mpointer model) :pointer) ((mpointer observations) :pointer)
(tolerance :double) (tolerance :double)
((pointer parameters) :pointer) (covariance :pointer) (chisq :double) ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
(workspace :pointer)) (workspace :pointer))
:documentation ; FDL :documentation ; FDL
"Compute the best-fit parameters c of the model "Compute the best-fit parameters c of the model
...@@ -168,10 +168,10 @@ ...@@ -168,10 +168,10 @@
(defmfun linear-mfit-svd (defmfun linear-mfit-svd
(model observations parameters covariance tolerance workspace) (model observations parameters covariance tolerance workspace)
"gsl_multifit_linear_svd" "gsl_multifit_linear_svd"
(((pointer model) :pointer) ((pointer observations) :pointer) (((mpointer model) :pointer) ((mpointer observations) :pointer)
(tolerance :double) (tolerance :double)
(rank size) (rank sizet)
((pointer parameters) :pointer) (covariance :pointer) (chisq :double) ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
(workspace :pointer)) (workspace :pointer))
:return ((dcref chisq) (scref rank)) :return ((dcref chisq) (scref rank))
:documentation ; FDL :documentation ; FDL
...@@ -195,11 +195,11 @@ ...@@ -195,11 +195,11 @@
(defmfun weighted-linear-mfit (defmfun weighted-linear-mfit
(model weight observations parameters covariance workspace) (model weight observations parameters covariance workspace)
"gsl_multifit_wlinear" "gsl_multifit_wlinear"
(((pointer model) :pointer) (((mpointer model) :pointer)
((pointer weight) :pointer) ((mpointer weight) :pointer)
((pointer observations) :pointer) ((mpointer observations) :pointer)
((pointer parameters) :pointer) ((mpointer parameters) :pointer)
((pointer covariance) :pointer) (chisq :double) ((mpointer covariance) :pointer) (chisq :double)
(workspace :pointer)) (workspace :pointer))
:documentation ; FDL :documentation ; FDL
"Compute the best-fit parameters c of the weighted "Compute the best-fit parameters c of the weighted
...@@ -217,12 +217,12 @@ ...@@ -217,12 +217,12 @@
(defmfun weighted-linear-mfit-svd (defmfun weighted-linear-mfit-svd
(model weight observations parameters covariance tolerance workspace) (model weight observations parameters covariance tolerance workspace)
"gsl_multifit_wlinear_svd" "gsl_multifit_wlinear_svd"
(((pointer model) :pointer) (((mpointer model) :pointer)
((pointer weight) :pointer) ((mpointer weight) :pointer)
((pointer observations) :pointer) ((mpointer observations) :pointer)
(tolerance :double) (tolerance :double)
(rank size) (rank sizet)
((pointer parameters) :pointer) (covariance :pointer) (chisq :double) ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
(workspace :pointer)) (workspace :pointer))
:return ((dcref chisq) (scref rank)) :return ((dcref chisq) (scref rank))
:documentation ; FDL :documentation ; FDL
...@@ -243,8 +243,8 @@ ...@@ -243,8 +243,8 @@
(defmfun multi-linear-estimate (x coefficients covariance) (defmfun multi-linear-estimate (x coefficients covariance)
"gsl_multifit_linear_est" "gsl_multifit_linear_est"
(((pointer x) :pointer) ((pointer coefficients) :pointer) (((mpointer x) :pointer) ((mpointer coefficients) :pointer)
((pointer covariance) :pointer) (y :double) (y-error :double)) ((mpointer covariance) :pointer) (y :double) (y-error :double))
:documentation ; FDL :documentation ; FDL
"Use the best-fit multilinear regression coefficients "Use the best-fit multilinear regression coefficients
and their covariance matrix to compute the fitted function value and their covariance matrix to compute the fitted function value
...@@ -258,9 +258,9 @@ ...@@ -258,9 +258,9 @@
(defun univariate-linear-least-squares-example () (defun univariate-linear-least-squares-example ()
"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.
(letm ((x (vector-double-float #(1970.0d0 1980.0d0 1990.0d0 2000.0d0))) (letm ((x (vector-double-float (a 1970.0d0 1980.0d0 1990.0d0 2000.0d0)))
(y (vector-double-float #(12.0d0 11.0d0 14.0d0 13.0d0))) (y (vector-double-float (a 12.0d0 11.0d0 14.0d0 13.0d0)))
(w (vector-double-float #(0.1d0 0.2d0 0.3d0 0.4d0)))) (w (vector-double-float (a 0.1d0 0.2d0 0.3d0 0.4d0))))
(multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq) (multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq)
(weighted-linear-fit x w y) (weighted-linear-fit x w y)
(format t "~&Best fit: Y = ~8,5f + ~8,5f X" c0 c1) (format t "~&Best fit: Y = ~8,5f + ~8,5f X" c0 c1)
...@@ -301,8 +301,8 @@ ...@@ -301,8 +301,8 @@
(defun mv-linear-least-squares-example (data) (defun mv-linear-least-squares-example (data)
"Second example in Section 36.5 of the GSL manual." "Second example in Section 36.5 of the GSL manual."
(letm ((n (length data)) chisq (letm ((n (length data)) chisq
(x (matrix-double-float n 3)) (x (matrix-double-float (list n 3)))
(cov (matrix-double-float 3 3)) (cov (matrix-double-float '(3 3)))
(y (vector-double-float n)) (y (vector-double-float n))
(w (vector-double-float n)) (w (vector-double-float n))
(c (vector-double-float 3))) (c (vector-double-float 3)))
......
;; Multivariate minimization. ;; Multivariate minimization.
;; Liam Healy <Tue Jan 8 2008 - 21:28> ;; 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$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -56,7 +56,7 @@ ...@@ -56,7 +56,7 @@
(defmfun allocate-mfminimizer (type dimension) (defmfun allocate-mfminimizer (type dimension)
"gsl_multimin_fminimizer_alloc" "gsl_multimin_fminimizer_alloc"
((type :pointer) (dimension size)) ((type :pointer) (dimension sizet))
:c-return :pointer :c-return :pointer
:export nil :export nil
:index (letm mfminimizer) :index (letm mfminimizer)
...@@ -66,7 +66,7 @@ ...@@ -66,7 +66,7 @@
(defmfun allocate-mfdfminimizer (type dimension) (defmfun allocate-mfdfminimizer (type dimension)
"gsl_multimin_fdfminimizer_alloc" "gsl_multimin_fdfminimizer_alloc"
((type :pointer) (dimension size)) ((type :pointer) (dimension sizet))
:c-return :pointer :c-return :pointer
:export nil :export nil
:index (letm mfdfminimizer) :index (letm mfdfminimizer)
...@@ -77,7 +77,7 @@ ...@@ -77,7 +77,7 @@
(defmfun set-mfminimizer (minimizer function initial step-size) (defmfun set-mfminimizer (minimizer function initial step-size)
"gsl_multimin_fminimizer_set" "gsl_multimin_fminimizer_set"
((minimizer :pointer) (function :pointer) ((minimizer :pointer) (function :pointer)
((pointer initial) :pointer) ((pointer step-size) :pointer)) ((mpointer initial) :pointer) ((mpointer step-size) :pointer))
:export nil :export nil
:index (letm mfminimizer) :index (letm mfminimizer)
:documentation ; FDL :documentation ; FDL
...@@ -91,7 +91,7 @@ ...@@ -91,7 +91,7 @@
(minimizer function-derivative initial step-size tolerance) (minimizer function-derivative initial step-size tolerance)
"gsl_multimin_fdfminimizer_set" "gsl_multimin_fdfminimizer_set"
((minimizer :pointer) (function-derivative :pointer) ((minimizer :pointer) (function-derivative :pointer)
((pointer initial) :pointer) (step-size :double) ((mpointer initial) :pointer) (step-size :double)
(tolerance :double)) (tolerance :double))
:export nil :export nil
:index (letm mfdfminimizer) :index (letm mfdfminimizer)
...@@ -160,16 +160,16 @@ ...@@ -160,16 +160,16 @@
(defmfun mfminimizer-x (minimizer) (defmfun mfminimizer-x (minimizer)
"gsl_multimin_fminimizer_x" "gsl_multimin_fminimizer_x"
((minimizer :pointer)) ((minimizer :pointer))
:c-return (canswer :pointer) :c-return :pointer
:return ((make-data-from-pointer canswer)) :return (:c-return)
:documentation ; FDL :documentation ; FDL
"The current best estimate of the location of the minimum.") "The current best estimate of the location of the minimum.")
(defmfun mfdfminimizer-x (minimizer) (defmfun mfdfminimizer-x (minimizer)
"gsl_multimin_fdfminimizer_x" "gsl_multimin_fdfminimizer_x"
((minimizer :pointer)) ((minimizer :pointer))
:c-return (canswer :pointer) :c-return :pointer
:return ((make-data-from-pointer canswer)) :return (:c-return)
:documentation ; FDL :documentation ; FDL
"The current best estimate of the location of the minimum.") "The current best estimate of the location of the minimum.")
...@@ -197,8 +197,8 @@ ...@@ -197,8 +197,8 @@
(defmfun mfdfminimizer-gradient (minimizer) (defmfun mfdfminimizer-gradient (minimizer)
"gsl_multimin_fdfminimizer_gradient" "gsl_multimin_fdfminimizer_gradient"
((minimizer :pointer)) ((minimizer :pointer))
:c-return (canswer :pointer) :c-return :pointer
:return ((make-data-from-pointer canswer)) :return (:c-return)
:documentation ; FDL :documentation ; FDL
"The current best estimate of the gradient for the minimizer.") "The current best estimate of the gradient for the minimizer.")
...@@ -215,7 +215,7 @@ ...@@ -215,7 +215,7 @@
(defmfun min-test-gradient (gradient absolute-error) (defmfun min-test-gradient (gradient absolute-error)
"gsl_multimin_test_gradient" "gsl_multimin_test_gradient"
(((pointer gradient) :pointer) (absolute-error :double)) ((gradient :pointer) (absolute-error :double))
:c-return :success-continue :c-return :success-continue
:documentation ; FDL :documentation ; FDL
"Test the norm of the gradient against the "Test the norm of the gradient against the
...@@ -316,8 +316,8 @@ ...@@ -316,8 +316,8 @@
(defun parabaloid (gsl-vector-pointer) (defun parabaloid (gsl-vector-pointer)
"A parabaloid function of two arguments, given in GSL manual Sec. 35.4." "A parabaloid function of two arguments, given in GSL manual Sec. 35.4."
(let ((x (vref gsl-vector-pointer 0)) (let ((x (maref gsl-vector-pointer 0))
(y (vref gsl-vector-pointer 1)) (y (maref gsl-vector-pointer 1))
(dp0 (aref *parabaloid-center* 0)) (dp0 (aref *parabaloid-center* 0))
(dp1 (aref *parabaloid-center* 1))) (dp1 (aref *parabaloid-center* 1)))
(+ (* 10 (expt (- x dp0) 2)) (+ (* 10 (expt (- x dp0) 2))
...@@ -326,13 +326,13 @@ ...@@ -326,13 +326,13 @@
(defun parabaloid-derivative (defun parabaloid-derivative
(arguments-gv-pointer derivative-gv-pointer) (arguments-gv-pointer derivative-gv-pointer)
(let ((x (vref arguments-gv-pointer 0)) (let ((x (maref arguments-gv-pointer 0))
(y (vref arguments-gv-pointer 1)) (y (maref arguments-gv-pointer 1))
(dp0 (aref *parabaloid-center* 0)) (dp0 (aref *parabaloid-center* 0))
(dp1 (aref *parabaloid-center* 1))) (dp1 (aref *parabaloid-center* 1)))
(setf (vref derivative-gv-pointer 0) (setf (maref derivative-gv-pointer 0)
(* 20 (- x dp0)) (* 20 (- x dp0))
(vref derivative-gv-pointer 1) (maref derivative-gv-pointer 1)
(* 40 (- y dp1))))) (* 40 (- y dp1)))))
(defun parabaloid-and-derivative (defun parabaloid-and-derivative
...@@ -345,7 +345,7 @@ ...@@ -345,7 +345,7 @@
parabaloid 2 parabaloid-derivative parabaloid-and-derivative) parabaloid 2 parabaloid-derivative parabaloid-and-derivative)
(defun multimin-example-fletcher-reeves () (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 (minimizer
(mfdfminimizer *conjugate-fletcher-reeves* 2 parabaloid (mfdfminimizer *conjugate-fletcher-reeves* 2 parabaloid
initial 0.01d0 1.0d-4))) initial 0.01d0 1.0d-4)))
...@@ -375,7 +375,7 @@ ...@@ -375,7 +375,7 @@
(def-minimization-functions parabaloid-f 2) (def-minimization-functions parabaloid-f 2)
(defun multimin-example-nelder-mead () (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))) (step-size (vector-double-float 2)))
(set-all step-size 1.0d0) (set-all step-size 1.0d0)
(letm ((minimizer (letm ((minimizer
......
;; 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: <2008-03-09 19:30:55EDT nonlinear-least-squares.lisp> ;; Time-stamp: <2008-08-31 15:13:29EDT nonlinear-least-squares.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
(defmfun allocate-ffit (solver-type number-of-observations number-of-parameters) (defmfun allocate-ffit (solver-type number-of-observations number-of-parameters)
"gsl_multifit_fsolver_alloc" "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 :c-return :pointer
:export nil :export nil
:index (letm nonlinear-ffit) :index (letm nonlinear-ffit)
...@@ -59,7 +59,7 @@ ...@@ -59,7 +59,7 @@
(defmfun allocate-fdffit (solver-type number-of-observations number-of-parameters) (defmfun allocate-fdffit (solver-type number-of-observations number-of-parameters)
"gsl_multifit_fdfsolver_alloc" "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 :c-return :pointer
:export nil :export nil
:index (letm nonlinear-fdffit) :index (letm nonlinear-fdffit)
...@@ -69,7 +69,7 @@ ...@@ -69,7 +69,7 @@
(defmfun set-fdffit (solver function initial-guess) (defmfun set-fdffit (solver function initial-guess)
"gsl_multifit_fdfsolver_set" "gsl_multifit_fdfsolver_set"
((solver :pointer) (function :pointer) ((pointer initial-guess) :pointer)) ((solver :pointer) (function :pointer) ((mpointer initial-guess) :pointer))
:documentation ; FDL :documentation ; FDL
"Initialize or reinitialize an existing solver "Initialize or reinitialize an existing solver
to use the function and the initial guess.") to use the function and the initial guess.")
...@@ -114,7 +114,7 @@ ...@@ -114,7 +114,7 @@
;; See /usr/include/gsl/gsl_multifit_nlin.h ;; 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."
(function :pointer) (function :pointer)
(dimensions size) (dimensions sizet)
(parameters :pointer)) (parameters :pointer))
(cffi:defcstruct gsl-fdffit-function (cffi:defcstruct gsl-fdffit-function
...@@ -124,8 +124,8 @@ ...@@ -124,8 +124,8 @@
(function :pointer) (function :pointer)
(df :pointer) (df :pointer)
(fdf :pointer) (fdf :pointer)
(number-of-observations size) (number-of-observations sizet)
(number-of-parameters size) (number-of-parameters sizet)
(parameters :pointer)) (parameters :pointer))
(export 'def-fitting-functions) (export 'def-fitting-functions)
...@@ -170,16 +170,16 @@ ...@@ -170,16 +170,16 @@
(defmfun current-parameters-ffit (solver) (defmfun current-parameters-ffit (solver)
"gsl_multifit_fsolver_position" "gsl_multifit_fsolver_position"
((solver :pointer)) ((solver :pointer))
:c-return (canswer :pointer) :c-return :pointer
:return ((make-data-from-pointer canswer)) :return (:c-return)
:documentation ; FDL :documentation ; FDL
"The current best-fit parameters.") "The current best-fit parameters.")
(defmfun current-parameters-fdffit (solver) (defmfun current-parameters-fdffit (solver)
"gsl_multifit_fdfsolver_position" "gsl_multifit_fdfsolver_position"
((solver :pointer)) ((solver :pointer))
:c-return (canswer :pointer) :c-return :pointer
:return ((make-data-from-pointer canswer)) :return (:c-return)
:documentation ; FDL :documentation ; FDL
"The current best-fit parameters.") "The current best-fit parameters.")
...@@ -216,7 +216,7 @@ ...@@ -216,7 +216,7 @@
(defmfun fit-gradient (jacobian function-values gradient) (defmfun fit-gradient (jacobian function-values gradient)
"gsl_multifit_gradient" "gsl_multifit_gradient"
((jacobian :pointer) ((pointer function-values) :pointer) (gradient :pointer)) ((jacobian :pointer) ((mpointer function-values) :pointer) (gradient :pointer))
:documentation ; FDL :documentation ; FDL
"Compute the gradient of \Phi(x) = (1/2) ||F(x)||^2 "Compute the gradient of \Phi(x) = (1/2) ||F(x)||^2
from the Jacobian matrix and the function values using from the Jacobian matrix and the function values using
...@@ -279,7 +279,7 @@ ...@@ -279,7 +279,7 @@
(defmfun ls-covariance (jacobian relative-error covariance) (defmfun ls-covariance (jacobian relative-error covariance)
"gsl_multifit_covar" "gsl_multifit_covar"
((jacobian :pointer) (relative-error :double) ((pointer covariance) :pointer)) ((jacobian :pointer) (relative-error :double) ((mpointer covariance) :pointer))
:return (covariance) :return (covariance)
:documentation ; FDL :documentation ; FDL
"Compute the covariance matrix of the best-fit parameters "Compute the covariance matrix of the best-fit parameters
...@@ -333,25 +333,25 @@ ...@@ -333,25 +333,25 @@
(make-exponent-fit-data (make-exponent-fit-data
:n *number-of-observations* :n *number-of-observations*
:y :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))) (letm ((rng (random-number-generator *mt19937* 0)))
(dotimes (i *number-of-observations* arr) (dotimes (i *number-of-observations* arr)
(setf (aref arr i) (setf (aref 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-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) (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
for the nonlinear least squares example." for the nonlinear least squares example."
(let ((A (vref x 0)) (let ((A (maref x 0))
(lambda (vref x 1)) (lambda (maref x 1))
(b (vref x 2))) (b (maref x 2)))
(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 *number-of-observations*)
(setf (vref f i) (setf (maref f i)
;; the difference model - observation = - residual ;; the difference model - observation = - residual
(/ (- (+ (* A (exp (* (- lambda) i))) b) (aref y i)) (/ (- (+ (* A (exp (* (- lambda) i))) b) (aref y i))
(aref sigma i))))))) (aref sigma i)))))))
...@@ -360,16 +360,16 @@ ...@@ -360,16 +360,16 @@
"Compute the partial derivatives of the negative of the "Compute the partial derivatives of the negative of the
residuals with the exponential model residuals with the exponential model
for the nonlinear least squares example." for the nonlinear least squares example."
(let ((A (vref x 0)) (let ((A (maref x 0))
(lambda (vref 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 *number-of-observations*)
(let ((e (exp (* (- lambda) i))) (let ((e (exp (* (- lambda) i)))
(s (aref sigma i))) (s (aref sigma i)))
(setf (mref jacobian i 0) (/ e s) (setf (maref jacobian i 0) (/ e s)
(mref jacobian i 1) (* -1 i A (/ e s)) (maref jacobian i 1) (* -1 i A (/ e s))
(mref jacobian i 2) (/ s))))))) (maref jacobian i 2) (/ s)))))))
(defun exponential-residual-fdf (x f jacobian) (defun exponential-residual-fdf (x f jacobian)
"Compute the function and partial derivatives of the negative of the "Compute the function and partial derivatives of the negative of the
...@@ -382,42 +382,47 @@ ...@@ -382,42 +382,47 @@
exponential-residual *number-of-observations* *number-of-parameters* exponential-residual *number-of-observations* *number-of-parameters*
exponential-residual-derivative exponential-residual-fdf) 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 () (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 (covariance
(matrix-double-float *number-of-parameters* *number-of-parameters*)) (matrix-double-float (list *number-of-parameters* *number-of-parameters*)))
(fit (nonlinear-fdffit (fit (nonlinear-fdffit
*levenberg-marquardt* *levenberg-marquardt*
*number-of-observations* *number-of-observations*
*number-of-parameters* *number-of-parameters*
exponential-residual exponential-residual
init))) 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)))) (err (i) `(sqrt (maref covariance ,i ,i))))
(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 (make-data-from-pointer (fdffit-slot fit 'f)))) (norm-f fit))
(loop for iter from 0 below 25 (loop for iter from 0 below 25
until until
(and (plusp iter) (and (plusp iter)
(fit-test-delta (fit-test-delta
(fdffit-slot fit 'dx) (fdffit-slot fit 'x) (fdffit-slot fit 'dx) (fdffit-slot fit 'x)
1.0d-4 1.0d-4)) 1.0d-4 1.0d-4))
do do
(iterate-fdffit fit) (iterate-fdffit fit)
(ls-covariance (fdffit-slot fit 'jacobian) 0.0d0 covariance) (ls-covariance (fdffit-slot fit 'jacobian) 0.0d0 covariance)
(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)