Commit 78b88074 authored by Liam Healy's avatar Liam Healy

Minimization will take scalar functions, added tests, optional printing

The minimization functions may now take and return scalars or marrays.
The examples take optional arguments with the solver method and
whether to print the steps.  They are now part of the #'examples and
test suite.  Use maref instead of cl-array in
#'roots-multi-example-no-derivative.  Add optional solver method to
roots one and roots multi, and add tests for all solvers.
parent 1060d1dc
...@@ -34,6 +34,7 @@ Fail in CCL but not in SBCL: ...@@ -34,6 +34,7 @@ Fail in CCL but not in SBCL:
series-acceleration series-acceleration
linear-least-squares linear-least-squares
nonlinear-least-squares nonlinear-least-squares
minimization-multi
=== GSL questions === === GSL questions ===
Sun Nov 9 2008: Sun Nov 9 2008:
......
;; Definition of GSLL system ;; Definition of GSLL system
;; Liam Healy ;; Liam Healy
;; Time-stamp: <2009-01-21 22:58:22EST gsll-tests.asd> ;; Time-stamp: <2009-01-22 21:37:49EST gsll-tests.asd>
;; $Id$ ;; $Id$
(asdf:defsystem "gsll-tests" (asdf:defsystem "gsll-tests"
...@@ -116,6 +116,8 @@ ...@@ -116,6 +116,8 @@
(:file "matrix-variance-with-fixed-mean") (:file "matrix-variance-with-fixed-mean")
(:file "matrix-variance-with-mean") (:file "matrix-variance-with-mean")
(:file "median-percentile") (:file "median-percentile")
(:file "minimization-one")
(:file "minimization-multi")
(:file "monte-carlo") (:file "monte-carlo")
(:file "multinomial") (:file "multinomial")
(:file "negative-binomial") (:file "negative-binomial")
......
;; Foreign callback functions. ;; Foreign callback functions.
;; Liam Healy ;; Liam Healy
;; Time-stamp: <2009-01-19 16:44:13EST callback.lisp> ;; Time-stamp: <2009-01-22 22:36:00EST callback.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -20,13 +20,17 @@ ...@@ -20,13 +20,17 @@
;;; Other GSL tasks make use of callback functions with different ;;; Other GSL tasks make use of callback functions with different
;;; characteristics. Since they are specific to each of the tasks, ;;; characteristics. Since they are specific to each of the tasks,
;;; they are defined with those tasks. A complexity encountered with ;;; they are defined with those tasks. The macro #'defmcallback can
;;; using a vector of double floats not present with scalars is that ;;; specify that the CL function is to expect in arglist and return as
;;; there is no portable way to make a C array available directly to ;;; multiple values scalar quantities that come from and will be bound
;;; Lisp, so callbacks using C arrays must read them using a macro. ;;; to either marrays or C vectors. This is done with a list of the
;;; Therefore it is necessary to define the function in a way that ;;; type and size, e.g. (:double 3), and for setting :set, type size,
;;; prevents its use in Lisp; to ameliorate this, the macro ;;; e.g. (:set :double 3). If the 'marray argument is nil, it will
;;; #'with-c-double is provided to give named access to the elements. ;;; expand to read or set a C vector; if it is T, it will expand to
;;; read or set a marray. This allows the user to define ordinary CL
;;; functions with scalars as input and output. However, it may be
;;; desirable to read and set marrays, in which case :pointer is the
;;; right specification.
(export (export
'(def-single-function undef-cbstruct defun-single '(def-single-function undef-cbstruct defun-single
...@@ -90,6 +94,7 @@ ...@@ -90,6 +94,7 @@
;;; (callback-args '(:double (:double 2) (:set :double 2))) ;;; (callback-args '(:double (:double 2) (:set :double 2)))
;;; ((#:ARG1193 :DOUBLE) (#:ARG1194 :POINTER) (#:ARG1195 :POINTER)) ;;; ((#:ARG1193 :DOUBLE) (#:ARG1194 :POINTER) (#:ARG1195 :POINTER))
(defun callback-args (types) (defun callback-args (types)
"The arguments passed by GSL to the callback function." "The arguments passed by GSL to the callback function."
(mapcar (lambda (type) (mapcar (lambda (type)
...@@ -102,12 +107,18 @@ ...@@ -102,12 +107,18 @@
;;; (embedded-clfunc-args '(:double (:double 2) (:set :double 2)) (callback-args '(:double (:double 2) (:set :double 2)))) ;;; (embedded-clfunc-args '(:double (:double 2) (:set :double 2)) (callback-args '(:double (:double 2) (:set :double 2))))
;;; (#:ARG1244 (MEM-AREF #:ARG1245 ':DOUBLE 0) (MEM-AREF #:ARG1245 ':DOUBLE 1)) ;;; (#:ARG1244 (MEM-AREF #:ARG1245 ':DOUBLE 0) (MEM-AREF #:ARG1245 ':DOUBLE 1))
(defvar *setting-spec* '(:set))
(defun embedded-clfunc-args (types callback-args &optional marray) (defun embedded-clfunc-args (types callback-args &optional marray)
"The arguments passed to the CL function call embedded in the callback." "The arguments passed to the CL function call embedded in the callback.
If 'marray is T, then reference GSL arrays; otherwise reference raw
C vectors. A specification (:set ...) means that the CL function
will define the array as multiple values; if the size is negative,
then the opposite value will be used for marray."
(loop for spec in types (loop for spec in types
for (symbol nil) in callback-args for (symbol nil) in callback-args
append append
(unless (and (listp spec) (eq (first spec) :set)) (unless (and (listp spec) (member (first spec) *setting-spec*))
(if (listp spec) (if (listp spec)
(if (third spec) (if (third spec)
;; matrix, marrays only ;; matrix, marrays only
...@@ -115,10 +126,10 @@ ...@@ -115,10 +126,10 @@
append append
(loop for j from 0 below (third spec) (loop for j from 0 below (third spec)
collect collect
`(maref ,symbol ,i ,j ',(cffi-cl (first spec))))) `(maref ,symbol ,i ,j ',(cffi-cl (first spec)))))
;; vector, marray or C array ;; vector, marray or C array
(loop for ind from 0 below (second spec) (loop for ind from 0 below (abs (second spec))
collect (if marray collect (if (if (minusp (second spec)) (not marray) marray)
`(maref ,symbol ,ind nil ',(cffi-cl (first spec))) `(maref ,symbol ,ind nil ',(cffi-cl (first spec)))
`(cffi:mem-aref ,symbol ',(first spec) ,ind)))) `(cffi:mem-aref ,symbol ',(first spec) ,ind))))
(list symbol))))) (list symbol)))))
...@@ -128,14 +139,17 @@ ...@@ -128,14 +139,17 @@
(multiple-value-bind (settype setcba) (multiple-value-bind (settype setcba)
(loop for cba in callback-args (loop for cba in callback-args
for type in types for type in types
for setting = (and (listp type) (eq (first type) :set)) for setting = (and (listp type) (member (first type) *setting-spec*))
when setting when setting
collect cba into setcba collect cba into setcba
when setting when setting
collect type into settype collect type into settype
finally (return (values (mapcar 'rest settype) setcba))) finally (return (values (mapcar 'rest settype) setcba)))
(let* ((setvbls (embedded-clfunc-args settype setcba marray)) (let* ((setvbls (embedded-clfunc-args settype setcba marray))
(count (apply '+ (mapcar (lambda (inds) (apply '* (rest inds))) settype))) (count
(apply
'+
(mapcar (lambda (inds) (abs (apply '* (rest inds)))) settype)))
(mvbvbls (loop repeat count collect (gensym "SETCB")))) (mvbvbls (loop repeat count collect (gensym "SETCB"))))
(if (zerop count) (if (zerop count)
form form
...@@ -145,9 +159,6 @@ ...@@ -145,9 +159,6 @@
for setvbl in setvbls for setvbl in setvbls
append (list setvbl mvbvbl)))))))) append (list setvbl mvbvbl))))))))
;;; (DEFMCALLBACK VANDERPOL :SUCCESS-FAILURE (:DOUBLE (:DOUBLE 2) (:SET :DOUBLE 2)))
;;; (DEFMCALLBACK VANDERPOL :SUCCESS-FAILURE (:DOUBLE (:DOUBLE 2)))
(defmacro defmcallback (defmacro defmcallback
(name &optional (return-type :double) (argument-types :double) (name &optional (return-type :double) (argument-types :double)
additional-argument-types marray) additional-argument-types marray)
......
;; Multivariate minimization. ;; Multivariate minimization.
;; Liam Healy <Tue Jan 8 2008 - 21:28> ;; Liam Healy <Tue Jan 8 2008 - 21:28>
;; Time-stamp: <2009-01-03 13:27:12EST minimization-multi.lisp> ;; Time-stamp: <2009-01-22 21:35:26EST minimization-multi.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
;;; /usr/include/gsl/gsl_multimin.h ;;; /usr/include/gsl/gsl_multimin.h
;; In the parabaloid example, I notice that the consruct ;;; In the parabaloid example, I notice that the consruct
;; (min-test-gradient (mfdfminimizer-gradient minimizer) 1.0d-3) ;;; (min-test-gradient (mfdfminimizer-gradient minimizer) 1.0d-3)
;; is constructing a CL vector-double-float (in mfdfminimizer-gradient) and ;;; is constructing a CL vector-double-float (in mfdfminimizer-gradient) and
;; then immediately pulling out the pointer (in min-test-gradient). It ;;; then immediately pulling out the pointer (in min-test-gradient). It
;; is easy enough to eliminate this, but then mfdfminimizer-gradient ;;; is easy enough to eliminate this, but then mfdfminimizer-gradient
;; would not be useful to a CL user. ;;; would not be useful to a CL user.
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Function definition ;;;; Function definition
...@@ -28,23 +28,24 @@ ...@@ -28,23 +28,24 @@
;;; they point to have different signatures. ;;; they point to have different signatures.
(export 'def-minimization-functions) (export 'def-minimization-functions)
(defmacro def-minimization-functions (function dimensions &optional df fdf) (defmacro def-minimization-functions (function dimensions &optional df fdf array)
"Setup functions for multivariate minimization. "Setup functions for multivariate minimization.
The CL functions name and derivative should be defined previously The CL functions name and derivative should be defined previously
with defuns." with defuns."
`(progn (let ((vdec (if array :pointer `((:double ,dimensions))))
(defmcallback ,function :double :pointer) (svdec (if array :pointer `(:set :double ,dimensions))))
,@(when df `(progn
`((defmcallback ,df :pointer :pointer :pointer) (defmcallback ,function :double ,vdec nil t)
(defmcallback ,fdf :pointer :pointer (:pointer :pointer)))) ,@(when df
,(if df `((defmcallback ,df :void ,vdec (,svdec) t)
`(defcbstruct (,function function ,df df ,fdf fdf) (defmcallback ,fdf :void ,vdec ((:set :double -1) ,svdec) t)))
gsl-mfunction-fdf ,(if df
((dimensions ,dimensions))) `(defcbstruct (,function function ,df df ,fdf fdf)
`(defcbstruct (,function function) gsl-mfunction-fdf
gsl-mfunction ((dimensions ,dimensions)))
((dimensions ,dimensions)))))) `(defcbstruct (,function function)
gsl-mfunction
((dimensions ,dimensions)))))))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Initialization ;;;; Initialization
...@@ -280,10 +281,13 @@ ...@@ -280,10 +281,13 @@
;;; Examples from Sec. 35.8. ;;; Examples from Sec. 35.8.
;;; Example using derivatives, taking a vector argument.
(defparameter *parabaloid-center* #(1.0d0 2.0d0)) (defparameter *parabaloid-center* #(1.0d0 2.0d0))
(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.
This version takes a vector-double-float argument."
(let ((x (maref gsl-vector-pointer 0)) (let ((x (maref gsl-vector-pointer 0))
(y (maref gsl-vector-pointer 1)) (y (maref gsl-vector-pointer 1))
(dp0 (aref *parabaloid-center* 0)) (dp0 (aref *parabaloid-center* 0))
...@@ -303,20 +307,21 @@ ...@@ -303,20 +307,21 @@
(maref derivative-gv-pointer 1) (maref derivative-gv-pointer 1)
(* 40 (- y dp1))))) (* 40 (- y dp1)))))
(defun parabaloid-and-derivative (defun parabaloid-and-derivative (arguments-gv-pointer derivative-gv-pointer)
(arguments-gv-pointer fnval derivative-gv-pointer) (prog1
(setf (dcref fnval) (parabaloid arguments-gv-pointer)) (parabaloid arguments-gv-pointer)
(parabaloid-derivative (parabaloid-derivative
arguments-gv-pointer derivative-gv-pointer)) arguments-gv-pointer derivative-gv-pointer)))
(def-minimization-functions (def-minimization-functions
parabaloid 2 parabaloid-derivative parabaloid-and-derivative) parabaloid 2 parabaloid-derivative parabaloid-and-derivative t)
(defun multimin-example-fletcher-reeves () (defun multimin-example-derivative
(&optional (method *conjugate-fletcher-reeves*) (print-steps t))
(let* ((initial #m(5.0d0 7.0d0)) (let* ((initial #m(5.0d0 7.0d0))
(minimizer (minimizer
(make-multi-dimensional-minimizer-fdf (make-multi-dimensional-minimizer-fdf
*conjugate-fletcher-reeves* 2 parabaloid method 2 parabaloid
initial 0.01d0 1.0d-4))) initial 0.01d0 1.0d-4)))
(loop with status = T (loop with status = T
for iter from 0 below 100 for iter from 0 below 100
...@@ -327,29 +332,36 @@ ...@@ -327,29 +332,36 @@
(not (min-test-gradient (not (min-test-gradient
(mfdfminimizer-gradient minimizer) (mfdfminimizer-gradient minimizer)
1.0d-3))) 1.0d-3)))
(let ((x (solution minimizer))) (when print-steps
(format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f" (let ((x (solution minimizer)))
iter (maref x 0) (maref x 1) (format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~&"
(function-value minimizer))) iter (maref x 0) (maref x 1)
finally (return (function-value minimizer))))
(let ((x (solution minimizer))) finally
(values (maref x 0) (maref x 1))))))) (return
(let ((x (solution minimizer)))
;;; Because def-minimization-functions bind a symbol (values (maref x 0) (maref x 1) (function-value minimizer)))))))
;;; of the same name as the first function, and we want both to run,
;;; we'll make an alias function so we can use both. ;;; Example without derivatives, same function but now defined using
(defun parabaloid-f (gsl-vector-pointer)
(parabaloid gsl-vector-pointer)) (defun parabaloid-f (x y)
"A parabaloid function of two arguments, given in GSL manual Sec. 35.4.
This version takes scalar arguments."
(let ((dp0 (aref *parabaloid-center* 0))
(dp1 (aref *parabaloid-center* 1)))
(+ (* 10 (expt (- x dp0) 2))
(* 20 (expt (- y dp1) 2))
30)))
(def-minimization-functions parabaloid-f 2) (def-minimization-functions parabaloid-f 2)
(defun multimin-example-nelder-mead () (defun multimin-example-no-derivative
(let ((initial #m(5.0d0 7.0d0)) (&optional (method *simplex-nelder-mead*) (print-steps t))
(step-size (make-marray 'double-float :dimensions 2))) (let ((step-size (make-marray 'double-float :dimensions 2)))
(set-all step-size 1.0d0) (set-all step-size 1.0d0)
(let ((minimizer (let ((minimizer
(make-multi-dimensional-minimizer-f (make-multi-dimensional-minimizer-f
*simplex-nelder-mead* 2 parabaloid-f initial step-size))) method 2 parabaloid-f #m(5.0d0 7.0d0) step-size)))
(loop with status = T and size (loop with status = T and size
for iter from 0 below 100 for iter from 0 below 100
while status while status
...@@ -358,11 +370,19 @@ ...@@ -358,11 +370,19 @@
(mfminimizer-size minimizer) (mfminimizer-size minimizer)
status status
(not (min-test-size size 1.0d-2))) (not (min-test-size size 1.0d-2)))
(let ((x (solution minimizer))) (when print-steps
(format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f~40t~8,3f" (let ((x (solution minimizer)))
iter (maref x 0) (maref x 1) (format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~40t~8,3f~&"
(function-value minimizer) iter (maref x 0) (maref x 1)
size)) (function-value minimizer)
finally (return size)))
(let ((x (solution minimizer))) finally
(values (maref x 0) (maref x 1)))))))) (return
(let ((x (solution minimizer)))
(values (maref x 0) (maref x 1) (function-value minimizer))))))))
(save-test minimization-multi
(multimin-example-no-derivative *simplex-nelder-mead* nil)
(multimin-example-derivative *conjugate-fletcher-reeves* nil)
(multimin-example-derivative *conjugate-polak-ribiere* nil)
(multimin-example-derivative *vector-bfgs* nil))
;; Univariate minimization ;; Univariate minimization
;; Liam Healy Tue Jan 8 2008 - 21:02 ;; Liam Healy Tue Jan 8 2008 - 21:02
;; Time-stamp: <2009-01-03 15:42:59EST minimization-one.lisp> ;; Time-stamp: <2009-01-22 18:59:23EST minimization-one.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -178,24 +178,34 @@ ...@@ -178,24 +178,34 @@
(defun-single minimization-one-fn (x) (defun-single minimization-one-fn (x)
(1+ (cos x))) (1+ (cos x)))
(defun minimization-one-example () (defun minimization-one-example
(&optional (minimizer-type *brent-fminimizer*) (print-steps t))
"Solving a minimum, the example given in Sec. 33.8 of the GSL manual." "Solving a minimum, the example given in Sec. 33.8 of the GSL manual."
(let ((max-iter 100) (let ((max-iter 100)
(minimizer (minimizer
(make-one-dimensional-minimizer (make-one-dimensional-minimizer
*brent-fminimizer* minimization-one-fn 2.0d0 0.0d0 6.0d0))) minimizer-type minimization-one-fn 2.0d0 0.0d0 6.0d0)))
(format t "~&iter ~6t [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower") (when print-steps
(format
t
"iter ~6t [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower~&"))
(loop for iter from 0 (loop for iter from 0
for min = (solution minimizer) for min = (solution minimizer)
for lower = (fminimizer-x-lower minimizer) for lower = (fminimizer-x-lower minimizer)
for upper = (fminimizer-x-upper minimizer) for upper = (fminimizer-x-upper minimizer)
do (iterate minimizer) do (iterate minimizer)
(when print-steps
(format t "~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g~&"
iter lower upper
min (- min pi)
(- upper lower)))
while (and (< iter max-iter) while (and (< iter max-iter)
;; abs and rel error swapped in example? ;; abs and rel error swapped in example?
(not (min-test-interval lower upper 0.001d0 0.0d0))) (not (min-test-interval lower upper 0.001d0 0.0d0)))
do finally
(format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g" (return (values iter lower upper min (- min pi) (- upper lower))))))
iter lower upper
min (- min pi)
(- upper lower)))))
(save-test minimization-one
(minimization-one-example *brent-fminimizer* nil)
(minimization-one-example *golden-section-fminimizer* 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-21 22:38:26EST roots-multi.lisp> ;;; Time-stamp: <2009-01-22 22:19:34EST roots-multi.lisp>
;;; $Id$ ;;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -385,32 +385,33 @@ ...@@ -385,32 +385,33 @@
(def-mfunction rosenbrock 2) (def-mfunction rosenbrock 2)
(defun roots-multi-example (&optional (print-steps t)) (defun roots-multi-example-no-derivative
(&optional (method *hybrid-scaled*) (print-steps t))
"Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual." "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual."
(let ((max-iter 1000)) (let ((max-iter 1000)
(let* ((vect #m(-10.0d0 -5.0d0)) (solver (make-multi-dimensional-root-solver-f
(solver (make-multi-dimensional-root-solver-f method rosenbrock
*hybrid-scaled* rosenbrock vect))) #m(-10.0d0 -5.0d0))))
(loop for iter from 0 (loop for iter from 0
with fnval and argval with fnval and argval
while (and (< iter max-iter) while (and (< iter max-iter)
(not (multiroot-test-residual solver 1.0d-7))) (not (multiroot-test-residual solver 1.0d-7)))
do do
(iterate solver) (iterate solver)
(setf fnval (cl-array (function-value solver)) (setf fnval (function-value solver)
argval (cl-array (solution solver))) argval (solution solver))
(when print-steps (when print-steps
(format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&" (format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
iter iter
(aref argval 0) (maref argval 0)
(aref argval 1) (maref argval 1)
(aref fnval 0) (maref fnval 0)
(aref fnval 1))) (maref fnval 1)))
finally (return finally (return
(values (aref argval 0) (values (maref argval 0)
(aref argval 1) (maref argval 1)
(aref fnval 0) (maref fnval 0)
(aref fnval 1))))))) (maref fnval 1))))))
(defun rosenbrock-df (arg0 arg1) (defun rosenbrock-df (arg0 arg1)
"The partial derivatives of the Rosenbrock functions." "The partial derivatives of the Rosenbrock functions."
...@@ -436,7 +437,7 @@ ...@@ -436,7 +437,7 @@
(def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2) (def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2)
(defun roots-multi-example-df (&optional (print-steps t)) (defun roots-multi-example-derivative (&optional (method *gnewton-mfdfsolver*) (print-steps t))
"Solving Rosenbrock with derivatives, the example given in Sec. 34.8 "Solving Rosenbrock with derivatives, the example given in Sec. 34.8
of the GSL manual." of the GSL manual."
(flet ((print-state (iter argval fnval) (flet ((print-state (iter argval fnval)
...@@ -447,27 +448,35 @@ ...@@ -447,27 +448,35 @@
(maref argval 1) (maref argval 1)
(maref fnval 0) (maref fnval 0)
(maref fnval 1))))) (maref fnval 1)))))
(let ((max-iter 1000)) (let ((max-iter 1000)
(let* ((vect #m(-10.0d0 -5.0d0)) (solver (make-multi-dimensional-root-solver-fdf
(solver (make-multi-dimensional-root-solver-fdf method rosenbrock-f
*gnewton-mfdfsolver* rosenbrock-f vect))) #m(-10.0d0 -5.0d0))))
(loop for iter from 0 (loop for iter from 0
with fnval = (function-value solver) with fnval = (function-value solver)
and argval = (solution solver) and argval = (solution solver)
while (and (< iter max-iter) while (and (< iter max-iter)
(not (multiroot-test-residual solver 1.0d-7))) (not (multiroot-test-residual solver 1.0d-7)))
initially (print-state iter argval fnval) initially (print-state iter argval fnval)
do do
(iterate solver) (iterate solver)
(setf fnval (function-value solver) (setf fnval (function-value solver)
argval (solution solver)) argval (solution solver))
(print-state iter argval fnval) (print-state iter argval fnval)
finally (return finally (return
(values (maref argval 0) (values (maref argval 0)
(maref argval 1) (maref argval 1)
(maref fnval 0) (maref fnval 0)
(maref fnval 1)))))))) (maref fnval 1)))))))
;; To see step-by-step information as the solution progresses, make
;; the last argument T.
(save-test roots-multi (save-test roots-multi
(roots-multi-example nil) (roots-multi-example-no-derivative *hybrid-unscaled* nil)
(roots-multi-example-df nil)) (roots-multi-example-no-derivative *hybrid-scaled* nil)
(roots-multi-example-no-derivative *discrete-newton* nil)
(roots-multi-example-no-derivative *broyden* nil)
(roots-multi-example-derivative *newton-mfdfsolver* nil)
(roots-multi-example-derivative *gnewton-mfdfsolver* nil)
(roots-multi-example-derivative *powells-hybrid* nil)
(roots-multi-example-derivative *powells-hybrid-unscaled* nil))
;; One-dimensional root solver. ;; One-dimensional root solver.
;; Liam Healy ;; Liam Healy
;; Time-stamp: <2009-01-19 22:04:19EST roots-one.lisp> ;; Time-stamp: <2009-01-22 22:15:50EST roots-one.lisp>
;; $Id$ ;; $Id$
(in-package :gsl) (in-package :gsl)
...@@ -330,12 +330,12 @@ ...@@ -330,12 +330,12 @@
(def-single-function quadratic) (def-single-function quadratic)
(defun roots-one-example (&optional (print-steps t)) (defun roots-one-example-no-derivative
(&optional (method *brent-fsolver*) (print-steps t))
"Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."
(let ((max-iter 50) (let ((max-iter 50)
(solver (solver
(make-one-dimensional-root-solver-f (make-one-dimensional-root-solver-f method quadratic 0.0d0 5.0d0)))
*brent-fsolver* quadratic 0.0d0 5.0d0)))
(when print-steps (when print-steps
(format t "iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)~&")) (format t "iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)~&"))
(loop for iter from 0 (loop for iter from 0
...@@ -361,12 +361,13 @@ ...@@ -361,12 +361,13 @@
(def-solver-functions (def-solver-functions
quadratic-df quadratic-derivative quadratic-and-derivative) quadratic-df quadratic-derivative quadratic-and-derivative)
(defun roots-one-fdf-example (&optional (print-steps t)) (defun roots-one-example-derivative
(&optional (method *newton-fdfsolver*) (print-steps t))
"Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."
(let* ((max-iter 100) (let* ((max-iter 100)
(initial 5.0d0) (initial 5.0d0)
(solver (make-one-dimensional-root-solver-fdf (solver (make-one-dimensional-root-solver-fdf
*newton-fdfsolver* quadratic-df initial))) method quadratic-df initial)))
(when print-steps (when print-steps
(format t "iter ~6t ~8troot ~22terr ~34terr(est)~&")) (format t "iter ~6t ~8troot ~22terr ~34terr(est)~&"))
(loop for iter from 0 (loop for iter from 0
...@@ -380,6 +381,12 @@ ...@@ -380,6 +381,12 @@
iter root (- root (sqrt 5.0d0)) (- root oldroot))) iter root (- root (sqrt 5.0d0)) (- root oldroot)))
finally (return root)))) finally (return root))))
;; To see step-by-step information as the solution progresses, make
;; the last argument T.
(save-test roots-one (save-test roots-one
(roots-one-example nil) (roots-one-example-no-derivative *bisection-fsolver* nil)
(roots-one-fdf-example nil)) (roots-one-example-no-derivative *false-position-fsolver* nil)
(roots-one-example-no-derivative *brent-fsolver* nil)
(roots-one-example-derivative *newton-fdfsolver* nil)
(roots-one-example-derivative *secant-fdfsolver* nil)
(roots-one-example-derivative *steffenson-fdfsolver* nil))
;; Regression test MINIMIZATION-MULTI for GSLL, automatically generated
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST MINIMIZATION-MULTI
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.9920430849306285d0 1.9969168063253164d0
30.000823246638923d0)
(MULTIPLE-VALUE-LIST
(MULTIMIN-EXAMPLE-NO-DERIVATIVE *SIMPLEX-NELDER-MEAD*
NIL)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.9999999999999997d0 2.0d0 30.0d0)
(MULTIPLE-VALUE-LIST
(MULTIMIN-EXAMPLE-DERIVATIVE
*CONJUGATE-FLETCHER-REEVES* NIL)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.9999999999999997d0 2.0d0 30.0d0)
(MULTIPLE-VALUE-LIST
(MULTIMIN-EXAMPLE-DERIVATIVE *CONJUGATE-POLAK-RIBIERE*
NIL)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.9999999999999997d0 2.0d0 30.0d0)
(MULTIPLE-VALUE-LIST
(MULTIMIN-EXAMPLE-DERIVATIVE *VECTOR-BFGS* NIL))))
;; Regression test MINIMIZATION-ONE for GSLL, automatically generated
(in-package :gsl)