Mercurial > octave-nkf

changeset 10296:035ac548a67e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement fminbnd
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 10 Feb 2010 13:33:54 +0100
parents 13d382fc758b
children ed88ea036716
files scripts/ChangeLog scripts/optimization/fminbnd.m scripts/optimization/fzero.m scripts/optimization/module.mk
diffstat 4 files changed, 228 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Wed Feb 10 08:41:49 2010 +0100
+++ b/scripts/ChangeLog	Wed Feb 10 13:33:54 2010 +0100
@@ -1,3 +1,9 @@
+2010-02-10  Jaroslav Hajek  <highegg@gmail.com>
+
+	* optimization/fminbnd.m: New function.
+	* optimization/module.mk: Add it here.
+	* optimization/fzero.m: Make default TolX consistent with fminbnd.
+
 2010-02-10  Jaroslav Hajek  <highegg@gmail.com>
 
 	* general/private/__isequal__.m: Remove list-related stub.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/optimization/fminbnd.m	Wed Feb 10 13:33:54 2010 +0100
@@ -0,0 +1,212 @@
+## Copyright (C) 2008, 2009, 2010 VZLU Prague, a.s.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+##
+## Author: Jaroslav Hajek <highegg@gmail.com>
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options})
+## Find a minimum point of a univariate function.  @var{fun} should be a function
+## handle or name.  @var{a}, @var{b} specify a starting interval. @var{options} is a
+## structure specifying additional options.  Currently, @code{fminbnd}
+## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"},
+## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}. 
+## For description of these options, see @ref{doc-optimset,,optimset}.
+## 
+## On exit, the function returns @var{x}, the approximate minimum point
+## and @var{fval}, the function value thereof.
+## @var{info} is an exit flag that can have these values:
+## @itemize
+## @item 1
+## The algorithm converged to a solution.
+## @item 0
+## Maximum number of iterations or function evaluations has been exhausted.
+## @item -1
+## The algorithm has been terminated from user output function.
+## @item -2 
+## A general unexpected error.
+## @item -3
+## A non-real value encountered.
+## @item -4
+## A NaN value encountered.
+## @end itemize
+## @seealso{optimset, fzero, fminunc} 
+## @end deftypefn
+
+## This is patterned after opt/fmin.f from Netlib, which in turn is taken from
+## Richard Brent: Algorithms For Minimization Without Derivatives, Prentice-Hall (1973)
+
+## PKG_ADD: __all_opts__ ("fminbnd");
+
+function [x, fval, info, output] = fminbnd (fun, xmin, xmax, options = struct ())
+
+  ## Get default options if requested.
+  if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults'))
+    x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, \
+    "OutputFcn", [], "FunValCheck", "off");
+    return;
+  endif
+
+  if (nargin < 2 || nargin > 4)
+    print_usage ();
+  endif
+
+  if (ischar (fun))
+    fun = str2func (fun, "global");
+  endif
+
+  ## TODO
+  ## displev = optimget (options, "Display", "notify");
+  funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
+  outfcn = optimget (options, "OutputFcn");
+  tolx = optimget (options, "TolX", 1e-8);
+  maxiter = optimget (options, "MaxIter", Inf);
+  maxfev = optimget (options, "MaxFunEvals", Inf);
+
+  persistent mu = 0.5;
+
+  if (funvalchk)
+    ## Replace fun with a guarded version.
+    fun = @(x) guarded_eval (fun, x);
+  endif
+
+  ## The default exit flag if exceeded number of iterations.
+  info = 0;
+  niter = 0;
+  nfev = 0;
+  eps = eps (class (xmin + xmax));
+
+  c = 0.5*(3-sqrt(5));
+  a = xmin; b = xmax;
+  v = a + c*(b-a);
+  w = x = v;
+  e = 0;
+  fv = fw = fval = fun (x);
+  nfev++;
+
+  while (niter < maxiter && nfev < maxfev)
+    xm = 0.5*(a+b);
+    tol = 2 * eps * abs (x) + tolx / 3;
+    if (abs (x - xm) <= (2*tol - 0.5*(b-a)))
+      info = 1;
+      break;
+    endif
+
+    if (abs (e) > tol)
+      dogs = false;
+      ## Try inverse parabolic step.
+      r = (x - w)*(fval - fv);
+      q = (x - v)*(fval - fw);
+      p = (x - v)*q - (x - w)*r;
+      q = 2*(q - r);
+      p *= -sign (q);
+      q = abs (q);
+      r = e;
+      e = d;
+
+      if (abs (p) < abs (0.5*q*r) && p > q*(a-x) && p < q*(b-x))
+        ## The parabolic step is acceptable.
+        d = p / q;
+        u = x + d;
+
+        ## f must not be evaluated too close to ax or bx.
+        if ((u-a) < 2*tol && (b-u) < 2*tol)
+          d = tol * (sign (xm - x) + (xm == x));
+        endif
+      else
+        dogs = true;
+      endif
+    else
+      dogs = true;
+    endif
+    if (dogs)
+      ## Default to golden section step.
+      e = ifelse (x >= xm, a - x, b - x);
+      d = c * e;
+    endif
+
+     ## f must not be evaluated too close to x.
+     u = x + max (abs (d), tol) * (sign (d) + (d == 0));
+     
+     fu = fun (u);
+     nfev++;
+     niter++;
+
+     ## update  a, b, v, w, and x
+
+     if (fu <= fval)
+       if (u < x)
+         b = x;
+       else
+         a = x;
+       endif
+       v = w; fv = fw;
+       w = x; fw = fval;
+       x = u; fval = fu;
+     else
+       ## The following if-statement was originally executed even if fu == fval.
+       if (u < x)
+         a = u;
+       else
+         b = u;
+       endif
+       if (fu <= fw || w == x)
+         v = w; fv = fw;
+         w = u; fw = fu;
+       elseif (fu <= fv || v == x || v == w)
+         v = u;
+         fv = fu;
+       endif
+     endif
+
+    ## If there's an output function, use it now.
+    if (outfcn)
+      optv.funccount = nfev;
+      optv.fval = fval;
+      optv.iteration = niter;
+      if (outfcn (x, optv, "iter"))
+        info = -1;
+        break;
+      endif
+    endif
+  endwhile
+
+  output.iterations = niter;
+  output.funcCount = nfev;
+  output.bracket = [a, b];
+  ## FIXME: bracketf possibly unavailable.
+
+endfunction
+
+## An assistant function that evaluates a function handle and checks for
+## bad results.
+function fx = guarded_eval (fun, x)
+  fx = fun (x);
+  fx = fx(1);
+  if (! isreal (fx))
+    error ("fminbnd:notreal", "fminbnd: non-real value encountered"); 
+  elseif (isnan (fx))
+    error ("fminbnd:isnan", "fminbnd: NaN value encountered"); 
+  endif
+endfunction
+
+%!shared opt0
+%! opt0 = optimset ("tolx", 0);
+%!assert (fminbnd (@cos, pi/2, 3*pi/2, opt0), pi, 10*eps)
+%!assert (fminbnd (@(x) (x - 1e-3)^4, -1, 1, opt0), 1e-3, 10e-3*eps)
+%!assert (fminbnd (@(x) abs(x-1e7), 0, 1e10, opt0), 1e7, 10e7*eps)
+%!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), 1e-7)
--- a/scripts/optimization/fzero.m	Wed Feb 10 08:41:49 2010 +0100
+++ b/scripts/optimization/fzero.m	Wed Feb 10 13:33:54 2010 +0100
@@ -65,7 +65,7 @@
 
   ## Get default options if requested.
   if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults'))
