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