# HG changeset patch
# User mmuetzel <markus.muetzel@gmx.de>
# Date 1454502734 -3600
#      Wed Feb 03 13:32:14 2016 +0100
# Node ID c06087af9987fc0cae283f75559c2709024f5b23
# Parent  a428c6ba63344734b9ea045928738aeb62af2740
fminsearch.m:
Add support for "exitflag" output variable.
Add support for "output" output variable.
For Display option, add 'final' and 'notify' options.
Add support for OutputFcn.
Add support for exiting based on TolFun.

diff -r a428c6ba6334 -r c06087af9987 scripts/optimization/fminsearch.m
--- a/scripts/optimization/fminsearch.m	Tue Feb 02 23:51:45 2016 -0500
+++ b/scripts/optimization/fminsearch.m	Wed Feb 03 13:32:14 2016 +0100
@@ -1,5 +1,6 @@
 ## Copyright (C) 2003,2012 Andy Adler
 ## Copyright (C) 2002, 2013 N.J.Higham
+## Copyright (C) 2016 M. Muetzel
 ##
 ## This file is part of Octave.
 ##
@@ -21,22 +22,39 @@
 ## @deftypefn  {} {@var{x} =} fminsearch (@var{fun}, @var{x0})
 ## @deftypefnx {} {@var{x} =} fminsearch (@var{fun}, @var{x0}, @var{options})
 ## @deftypefnx {} {[@var{x}, @var{fval}] =} fminsearch (@dots{})
+## @deftypefnx {} {[@var{x}, @var{fval}, @var{exitflag}] =} fminsearch (@dots{})
+## @deftypefnx {} {[@var{x}, @var{fval}, @var{exitflag}, @var{output}] =} fminsearch (@dots{})
 ##
 ## Find a value of @var{x} which minimizes the function @var{fun}.
 ##
 ## The search begins at the point @var{x0} and iterates using the
 ## @nospell{Nelder & Mead} Simplex algorithm (a derivative-free method).  This
 ## algorithm is better-suited to functions which have discontinuities or for
-## which a gradient-based search such as @code{fminunc} fails.
+## which a gradient-based search such as @command{fminunc} fails.
 ##
 ## Options for the search are provided in the parameter @var{options} using the
-## function @code{optimset}.  Currently, @code{fminsearch} accepts the options:
-## @qcode{"TolX"}, @qcode{"MaxFunEvals"}, @qcode{"MaxIter"}, @qcode{"Display"}.
-## For a description of these options, see @code{optimset}.
+## function @command{optimset}.  Currently, @code{fminsearch} accepts the
+## options: @qcode{"TolX"}, @qcode{"TolFun"}, @qcode{"MaxFunEvals"},
+## @qcode{"MaxIter"}, @qcode{"Display"} and @qcode{"OutputFcn"}.
+## For a description of these options, see @command{optimset}.
 ##
 ## On exit, the function returns @var{x}, the minimum point, and @var{fval},
 ## the function value thereof.
 ##
+## Optionally, the function returns @var{exitflag}. @var{exitflag} is 1 if the
+## size of the simplex is smaller than @var{options}.@qcode{TolX} AND the step
+## in the function value between iterations is smaller than
+## @var{options}.@qcode{TolFun}. It is -1 if the iteration is stopped by the
+## @qcode{"OutputFcn"}, and it is 0 if the maximum number of iterations or the
+## maximum number of function evaluations are exceeded.
+##
+## Optionally, the structure @var{output} is returned with the fields,
+## @code{funcCount} containing the number of function calls to @var{fun},
+## @code{iterations} containing the number of iteration steps,
+## @code{algorithm} with the name of the used search algorithm (always: 
+## @qcode{"Nelder-Mead simplex direct search"}) and @code{message} with the
+## exit message.
+##
 ## Example usages:
 ##
 ## @example
@@ -52,13 +70,9 @@
 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
 ## PKG_ADD: [~] = __all_opts__ ("fminsearch");
 
-## FIXME: Add support for "exitflag" output variable.
-## FIXME: Add support for "output" output variable.
-## FIXME: For Display option, add 'final' and 'notify' options.  Not too hard.
-## FIXME: Add support for OutputFcn.  See fminunc for a template.
-## FIXME: Add support for exiting based on TolFun.  See fminunc for an idea.
+## FIXME: Add support for output function with "state" set to "interrupt".
 
