1 | /*************************************************************************
|
---|
2 | Copyright (c) 2007-2008, Sergey Bochkanov (ALGLIB project).
|
---|
3 |
|
---|
4 | >>> SOURCE LICENSE >>>
|
---|
5 | This program is free software; you can redistribute it and/or modify
|
---|
6 | it under the terms of the GNU General Public License as published by
|
---|
7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
8 | License, or (at your option) any later version.
|
---|
9 |
|
---|
10 | This program is distributed in the hope that it will be useful,
|
---|
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
13 | GNU General Public License for more details.
|
---|
14 |
|
---|
15 | A copy of the GNU General Public License is available at
|
---|
16 | http://www.fsf.org/licensing/licenses
|
---|
17 |
|
---|
18 | >>> END OF LICENSE >>>
|
---|
19 | *************************************************************************/
|
---|
20 |
|
---|
21 | using System;
|
---|
22 |
|
---|
23 | namespace alglib
|
---|
24 | {
|
---|
25 | public class lbfgs
|
---|
26 | {
|
---|
27 | public struct lbfgsstate
|
---|
28 | {
|
---|
29 | public int n;
|
---|
30 | public int m;
|
---|
31 | public double epsg;
|
---|
32 | public double epsf;
|
---|
33 | public double epsx;
|
---|
34 | public int maxits;
|
---|
35 | public int flags;
|
---|
36 | public int nfev;
|
---|
37 | public int mcstage;
|
---|
38 | public int k;
|
---|
39 | public int q;
|
---|
40 | public int p;
|
---|
41 | public double[] rho;
|
---|
42 | public double[,] y;
|
---|
43 | public double[,] s;
|
---|
44 | public double[] theta;
|
---|
45 | public double[] d;
|
---|
46 | public double stp;
|
---|
47 | public double[] work;
|
---|
48 | public double fold;
|
---|
49 | public double gammak;
|
---|
50 | public double[] x;
|
---|
51 | public double f;
|
---|
52 | public double[] g;
|
---|
53 | public bool xupdated;
|
---|
54 | public AP.rcommstate rstate;
|
---|
55 | public int repiterationscount;
|
---|
56 | public int repnfev;
|
---|
57 | public int repterminationtype;
|
---|
58 | public bool brackt;
|
---|
59 | public bool stage1;
|
---|
60 | public int infoc;
|
---|
61 | public double dg;
|
---|
62 | public double dgm;
|
---|
63 | public double dginit;
|
---|
64 | public double dgtest;
|
---|
65 | public double dgx;
|
---|
66 | public double dgxm;
|
---|
67 | public double dgy;
|
---|
68 | public double dgym;
|
---|
69 | public double finit;
|
---|
70 | public double ftest1;
|
---|
71 | public double fm;
|
---|
72 | public double fx;
|
---|
73 | public double fxm;
|
---|
74 | public double fy;
|
---|
75 | public double fym;
|
---|
76 | public double stx;
|
---|
77 | public double sty;
|
---|
78 | public double stmin;
|
---|
79 | public double stmax;
|
---|
80 | public double width;
|
---|
81 | public double width1;
|
---|
82 | public double xtrapf;
|
---|
83 | };
|
---|
84 |
|
---|
85 |
|
---|
86 | public struct lbfgsreport
|
---|
87 | {
|
---|
88 | public int iterationscount;
|
---|
89 | public int nfev;
|
---|
90 | public int terminationtype;
|
---|
91 | };
|
---|
92 |
|
---|
93 |
|
---|
94 |
|
---|
95 |
|
---|
96 | public const double ftol = 0.0001;
|
---|
97 | public const double xtol = 100*AP.Math.MachineEpsilon;
|
---|
98 | public const double gtol = 0.9;
|
---|
99 | public const int maxfev = 20;
|
---|
100 | public const double stpmin = 1.0E-20;
|
---|
101 | public const double stpmax = 1.0E20;
|
---|
102 |
|
---|
103 |
|
---|
104 | /*************************************************************************
|
---|
105 | LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
|
---|
106 |
|
---|
107 | The subroutine minimizes function F(x) of N arguments by using a quasi-
|
---|
108 | Newton method (LBFGS scheme) which is optimized to use a minimum amount
|
---|
109 | of memory.
|
---|
110 |
|
---|
111 | The subroutine generates the approximation of an inverse Hessian matrix by
|
---|
112 | using information about the last M steps of the algorithm (instead of N).
|
---|
113 | It lessens a required amount of memory from a value of order N^2 to a
|
---|
114 | value of order 2*N*M.
|
---|
115 |
|
---|
116 | Input parameters:
|
---|
117 | N - problem dimension. N>0
|
---|
118 | M - number of corrections in the BFGS scheme of Hessian
|
---|
119 | approximation update. Recommended value: 3<=M<=7. The smaller
|
---|
120 | value causes worse convergence, the bigger will not cause a
|
---|
121 | considerably better convergence, but will cause a fall in the
|
---|
122 | performance. M<=N.
|
---|
123 | X - initial solution approximation, array[0..N-1].
|
---|
124 | EpsG - positive number which defines a precision of search. The
|
---|
125 | subroutine finishes its work if the condition ||G|| < EpsG is
|
---|
126 | satisfied, where ||.|| means Euclidian norm, G - gradient, X -
|
---|
127 | current approximation.
|
---|
128 | EpsF - positive number which defines a precision of search. The
|
---|
129 | subroutine finishes its work if on iteration number k+1 the
|
---|
130 | condition |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1} is
|
---|
131 | satisfied.
|
---|
132 | EpsX - positive number which defines a precision of search. The
|
---|
133 | subroutine finishes its work if on iteration number k+1 the
|
---|
134 | condition |X(k+1)-X(k)| <= EpsX is fulfilled.
|
---|
135 | MaxIts- maximum number of iterations. If MaxIts=0, the number of
|
---|
136 | iterations is unlimited.
|
---|
137 | Flags - additional settings:
|
---|
138 | * Flags = 0 means no additional settings
|
---|
139 | * Flags = 1 "do not allocate memory". used when solving
|
---|
140 | a many subsequent tasks with same N/M values.
|
---|
141 | First call MUST be without this flag bit set,
|
---|
142 | subsequent calls of MinLBFGS with same LBFGSState
|
---|
143 | structure can set Flags to 1.
|
---|
144 |
|
---|
145 | Output parameters:
|
---|
146 | State - structure used for reverse communication.
|
---|
147 |
|
---|
148 | See also MinLBFGSIteration, MinLBFGSResults
|
---|
149 |
|
---|
150 | -- ALGLIB --
|
---|
151 | Copyright 14.11.2007 by Bochkanov Sergey
|
---|
152 | *************************************************************************/
|
---|
153 | public static void minlbfgs(int n,
|
---|
154 | int m,
|
---|
155 | ref double[] x,
|
---|
156 | double epsg,
|
---|
157 | double epsf,
|
---|
158 | double epsx,
|
---|
159 | int maxits,
|
---|
160 | int flags,
|
---|
161 | ref lbfgsstate state)
|
---|
162 | {
|
---|
163 | bool allocatemem = new bool();
|
---|
164 | int i_ = 0;
|
---|
165 |
|
---|
166 | System.Diagnostics.Debug.Assert(n>=1, "MinLBFGS: N too small!");
|
---|
167 | System.Diagnostics.Debug.Assert(m>=1, "MinLBFGS: M too small!");
|
---|
168 | System.Diagnostics.Debug.Assert(m<=n, "MinLBFGS: M too large!");
|
---|
169 | System.Diagnostics.Debug.Assert((double)(epsg)>=(double)(0), "MinLBFGS: negative EpsG!");
|
---|
170 | System.Diagnostics.Debug.Assert((double)(epsf)>=(double)(0), "MinLBFGS: negative EpsF!");
|
---|
171 | System.Diagnostics.Debug.Assert((double)(epsx)>=(double)(0), "MinLBFGS: negative EpsX!");
|
---|
172 | System.Diagnostics.Debug.Assert(maxits>=0, "MinLBFGS: negative MaxIts!");
|
---|
173 |
|
---|
174 | //
|
---|
175 | // Initialize
|
---|
176 | //
|
---|
177 | state.n = n;
|
---|
178 | state.m = m;
|
---|
179 | state.epsg = epsg;
|
---|
180 | state.epsf = epsf;
|
---|
181 | state.epsx = epsx;
|
---|
182 | state.maxits = maxits;
|
---|
183 | state.flags = flags;
|
---|
184 | allocatemem = flags%2==0;
|
---|
185 | flags = flags/2;
|
---|
186 | if( allocatemem )
|
---|
187 | {
|
---|
188 | state.rho = new double[m-1+1];
|
---|
189 | state.theta = new double[m-1+1];
|
---|
190 | state.y = new double[m-1+1, n-1+1];
|
---|
191 | state.s = new double[m-1+1, n-1+1];
|
---|
192 | state.d = new double[n-1+1];
|
---|
193 | state.x = new double[n-1+1];
|
---|
194 | state.g = new double[n-1+1];
|
---|
195 | state.work = new double[n-1+1];
|
---|
196 | }
|
---|
197 |
|
---|
198 | //
|
---|
199 | // Initialize Rep structure
|
---|
200 | //
|
---|
201 | state.xupdated = false;
|
---|
202 |
|
---|
203 | //
|
---|
204 | // Prepare first run
|
---|
205 | //
|
---|
206 | state.k = 0;
|
---|
207 | for(i_=0; i_<=n-1;i_++)
|
---|
208 | {
|
---|
209 | state.x[i_] = x[i_];
|
---|
210 | }
|
---|
211 | state.rstate.ia = new int[6+1];
|
---|
212 | state.rstate.ra = new double[4+1];
|
---|
213 | state.rstate.stage = -1;
|
---|
214 | }
|
---|
215 |
|
---|
216 |
|
---|
217 | /*************************************************************************
|
---|
218 | One L-BFGS iteration
|
---|
219 |
|
---|
220 | Called after initialization with MinLBFGS.
|
---|
221 | See HTML documentation for examples.
|
---|
222 |
|
---|
223 | Input parameters:
|
---|
224 | State - structure which stores algorithm state between calls and
|
---|
225 | which is used for reverse communication. Must be initialized
|
---|
226 | with MinLBFGS.
|
---|
227 |
|
---|
228 | If suborutine returned False, iterative proces has converged.
|
---|
229 |
|
---|
230 | If subroutine returned True, caller should calculate function value
|
---|
231 | State.F an gradient State.G[0..N-1] at State.X[0..N-1] and call
|
---|
232 | MinLBFGSIteration again.
|
---|
233 |
|
---|
234 | -- ALGLIB --
|
---|
235 | Copyright 20.04.2009 by Bochkanov Sergey
|
---|
236 | *************************************************************************/
|
---|
237 | public static bool minlbfgsiteration(ref lbfgsstate state)
|
---|
238 | {
|
---|
239 | bool result = new bool();
|
---|
240 | int n = 0;
|
---|
241 | int m = 0;
|
---|
242 | int maxits = 0;
|
---|
243 | double epsf = 0;
|
---|
244 | double epsg = 0;
|
---|
245 | double epsx = 0;
|
---|
246 | int i = 0;
|
---|
247 | int j = 0;
|
---|
248 | int ic = 0;
|
---|
249 | int mcinfo = 0;
|
---|
250 | double v = 0;
|
---|
251 | double vv = 0;
|
---|
252 | int i_ = 0;
|
---|
253 |
|
---|
254 |
|
---|
255 | //
|
---|
256 | // Reverse communication preparations
|
---|
257 | // I know it looks ugly, but it works the same way
|
---|
258 | // anywhere from C++ to Python.
|
---|
259 | //
|
---|
260 | // This code initializes locals by:
|
---|
261 | // * random values determined during code
|
---|
262 | // generation - on first subroutine call
|
---|
263 | // * values from previous call - on subsequent calls
|
---|
264 | //
|
---|
265 | if( state.rstate.stage>=0 )
|
---|
266 | {
|
---|
267 | n = state.rstate.ia[0];
|
---|
268 | m = state.rstate.ia[1];
|
---|
269 | maxits = state.rstate.ia[2];
|
---|
270 | i = state.rstate.ia[3];
|
---|
271 | j = state.rstate.ia[4];
|
---|
272 | ic = state.rstate.ia[5];
|
---|
273 | mcinfo = state.rstate.ia[6];
|
---|
274 | epsf = state.rstate.ra[0];
|
---|
275 | epsg = state.rstate.ra[1];
|
---|
276 | epsx = state.rstate.ra[2];
|
---|
277 | v = state.rstate.ra[3];
|
---|
278 | vv = state.rstate.ra[4];
|
---|
279 | }
|
---|
280 | else
|
---|
281 | {
|
---|
282 | n = -983;
|
---|
283 | m = -989;
|
---|
284 | maxits = -834;
|
---|
285 | i = 900;
|
---|
286 | j = -287;
|
---|
287 | ic = 364;
|
---|
288 | mcinfo = 214;
|
---|
289 | epsf = -338;
|
---|
290 | epsg = -686;
|
---|
291 | epsx = 912;
|
---|
292 | v = 585;
|
---|
293 | vv = 497;
|
---|
294 | }
|
---|
295 | if( state.rstate.stage==0 )
|
---|
296 | {
|
---|
297 | goto lbl_0;
|
---|
298 | }
|
---|
299 | if( state.rstate.stage==1 )
|
---|
300 | {
|
---|
301 | goto lbl_1;
|
---|
302 | }
|
---|
303 |
|
---|
304 | //
|
---|
305 | // Routine body
|
---|
306 | //
|
---|
307 |
|
---|
308 | //
|
---|
309 | // Unload frequently used variables from State structure
|
---|
310 | // (just for typing convinience)
|
---|
311 | //
|
---|
312 | n = state.n;
|
---|
313 | m = state.m;
|
---|
314 | epsg = state.epsg;
|
---|
315 | epsf = state.epsf;
|
---|
316 | epsx = state.epsx;
|
---|
317 | maxits = state.maxits;
|
---|
318 | state.repterminationtype = 0;
|
---|
319 | state.repiterationscount = 0;
|
---|
320 | state.repnfev = 0;
|
---|
321 |
|
---|
322 | //
|
---|
323 | // Update info
|
---|
324 | //
|
---|
325 | state.xupdated = false;
|
---|
326 |
|
---|
327 | //
|
---|
328 | // Calculate F/G
|
---|
329 | //
|
---|
330 | state.rstate.stage = 0;
|
---|
331 | goto lbl_rcomm;
|
---|
332 | lbl_0:
|
---|
333 | state.repnfev = 1;
|
---|
334 |
|
---|
335 | //
|
---|
336 | // Preparations
|
---|
337 | //
|
---|
338 | state.fold = state.f;
|
---|
339 | v = 0.0;
|
---|
340 | for(i_=0; i_<=n-1;i_++)
|
---|
341 | {
|
---|
342 | v += state.g[i_]*state.g[i_];
|
---|
343 | }
|
---|
344 | v = Math.Sqrt(v);
|
---|
345 | if( (double)(v)==(double)(0) )
|
---|
346 | {
|
---|
347 | state.repterminationtype = 4;
|
---|
348 | result = false;
|
---|
349 | return result;
|
---|
350 | }
|
---|
351 | state.stp = 1.0/v;
|
---|
352 | for(i_=0; i_<=n-1;i_++)
|
---|
353 | {
|
---|
354 | state.d[i_] = -state.g[i_];
|
---|
355 | }
|
---|
356 |
|
---|
357 | //
|
---|
358 | // Main cycle
|
---|
359 | //
|
---|
360 | lbl_2:
|
---|
361 | if( false )
|
---|
362 | {
|
---|
363 | goto lbl_3;
|
---|
364 | }
|
---|
365 |
|
---|
366 | //
|
---|
367 | // Main cycle: prepare to 1-D line search
|
---|
368 | //
|
---|
369 | state.p = state.k%m;
|
---|
370 | state.q = Math.Min(state.k, m-1);
|
---|
371 |
|
---|
372 | //
|
---|
373 | // Store X[k], G[k]
|
---|
374 | //
|
---|
375 | for(i_=0; i_<=n-1;i_++)
|
---|
376 | {
|
---|
377 | state.s[state.p,i_] = -state.x[i_];
|
---|
378 | }
|
---|
379 | for(i_=0; i_<=n-1;i_++)
|
---|
380 | {
|
---|
381 | state.y[state.p,i_] = -state.g[i_];
|
---|
382 | }
|
---|
383 |
|
---|
384 | //
|
---|
385 | // Minimize F(x+alpha*d)
|
---|
386 | //
|
---|
387 | state.mcstage = 0;
|
---|
388 | if( state.k!=0 )
|
---|
389 | {
|
---|
390 | state.stp = 1.0;
|
---|
391 | }
|
---|
392 | mcsrch(n, ref state.x, ref state.f, ref state.g, ref state.d, ref state.stp, ref mcinfo, ref state.nfev, ref state.work, ref state, ref state.mcstage);
|
---|
393 | lbl_4:
|
---|
394 | if( state.mcstage==0 )
|
---|
395 | {
|
---|
396 | goto lbl_5;
|
---|
397 | }
|
---|
398 | state.rstate.stage = 1;
|
---|
399 | goto lbl_rcomm;
|
---|
400 | lbl_1:
|
---|
401 | mcsrch(n, ref state.x, ref state.f, ref state.g, ref state.d, ref state.stp, ref mcinfo, ref state.nfev, ref state.work, ref state, ref state.mcstage);
|
---|
402 | goto lbl_4;
|
---|
403 | lbl_5:
|
---|
404 |
|
---|
405 | //
|
---|
406 | // Main cycle: update information and Hessian.
|
---|
407 | // Check stopping conditions.
|
---|
408 | //
|
---|
409 | state.repnfev = state.repnfev+state.nfev;
|
---|
410 | state.repiterationscount = state.repiterationscount+1;
|
---|
411 |
|
---|
412 | //
|
---|
413 | // Calculate S[k], Y[k], Rho[k], GammaK
|
---|
414 | //
|
---|
415 | for(i_=0; i_<=n-1;i_++)
|
---|
416 | {
|
---|
417 | state.s[state.p,i_] = state.s[state.p,i_] + state.x[i_];
|
---|
418 | }
|
---|
419 | for(i_=0; i_<=n-1;i_++)
|
---|
420 | {
|
---|
421 | state.y[state.p,i_] = state.y[state.p,i_] + state.g[i_];
|
---|
422 | }
|
---|
423 |
|
---|
424 | //
|
---|
425 | // Stopping conditions
|
---|
426 | //
|
---|
427 | if( state.repiterationscount>=maxits & maxits>0 )
|
---|
428 | {
|
---|
429 |
|
---|
430 | //
|
---|
431 | // Too many iterations
|
---|
432 | //
|
---|
433 | state.repterminationtype = 5;
|
---|
434 | result = false;
|
---|
435 | return result;
|
---|
436 | }
|
---|
437 | v = 0.0;
|
---|
438 | for(i_=0; i_<=n-1;i_++)
|
---|
439 | {
|
---|
440 | v += state.g[i_]*state.g[i_];
|
---|
441 | }
|
---|
442 | if( (double)(Math.Sqrt(v))<=(double)(epsg) )
|
---|
443 | {
|
---|
444 |
|
---|
445 | //
|
---|
446 | // Gradient is small enough
|
---|
447 | //
|
---|
448 | state.repterminationtype = 4;
|
---|
449 | result = false;
|
---|
450 | return result;
|
---|
451 | }
|
---|
452 | if( (double)(state.fold-state.f)<=(double)(epsf*Math.Max(Math.Abs(state.fold), Math.Max(Math.Abs(state.f), 1.0))) )
|
---|
453 | {
|
---|
454 |
|
---|
455 | //
|
---|
456 | // F(k+1)-F(k) is small enough
|
---|
457 | //
|
---|
458 | state.repterminationtype = 1;
|
---|
459 | result = false;
|
---|
460 | return result;
|
---|
461 | }
|
---|
462 | v = 0.0;
|
---|
463 | for(i_=0; i_<=n-1;i_++)
|
---|
464 | {
|
---|
465 | v += state.s[state.p,i_]*state.s[state.p,i_];
|
---|
466 | }
|
---|
467 | if( (double)(Math.Sqrt(v))<=(double)(epsx) )
|
---|
468 | {
|
---|
469 |
|
---|
470 | //
|
---|
471 | // X(k+1)-X(k) is small enough
|
---|
472 | //
|
---|
473 | state.repterminationtype = 2;
|
---|
474 | result = false;
|
---|
475 | return result;
|
---|
476 | }
|
---|
477 |
|
---|
478 | //
|
---|
479 | // Calculate Rho[k], GammaK
|
---|
480 | //
|
---|
481 | v = 0.0;
|
---|
482 | for(i_=0; i_<=n-1;i_++)
|
---|
483 | {
|
---|
484 | v += state.y[state.p,i_]*state.s[state.p,i_];
|
---|
485 | }
|
---|
486 | vv = 0.0;
|
---|
487 | for(i_=0; i_<=n-1;i_++)
|
---|
488 | {
|
---|
489 | vv += state.y[state.p,i_]*state.y[state.p,i_];
|
---|
490 | }
|
---|
491 | if( (double)(v)==(double)(0) | (double)(vv)==(double)(0) )
|
---|
492 | {
|
---|
493 |
|
---|
494 | //
|
---|
495 | // Rounding errors make further iterations impossible.
|
---|
496 | //
|
---|
497 | state.repterminationtype = -2;
|
---|
498 | result = false;
|
---|
499 | return result;
|
---|
500 | }
|
---|
501 | state.rho[state.p] = 1/v;
|
---|
502 | state.gammak = v/vv;
|
---|
503 |
|
---|
504 | //
|
---|
505 | // Calculate d(k+1) = H(k+1)*g(k+1)
|
---|
506 | //
|
---|
507 | // for I:=K downto K-Q do
|
---|
508 | // V = s(i)^T * work(iteration:I)
|
---|
509 | // theta(i) = V
|
---|
510 | // work(iteration:I+1) = work(iteration:I) - V*Rho(i)*y(i)
|
---|
511 | // work(last iteration) = H0*work(last iteration)
|
---|
512 | // for I:=K-Q to K do
|
---|
513 | // V = y(i)^T*work(iteration:I)
|
---|
514 | // work(iteration:I+1) = work(iteration:I) +(-V+theta(i))*Rho(i)*s(i)
|
---|
515 | //
|
---|
516 | // NOW WORK CONTAINS d(k+1)
|
---|
517 | //
|
---|
518 | for(i_=0; i_<=n-1;i_++)
|
---|
519 | {
|
---|
520 | state.work[i_] = state.g[i_];
|
---|
521 | }
|
---|
522 | for(i=state.k; i>=state.k-state.q; i--)
|
---|
523 | {
|
---|
524 | ic = i%m;
|
---|
525 | v = 0.0;
|
---|
526 | for(i_=0; i_<=n-1;i_++)
|
---|
527 | {
|
---|
528 | v += state.s[ic,i_]*state.work[i_];
|
---|
529 | }
|
---|
530 | state.theta[ic] = v;
|
---|
531 | vv = v*state.rho[ic];
|
---|
532 | for(i_=0; i_<=n-1;i_++)
|
---|
533 | {
|
---|
534 | state.work[i_] = state.work[i_] - vv*state.y[ic,i_];
|
---|
535 | }
|
---|
536 | }
|
---|
537 | v = state.gammak;
|
---|
538 | for(i_=0; i_<=n-1;i_++)
|
---|
539 | {
|
---|
540 | state.work[i_] = v*state.work[i_];
|
---|
541 | }
|
---|
542 | for(i=state.k-state.q; i<=state.k; i++)
|
---|
543 | {
|
---|
544 | ic = i%m;
|
---|
545 | v = 0.0;
|
---|
546 | for(i_=0; i_<=n-1;i_++)
|
---|
547 | {
|
---|
548 | v += state.y[ic,i_]*state.work[i_];
|
---|
549 | }
|
---|
550 | vv = state.rho[ic]*(-v+state.theta[ic]);
|
---|
551 | for(i_=0; i_<=n-1;i_++)
|
---|
552 | {
|
---|
553 | state.work[i_] = state.work[i_] + vv*state.s[ic,i_];
|
---|
554 | }
|
---|
555 | }
|
---|
556 | for(i_=0; i_<=n-1;i_++)
|
---|
557 | {
|
---|
558 | state.d[i_] = -state.work[i_];
|
---|
559 | }
|
---|
560 |
|
---|
561 | //
|
---|
562 | // Next step
|
---|
563 | //
|
---|
564 | state.fold = state.f;
|
---|
565 | state.k = state.k+1;
|
---|
566 | state.xupdated = true;
|
---|
567 | goto lbl_2;
|
---|
568 | lbl_3:
|
---|
569 | result = false;
|
---|
570 | return result;
|
---|
571 |
|
---|
572 | //
|
---|
573 | // Saving state
|
---|
574 | //
|
---|
575 | lbl_rcomm:
|
---|
576 | result = true;
|
---|
577 | state.rstate.ia[0] = n;
|
---|
578 | state.rstate.ia[1] = m;
|
---|
579 | state.rstate.ia[2] = maxits;
|
---|
580 | state.rstate.ia[3] = i;
|
---|
581 | state.rstate.ia[4] = j;
|
---|
582 | state.rstate.ia[5] = ic;
|
---|
583 | state.rstate.ia[6] = mcinfo;
|
---|
584 | state.rstate.ra[0] = epsf;
|
---|
585 | state.rstate.ra[1] = epsg;
|
---|
586 | state.rstate.ra[2] = epsx;
|
---|
587 | state.rstate.ra[3] = v;
|
---|
588 | state.rstate.ra[4] = vv;
|
---|
589 | return result;
|
---|
590 | }
|
---|
591 |
|
---|
592 |
|
---|
593 | /*************************************************************************
|
---|
594 | L-BFGS algorithm results
|
---|
595 |
|
---|
596 | Called after MinLBFGSIteration returned False.
|
---|
597 |
|
---|
598 | Input parameters:
|
---|
599 | State - algorithm state (used by MinLBFGSIteration).
|
---|
600 |
|
---|
601 | Output parameters:
|
---|
602 | X - array[0..N-1], solution
|
---|
603 | Rep - optimization report:
|
---|
604 | * Rep.TerminationType completetion code:
|
---|
605 | * -2 rounding errors prevent further improvement.
|
---|
606 | X contains best point found.
|
---|
607 | * -1 incorrect parameters were specified
|
---|
608 | * 1 relative function improvement is no more than
|
---|
609 | EpsF.
|
---|
610 | * 2 relative step is no more than EpsX.
|
---|
611 | * 4 gradient norm is no more than EpsG
|
---|
612 | * 5 MaxIts steps was taken
|
---|
613 | * Rep.IterationsCount contains iterations count
|
---|
614 | * NFEV countains number of function calculations
|
---|
615 |
|
---|
616 | -- ALGLIB --
|
---|
617 | Copyright 14.11.2007 by Bochkanov Sergey
|
---|
618 | *************************************************************************/
|
---|
619 | public static void minlbfgsresults(ref lbfgsstate state,
|
---|
620 | ref double[] x,
|
---|
621 | ref lbfgsreport rep)
|
---|
622 | {
|
---|
623 | int i_ = 0;
|
---|
624 |
|
---|
625 | x = new double[state.n-1+1];
|
---|
626 | for(i_=0; i_<=state.n-1;i_++)
|
---|
627 | {
|
---|
628 | x[i_] = state.x[i_];
|
---|
629 | }
|
---|
630 | rep.iterationscount = state.repiterationscount;
|
---|
631 | rep.nfev = state.repnfev;
|
---|
632 | rep.terminationtype = state.repterminationtype;
|
---|
633 | }
|
---|
634 |
|
---|
635 |
|
---|
636 | /*************************************************************************
|
---|
637 | THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES A SUFFICIENT
|
---|
638 | DECREASE CONDITION AND A CURVATURE CONDITION.
|
---|
639 |
|
---|
640 | AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF UNCERTAINTY WITH
|
---|
641 | ENDPOINTS STX AND STY. THE INTERVAL OF UNCERTAINTY IS INITIALLY CHOSEN
|
---|
642 | SO THAT IT CONTAINS A MINIMIZER OF THE MODIFIED FUNCTION
|
---|
643 |
|
---|
644 | F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
|
---|
645 |
|
---|
646 | IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE
|
---|
647 | FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, THEN THE INTERVAL OF
|
---|
648 | UNCERTAINTY IS CHOSEN SO THAT IT CONTAINS A MINIMIZER OF F(X+STP*S).
|
---|
649 |
|
---|
650 | THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES THE SUFFICIENT
|
---|
651 | DECREASE CONDITION
|
---|
652 |
|
---|
653 | F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
|
---|
654 |
|
---|
655 | AND THE CURVATURE CONDITION
|
---|
656 |
|
---|
657 | ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
|
---|
658 |
|
---|
659 | IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION IS BOUNDED
|
---|
660 | BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES BOTH CONDITIONS.
|
---|
661 | IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH CONDITIONS, THEN THE
|
---|
662 | ALGORITHM USUALLY STOPS WHEN ROUNDING ERRORS PREVENT FURTHER PROGRESS.
|
---|
663 | IN THIS CASE STP ONLY SATISFIES THE SUFFICIENT DECREASE CONDITION.
|
---|
664 |
|
---|
665 | PARAMETERS DESCRIPRION
|
---|
666 |
|
---|
667 | N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER OF VARIABLES.
|
---|
668 |
|
---|
669 | X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE BASE POINT FOR
|
---|
670 | THE LINE SEARCH. ON OUTPUT IT CONTAINS X+STP*S.
|
---|
671 |
|
---|
672 | F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F AT X. ON OUTPUT
|
---|
673 | IT CONTAINS THE VALUE OF F AT X + STP*S.
|
---|
674 |
|
---|
675 | G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE GRADIENT OF F AT X.
|
---|
676 | ON OUTPUT IT CONTAINS THE GRADIENT OF F AT X + STP*S.
|
---|
677 |
|
---|
678 | S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE SEARCH DIRECTION.
|
---|
679 |
|
---|
680 | STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN INITIAL ESTIMATE
|
---|
681 | OF A SATISFACTORY STEP. ON OUTPUT STP CONTAINS THE FINAL ESTIMATE.
|
---|
682 |
|
---|
683 | FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION OCCURS WHEN THE
|
---|
684 | SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
|
---|
685 | SATISFIED.
|
---|
686 |
|
---|
687 | XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS WHEN THE RELATIVE
|
---|
688 | WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL.
|
---|
689 |
|
---|
690 | STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH SPECIFY LOWER AND
|
---|
691 | UPPER BOUNDS FOR THE STEP.
|
---|
692 |
|
---|
693 | MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION OCCURS WHEN THE
|
---|
694 | NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV BY THE END OF AN ITERATION.
|
---|
695 |
|
---|
696 | INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
|
---|
697 | INFO = 0 IMPROPER INPUT PARAMETERS.
|
---|
698 |
|
---|
699 | INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE
|
---|
700 | DIRECTIONAL DERIVATIVE CONDITION HOLD.
|
---|
701 |
|
---|
702 | INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
|
---|
703 | IS AT MOST XTOL.
|
---|
704 |
|
---|
705 | INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
|
---|
706 |
|
---|
707 | INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN.
|
---|
708 |
|
---|
709 | INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX.
|
---|
710 |
|
---|
711 | INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
|
---|
712 | THERE MAY NOT BE A STEP WHICH SATISFIES THE
|
---|
713 | SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
|
---|
714 | TOLERANCES MAY BE TOO SMALL.
|
---|
715 |
|
---|
716 | NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF CALLS TO FCN.
|
---|
717 |
|
---|
718 | WA IS A WORK ARRAY OF LENGTH N.
|
---|
719 |
|
---|
720 | ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
|
---|
721 | JORGE J. MORE', DAVID J. THUENTE
|
---|
722 | *************************************************************************/
|
---|
723 | private static void mcsrch(int n,
|
---|
724 | ref double[] x,
|
---|
725 | ref double f,
|
---|
726 | ref double[] g,
|
---|
727 | ref double[] s,
|
---|
728 | ref double stp,
|
---|
729 | ref int info,
|
---|
730 | ref int nfev,
|
---|
731 | ref double[] wa,
|
---|
732 | ref lbfgsstate state,
|
---|
733 | ref int stage)
|
---|
734 | {
|
---|
735 | double v = 0;
|
---|
736 | double p5 = 0;
|
---|
737 | double p66 = 0;
|
---|
738 | double zero = 0;
|
---|
739 | int i_ = 0;
|
---|
740 |
|
---|
741 |
|
---|
742 | //
|
---|
743 | // init
|
---|
744 | //
|
---|
745 | p5 = 0.5;
|
---|
746 | p66 = 0.66;
|
---|
747 | state.xtrapf = 4.0;
|
---|
748 | zero = 0;
|
---|
749 |
|
---|
750 | //
|
---|
751 | // Main cycle
|
---|
752 | //
|
---|
753 | while( true )
|
---|
754 | {
|
---|
755 | if( stage==0 )
|
---|
756 | {
|
---|
757 |
|
---|
758 | //
|
---|
759 | // NEXT
|
---|
760 | //
|
---|
761 | stage = 2;
|
---|
762 | continue;
|
---|
763 | }
|
---|
764 | if( stage==2 )
|
---|
765 | {
|
---|
766 | state.infoc = 1;
|
---|
767 | info = 0;
|
---|
768 |
|
---|
769 | //
|
---|
770 | // CHECK THE INPUT PARAMETERS FOR ERRORS.
|
---|
771 | //
|
---|
772 | if( n<=0 | (double)(stp)<=(double)(0) | (double)(ftol)<(double)(0) | (double)(gtol)<(double)(zero) | (double)(xtol)<(double)(zero) | (double)(stpmin)<(double)(zero) | (double)(stpmax)<(double)(stpmin) | maxfev<=0 )
|
---|
773 | {
|
---|
774 | stage = 0;
|
---|
775 | return;
|
---|
776 | }
|
---|
777 |
|
---|
778 | //
|
---|
779 | // COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
|
---|
780 | // AND CHECK THAT S IS A DESCENT DIRECTION.
|
---|
781 | //
|
---|
782 | v = 0.0;
|
---|
783 | for(i_=0; i_<=n-1;i_++)
|
---|
784 | {
|
---|
785 | v += g[i_]*s[i_];
|
---|
786 | }
|
---|
787 | state.dginit = v;
|
---|
788 | if( (double)(state.dginit)>=(double)(0) )
|
---|
789 | {
|
---|
790 | stage = 0;
|
---|
791 | return;
|
---|
792 | }
|
---|
793 |
|
---|
794 | //
|
---|
795 | // INITIALIZE LOCAL VARIABLES.
|
---|
796 | //
|
---|
797 | state.brackt = false;
|
---|
798 | state.stage1 = true;
|
---|
799 | nfev = 0;
|
---|
800 | state.finit = f;
|
---|
801 | state.dgtest = ftol*state.dginit;
|
---|
802 | state.width = stpmax-stpmin;
|
---|
803 | state.width1 = state.width/p5;
|
---|
804 | for(i_=0; i_<=n-1;i_++)
|
---|
805 | {
|
---|
806 | wa[i_] = x[i_];
|
---|
807 | }
|
---|
808 |
|
---|
809 | //
|
---|
810 | // THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
|
---|
811 | // FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
|
---|
812 | // THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
|
---|
813 | // FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
|
---|
814 | // THE INTERVAL OF UNCERTAINTY.
|
---|
815 | // THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
|
---|
816 | // FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
|
---|
817 | //
|
---|
818 | state.stx = 0;
|
---|
819 | state.fx = state.finit;
|
---|
820 | state.dgx = state.dginit;
|
---|
821 | state.sty = 0;
|
---|
822 | state.fy = state.finit;
|
---|
823 | state.dgy = state.dginit;
|
---|
824 |
|
---|
825 | //
|
---|
826 | // NEXT
|
---|
827 | //
|
---|
828 | stage = 3;
|
---|
829 | continue;
|
---|
830 | }
|
---|
831 | if( stage==3 )
|
---|
832 | {
|
---|
833 |
|
---|
834 | //
|
---|
835 | // START OF ITERATION.
|
---|
836 | //
|
---|
837 | // SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
|
---|
838 | // TO THE PRESENT INTERVAL OF UNCERTAINTY.
|
---|
839 | //
|
---|
840 | if( state.brackt )
|
---|
841 | {
|
---|
842 | if( (double)(state.stx)<(double)(state.sty) )
|
---|
843 | {
|
---|
844 | state.stmin = state.stx;
|
---|
845 | state.stmax = state.sty;
|
---|
846 | }
|
---|
847 | else
|
---|
848 | {
|
---|
849 | state.stmin = state.sty;
|
---|
850 | state.stmax = state.stx;
|
---|
851 | }
|
---|
852 | }
|
---|
853 | else
|
---|
854 | {
|
---|
855 | state.stmin = state.stx;
|
---|
856 | state.stmax = stp+state.xtrapf*(stp-state.stx);
|
---|
857 | }
|
---|
858 |
|
---|
859 | //
|
---|
860 | // FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
|
---|
861 | //
|
---|
862 | if( (double)(stp)>(double)(stpmax) )
|
---|
863 | {
|
---|
864 | stp = stpmax;
|
---|
865 | }
|
---|
866 | if( (double)(stp)<(double)(stpmin) )
|
---|
867 | {
|
---|
868 | stp = stpmin;
|
---|
869 | }
|
---|
870 |
|
---|
871 | //
|
---|
872 | // IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
|
---|
873 | // STP BE THE LOWEST POINT OBTAINED SO FAR.
|
---|
874 | //
|
---|
875 | if( state.brackt & ((double)(stp)<=(double)(state.stmin) | (double)(stp)>=(double)(state.stmax)) | nfev>=maxfev-1 | state.infoc==0 | state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
|
---|
876 | {
|
---|
877 | stp = state.stx;
|
---|
878 | }
|
---|
879 |
|
---|
880 | //
|
---|
881 | // EVALUATE THE FUNCTION AND GRADIENT AT STP
|
---|
882 | // AND COMPUTE THE DIRECTIONAL DERIVATIVE.
|
---|
883 | //
|
---|
884 | for(i_=0; i_<=n-1;i_++)
|
---|
885 | {
|
---|
886 | x[i_] = wa[i_];
|
---|
887 | }
|
---|
888 | for(i_=0; i_<=n-1;i_++)
|
---|
889 | {
|
---|
890 | x[i_] = x[i_] + stp*s[i_];
|
---|
891 | }
|
---|
892 |
|
---|
893 | //
|
---|
894 | // NEXT
|
---|
895 | //
|
---|
896 | stage = 4;
|
---|
897 | return;
|
---|
898 | }
|
---|
899 | if( stage==4 )
|
---|
900 | {
|
---|
901 | info = 0;
|
---|
902 | nfev = nfev+1;
|
---|
903 | v = 0.0;
|
---|
904 | for(i_=0; i_<=n-1;i_++)
|
---|
905 | {
|
---|
906 | v += g[i_]*s[i_];
|
---|
907 | }
|
---|
908 | state.dg = v;
|
---|
909 | state.ftest1 = state.finit+stp*state.dgtest;
|
---|
910 |
|
---|
911 | //
|
---|
912 | // TEST FOR CONVERGENCE.
|
---|
913 | //
|
---|
914 | if( state.brackt & ((double)(stp)<=(double)(state.stmin) | (double)(stp)>=(double)(state.stmax)) | state.infoc==0 )
|
---|
915 | {
|
---|
916 | info = 6;
|
---|
917 | }
|
---|
918 | if( (double)(stp)==(double)(stpmax) & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)<=(double)(state.dgtest) )
|
---|
919 | {
|
---|
920 | info = 5;
|
---|
921 | }
|
---|
922 | if( (double)(stp)==(double)(stpmin) & ((double)(f)>(double)(state.ftest1) | (double)(state.dg)>=(double)(state.dgtest)) )
|
---|
923 | {
|
---|
924 | info = 4;
|
---|
925 | }
|
---|
926 | if( nfev>=maxfev )
|
---|
927 | {
|
---|
928 | info = 3;
|
---|
929 | }
|
---|
930 | if( state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
|
---|
931 | {
|
---|
932 | info = 2;
|
---|
933 | }
|
---|
934 | if( (double)(f)<=(double)(state.ftest1) & (double)(Math.Abs(state.dg))<=(double)(-(gtol*state.dginit)) )
|
---|
935 | {
|
---|
936 | info = 1;
|
---|
937 | }
|
---|
938 |
|
---|
939 | //
|
---|
940 | // CHECK FOR TERMINATION.
|
---|
941 | //
|
---|
942 | if( info!=0 )
|
---|
943 | {
|
---|
944 | stage = 0;
|
---|
945 | return;
|
---|
946 | }
|
---|
947 |
|
---|
948 | //
|
---|
949 | // IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
|
---|
950 | // FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
|
---|
951 | //
|
---|
952 | if( state.stage1 & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)>=(double)(Math.Min(ftol, gtol)*state.dginit) )
|
---|
953 | {
|
---|
954 | state.stage1 = false;
|
---|
955 | }
|
---|
956 |
|
---|
957 | //
|
---|
958 | // A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
|
---|
959 | // WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
|
---|
960 | // FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
|
---|
961 | // DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
|
---|
962 | // OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
|
---|
963 | //
|
---|
964 | if( state.stage1 & (double)(f)<=(double)(state.fx) & (double)(f)>(double)(state.ftest1) )
|
---|
965 | {
|
---|
966 |
|
---|
967 | //
|
---|
968 | // DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
|
---|
969 | //
|
---|
970 | state.fm = f-stp*state.dgtest;
|
---|
971 | state.fxm = state.fx-state.stx*state.dgtest;
|
---|
972 | state.fym = state.fy-state.sty*state.dgtest;
|
---|
973 | state.dgm = state.dg-state.dgtest;
|
---|
974 | state.dgxm = state.dgx-state.dgtest;
|
---|
975 | state.dgym = state.dgy-state.dgtest;
|
---|
976 |
|
---|
977 | //
|
---|
978 | // CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
|
---|
979 | // AND TO COMPUTE THE NEW STEP.
|
---|
980 | //
|
---|
981 | mcstep(ref state.stx, ref state.fxm, ref state.dgxm, ref state.sty, ref state.fym, ref state.dgym, ref stp, state.fm, state.dgm, ref state.brackt, state.stmin, state.stmax, ref state.infoc);
|
---|
982 |
|
---|
983 | //
|
---|
984 | // RESET THE FUNCTION AND GRADIENT VALUES FOR F.
|
---|
985 | //
|
---|
986 | state.fx = state.fxm+state.stx*state.dgtest;
|
---|
987 | state.fy = state.fym+state.sty*state.dgtest;
|
---|
988 | state.dgx = state.dgxm+state.dgtest;
|
---|
989 | state.dgy = state.dgym+state.dgtest;
|
---|
990 | }
|
---|
991 | else
|
---|
992 | {
|
---|
993 |
|
---|
994 | //
|
---|
995 | // CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
|
---|
996 | // AND TO COMPUTE THE NEW STEP.
|
---|
997 | //
|
---|
998 | mcstep(ref state.stx, ref state.fx, ref state.dgx, ref state.sty, ref state.fy, ref state.dgy, ref stp, f, state.dg, ref state.brackt, state.stmin, state.stmax, ref state.infoc);
|
---|
999 | }
|
---|
1000 |
|
---|
1001 | //
|
---|
1002 | // FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
|
---|
1003 | // INTERVAL OF UNCERTAINTY.
|
---|
1004 | //
|
---|
1005 | if( state.brackt )
|
---|
1006 | {
|
---|
1007 | if( (double)(Math.Abs(state.sty-state.stx))>=(double)(p66*state.width1) )
|
---|
1008 | {
|
---|
1009 | stp = state.stx+p5*(state.sty-state.stx);
|
---|
1010 | }
|
---|
1011 | state.width1 = state.width;
|
---|
1012 | state.width = Math.Abs(state.sty-state.stx);
|
---|
1013 | }
|
---|
1014 |
|
---|
1015 | //
|
---|
1016 | // NEXT.
|
---|
1017 | //
|
---|
1018 | stage = 3;
|
---|
1019 | continue;
|
---|
1020 | }
|
---|
1021 | }
|
---|
1022 | }
|
---|
1023 |
|
---|
1024 |
|
---|
1025 | private static void mcstep(ref double stx,
|
---|
1026 | ref double fx,
|
---|
1027 | ref double dx,
|
---|
1028 | ref double sty,
|
---|
1029 | ref double fy,
|
---|
1030 | ref double dy,
|
---|
1031 | ref double stp,
|
---|
1032 | double fp,
|
---|
1033 | double dp,
|
---|
1034 | ref bool brackt,
|
---|
1035 | double stmin,
|
---|
1036 | double stmax,
|
---|
1037 | ref int info)
|
---|
1038 | {
|
---|
1039 | bool bound = new bool();
|
---|
1040 | double gamma = 0;
|
---|
1041 | double p = 0;
|
---|
1042 | double q = 0;
|
---|
1043 | double r = 0;
|
---|
1044 | double s = 0;
|
---|
1045 | double sgnd = 0;
|
---|
1046 | double stpc = 0;
|
---|
1047 | double stpf = 0;
|
---|
1048 | double stpq = 0;
|
---|
1049 | double theta = 0;
|
---|
1050 |
|
---|
1051 | info = 0;
|
---|
1052 |
|
---|
1053 | //
|
---|
1054 | // CHECK THE INPUT PARAMETERS FOR ERRORS.
|
---|
1055 | //
|
---|
1056 | if( brackt & ((double)(stp)<=(double)(Math.Min(stx, sty)) | (double)(stp)>=(double)(Math.Max(stx, sty))) | (double)(dx*(stp-stx))>=(double)(0) | (double)(stmax)<(double)(stmin) )
|
---|
1057 | {
|
---|
1058 | return;
|
---|
1059 | }
|
---|
1060 |
|
---|
1061 | //
|
---|
1062 | // DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
|
---|
1063 | //
|
---|
1064 | sgnd = dp*(dx/Math.Abs(dx));
|
---|
1065 |
|
---|
1066 | //
|
---|
1067 | // FIRST CASE. A HIGHER FUNCTION VALUE.
|
---|
1068 | // THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
|
---|
1069 | // TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
|
---|
1070 | // ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
|
---|
1071 | //
|
---|
1072 | if( (double)(fp)>(double)(fx) )
|
---|
1073 | {
|
---|
1074 | info = 1;
|
---|
1075 | bound = true;
|
---|
1076 | theta = 3*(fx-fp)/(stp-stx)+dx+dp;
|
---|
1077 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
|
---|
1078 | gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dx/s*(dp/s));
|
---|
1079 | if( (double)(stp)<(double)(stx) )
|
---|
1080 | {
|
---|
1081 | gamma = -gamma;
|
---|
1082 | }
|
---|
1083 | p = gamma-dx+theta;
|
---|
1084 | q = gamma-dx+gamma+dp;
|
---|
1085 | r = p/q;
|
---|
1086 | stpc = stx+r*(stp-stx);
|
---|
1087 | stpq = stx+dx/((fx-fp)/(stp-stx)+dx)/2*(stp-stx);
|
---|
1088 | if( (double)(Math.Abs(stpc-stx))<(double)(Math.Abs(stpq-stx)) )
|
---|
1089 | {
|
---|
1090 | stpf = stpc;
|
---|
1091 | }
|
---|
1092 | else
|
---|
1093 | {
|
---|
1094 | stpf = stpc+(stpq-stpc)/2;
|
---|
1095 | }
|
---|
1096 | brackt = true;
|
---|
1097 | }
|
---|
1098 | else
|
---|
1099 | {
|
---|
1100 | if( (double)(sgnd)<(double)(0) )
|
---|
1101 | {
|
---|
1102 |
|
---|
1103 | //
|
---|
1104 | // SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
|
---|
1105 | // OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
|
---|
1106 | // STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
|
---|
1107 | // THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
|
---|
1108 | //
|
---|
1109 | info = 2;
|
---|
1110 | bound = false;
|
---|
1111 | theta = 3*(fx-fp)/(stp-stx)+dx+dp;
|
---|
1112 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
|
---|
1113 | gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dx/s*(dp/s));
|
---|
1114 | if( (double)(stp)>(double)(stx) )
|
---|
1115 | {
|
---|
1116 | gamma = -gamma;
|
---|
1117 | }
|
---|
1118 | p = gamma-dp+theta;
|
---|
1119 | q = gamma-dp+gamma+dx;
|
---|
1120 | r = p/q;
|
---|
1121 | stpc = stp+r*(stx-stp);
|
---|
1122 | stpq = stp+dp/(dp-dx)*(stx-stp);
|
---|
1123 | if( (double)(Math.Abs(stpc-stp))>(double)(Math.Abs(stpq-stp)) )
|
---|
1124 | {
|
---|
1125 | stpf = stpc;
|
---|
1126 | }
|
---|
1127 | else
|
---|
1128 | {
|
---|
1129 | stpf = stpq;
|
---|
1130 | }
|
---|
1131 | brackt = true;
|
---|
1132 | }
|
---|
1133 | else
|
---|
1134 | {
|
---|
1135 | if( (double)(Math.Abs(dp))<(double)(Math.Abs(dx)) )
|
---|
1136 | {
|
---|
1137 |
|
---|
1138 | //
|
---|
1139 | // THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
|
---|
1140 | // SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
|
---|
1141 | // THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
|
---|
1142 | // IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
|
---|
1143 | // IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
|
---|
1144 | // EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
|
---|
1145 | // COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
|
---|
1146 | // CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
|
---|
1147 | //
|
---|
1148 | info = 3;
|
---|
1149 | bound = true;
|
---|
1150 | theta = 3*(fx-fp)/(stp-stx)+dx+dp;
|
---|
1151 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
|
---|
1152 |
|
---|
1153 | //
|
---|
1154 | // THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
|
---|
1155 | // TO INFINITY IN THE DIRECTION OF THE STEP.
|
---|
1156 | //
|
---|
1157 | gamma = s*Math.Sqrt(Math.Max(0, AP.Math.Sqr(theta/s)-dx/s*(dp/s)));
|
---|
1158 | if( (double)(stp)>(double)(stx) )
|
---|
1159 | {
|
---|
1160 | gamma = -gamma;
|
---|
1161 | }
|
---|
1162 | p = gamma-dp+theta;
|
---|
1163 | q = gamma+(dx-dp)+gamma;
|
---|
1164 | r = p/q;
|
---|
1165 | if( (double)(r)<(double)(0) & (double)(gamma)!=(double)(0) )
|
---|
1166 | {
|
---|
1167 | stpc = stp+r*(stx-stp);
|
---|
1168 | }
|
---|
1169 | else
|
---|
1170 | {
|
---|
1171 | if( (double)(stp)>(double)(stx) )
|
---|
1172 | {
|
---|
1173 | stpc = stmax;
|
---|
1174 | }
|
---|
1175 | else
|
---|
1176 | {
|
---|
1177 | stpc = stmin;
|
---|
1178 | }
|
---|
1179 | }
|
---|
1180 | stpq = stp+dp/(dp-dx)*(stx-stp);
|
---|
1181 | if( brackt )
|
---|
1182 | {
|
---|
1183 | if( (double)(Math.Abs(stp-stpc))<(double)(Math.Abs(stp-stpq)) )
|
---|
1184 | {
|
---|
1185 | stpf = stpc;
|
---|
1186 | }
|
---|
1187 | else
|
---|
1188 | {
|
---|
1189 | stpf = stpq;
|
---|
1190 | }
|
---|
1191 | }
|
---|
1192 | else
|
---|
1193 | {
|
---|
1194 | if( (double)(Math.Abs(stp-stpc))>(double)(Math.Abs(stp-stpq)) )
|
---|
1195 | {
|
---|
1196 | stpf = stpc;
|
---|
1197 | }
|
---|
1198 | else
|
---|
1199 | {
|
---|
1200 | stpf = stpq;
|
---|
1201 | }
|
---|
1202 | }
|
---|
1203 | }
|
---|
1204 | else
|
---|
1205 | {
|
---|
1206 |
|
---|
1207 | //
|
---|
1208 | // FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
|
---|
1209 | // SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
|
---|
1210 | // NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
|
---|
1211 | // IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
|
---|
1212 | //
|
---|
1213 | info = 4;
|
---|
1214 | bound = false;
|
---|
1215 | if( brackt )
|
---|
1216 | {
|
---|
1217 | theta = 3*(fp-fy)/(sty-stp)+dy+dp;
|
---|
1218 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dy), Math.Abs(dp)));
|
---|
1219 | gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dy/s*(dp/s));
|
---|
1220 | if( (double)(stp)>(double)(sty) )
|
---|
1221 | {
|
---|
1222 | gamma = -gamma;
|
---|
1223 | }
|
---|
1224 | p = gamma-dp+theta;
|
---|
1225 | q = gamma-dp+gamma+dy;
|
---|
1226 | r = p/q;
|
---|
1227 | stpc = stp+r*(sty-stp);
|
---|
1228 | stpf = stpc;
|
---|
1229 | }
|
---|
1230 | else
|
---|
1231 | {
|
---|
1232 | if( (double)(stp)>(double)(stx) )
|
---|
1233 | {
|
---|
1234 | stpf = stmax;
|
---|
1235 | }
|
---|
1236 | else
|
---|
1237 | {
|
---|
1238 | stpf = stmin;
|
---|
1239 | }
|
---|
1240 | }
|
---|
1241 | }
|
---|
1242 | }
|
---|
1243 | }
|
---|
1244 |
|
---|
1245 | //
|
---|
1246 | // UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
|
---|
1247 | // DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
|
---|
1248 | //
|
---|
1249 | if( (double)(fp)>(double)(fx) )
|
---|
1250 | {
|
---|
1251 | sty = stp;
|
---|
1252 | fy = fp;
|
---|
1253 | dy = dp;
|
---|
1254 | }
|
---|
1255 | else
|
---|
1256 | {
|
---|
1257 | if( (double)(sgnd)<(double)(0.0) )
|
---|
1258 | {
|
---|
1259 | sty = stx;
|
---|
1260 | fy = fx;
|
---|
1261 | dy = dx;
|
---|
1262 | }
|
---|
1263 | stx = stp;
|
---|
1264 | fx = fp;
|
---|
1265 | dx = dp;
|
---|
1266 | }
|
---|
1267 |
|
---|
1268 | //
|
---|
1269 | // COMPUTE THE NEW STEP AND SAFEGUARD IT.
|
---|
1270 | //
|
---|
1271 | stpf = Math.Min(stmax, stpf);
|
---|
1272 | stpf = Math.Max(stmin, stpf);
|
---|
1273 | stp = stpf;
|
---|
1274 | if( brackt & bound )
|
---|
1275 | {
|
---|
1276 | if( (double)(sty)>(double)(stx) )
|
---|
1277 | {
|
---|
1278 | stp = Math.Min(stx+0.66*(sty-stx), stp);
|
---|
1279 | }
|
---|
1280 | else
|
---|
1281 | {
|
---|
1282 | stp = Math.Max(stx+0.66*(sty-stx), stp);
|
---|
1283 | }
|
---|
1284 | }
|
---|
1285 | }
|
---|
1286 | }
|
---|
1287 | }
|
---|