[2563] | 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 | }
|
---|