-function [x, fval] = fminsearch (fun, x0, options = struct ())
+function [x, fval, exitflag, output] = fminsearch (fun, x0, options, varargin)
 
   ## Get default options if requested.
   if (nargin == 1 && ischar (fun) && strcmp (fun, "defaults"))
@@ -69,11 +83,15 @@
     return;
   endif
 
-  if (nargin < 2 || nargin > 3)
+  if (nargin < 2)
     print_usage ();
   endif
+  
+  if (nargin < 3 || isempty (options))
+    options = struct ();
+  endif
 
-  x = nmsmax (fun, x0, options);
+  [x, exitflag, output] = nmsmax (fun, x0, options, varargin{:});
 
   if (isargout (2))
     fval = feval (fun, x);
@@ -121,11 +139,14 @@
 ##
 ## Modifications for Octave by A.Adler 2003
 
-function [stopit, savit, dirn, trace, tol, maxiter] = parse_options (options, x );
+function [stopit, savit, dirn, trace, tol, maxiter, tol_f, outfcn] = parse_options (options, x);
 
   ## Tolerance for cgce test based on relative size of simplex.
   stopit(1) = tol = optimget (options, "TolX", 1e-4);
 
+  ## Tolerance for cgce test based on step in function value
+  tol_f = optimget (options, "TolFun", 1e-7);
+
   ## Max no. of f-evaluations.
   stopit(2) = optimget (options, "MaxFunEvals", length (x) * 200);
 
@@ -133,6 +154,7 @@
   maxiter = optimget (options, "MaxIter", length (x) * 200);
 
   ## Default target for f-values.
+  ## FIXME: This can be handled by an OutputFcn. Remove stopit(3)?
   stopit(3) = Inf;  # FIXME: expose this parameter to the outside
 
   ## Default initial simplex.
@@ -140,11 +162,16 @@
 
   ## Default: show progress.
   display = optimget (options, "Display", "notify");
-  if (strcmp (display, "iter"))
-    stopit(5) = 1;
-  else
-    stopit(5) = 0;
-  endif
+  switch (display)
+    case "iter"
+      stopit(5) = 1;
+    case "final"
+      stopit(5) = 2;
+    case "notify"
+      stopit(5) = 3;
+    otherwise ## "none"
+      stopit(5) = 0;
+  endswitch
   trace = stopit(5);
 
   ## Use function to minimize, not maximize
@@ -153,11 +180,14 @@
   ## Filename for snapshots.
   savit = [];  # FIXME: expose this parameter to the outside
 
+  ## OutputFcn
+  outfcn = optimget (options, "OutputFcn");
+
 endfunction
 
-function [x, fmax, nf] = nmsmax (fun, x, options, savit, varargin)
+function [x, exitflag, output] = nmsmax (fun, x, options, varargin)
 
-  [stopit, savit, dirn, trace, tol, maxiter] = parse_options (options, x);
+  [stopit, savit, dirn, trace, tol, maxiter, tol_f, outfcn] = parse_options (options, x);
 
   if (strcmpi (optimget (options, "FunValCheck", "off"), "on"))
     ## Replace fcn with a guarded version.
@@ -170,47 +200,61 @@
   V = [zeros(n,1) eye(n)];
   f = zeros (n+1,1);
   V(:,1) = x0;
-  f(1) = dirn * feval (fun,x,varargin{:});
+  f(1) = dirn * feval (fun, x, varargin{:});
   fmax_old = f(1);
 
-  if (trace)
+  if (trace == 1)
     fprintf ("f(x0) = %9.4e\n", f(1));
   endif
 
   k = 0; m = 0;
 
   ## Set up initial simplex.
-  scale = max (norm (x0,Inf), 1);
+  scale = max (norm (x0, Inf), 1);
   if (stopit(4) == 0)
     ## Regular simplex - all edges have same length.
     ## Generated from construction given in reference [18, pp. 80-81] of [1].
     alpha = scale / (n*sqrt (2)) * [sqrt(n+1)-1+n, sqrt(n+1)-1];
