Mercurial > octave-nkf

comparison scripts/optimization/fminbnd.m @ 15706:242e9efd4315

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added Display option for fminbnd()
author Júlio Hoffimann <julio.hoffimann@gmail.com>
date Sun, 02 Sep 2012 16:21:48 -0300
parents e0525ecf156e
children 7eff3032d144
comparison
equal deleted inserted replaced
15705:a3189d329906 15706:242e9efd4315
19 ## Author: Jaroslav Hajek <highegg@gmail.com> 19 ## Author: Jaroslav Hajek <highegg@gmail.com>
20 20
21 ## -*- texinfo -*- 21 ## -*- texinfo -*-
22 ## @deftypefn {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options}) 22 ## @deftypefn {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options})
23 ## Find a minimum point of a univariate function. 23 ## Find a minimum point of a univariate function.
24 ## 24 ##
25 ## @var{fun} should be a function handle or name. @var{a}, @var{b} specify a 25 ## @var{fun} should be a function handle or name. @var{a}, @var{b} specify a
26 ## starting interval. @var{options} is a structure specifying additional 26 ## starting interval. @var{options} is a structure specifying additional
27 ## options. Currently, @code{fminbnd} recognizes these options: 27 ## options. Currently, @code{fminbnd} recognizes these options:
28 ## "FunValCheck", "OutputFcn", "TolX", "MaxIter", "MaxFunEvals". For a 28 ## "FunValCheck", "OutputFcn", "TolX", "MaxIter", "MaxFunEvals". For a
29 ## description of these options, see @ref{doc-optimset,,optimset}. 29 ## description of these options, see @ref{doc-optimset,,optimset}.
72 72
73 if (ischar (fun)) 73 if (ischar (fun))
74 fun = str2func (fun, "global"); 74 fun = str2func (fun, "global");
75 endif 75 endif
76 76
77 ## TODO 77 displ = optimget (options, "Display", "notify");
78 ## displev = optimget (options, "Display", "notify");
79 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on"); 78 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
80 outfcn = optimget (options, "OutputFcn"); 79 outfcn = optimget (options, "OutputFcn");
81 tolx = optimget (options, "TolX", 1e-8); 80 tolx = optimget (options, "TolX", 1e-8);
82 maxiter = optimget (options, "MaxIter", Inf); 81 maxiter = optimget (options, "MaxIter", Inf);
83 maxfev = optimget (options, "MaxFunEvals", Inf); 82 maxfev = optimget (options, "MaxFunEvals", Inf);
98 v = a + c*(b-a); 97 v = a + c*(b-a);
99 w = x = v; 98 w = x = v;
100 e = 0; 99 e = 0;
101 fv = fw = fval = fun (x); 100 fv = fw = fval = fun (x);
102 nfev++; 101 nfev++;
102
103 ## Only for display purposes.
104 iter(1).funccount = nfev;
105 iter(1).x = x;
106 iter(1).fx = fval;
103 107
104 while (niter < maxiter && nfev < maxfev) 108 while (niter < maxiter && nfev < maxfev)
105 xm = 0.5*(a+b); 109 xm = 0.5*(a+b);
106 ## FIXME: the golden section search can actually get closer than sqrt(eps) 110 ## FIXME: the golden section search can actually get closer than sqrt(eps)
107 ## sometimes. Sometimes not, it depends on the function. This is the 111 ## sometimes. Sometimes not, it depends on the function. This is the
113 endif 117 endif
114 118
115 if (abs (e) > tol) 119 if (abs (e) > tol)
116 dogs = false; 120 dogs = false;
117 ## Try inverse parabolic step. 121 ## Try inverse parabolic step.
122 iter(niter+1).procedure = "parabolic";
123
118 r = (x - w)*(fval - fv); 124 r = (x - w)*(fval - fv);
119 q = (x - v)*(fval - fw); 125 q = (x - v)*(fval - fw);
120 p = (x - v)*q - (x - w)*r; 126 p = (x - v)*q - (x - w)*r;
121 q = 2*(q - r); 127 q = 2*(q - r);
122 p *= -sign (q); 128 p *= -sign (q);
139 else 145 else
140 dogs = true; 146 dogs = true;
141 endif 147 endif
142 if (dogs) 148 if (dogs)
143 ## Default to golden section step. 149 ## Default to golden section step.
150
151 ## WARNING: This is also the "initial" procedure following
152 ## MATLAB nomenclature. After the loop we'll fix the string
153 ## for the first step.
154 iter(niter+1).procedure = "golden";
155
144 e = ifelse (x >= xm, a - x, b - x); 156 e = ifelse (x >= xm, a - x, b - x);
145 d = c * e; 157 d = c * e;
146 endif 158 endif
147 159
148 ## f must not be evaluated too close to x. 160 ## f must not be evaluated too close to x.
149 u = x + max (abs (d), tol) * (sign (d) + (d == 0)); 161 u = x + max (abs (d), tol) * (sign (d) + (d == 0));
150 162 fu = fun (u);
151 fu = fun (u); 163
152 nfev++; 164 niter++;
153 niter++; 165
154 166 iter(niter).funccount = nfev++;
155 ## update a, b, v, w, and x 167 iter(niter).x = u;
156 168 iter(niter).fx = fu;
157 if (fu <= fval) 169
158 if (u < x) 170 ## update a, b, v, w, and x
159 b = x; 171
160 else 172 if (fu <= fval)
161 a = x; 173 if (u < x)
162 endif 174 b = x;
163 v = w; fv = fw; 175 else
164 w = x; fw = fval; 176 a = x;
165 x = u; fval = fu; 177 endif
166 else 178 v = w; fv = fw;
167 ## The following if-statement was originally executed even if fu == fval. 179 w = x; fw = fval;
168 if (u < x) 180 x = u; fval = fu;
169 a = u; 181 else
170 else 182 ## The following if-statement was originally executed even if fu == fval.
171 b = u; 183 if (u < x)
172 endif 184 a = u;
173 if (fu <= fw || w == x) 185 else
174 v = w; fv = fw; 186 b = u;
175 w = u; fw = fu; 187 endif
176 elseif (fu <= fv || v == x || v == w) 188 if (fu <= fw || w == x)
177 v = u; 189 v = w; fv = fw;
178 fv = fu; 190 w = u; fw = fu;
179 endif 191 elseif (fu <= fv || v == x || v == w)
180 endif 192 v = u;
193 fv = fu;
194 endif
195 endif
181 196
182 ## If there's an output function, use it now. 197 ## If there's an output function, use it now.
183 if (outfcn) 198 if (outfcn)
184 optv.funccount = nfev; 199 optv.funccount = nfev;
185 optv.fval = fval; 200 optv.fval = fval;
188 info = -1; 203 info = -1;
189 break; 204 break;
190 endif 205 endif
191 endif 206 endif
192 endwhile 207 endwhile
208
209 ## Fix the first step procedure.
210 iter(1).procedure = "initial";
211
212 ## Handle the "Display" option
213 switch displ
214 case "iter"
215 print_formatted_table (iter);
216 print_exit_msg (info, struct("TolX", tolx, "fx", fval));
217 case "notify"
218 if (info == 0)
219 print_exit_msg (info, struct("fx",fval));
220 endif
221 case "final"
222 print_exit_msg (info, struct("TolX", tolx, "fx", fval));
223 case "off"
224 "skip";
225 otherwise
226 warning ("unknown option for Display: '%s'", displ);
227 endswitch
193 228
194 output.iterations = niter; 229 output.iterations = niter;
195 output.funcCount = nfev; 230 output.funcCount = nfev;
196 output.bracket = [a, b]; 231 output.bracket = [a, b];
197 ## FIXME: bracketf possibly unavailable. 232 ## FIXME: bracketf possibly unavailable.
208 elseif (isnan (fx)) 243 elseif (isnan (fx))
209 error ("fminbnd:isnan", "fminbnd: NaN value encountered"); 244 error ("fminbnd:isnan", "fminbnd: NaN value encountered");
210 endif 245 endif
211 endfunction 246 endfunction
212 247
248 ## A hack for printing a formatted table
249 function print_formatted_table (table)
250 printf ("\n Func-count x f(x) Procedure\n");
251 for row=table
252 printf("%5.5s %7.7s %8.8s\t%s\n",
253 int2str (row.funccount), num2str (row.x,"%.5f"),
254 num2str (row.fx,"%.6f"), row.procedure);
255 endfor
256 printf ("\n");
257 endfunction
258
259 ## Print either a success termination message or bad news
260 function print_exit_msg (info, opt=struct())
261 printf ("");
262 switch info
263 case 1
264 printf ("Optimization terminated:\n");
265 printf (" the current x satisfies the termination criteria using OPTIONS.TolX of %e\n", opt.TolX);
266 case 0
267 printf ("Exiting: Maximum number of iterations has been exceeded\n");
268 printf (" - increase MaxIter option.\n");
269 printf (" Current function value: %.6f\n", opt.fx);
270 case -1
271 "FIXME"; ## FIXME: what's the message MATLAB prints for this case?
272 otherwise
273 error ("internal error - fminbnd() is bug, sorry!");
274 endswitch
275 printf ("\n");
276 endfunction
277
213 278
214 %!shared opt0 279 %!shared opt0
215 %! opt0 = optimset ("tolx", 0); 280 %! opt0 = optimset ("tolx", 0);
216 %!assert (fminbnd (@cos, pi/2, 3*pi/2, opt0), pi, 10*sqrt (eps)) 281 %!assert (fminbnd (@cos, pi/2, 3*pi/2, opt0), pi, 10*sqrt (eps))
217 %!assert (fminbnd (@(x) (x - 1e-3)^4, -1, 1, opt0), 1e-3, 10e-3*sqrt (eps)) 282 %!assert (fminbnd (@(x) (x - 1e-3)^4, -1, 1, opt0), 1e-3, 10e-3*sqrt (eps))