0001 ########################################################################
0002 ##
0003 ## Copyright (C) 2008-2022 The Octave Project Developers
0004 ##
0005 ## See the file COPYRIGHT.md in the top-level directory of this
0006 ## distribution or <https://octave.org/copyright/>.
0007 ##
0008 ## This file is part of Octave.
0009 ##
0010 ## Octave is free software: you can redistribute it and/or modify it
0011 ## under the terms of the GNU General Public License as published by
0012 ## the Free Software Foundation, either version 3 of the License, or
0013 ## (at your option) any later version.
0014 ##
0015 ## Octave is distributed in the hope that it will be useful, but
0016 ## WITHOUT ANY WARRANTY; without even the implied warranty of
0017 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0018 ## GNU General Public License for more details.
0019 ##
0020 ## You should have received a copy of the GNU General Public License
0021 ## along with Octave; see the file COPYING. If not, see
0022 ## <https://www.gnu.org/licenses/>.
0023 ##
0024 ########################################################################
0025
0026 ## -*- texinfo -*-
0027 ## @deftypefn {} {@var{x} =} fminbnd (@var{fun}, @var{a}, @var{b})
0028 ## @deftypefnx {} {@var{x} =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options})
0029 ## @deftypefnx {} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@dots{})
0030 ## Find a minimum point of a univariate function.
0031 ##
0032 ## @var{fun} is a function handle, inline function, or string containing the
0033 ## name of the function to evaluate.
0034 ##
0035 ## The starting interval is specified by @var{a} (left boundary) and @var{b}
0036 ## (right boundary). The endpoints must be finite.
0037 ##
0038 ## @var{options} is a structure specifying additional parameters which
0039 ## control the algorithm. Currently, @code{fminbnd} recognizes these options:
0040 ## @qcode{"Display"}, @qcode{"FunValCheck"}, @qcode{"MaxFunEvals"},
0041 ## @qcode{"MaxIter"}, @qcode{"OutputFcn"}, @qcode{"TolX"}.
0042 ##
0043 ## @qcode{"MaxFunEvals"} proscribes the maximum number of function evaluations
0044 ## before optimization is halted. The default value is 500.
0045 ## The value must be a positive integer.
0046 ##
0047 ## @qcode{"MaxIter"} proscribes the maximum number of algorithm iterations
0048 ## before optimization is halted. The default value is 500.
0049 ## The value must be a positive integer.
0050 ##
0051 ## @qcode{"TolX"} specifies the termination tolerance for the solution @var{x}.
0052 ## The default is @code{1e-4}.
0053 ##
0054 ## For a description of the other options,
0055 ## @pxref{XREFoptimset,,@code{optimset}}.
0056 ## To initialize an options structure with default values for @code{fminbnd}
0057 ## use @code{options = optimset ("fminbnd")}.
0058 ##
0059 ## On exit, the function returns @var{x}, the approximate minimum point, and
0060 ## @var{fval}, the function evaluated @var{x}.
0061 ##
0062 ## The third output @var{info} reports whether the algorithm succeeded and may
0063 ## take one of the following values:
0064 ##
0065 ## @itemize
0066 ## @item 1
0067 ## The algorithm converged to a solution.
0068 ##
0069 ## @item 0
0070 ## Iteration limit (either @code{MaxIter} or @code{MaxFunEvals}) exceeded.
0071 ##
0072 ## @item -1
0073 ## The algorithm was terminated by a user @code{OutputFcn}.
0074 ## @end itemize
0075 ##
0076 ## Programming Notes: The search for a minimum is restricted to be in the
0077 ## finite interval bound by @var{a} and @var{b}. If you have only one initial
0078 ## point to begin searching from then you will need to use an unconstrained
0079 ## minimization algorithm such as @code{fminunc} or @code{fminsearch}.
0080 ## @code{fminbnd} internally uses a Golden Section search strategy.
0081 ## @seealso{fzero, fminunc, fminsearch, optimset}
0082 ## @end deftypefn
0083
0084 ## This is patterned after opt/fmin.f from Netlib, which in turn is taken from
0085 ## Richard Brent: Algorithms For Minimization Without Derivatives,
0086 ## Prentice-Hall (1973)
0087
0088 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
0089 ## PKG_ADD: [~] = __all_opts__ ("fminbnd");
0090
0091 ## Added ability to pass extra params to fn
0092 ## A Adler, Dec 2022
0093
0094 function [x, fval, info, output] = fminbnd (fun, a, b, options = struct (), varargin = {})
0095
0096 ## Get default options if requested.
0097 if (nargin == 1 && ischar (fun) && strcmp (fun, "defaults"))
0098 x = struct ("Display", "notify", "FunValCheck", "off",
0099 "MaxFunEvals", 500, "MaxIter", 500,
0100 "OutputFcn", [], "TolX", 1e-4);
0101 return;
0102 endif
0103
0104 if (nargin < 2)
0105 print_usage ();
0106 endif
0107
0108 if (a > b)
0109 error ("Octave:invalid-input-arg",
0110 "fminbnd: the lower bound cannot be greater than the upper one");
0111 endif
0112
0113 if (ischar (fun))
0114 fun = str2func (fun);
0115 endif
0116
0117 displ = optimget (options, "Display", "notify");
0118 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
0119 outfcn = optimget (options, "OutputFcn");
0120 tolx = optimget (options, "TolX", 1e-4);
0121 maxiter = optimget (options, "MaxIter", 500);
0122 maxfev = optimget (options, "MaxFunEvals", 500);
0123
0124 if (funvalchk)
0125 ## Replace fun with a guarded version.
0126 fun = @(x) guarded_eval (fun, x, varargin{:});
0127 endif
0128
0129 ## The default exit flag if exceeded number of iterations.
0130 info = 0;
0131 niter = 0;
0132 nfev = 0;
0133
0134 c = 0.5*(3 - sqrt (5));
0135 v = a + c*(b-a);
0136 w = x = v;
0137 e = 0;
0138 fv = fw = fval = fun (x, varargin{:});
0139 nfev += 1;
0140
0141 if (isa (a, "single") || isa (b, "single") || isa (fval, "single"))
0142 sqrteps = eps ("single");
0143 else
0144 sqrteps = eps ("double");
0145 endif
0146
0147 ## Only for display purposes.
0148 iter(1).funccount = nfev;
0149 iter(1).x = x;
0150 iter(1).fx = fval;
0151
0152 while (niter < maxiter && nfev < maxfev)
0153 xm = 0.5*(a+b);
0154 ## FIXME: the golden section search can actually get closer than sqrt(eps)
0155 ## sometimes. Sometimes not, it depends on the function. This is the
0156 ## strategy from the Netlib code. Something smarter would be good.
0157 tol = 2 * sqrteps * abs (x) + tolx / 3;
0158 if (abs (x - xm) <= (2*tol - 0.5*(b-a)))
0159 info = 1;
0160 break;
0161 endif
0162
0163 if (abs (e) > tol)
0164 dogs = false;
0165 ## Try inverse parabolic step.
0166 iter(niter+1).procedure = "parabolic";
0167
0168 r = (x - w)*(fval - fv);
0169 q = (x - v)*(fval - fw);
0170 p = (x - v)*q - (x - w)*r;
0171 q = 2*(q - r);
0172 p *= -sign (q);
0173 q = abs (q);
0174 r = e;
0175 e = d;
0176
0177 if (abs (p) < abs (0.5*q*r) && p > q*(a-x) && p < q*(b-x))
0178 ## The parabolic step is acceptable.
0179 d = p / q;
0180 u = x + d;
0181
0182 ## f must not be evaluated too close to ax or bx.
0183 if (min (u-a, b-u) < 2*tol)
0184 d = tol * (sign (xm - x) + (xm == x));
0185 endif
0186 else
0187 dogs = true;
0188 endif
0189 else
0190 dogs = true;
0191 endif
0192 if (dogs)
0193 ## Default to golden section step.
0194
0195 ## WARNING: This is also the "initial" procedure following MATLAB
0196 ## nomenclature. After the loop we'll fix the string for the first step.
0197 iter(niter+1).procedure = "golden";
0198
0199 e = ifelse (x >= xm, a - x, b - x);
0200 d = c * e;
0201 endif
0202
0203 ## f must not be evaluated too close to x.
0204 u = x + max (abs (d), tol) * (sign (d) + (d == 0));
0205 fu = fun (u, varargin{:});
0206
0207 niter += 1;
0208
0209 iter(niter).funccount = nfev++;
0210 iter(niter).x = u;
0211 iter(niter).fx = fu;
0212
0213 ## update a, b, v, w, and x
0214
0215 if (fu < fval)
0216 if (u < x)
0217 b = x;
0218 else
0219 a = x;
0220 endif
0221 v = w; fv = fw;
0222 w = x; fw = fval;
0223 x = u; fval = fu;
0224 else
0225 ## The following if-statement was originally executed even if fu == fval.
0226 if (u < x)
0227 a = u;
0228 else
0229 b = u;
0230 endif
0231 if (fu <= fw || w == x)
0232 v = w; fv = fw;
0233 w = u; fw = fu;
0234 elseif (fu <= fv || v == x || v == w)
0235 v = u;
0236 fv = fu;
0237 endif
0238 endif
0239
0240 ## If there's an output function, use it now.
0241 if (! isempty (outfcn))
0242 optv.funccount = nfev;
0243 optv.fval = fval;
0244 optv.iteration = niter;
0245 if (outfcn (x, optv, "iter"))
0246 info = -1;
0247 break;
0248 endif
0249 endif
0250 endwhile
0251
0252 ## Fix the first step procedure.
0253 iter(1).procedure = "initial";
0254
0255 ## Handle the "Display" option
0256 switch (displ)
0257 case "iter"
0258 print_formatted_table (iter);
0259 print_exit_msg (info, struct ("TolX", tolx, "fx", fval));
0260 case "notify"
0261 if (info == 0)
0262 print_exit_msg (info, struct ("fx",fval));
0263 endif
0264 case "final"
0265 print_exit_msg (info, struct ("TolX", tolx, "fx", fval));
0266 case "off"
0267 "skip";
0268 otherwise
0269 warning ("fminbnd: unknown option for Display: '%s'", displ);
0270 endswitch
0271
0272 output.iterations = niter;
0273 output.funcCount = nfev;
0274 output.algorithm = "golden section search, parabolic interpolation";
0275 output.bracket = [a, b];
0276 ## FIXME: bracketf possibly unavailable.
0277
0278 endfunction
0279
0280 ## A helper function that evaluates a function and checks for bad results.
0281 function fx = guarded_eval (fun, x)
0282
0283 fx = fun (x);
0284 fx = fx(1);
0285 if (! isreal (fx))
0286 error ("Octave:fmindbnd:notreal", "fminbnd: non-real value encountered");
0287 elseif (isnan (fx))
0288 error ("Octave:fmindbnd:isnan", "fminbnd: NaN value encountered");
0289 endif
0290
0291 endfunction
0292
0293 ## A hack for printing a formatted table
0294 function print_formatted_table (table)
0295 printf ("\n Func-count x f(x) Procedure\n");
0296 for row=table
0297 printf ("
0298 int2str (row.funccount), num2str (row.x,"
0299 num2str (row.fx,"
0300 endfor
0301 printf ("\n");
0302 endfunction
0303
0304 ## Print either a success termination message or bad news
0305 function print_exit_msg (info, opt=struct ())
0306
0307 printf ("");
0308 switch (info)
0309 case 1
0310 printf ("Optimization terminated:\n");
0311 printf (" the current x satisfies the termination criteria using OPTIONS.TolX of
0312 case 0
0313 printf ("Exiting: Maximum number of iterations has been exceeded\n");
0314 printf (" - increase MaxIter option.\n");
0315 printf (" Current function value:
0316 case -1
0317 "FIXME"; # FIXME: what's the message MATLAB prints for this case?
0318 otherwise
0319 error ("fminbnd: internal error, info return code was
0320 endswitch
0321 printf ("\n");
0322
0323 endfunction
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335 