-    x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 0, \
+    x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, \
     "OutputFcn", [], "FunValCheck", "off");
     return;
   endif
@@ -82,7 +82,7 @@
   ## displev = optimget (options, "Display", "notify");
   funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
   outfcn = optimget (options, "OutputFcn");
-  tolx = optimget (options, "TolX", 0);
+  tolx = optimget (options, "TolX", 1e-8);
   maxiter = optimget (options, "MaxIter", Inf);
   maxfev = optimget (options, "MaxFunEvals", Inf);
 
@@ -99,6 +99,7 @@
   nfev = 0;
 
   x = fval = a = fa = b = fb = NaN;
+  eps = eps (class (x0));
 
   ## Prepare...
   a = x0(1);
@@ -257,7 +258,7 @@
 
     ## If there's an output function, use it now.
     if (outfcn)
-      optv.funccount = niter + 2;
+      optv.funccount = nfev;
       optv.fval = fval;
       optv.iteration = niter;
       if (outfcn (x, optv, "iter"))
@@ -286,7 +287,7 @@
   endwhile
 
   output.iterations = niter;
-  output.funcCount = niter + 2;
+  output.funcCount = nfev;
   output.bracket = [a, b];
   output.bracketf = [fa, fb];
 
@@ -304,5 +305,7 @@
   endif
 endfunction
 
-%!assert(fzero(@cos, [0, 3]), pi/2, 10*eps)
-%!assert(fzero(@(x) x^(1/3) - 1e-8, [0,1]), 1e-24, 1e-22*eps)
+%!shared opt0
+%! opt0 = optimset ("tolx", 0);
+%!assert(fzero(@cos, [0, 3], opt0), pi/2, 10*eps)
+%!assert(fzero(@(x) x^(1/3) - 1e-8, [0,1], opt0), 1e-24, 1e-22*eps)
--- a/scripts/optimization/module.mk	Wed Feb 10 08:41:49 2010 +0100
+++ b/scripts/optimization/module.mk	Wed Feb 10 13:33:54 2010 +0100
@@ -5,6 +5,7 @@
 
 optimization_FCN_FILES = \
   optimization/__all_opts__.m \
+  optimization/fminbnd.m \
   optimization/fminunc.m \
   optimization/fsolve.m \
   optimization/fzero.m \