octave_packages/m/optimization/fminbnd.m

   1 ## Copyright (C) 2008-2012 VZLU Prague, a.s.
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18 ##
  19 ## Author: Jaroslav Hajek <highegg@gmail.com>
  20
  21 ## -*- texinfo -*-
  22 ## @deftypefn {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options})
  23 ## Find a minimum point of a univariate function.  @var{fun} should be a
  24 ## function
  25 ## handle or name.  @var{a}, @var{b} specify a starting interval.  @var{options}
  26 ## is a
  27 ## structure specifying additional options.  Currently, @code{fminbnd}
  28 ## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"},
  29 ## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}.
  30 ## For description of these options, see @ref{doc-optimset,,optimset}.
  31 ##
  32 ## On exit, the function returns @var{x}, the approximate minimum point
  33 ## and @var{fval}, the function value thereof.
  34 ## @var{info} is an exit flag that can have these values:
  35 ##
  36 ## @itemize
  37 ## @item 1
  38 ## The algorithm converged to a solution.
  39 ##
  40 ## @item 0
  41 ## Maximum number of iterations or function evaluations has been exhausted.
  42 ##
  43 ## @item -1
  44 ## The algorithm has been terminated from user output function.
  45 ## @end itemize
  46 ## @seealso{optimset, fzero, fminunc}
  47 ## @end deftypefn
  48
  49 ## This is patterned after opt/fmin.f from Netlib, which in turn is taken from
  50 ## Richard Brent: Algorithms For Minimization Without Derivatives, Prentice-Hall (1973)
  51
  52 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
  53 ## PKG_ADD: [~] = __all_opts__ ("fminbnd");
  54
  55 function [x, fval, info, output] = fminbnd (fun, xmin, xmax, options = struct ())
  56
  57   ## Get default options if requested.
  58   if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults'))
  59     x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, \
  60     "OutputFcn", [], "FunValCheck", "off");
  61     return;
  62   endif
  63
  64   if (nargin < 2 || nargin > 4)
  65     print_usage ();
  66   endif
  67
  68   if (ischar (fun))
  69     fun = str2func (fun, "global");
  70   endif
  71
  72   ## TODO
  73   ## displev = optimget (options, "Display", "notify");
  74   funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
  75   outfcn = optimget (options, "OutputFcn");
  76   tolx = optimget (options, "TolX", 1e-8);
  77   maxiter = optimget (options, "MaxIter", Inf);
  78   maxfev = optimget (options, "MaxFunEvals", Inf);
  79
  80   if (funvalchk)
  81     ## Replace fun with a guarded version.
  82     fun = @(x) guarded_eval (fun, x);
  83   endif
  84
  85   ## The default exit flag if exceeded number of iterations.
  86   info = 0;
  87   niter = 0;
  88   nfev = 0;
  89   sqrteps = eps (class (xmin + xmax));
  90
  91   c = 0.5*(3-sqrt(5));
  92   a = xmin; b = xmax;
  93   v = a + c*(b-a);
  94   w = x = v;
  95   e = 0;
  96   fv = fw = fval = fun (x);
  97   nfev++;
  98
  99   while (niter < maxiter && nfev < maxfev)
 100     xm = 0.5*(a+b);
 101     ## FIXME: the golden section search can actually get closer than sqrt(eps)...
 102     ## sometimes. Sometimes not, it depends on the function. This is the strategy
 103     ## from the Netlib code. Something yet smarter would be good.
 104     tol = 2 * sqrteps * abs (x) + tolx / 3;
 105     if (abs (x - xm) <= (2*tol - 0.5*(b-a)))
 106       info = 1;
 107       break;
 108     endif
 109
 110     if (abs (e) > tol)
 111       dogs = false;
 112       ## Try inverse parabolic step.
 113       r = (x - w)*(fval - fv);
 114       q = (x - v)*(fval - fw);
 115       p = (x - v)*q - (x - w)*r;
 116       q = 2*(q - r);
 117       p *= -sign (q);
 118       q = abs (q);
 119       r = e;
 120       e = d;
 121
 122       if (abs (p) < abs (0.5*q*r) && p > q*(a-x) && p < q*(b-x))
 123         ## The parabolic step is acceptable.
 124         d = p / q;
 125         u = x + d;
 126
 127         ## f must not be evaluated too close to ax or bx.
 128         if (min (u-a, b-u) < 2*tol)
 129           d = tol * (sign (xm - x) + (xm == x));
 130         endif
 131       else
 132         dogs = true;
 133       endif
 134     else
 135       dogs = true;
 136     endif
 137     if (dogs)
 138       ## Default to golden section step.
 139       e = ifelse (x >= xm, a - x, b - x);
 140       d = c * e;
 141     endif
 142
 143      ## f must not be evaluated too close to x.
 144      u = x + max (abs (d), tol) * (sign (d) + (d == 0));
 145
 146      fu = fun (u);
 147      nfev++;
 148      niter++;
 149
 150      ## update  a, b, v, w, and x
 151
 152      if (fu <= fval)
 153        if (u < x)
 154          b = x;
 155        else
 156          a = x;
 157        endif
 158        v = w; fv = fw;
 159        w = x; fw = fval;
 160        x = u; fval = fu;
 161      else
 162        ## The following if-statement was originally executed even if fu == fval.
 163        if (u < x)
 164          a = u;
 165        else
 166          b = u;
 167        endif
 168        if (fu <= fw || w == x)
 169          v = w; fv = fw;
 170          w = u; fw = fu;
 171        elseif (fu <= fv || v == x || v == w)
 172          v = u;
 173          fv = fu;
 174        endif
 175      endif
 176
 177     ## If there's an output function, use it now.
 178     if (outfcn)
 179       optv.funccount = nfev;
 180       optv.fval = fval;
 181       optv.iteration = niter;
 182       if (outfcn (x, optv, "iter"))
 183         info = -1;
 184         break;
 185       endif
 186     endif
 187   endwhile
 188
 189   output.iterations = niter;
 190   output.funcCount = nfev;
 191   output.bracket = [a, b];
 192   ## FIXME: bracketf possibly unavailable.
 193
 194 endfunction
 195
 196 ## An assistant function that evaluates a function handle and checks for
 197 ## bad results.
 198 function fx = guarded_eval (fun, x)
 199   fx = fun (x);
 200   fx = fx(1);
 201   if (! isreal (fx))
 202     error ("fminbnd:notreal", "fminbnd: non-real value encountered");
 203   elseif (isnan (fx))
 204     error ("fminbnd:isnan", "fminbnd: NaN value encountered");
 205   endif
 206 endfunction
 207
 208 %!shared opt0
 209 %! opt0 = optimset ("tolx", 0);
 210 %!assert (fminbnd (@cos, pi/2, 3*pi/2, opt0), pi, 10*sqrt(eps))
 211 %!assert (fminbnd (@(x) (x - 1e-3)^4, -1, 1, opt0), 1e-3, 10e-3*sqrt(eps))
 212 %!assert (fminbnd (@(x) abs(x-1e7), 0, 1e10, opt0), 1e7, 10e7*sqrt(eps))
 213 %!assert (fminbnd (@(x) x^2 + sin(2*pi*x), 0.4, 1, opt0), fzero (@(x) 2*x + 2*pi*cos(2*pi*x), [0.4, 1], opt0), sqrt(eps))