-    V(:,2:n+1) = (x0 + alpha(2)*ones (n,1)) * ones (1,n);
+    V(:,2:n+1) = (x0 + alpha(2)*ones (n, 1)) * ones (1, n);
     for j = 2:n+1
       V(j-1,j) = x0(j-1) + alpha(1);
       x(:) = V(:,j);
-      f(j) = dirn * feval (fun,x,varargin{:});
+      f(j) = dirn * feval (fun, x, varargin{:});
     endfor
   else
     ## Right-angled simplex based on co-ordinate axes.
-    alpha = scale * ones(n+1,1);
+    alpha = scale * ones(n+1, 1);
     for j=2:n+1
       V(:,j) = x0 + alpha(j)*V(:,j);
       x(:) = V(:,j);
-      f(j) = dirn * feval (fun,x,varargin{:});
+      f(j) = dirn * feval (fun, x, varargin{:});
     endfor
   endif
-  nf = n+1;
-  how = "initial  ";
+  nf = n+j;
+  how = "initial";
 
   [~,j] = sort (f);
   j = j(n+1:-1:1);
   f = f(j);
   V = V(:,j);
 
+  exitflag = 0;
+
+  if (! isempty (outfcn))
+    optimvalues.iteration = 0;
+    optimvalues.funccount = nf;
+    optimvalues.fval = f(1);
+    optimvalues.procedure = how;
+    state = "init";
+    stop = outfcn (x, optimvalues, state);
+    if (stop)
+      exitflag = -1;
+    endif
+  endif
+
   alpha = 1;  beta = 1/2;  gamma = 2;
 
-  while (1)   # Outer (and only) loop.
+  while (exitflag != -1)   # Outer (and only) loop.
     k += 1;
 
     if (k > maxiter)
@@ -225,9 +269,9 @@
         eval (["save " savit " x fmax nf"]);
       endif
     endif
-    if (trace)
+    if (trace == 1)
       fprintf ("Iter. %2.0f,", k);
-      fprintf (["  how = " how "  "]);
+      fprintf ("  how = %-11s", [how ","]);
       fprintf ("nf = %3.0f,  f = %9.4e  (%2.1f%%)\n", nf, fmax, ...
                100*(fmax-fmax_old)/(abs(fmax_old)+eps));
     endif
@@ -238,6 +282,7 @@
     ## Stopping Test 1 - f reached target value?
     if (fmax >= stopit(3))
       msg = "Exceeded target...quitting\n";
+      exitflag = -1; ## FIXME: Add docu when stopit(3) gets exposed to the outside
       break;
     endif
 
@@ -247,15 +292,31 @@
       break;
     endif
 
-    ## Stopping Test 3 - converged?   This is test (4.3) in [1].
+    ## Stopping Test 3 - converged?   The first part is test (4.3) in [1].
     v1 = V(:,1);
-    size_simplex = norm (V(:,2:n+1)-v1(:,ones (1,n)),1) / max (1, norm (v1,1));
-    if (size_simplex <= tol)
-      msg = sprintf ("Simplex size %9.4e <= %9.4e...quitting\n", ...
-                      size_simplex, tol);
+    size_simplex = norm (V(:,2:n+1)-v1(:,ones (1, n)),1) / max (1, norm (v1, 1));
+    step_f = max (abs (f(1) - f(2:n+1)));
+    if (size_simplex <= tol && step_f <= tol_f )
+      msg = sprintf ("Simplex size %9.4e <= %9.4e and step in function value %9.4e <= %9.4e...quitting\n", ...
+                      size_simplex, tol, step_f, tol_f);
+      exitflag = 1;
       break;
     endif
 
+    ## Call OutputFcn
+    if (! isempty (outfcn))
+      optimvalues.funccount = nf;
+      optimvalues.fval = f(1);
+      optimvalues.iteration = k;
+      optimvalues.procedure = how;
+      state = "iter";
+      stop = outfcn (x, optimvalues, state);
+      if (stop)
+        exitflag = -1;
+        break;
+      endif
+    endif
+
     ##  One step of the Nelder-Mead simplex algorithm
     ##  NJH: Altered function calls and changed CNT to NF.
     ##       Changed each 'fr < f(1)' type test to '>' for maximization
