minimization-one.lisp 8.74 KB
Newer Older
1 2
;; Univariate minimization
;; Liam Healy Tue Jan  8 2008 - 21:02
3
;; Time-stamp: <2008-12-25 10:22:13EST minimization-one.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7 8 9 10 11

(in-package :gsl)

;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************

12 13 14 15 16 17 18 19 20 21 22 23 24
#|
;; The one-dimensional minimizer uses a function (no derivative) only.
(defmobject one-dimensional-minimizer
    "gsl_min_fminimizer"
  ((type :pointer))
  "one-dimensional minimizer" ; FDL
  "Make an instance of a minimizer of the given type.  Optionally
   set to use the function and the initial search interval [lower,
   upper], with a guess for the location of the minimum."
  "set"
  ((function :pointer) (minimum :double) (lower :double) (upper :double)))
|#

25
(defgo-s (fminimizer type function minimum lower upper)
26
	 allocate-fminimizer free-fminimizer set-fminimizer)
27

liam's avatar
liam committed
28
(defmfun allocate-fminimizer (type)
29 30 31
  "gsl_min_fminimizer_alloc"
  ((type :pointer))
  :c-return :pointer
32
  :export nil
33
  :index (letm fminimizer)
liam's avatar
liam committed
34
  :documentation			; FDL
35 36
  "Allocate an instance of a minimizer of the given type.")

