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