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 ## This was submitted to octave as https://savannah.gnu.org/bugs/?65551
0094 ## However, it is marked as "Won't fix". We thus need to keep overloads
0095
0096 function [x, fval, info, output] = fminbnd (fun, a, b, options = struct (), varargin = {})
0097
0098 ## Get default options if requested.
0099 if (nargin == 1 && ischar (fun) && strcmp (fun, "defaults"))
0100 x = struct ("Display", "notify", "FunValCheck", "off",
0101 "MaxFunEvals", 500, "MaxIter", 500,
0102 "OutputFcn", [], "TolX", 1e-4);
0103 return;
0104 endif
0105
0106 if (nargin < 2)
0107 print_usage ();
0108 endif
0109
0110 if (a > b)
0111 error ("Octave:invalid-input-arg",
0112 "fminbnd: the lower bound cannot be greater than the upper one");
0113 endif
0114
0115 if (ischar (fun))
0116 fun = str2func (fun);
0117 endif
0118
0119 displ = optimget (options, "Display", "notify");
0120 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
0121 outfcn = optimget (options, "OutputFcn");
0122 tolx = optimget (options, "TolX", 1e-4);
0123 maxiter = optimget (options, "MaxIter", 500);
0124 maxfev = optimget (options, "MaxFunEvals", 500);
0125
0126 if (funvalchk)
0127 ## Replace fun with a guarded version.
0128 fun = @(x) guarded_eval (fun, x, varargin{:});
0129 endif
0130
0131 ## The default exit flag if exceeded number of iterations.
0132 info = 0;
0133 niter = 0;
0134 nfev = 0;
0135
0136 c = 0.5*(3 - sqrt (5));
0137 v = a + c*(b-a);
0138 w = x = v;
0139 e = 0;
0140 fv = fw = fval = fun (x, varargin{:});
0141 nfev += 1;
0142
0143 if (isa (a, "single") || isa (b, "single") || isa (fval, "single"))
0144 sqrteps = eps ("single");
0145 else
0146 sqrteps = eps ("double");
0147 endif
0148
0149 ## Only for display purposes.
0150 iter(1).funccount = nfev;
0151 iter(1).x = x;
0152 iter(1).fx = fval;
0153
0154 while (niter < maxiter && nfev < maxfev)
0155 xm = 0.5*(a+b);
0156 ## FIXME: the golden section search can actually get closer than sqrt(eps)
0157 ## sometimes. Sometimes not, it depends on the function. This is the
0158 ## strategy from the Netlib code. Something smarter would be good.
0159 tol = 2 * sqrteps * abs (x) + tolx / 3;
0160 if (abs (x - xm) <= (2*tol - 0.5*(b-a)))
0161 info = 1;
0162 break;
0163 endif
0164
0165 if (abs (e) > tol)
0166 dogs = false;
0167 ## Try inverse parabolic step.
0168 iter(niter+1).procedure = "parabolic";
0169
0170 r = (x - w)*(fval - fv);
0171 q = (x - v)*(fval - fw);
0172 p = (x - v)*q - (x - w)*r;
0173 q = 2*(q - r);
0174 p *= -sign (q);
0175 q = abs (q);
0176 r = e;
0177 e = d;
0178
0179 if (abs (p) < abs (0.5*q*r) && p > q*(a-x) && p < q*(b-x))
0180 ## The parabolic step is acceptable.
0181 d = p / q;
0182 u = x + d;
0183
0184 ## f must not be evaluated too close to ax or bx.
0185 if (min (u-a, b-u) < 2*tol)
0186 d = tol * (sign (xm - x) + (xm == x));
0187 endif
0188 else
0189 dogs = true;
0190 endif
0191 else
0192 dogs = true;
0193 endif
0194 if (dogs)
0195 ## Default to golden section step.
0196
0197 ## WARNING: This is also the "initial" procedure following MATLAB
0198 ## nomenclature. After the loop we'll fix the string for the first step.
0199 iter(niter+1).procedure = "golden";
0200
0201 e = ifelse (x >= xm, a - x, b - x);
0202 d = c * e;
0203 endif
0204
0205 ## f must not be evaluated too close to x.
0206 u = x + max (abs (d), tol) * (sign (d) + (d == 0));
0207 fu = fun (u, varargin{:});
0208
0209 niter += 1;
0210
0211 iter(niter).funccount = nfev++;
0212 iter(niter).x = u;
0213 iter(niter).fx = fu;
0214
0215 ## update a, b, v, w, and x
0216
0217 if (fu < fval)
0218 if (u < x)
0219 b = x;
0220 else
0221 a = x;
0222 endif
0223 v = w; fv = fw;
0224 w = x; fw = fval;
0225 x = u; fval = fu;
0226 else
0227 ## The following if-statement was originally executed even if fu == fval.
0228 if (u < x)
0229 a = u;
0230 else
0231 b = u;
0232 endif
0233 if (fu <= fw || w == x)
0234 v = w; fv = fw;
0235 w = u; fw = fu;
0236 elseif (fu <= fv || v == x || v == w)
0237 v = u;
0238 fv = fu;
0239 endif
0240 endif
0241
0242 ## If there's an output function, use it now.
0243 if (! isempty (outfcn))
0244 optv.funccount = nfev;
0245 optv.fval = fval;
0246 optv.iteration = niter;
0247 if (outfcn (x, optv, "iter"))
0248 info = -1;
0249 break;
0250 endif
0251 endif
0252 endwhile
0253
0254 ## Fix the first step procedure.
0255 iter(1).procedure = "initial";
0256
0257 ## Handle the "Display" option
0258 switch (displ)
0259 case "iter"
0260 print_formatted_table (iter);
0261 print_exit_msg (info, struct ("TolX", tolx, "fx", fval));
0262 case "notify"
0263 if (info == 0)
0264 print_exit_msg (info, struct ("fx",fval));
0265 endif
0266 case "final"
0267 print_exit_msg (info, struct ("TolX", tolx, "fx", fval));
0268 case "off"
0269 "skip";
0270 otherwise
0271 warning ("fminbnd: unknown option for Display: '%s'", displ);
0272 endswitch
0273
0274 output.iterations = niter;
0275 output.funcCount = nfev;
0276 output.algorithm = "golden section search, parabolic interpolation";
0277 output.bracket = [a, b];
0278 ## FIXME: bracketf possibly unavailable.
0279
0280 endfunction
0281
0282 ## A helper function that evaluates a function and checks for bad results.
0283 function fx = guarded_eval (fun, x)
0284
0285 fx = fun (x);
0286 fx = fx(1);
0287 if (! isreal (fx))
0288 error ("Octave:fmindbnd:notreal", "fminbnd: non-real value encountered");
0289 elseif (isnan (fx))
0290 error ("Octave:fmindbnd:isnan", "fminbnd: NaN value encountered");
0291 endif
0292
0293 endfunction
0294
0295 ## A hack for printing a formatted table
0296 function print_formatted_table (table)
0297 printf ("\n Func-count x f(x) Procedure\n");
0298 for row=table
0299 printf ("
0300 int2str (row.funccount), num2str (row.x,"
0301 num2str (row.fx,"
0302 endfor
0303 printf ("\n");
0304 endfunction
0305
0306 ## Print either a success termination message or bad news
0307 function print_exit_msg (info, opt=struct ())
0308
0309 printf ("");
0310 switch (info)
0311 case 1
0312 printf ("Optimization terminated:\n");
0313 printf (" the current x satisfies the termination criteria using OPTIONS.TolX of
0314 case 0
0315 printf ("Exiting: Maximum number of iterations has been exceeded\n");
0316 printf (" - increase MaxIter option.\n");
0317 printf (" Current function value:
0318 case -1
0319 "FIXME"; # FIXME: what's the message MATLAB prints for this case?
0320 otherwise
0321 error ("fminbnd: internal error, info return code was
0322 endswitch
0323 printf ("\n");
0324
0325 endfunction
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337 