minimization-one.lisp 8.41 KB
Newer Older
1 2
;; Univariate minimization
;; Liam Healy Tue Jan  8 2008 - 21:02
3
;; Time-stamp: <2009-01-24 14:46:49EST minimization-one.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_min.h

10 11 12 13
;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************

14 15 16 17 18 19 20 21
;; 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."
22
  "set"					; should use set_with_values?
23 24
  ((function :pointer) (minimum :double) (lower :double) (upper :double)))

liam's avatar
liam committed
25
(defmfun set-fminimizer-with-values
26 27 28
    (minimizer function x-minimum x-lower x-upper
	       f-minimum f-lower f-upper)
  "gsl_min_fminimizer_set_with_values"
29
  (((mpointer minimizer) :pointer) (function :pointer)
30 31
   (x-minimum :double) (x-lower :double) (x-upper :double)
   (f-minimum :double) (f-lower :double) (f-upper :double))
liam's avatar
liam committed
32
  :documentation			; FDL
33 34 35 36 37
  "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.")

38
(defmfun name ((minimizer one-dimensional-minimizer))
39
  "gsl_min_fminimizer_name"
40 41
  (((mpointer minimizer) :pointer))
  :definition :method
42
  :c-return :string
liam's avatar
liam committed
43
  :documentation			; FDL
44 45 46 47 48 49
  "The name of the minimizer.")

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

50
(defmfun iterate ((minimizer one-dimensional-minimizer))
51
  "gsl_min_fminimizer_iterate"
52
  (((mpointer minimizer) :pointer))
53
  :definition :method
54
  :c-return :success-continue
liam's avatar
liam committed
55
  :documentation			; FDL
56
  "Perform a single iteration of the minimizer.  The following
57
   errors may be signalled: 'bad-function-supplied,
58 59 60 61 62
   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.")

63
(defmfun solution ((minimizer one-dimensional-minimizer))
64
  "gsl_min_fminimizer_x_minimum"
65
  (((mpointer minimizer) :pointer))
66
  :definition :method
67
  :c-return :double
liam's avatar
liam committed
68
  :documentation			; FDL
69 70
  "The current estimate of the position of the minimum for the minimizer.")

liam's avatar
liam committed
71
(defmfun fminimizer-x-lower (minimizer)
72
  "gsl_min_fminimizer_x_lower"
73
  (((mpointer minimizer) :pointer))
74
  :c-return :double
liam's avatar
liam committed
75
  :documentation			; FDL
76 77
  "The current lower bound of the interval for the minimizer.")

liam's avatar
liam committed
78
(defmfun fminimizer-x-upper (minimizer)
79
  "gsl_min_fminimizer_x_upper"
80
  (((mpointer minimizer) :pointer))
81
  :c-return :double
liam's avatar
liam committed
82
  :documentation			; FDL
83 84
  "The current upper bound of the interval for the minimizer.")

85
(defmfun function-value ((minimizer one-dimensional-minimizer))
86
  "gsl_min_fminimizer_f_minimum"
87
  (((mpointer minimizer) :pointer))
88
  :definition :method
89
  :c-return :double
liam's avatar
liam committed
90
  :documentation			; FDL
91 92 93
  "The value of the function at the current estimate of the minimum for the
   minimizer.")

liam's avatar
liam committed
94
(defmfun fminimizer-f-lower (minimizer)
95
  "gsl_min_fminimizer_f_lower"
96
  (((mpointer minimizer) :pointer))
97
  :c-return :double
liam's avatar
liam committed
98
  :documentation			; FDL
99 100 101
  "The value of the function at the current estimate of the lower bound
   for the minimizer.")

liam's avatar
liam committed
102
(defmfun fminimizer-f-upper (minimizer)
103
  "gsl_min_fminimizer_f_upper"
104
  (((mpointer minimizer) :pointer))
105
  :c-return :double
liam's avatar
liam committed
106
  :documentation			; FDL
107 108 109 110 111 112 113
  "The value of the function at the current estimate of the upper bound
   for the minimizer.")

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

114
(defmfun min-test-interval (lower upper absolute-error relative-error)
115 116 117
  "gsl_min_test_interval"
  ((lower :double) (upper :double)
   (absolute-error :double) (relative-error :double))
118
  :c-return :success-continue		; guess that this is s-c, not s-f
liam's avatar
liam committed
119
  :documentation			; FDL
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
  "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
140
(defmpar *golden-section-fminimizer* "gsl_min_fminimizer_goldensection"
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
  "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
156
(defmpar *brent-fminimizer* "gsl_min_fminimizer_brent"
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
  "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.

178
(defun minimization-one-fn (x)
179 180
  (1+ (cos x)))

181 182
(defun minimization-one-example
    (&optional (minimizer-type *brent-fminimizer*) (print-steps t))
183
  "Solving a minimum, the example given in Sec. 33.8 of the GSL manual."
184 185 186
  (let ((max-iter 100)
	(minimizer
	 (make-one-dimensional-minimizer
187 188 189
	  minimizer-type
	  (make-single-function minimization-one-fn)
	  2.0d0 0.0d0 6.0d0)))
190 191 192 193
    (when print-steps
      (format
       t
       "iter ~6t   [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower~&"))
194
    (loop for iter from 0
195
       for min = (solution minimizer)
196 197
       for lower = (fminimizer-x-lower minimizer)
       for upper = (fminimizer-x-upper minimizer)
198
       do (iterate minimizer)
199 200 201 202 203 204
       (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)))

205 206 207
       while  (and (< iter max-iter)
		   ;; abs and rel error swapped in example?
		   (not (min-test-interval lower upper 0.001d0 0.0d0)))
208 209
       finally
       (return (values iter lower upper min (- min pi) (- upper lower))))))
210

211 212 213
(save-test minimization-one
 (minimization-one-example *brent-fminimizer* nil)
 (minimization-one-example *golden-section-fminimizer* nil))