1########################################################################
2##
3## Copyright (C) 2008-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {@var{x} =} fminbnd (@var{fun}, @var{a}, @var{b})
28## @deftypefnx {} {@var{x} =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options})
29## @deftypefnx {} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@dots{})
30## Find a minimum point of a univariate function.
31##
32## @var{fun} is a function handle, inline function, or string containing the
33## name of the function to evaluate.
34##
35## The starting interval is specified by @var{a} (left boundary) and @var{b}
36## (right boundary).  The endpoints must be finite.
37##
38## @var{options} is a structure specifying additional parameters which
39## control the algorithm.  Currently, @code{fminbnd} recognizes these options:
40## @qcode{"Display"}, @qcode{"FunValCheck"}, @qcode{"MaxFunEvals"},
41## @qcode{"MaxIter"}, @qcode{"OutputFcn"}, @qcode{"TolX"}.
42##
43## @qcode{"MaxFunEvals"} proscribes the maximum number of function evaluations
44## before optimization is halted.  The default value is 500.
45## The value must be a positive integer.
46##
47## @qcode{"MaxIter"} proscribes the maximum number of algorithm iterations
48## before optimization is halted.  The default value is 500.
49## The value must be a positive integer.
50##
51## @qcode{"TolX"} specifies the termination tolerance for the solution @var{x}.
52## The default is @code{1e-4}.
53##
54## For a description of the other options, see @ref{XREFoptimset,,optimset}.
55## To initialize an options structure with default values for @code{fminbnd}
56## use @code{options = optimset ("fminbnd")}.
57##
58## On exit, the function returns @var{x}, the approximate minimum point, and
59## @var{fval}, the function evaluated @var{x}.
60##
61## The third output @var{info} reports whether the algorithm succeeded and may
62## take one of the following values:
63##
64## @itemize
65## @item 1
66## The algorithm converged to a solution.
67##
68## @item 0
69## Iteration limit (either @code{MaxIter} or @code{MaxFunEvals}) exceeded.
70##
71## @item -1
72## The algorithm was terminated by a user @code{OutputFcn}.
73## @end itemize
74##
75## Programming Notes: The search for a minimum is restricted to be in the
76## finite interval bound by @var{a} and @var{b}.  If you have only one initial
77## point to begin searching from then you will need to use an unconstrained
78## minimization algorithm such as @code{fminunc} or @code{fminsearch}.
79## @code{fminbnd} internally uses a Golden Section search strategy.
80## @seealso{fzero, fminunc, fminsearch, optimset}
81## @end deftypefn
82
83## This is patterned after opt/fmin.f from Netlib, which in turn is taken from
84## Richard Brent: Algorithms For Minimization Without Derivatives,
85## Prentice-Hall (1973)
86
87## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
88## PKG_ADD: [~] = __all_opts__ ("fminbnd");
89
90function [x, fval, info, output] = fminbnd (fun, a, b, options = struct ())
91
92  ## Get default options if requested.
93  if (nargin == 1 && ischar (fun) && strcmp (fun, "defaults"))
94    x = struct ("Display", "notify", "FunValCheck", "off",
95                "MaxFunEvals", 500, "MaxIter", 500,
96                "OutputFcn", [], "TolX", 1e-4);
97    return;
98  endif
99
100  if (nargin < 2 || nargin > 4)
101    print_usage ();
102  endif
103
104  if (a > b)
105    error ("Octave:invalid-input-arg",
106           "fminbnd: the lower bound cannot be greater than the upper one");
107  endif
108
109  if (ischar (fun))
110    fun = str2func (fun);
111  endif
112
113  displ = optimget (options, "Display", "notify");
114  funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
115  outfcn = optimget (options, "OutputFcn");
116  tolx = optimget (options, "TolX", 1e-4);
117  maxiter = optimget (options, "MaxIter", 500);
118  maxfev = optimget (options, "MaxFunEvals", 500);
119
120  if (funvalchk)
121    ## Replace fun with a guarded version.
122    fun = @(x) guarded_eval (fun, x);
123  endif
124
125  ## The default exit flag if exceeded number of iterations.
126  info = 0;
127  niter = 0;
128  nfev = 0;
129
130  c = 0.5*(3 - sqrt (5));
131  v = a + c*(b-a);
132  w = x = v;
133  e = 0;
134  fv = fw = fval = fun (x);
135  nfev += 1;
136
137  if (isa (a, "single") || isa (b, "single") || isa (fval, "single"))
138    sqrteps = eps ("single");
139  else
140    sqrteps = eps ("double");
141  endif
142
143  ## Only for display purposes.
144  iter(1).funccount = nfev;
145  iter(1).x = x;
146  iter(1).fx = fval;
147
148  while (niter < maxiter && nfev < maxfev)
149    xm = 0.5*(a+b);
150    ## FIXME: the golden section search can actually get closer than sqrt(eps)
151    ## sometimes.  Sometimes not, it depends on the function.  This is the
152    ## strategy from the Netlib code.  Something smarter would be good.
153    tol = 2 * sqrteps * abs (x) + tolx / 3;
154    if (abs (x - xm) <= (2*tol - 0.5*(b-a)))
155      info = 1;
156      break;
157    endif
158
159    if (abs (e) > tol)
160      dogs = false;
161      ## Try inverse parabolic step.
162      iter(niter+1).procedure = "parabolic";
163
164      r = (x - w)*(fval - fv);
165      q = (x - v)*(fval - fw);
166      p = (x - v)*q - (x - w)*r;
167      q = 2*(q - r);
168      p *= -sign (q);
169      q = abs (q);
170      r = e;
171      e = d;
172
173      if (abs (p) < abs (0.5*q*r) && p > q*(a-x) && p < q*(b-x))
174        ## The parabolic step is acceptable.
175        d = p / q;
176        u = x + d;
177
178        ## f must not be evaluated too close to ax or bx.
179        if (min (u-a, b-u) < 2*tol)
180          d = tol * (sign (xm - x) + (xm == x));
181        endif
182      else
183        dogs = true;
184      endif
185    else
186      dogs = true;
187    endif
188    if (dogs)
189      ## Default to golden section step.
190
191      ## WARNING: This is also the "initial" procedure following MATLAB
192      ## nomenclature.  After the loop we'll fix the string for the first step.
193      iter(niter+1).procedure = "golden";
194
195      e = ifelse (x >= xm, a - x, b - x);
196      d = c * e;
197    endif
198
199    ## f must not be evaluated too close to x.
200    u = x + max (abs (d), tol) * (sign (d) + (d == 0));
201    fu = fun (u);
202
203    niter += 1;
204
205    iter(niter).funccount = nfev++;
206    iter(niter).x = u;
207    iter(niter).fx = fu;
208
209    ## update a, b, v, w, and x
210
211    if (fu < fval)
212      if (u < x)
213        b = x;
214      else
215        a = x;
216      endif
217      v = w; fv = fw;
218      w = x; fw = fval;
219      x = u; fval = fu;
220    else
221      ## The following if-statement was originally executed even if fu == fval.
222      if (u < x)
223        a = u;
224      else
225        b = u;
226      endif
227      if (fu <= fw || w == x)
228        v = w; fv = fw;
229        w = u; fw = fu;
230      elseif (fu <= fv || v == x || v == w)
231        v = u;
232        fv = fu;
233      endif
234    endif
235
236    ## If there's an output function, use it now.
237    if (! isempty (outfcn))
238      optv.funccount = nfev;
239      optv.fval = fval;
240      optv.iteration = niter;
241      if (outfcn (x, optv, "iter"))
242        info = -1;
243        break;
244      endif
245    endif
246  endwhile
247
248  ## Fix the first step procedure.
249  iter(1).procedure = "initial";
250
251  ## Handle the "Display" option
252  switch (displ)
253    case "iter"
254      print_formatted_table (iter);
255      print_exit_msg (info, struct ("TolX", tolx, "fx", fval));
256    case "notify"
257      if (info == 0)
258        print_exit_msg (info, struct ("fx",fval));
259      endif
260    case "final"
261      print_exit_msg (info, struct ("TolX", tolx, "fx", fval));
262    case "off"
263      "skip";
264    otherwise
265      warning ("fminbnd: unknown option for Display: '%s'", displ);
266  endswitch
267
268  output.iterations = niter;
269  output.funcCount = nfev;
270  output.algorithm = "golden section search, parabolic interpolation";
271  output.bracket = [a, b];
272  ## FIXME: bracketf possibly unavailable.
273
274endfunction
275
276## A helper function that evaluates a function and checks for bad results.
277function fx = guarded_eval (fun, x)
278  fx = fun (x);
279  fx = fx(1);
280  if (! isreal (fx))
281    error ("Octave:fmindbnd:notreal", "fminbnd: non-real value encountered");
282  elseif (isnan (fx))
283    error ("Octave:fmindbnd:isnan", "fminbnd: NaN value encountered");
284  endif
285endfunction
286
287## A hack for printing a formatted table
288function print_formatted_table (table)
289  printf ("\n Func-count     x          f(x)         Procedure\n");
290  for row=table
291    printf("%5.5s        %7.7s    %8.8s\t%s\n",
292           int2str (row.funccount), num2str (row.x,"%.5f"),
293           num2str (row.fx,"%.6f"), row.procedure);
294  endfor
295  printf ("\n");
296endfunction
297
298## Print either a success termination message or bad news
299function print_exit_msg (info, opt=struct())
300  printf ("");
301  switch (info)
302    case 1
303      printf ("Optimization terminated:\n");
304      printf (" the current x satisfies the termination criteria using OPTIONS.TolX of %e\n", opt.TolX);
305    case 0
306      printf ("Exiting: Maximum number of iterations has been exceeded\n");
307      printf ("         - increase MaxIter option.\n");
308      printf ("         Current function value: %.6f\n", opt.fx);
309    case -1
310      "FIXME"; # FIXME: what's the message MATLAB prints for this case?
311    otherwise
312      error ("fminbnd: internal error, info return code was %d", info);
313  endswitch
314  printf ("\n");
315endfunction
316
317
318%!shared opt0
319%! opt0 = optimset ("tolx", 0);
320%!assert (fminbnd (@cos, pi/2, 3*pi/2, opt0), pi, 10*sqrt (eps))
321%!assert (fminbnd (@(x) (x - 1e-3)^4, -1, 1, opt0), 1e-3, 10e-3*sqrt (eps))
322%!assert (fminbnd (@(x) abs (x-1e7), 0, 1e10, opt0), 1e7, 10e7*sqrt (eps))
323%!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))
324%!assert (fminbnd (@(x) x > 0.3, 0, 1) < 0.3)
325%!assert (fminbnd (@(x) sin (x), 0, 0), 0, eps)
326
327%!error <lower bound cannot be greater> fminbnd (@(x) sin (x), 0, -pi)
328