liam's avatar
liam committed
37
(defmfun set-fminimizer (minimizer function minimum lower upper)
38 39 40
  "gsl_min_fminimizer_set"
  ((minimizer :pointer) (function :pointer)
   (minimum :double) (lower :double) (upper :double))
41
  :export nil
42
  :index (letm fminimizer)
liam's avatar
liam committed
43
  :documentation			; FDL
44 45 46 47
  "Set, or reset, an existing minimizer to use the
   function and the initial search interval [lower,
   upper], with a guess for the location of the minimum.")

48
;;; Use this in letm macro in any way?
liam's avatar
liam committed
49
(defmfun set-fminimizer-with-values
50 51 52 53 54 55
    (minimizer function x-minimum x-lower x-upper
	       f-minimum f-lower f-upper)
  "gsl_min_fminimizer_set_with_values"
  ((minimizer :pointer) (function :pointer)
   (x-minimum :double) (x-lower :double) (x-upper :double)
   (f-minimum :double) (f-lower :double) (f-upper :double))
liam's avatar
liam committed
56
  :documentation			; FDL
57 58 59 60 61
  "Set, or reset, an existing minimizer to use the
   function and the initial search interval [lower,
   upper], with a guess for the location of the minimum, using
   supplied rather than computed values of the function.")

liam's avatar
liam committed
62
(defmfun free-fminimizer (minimizer)
63 64 65
  "gsl_min_fminimizer_free"
  ((minimizer :pointer))
  :c-return :void
66
  :export nil
67
  :index (letm fminimizer)
liam's avatar
liam committed
68
  :documentation			; FDL
69 70
  "Free all the memory associated with the minimizer.")

liam's avatar
liam committed
71
(defmfun fminimizer-name (minimizer)
72 73 74
  "gsl_min_fminimizer_name"
  ((minimizer :pointer))
  :c-return :string
liam's avatar
liam committed
75
  :documentation			; FDL
76 77 78 79 80 81
  "The name of the minimizer.")

;;;;****************************************************************************
;;;; Iteration
;;;;****************************************************************************

liam's avatar
liam committed
82
(defmfun iterate-fminimizer (minimizer)
83 84 85
  "gsl_min_fminimizer_iterate"
  ((minimizer :pointer))
  :c-return :success-continue
liam's avatar
liam committed
86
  :documentation			; FDL
87
  "Perform a single iteration of the minimizer.  The following
88
   errors may be signalled: 'bad-function-supplied,
89 90 91 92 93
   the iteration encountered a singular point where the function or its
   derivative evaluated to infinity or NaN, or
   :FAILURE, the algorithm could not improve the current best approximation or
   bounding interval.")

liam's avatar
liam committed
94
(defmfun fminimizer-x-minimum (minimizer)
95
  "gsl_min_fminimizer_x_minimum"
96 97
  ((minimizer :pointer))
  :c-return :double
liam's avatar
liam committed
98
  :documentation			; FDL
99 100
  "The current estimate of the position of the minimum for the minimizer.")

liam's avatar
liam committed
101
(defmfun fminimizer-x-lower (minimizer)
102 103 104
  "gsl_min_fminimizer_x_lower"
  ((minimizer :pointer))
  :c-return :double
liam's avatar
liam committed
105
  :documentation			; FDL
106 107
  "The current lower bound of the interval for the minimizer.")

liam's avatar
liam committed
108
(defmfun fminimizer-x-upper (minimizer)
109 110 111
  "gsl_min_fminimizer_x_upper"
  ((minimizer :pointer))
  :c-return :double
liam's avatar
liam committed
112
  :documentation			; FDL
113 114
  "The current upper bound of the interval for the minimizer.")

liam's avatar
liam committed
115
(defmfun fminimizer-f-minimum (minimizer)
116 117 118
  "gsl_min_fminimizer_f_minimum"
  ((minimizer :pointer))
  :c-return :double
liam's avatar
liam committed
119
  :documentation			; FDL
120 121 122
  "The value of the function at the current estimate of the minimum for the
   minimizer.")

liam's avatar
liam committed
123
(defmfun fminimizer-f-lower (minimizer)
124 125 126
  "gsl_min_fminimizer_f_lower"
  ((minimizer :pointer))
  :c-return :double
liam's avatar
liam committed
127
  :documentation			; FDL
128 129 130
  "The value of the function at the current estimate of the lower bound
   for the minimizer.")

liam's avatar
liam committed
131
(defmfun fminimizer-f-upper (minimizer)
132 133 134
  "gsl_min_fminimizer_f_upper"
  ((minimizer :pointer))
  :c-return :double
liam's avatar
liam committed
135
  :documentation			; FDL
136 137 138 139 140 141 142
  "The value of the function at the current estimate of the upper bound
   for the minimizer.")

;;;;****************************************************************************
;;;; Stopping parameters
;;;;****************************************************************************

liam's avatar
liam committed
143 144
(defmfun
 min-test-interval (lower upper absolute-error relative-error)
145 146 147
  "gsl_min_test_interval"
  ((lower :double) (upper :double)
   (absolute-error :double) (relative-error :double))
148
  :c-return :success-continue		; guess that this is s-c, not s-f
liam's avatar
liam committed
149
  :documentation			; FDL
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
  "Test for the convergence of the interval [lower,upper]
   with absolute error and relative error specified.
   The test returns T if the following condition is achieved:
   |a - b| < epsabs + epsrel min(|a|,|b|) 
   when the interval x = [a,b] does not include the origin.  If the
   interval includes the origin then min(|a|,|b|) is replaced by
   zero (which is the minimum value of |x| over the interval).  This
   ensures that the relative error is accurately estimated for minima close
   to the origin.

   This condition on the interval also implies that any estimate of the
   minimum x_m in the interval satisfies the same condition with respect
   to the true minimum x_m^*,
   |x_m - x_m^*| < epsabs + epsrel x_m^*
   assuming that the true minimum x_m^* is contained within the interval.")

;;;;****************************************************************************
;;;; Minimization algorithms
;;;;****************************************************************************

liam's avatar
liam committed
170
(defmpar *golden-section-fminimizer* "gsl_min_fminimizer_goldensection"
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
  "The golden section algorithm is the simplest method of bracketing
   the minimum of a function.  It is the slowest algorithm provided by the
   library, with linear convergence.

   On each iteration, the algorithm first compares the subintervals from
   the endpoints to the current minimum.  The larger subinterval is divided
   in a golden section (using the famous ratio (3-sqrt 5)/2 =
   0.3189660...) and the value of the function at this new point is
   calculated.  The new value is used with the constraint f(a') >
   f(x') < f(b') to a select new interval containing the minimum, by
   discarding the least useful point.  This procedure can be continued
   indefinitely until the interval is sufficiently small.  Choosing the
   golden section as the bisection ratio can be shown to provide the
   fastest convergence for this type of algorithm.")

liam's avatar
liam committed
186
(defmpar *brent-fminimizer* "gsl_min_fminimizer_brent"
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
  "The Brent minimization algorithm combines a parabolic
   interpolation with the golden section algorithm.  This produces a fast
   algorithm which is still robust.

   The outline of the algorithm can be summarized as follows: on each
   iteration Brent's method approximates the function using an
   interpolating parabola through three existing points.  The minimum of the
   parabola is taken as a guess for the minimum.  If it lies within the
   bounds of the current interval then the interpolating point is accepted,
   and used to generate a smaller interval.  If the interpolating point is
   not accepted then the algorithm falls back to an ordinary golden section
   step.  The full details of Brent's method include some additional checks
   to improve convergence.")

;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************

;;; This is the example given in Sec. 33.8.  The results are different
;;; than given there.

208
(defun-single minimization-one-fn (x)
209 210 211 212
  (1+ (cos x)))

(defun minimization-one-example ()
  "Solving a minimum, the example given in Sec. 33.8 of the GSL manual."
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
  (letm ((max-iter 100)
	 (minimizer
	  (fminimizer *brent-fminimizer* minimization-one-fn 2.0d0 0.0d0 6.0d0)))
    (format t "~&iter ~6t   [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower")
    (loop for iter from 0
	  for min = (fminimizer-x-minimum minimizer)
	  for lower = (fminimizer-x-lower minimizer)
	  for upper = (fminimizer-x-upper minimizer)
	  do (iterate-fminimizer minimizer)
	  while  (and (< iter max-iter)
		      ;; abs and rel error swapped in example?
		      (not (min-test-interval lower upper 0.001d0 0.0d0)))
	  do
	  (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)))))
230