@@ -264,19 +325,19 @@
     vbar = (sum (V(:,1:n)')/n)';  # Mean value
     vr = (1 + alpha)*vbar - alpha*V(:,n+1);
     x(:) = vr;
-    fr = dirn * feval (fun,x,varargin{:});
+    fr = dirn * feval (fun, x, varargin{:});
     nf += 1;
-    vk = vr;  fk = fr; how = "reflect, ";
+    vk = vr;  fk = fr; how = "reflect";
     if (fr > f(n))
       if (fr > f(1))
         ve = gamma*vr + (1-gamma)*vbar;
         x(:) = ve;
-        fe = dirn * feval (fun,x,varargin{:});
+        fe = dirn * feval (fun, x, varargin{:});
         nf += 1;
         if (fe > f(1))
           vk = ve;
           fk = fe;
-          how = "expand,  ";
+          how = "expand";
         endif
       endif
     else
@@ -288,23 +349,23 @@
       endif
       vc = beta*vt + (1-beta)*vbar;
       x(:) = vc;
-      fc = dirn * feval (fun,x,varargin{:});
+      fc = dirn * feval (fun, x, varargin{:});
       nf += 1;
       if (fc > f(n))
         vk = vc; fk = fc;
-        how = "contract,";
+        how = "contract";
       else
         for j = 2:n
           V(:,j) = (V(:,1) + V(:,j))/2;
           x(:) = V(:,j);
-          f(j) = dirn * feval (fun,x,varargin{:});
+          f(j) = dirn * feval (fun, x, varargin{:});
         endfor
         nf += n-1;
         vk = (V(:,1) + V(:,n+1))/2;
         x(:) = vk;
-        fk = dirn * feval (fun,x,varargin{:});
+        fk = dirn * feval (fun, x, varargin{:});
         nf += 1;
-        how = "shrink,  ";
+        how = "shrink";
       endif
     endif
     V(:,n+1) = vk;
@@ -317,11 +378,27 @@
   endwhile   # End of outer (and only) loop.
 
   ## Finished.
-  if (trace)
+  if ( (trace == 1) || (trace == 2) || (trace == 3 && exitflag != 1) )
     fprintf (msg);
   endif
   x(:) = V(:,1);
 
+  ## FIXME: should outputfcn be called only after a successful iteration?
+  if (! isempty (outfcn))
+    optimvalues.funccount = nf;
+    optimvalues.fval = f(1);
+    optimvalues.iteration = k;
+    optimvalues.procedure = how;
+    state = "done";
+    outfcn (x, optimvalues, state);
+  endif
+
+  ## output
+  output.iterations = k;
+  output.funcCount = nf;
+  output.algorithm = "Nelder-Mead simplex direct search";
+  output.message = sprintf (msg);
+
 endfunction
 
 ## A helper function that evaluates a function and checks for bad results.
@@ -347,7 +424,11 @@
 
 %!assert (fminsearch (@sin, 3, optimset ("MaxIter", 3)), 4.8750, 1e-4)
 %!assert (fminsearch (@sin, 3, optimset ("MaxIter", 30)), 4.7124, 1e-4)
+%!assert (fminsearch (@sin, 3, optimset ("MaxFunEval", 10)), 4.8750, 1e-4)
 %!shared c
 %! c = 1.5;
-%!assert (fminsearch (@(x) x(1).^2+c*x(2).^2,[1;1]), [0;0], 1e-4)
+%!assert (fminsearch (@(x) x(1).^2+c*x(2).^2, [1;1]), [0;0], 1e-4)
+%!assert (fminsearch (@(x, a) x(1).^2+a*x(2).^2, [1;1], [], 1.5), [0;0], 1e-4)
+%!assert (fminsearch (@(x, a) x(1).^2+a*x(2).^2, [1;1], optimset ("TolX", 1e-5), 1.5), [0;0], 1e-5)
+%!assert (fminsearch (@(x, a) x(1).^2+a*x(2).^2, [1;1], optimset ("TolFun", 1e-9), 1.5), [0;0], 1e-5)
 
