[2806] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 2009, 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 minlm
|
---|
| 26 | {
|
---|
| 27 | public struct lmstate
|
---|
| 28 | {
|
---|
| 29 | public bool wrongparams;
|
---|
| 30 | public int n;
|
---|
| 31 | public int m;
|
---|
| 32 | public double epsf;
|
---|
| 33 | public double epsx;
|
---|
| 34 | public int maxits;
|
---|
| 35 | public int flags;
|
---|
| 36 | public int usermode;
|
---|
| 37 | public double[] x;
|
---|
| 38 | public double f;
|
---|
| 39 | public double[] fi;
|
---|
| 40 | public double[,] j;
|
---|
| 41 | public double[,] h;
|
---|
| 42 | public double[] g;
|
---|
| 43 | public bool needf;
|
---|
| 44 | public bool needfg;
|
---|
| 45 | public bool needfgh;
|
---|
| 46 | public bool needfij;
|
---|
| 47 | public bool xupdated;
|
---|
| 48 | public lbfgs.lbfgsstate internalstate;
|
---|
| 49 | public lbfgs.lbfgsreport internalrep;
|
---|
| 50 | public double[] xprec;
|
---|
| 51 | public double[] xbase;
|
---|
| 52 | public double[] xdir;
|
---|
| 53 | public double[] gbase;
|
---|
| 54 | public double[,] rawmodel;
|
---|
| 55 | public double[,] model;
|
---|
| 56 | public double[] work;
|
---|
| 57 | public AP.rcommstate rstate;
|
---|
| 58 | public int repiterationscount;
|
---|
| 59 | public int repterminationtype;
|
---|
| 60 | public int repnfunc;
|
---|
| 61 | public int repnjac;
|
---|
| 62 | public int repngrad;
|
---|
| 63 | public int repnhess;
|
---|
| 64 | public int repncholesky;
|
---|
| 65 | public int solverinfo;
|
---|
| 66 | public densesolver.densesolverreport solverrep;
|
---|
| 67 | };
|
---|
| 68 |
|
---|
| 69 |
|
---|
| 70 | public struct lmreport
|
---|
| 71 | {
|
---|
| 72 | public int iterationscount;
|
---|
| 73 | public int terminationtype;
|
---|
| 74 | public int nfunc;
|
---|
| 75 | public int njac;
|
---|
| 76 | public int ngrad;
|
---|
| 77 | public int nhess;
|
---|
| 78 | public int ncholesky;
|
---|
| 79 | };
|
---|
| 80 |
|
---|
| 81 |
|
---|
| 82 |
|
---|
| 83 |
|
---|
| 84 | public const int lmmodefj = 0;
|
---|
| 85 | public const int lmmodefgj = 1;
|
---|
| 86 | public const int lmmodefgh = 2;
|
---|
| 87 | public const int lmflagnoprelbfgs = 1;
|
---|
| 88 | public const int lmflagnointlbfgs = 2;
|
---|
| 89 | public const int lmprelbfgsm = 5;
|
---|
| 90 | public const int lmintlbfgsits = 5;
|
---|
| 91 | public const int lbfgsnorealloc = 1;
|
---|
| 92 |
|
---|
| 93 |
|
---|
| 94 | /*************************************************************************
|
---|
| 95 | LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
|
---|
| 96 |
|
---|
| 97 | Optimization using function gradient and Hessian. Algorithm - Levenberg-
|
---|
| 98 | Marquardt modification with L-BFGS pre-optimization and internal
|
---|
| 99 | pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
|
---|
| 100 |
|
---|
| 101 | Function F has general form (not "sum-of-squares"):
|
---|
| 102 |
|
---|
| 103 | F = F(x[0], ..., x[n-1])
|
---|
| 104 |
|
---|
| 105 | EXAMPLE
|
---|
| 106 |
|
---|
| 107 | See HTML-documentation.
|
---|
| 108 |
|
---|
| 109 | INPUT PARAMETERS:
|
---|
| 110 | N - dimension, N>1
|
---|
| 111 | X - initial solution, array[0..N-1]
|
---|
| 112 | EpsF - stopping criterion. Algorithm stops if
|
---|
| 113 | |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
|
---|
| 114 | EpsX - stopping criterion. Algorithm stops if
|
---|
| 115 | |X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
|
---|
| 116 | MaxIts - stopping criterion. Algorithm stops after MaxIts iterations.
|
---|
| 117 | MaxIts=0 means no stopping criterion.
|
---|
| 118 |
|
---|
| 119 | OUTPUT PARAMETERS:
|
---|
| 120 | State - structure which stores algorithm state between subsequent
|
---|
| 121 | calls of MinLMIteration. Used for reverse communication.
|
---|
| 122 | This structure should be passed to MinLMIteration subroutine.
|
---|
| 123 |
|
---|
| 124 | See also MinLMIteration, MinLMResults.
|
---|
| 125 |
|
---|
| 126 | NOTE
|
---|
| 127 |
|
---|
| 128 | Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
|
---|
| 129 | stopping criterion selection (small EpsX).
|
---|
| 130 |
|
---|
| 131 | -- ALGLIB --
|
---|
| 132 | Copyright 30.03.2009 by Bochkanov Sergey
|
---|
| 133 | *************************************************************************/
|
---|
| 134 | public static void minlmfgh(int n,
|
---|
| 135 | ref double[] x,
|
---|
| 136 | double epsf,
|
---|
| 137 | double epsx,
|
---|
| 138 | int maxits,
|
---|
| 139 | ref lmstate state)
|
---|
| 140 | {
|
---|
| 141 | int i_ = 0;
|
---|
| 142 |
|
---|
| 143 |
|
---|
| 144 | //
|
---|
| 145 | // Prepare RComm
|
---|
| 146 | //
|
---|
| 147 | state.rstate.ia = new int[3+1];
|
---|
| 148 | state.rstate.ba = new bool[0+1];
|
---|
| 149 | state.rstate.ra = new double[8+1];
|
---|
| 150 | state.rstate.stage = -1;
|
---|
| 151 |
|
---|
| 152 | //
|
---|
| 153 | // prepare internal structures
|
---|
| 154 | //
|
---|
| 155 | lmprepare(n, 0, true, ref state);
|
---|
| 156 |
|
---|
| 157 | //
|
---|
| 158 | // initialize, check parameters
|
---|
| 159 | //
|
---|
| 160 | state.xupdated = false;
|
---|
| 161 | state.n = n;
|
---|
| 162 | state.m = 0;
|
---|
| 163 | state.epsf = epsf;
|
---|
| 164 | state.epsx = epsx;
|
---|
| 165 | state.maxits = maxits;
|
---|
| 166 | state.flags = 0;
|
---|
| 167 | if( (double)(state.epsf)==(double)(0) & (double)(state.epsx)==(double)(0) & state.maxits==0 )
|
---|
| 168 | {
|
---|
| 169 | state.epsx = 1.0E-6;
|
---|
| 170 | }
|
---|
| 171 | state.usermode = lmmodefgh;
|
---|
| 172 | state.wrongparams = false;
|
---|
| 173 | if( n<1 | (double)(epsf)<(double)(0) | (double)(epsx)<(double)(0) | maxits<0 )
|
---|
| 174 | {
|
---|
| 175 | state.wrongparams = true;
|
---|
| 176 | return;
|
---|
| 177 | }
|
---|
| 178 | for(i_=0; i_<=n-1;i_++)
|
---|
| 179 | {
|
---|
| 180 | state.x[i_] = x[i_];
|
---|
| 181 | }
|
---|
| 182 | }
|
---|
| 183 |
|
---|
| 184 |
|
---|
| 185 | /*************************************************************************
|
---|
| 186 | LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
|
---|
| 187 |
|
---|
| 188 | Optimization using function gradient and Jacobian. Algorithm - Levenberg-
|
---|
| 189 | Marquardt modification with L-BFGS pre-optimization and internal
|
---|
| 190 | pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
|
---|
| 191 |
|
---|
| 192 | Function F is represented as sum of squares:
|
---|
| 193 |
|
---|
| 194 | F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
|
---|
| 195 |
|
---|
| 196 | EXAMPLE
|
---|
| 197 |
|
---|
| 198 | See HTML-documentation.
|
---|
| 199 |
|
---|
| 200 | INPUT PARAMETERS:
|
---|
| 201 | N - dimension, N>1
|
---|
| 202 | M - number of functions f[i]
|
---|
| 203 | X - initial solution, array[0..N-1]
|
---|
| 204 | EpsF - stopping criterion. Algorithm stops if
|
---|
| 205 | |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
|
---|
| 206 | EpsX - stopping criterion. Algorithm stops if
|
---|
| 207 | |X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
|
---|
| 208 | MaxIts - stopping criterion. Algorithm stops after MaxIts iterations.
|
---|
| 209 | MaxIts=0 means no stopping criterion.
|
---|
| 210 |
|
---|
| 211 | OUTPUT PARAMETERS:
|
---|
| 212 | State - structure which stores algorithm state between subsequent
|
---|
| 213 | calls of MinLMIteration. Used for reverse communication.
|
---|
| 214 | This structure should be passed to MinLMIteration subroutine.
|
---|
| 215 |
|
---|
| 216 | See also MinLMIteration, MinLMResults.
|
---|
| 217 |
|
---|
| 218 | NOTE
|
---|
| 219 |
|
---|
| 220 | Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
|
---|
| 221 | stopping criterion selection (small EpsX).
|
---|
| 222 |
|
---|
| 223 | -- ALGLIB --
|
---|
| 224 | Copyright 30.03.2009 by Bochkanov Sergey
|
---|
| 225 | *************************************************************************/
|
---|
| 226 | public static void minlmfgj(int n,
|
---|
| 227 | int m,
|
---|
| 228 | ref double[] x,
|
---|
| 229 | double epsf,
|
---|
| 230 | double epsx,
|
---|
| 231 | int maxits,
|
---|
| 232 | ref lmstate state)
|
---|
| 233 | {
|
---|
| 234 | int i_ = 0;
|
---|
| 235 |
|
---|
| 236 |
|
---|
| 237 | //
|
---|
| 238 | // Prepare RComm
|
---|
| 239 | //
|
---|
| 240 | state.rstate.ia = new int[3+1];
|
---|
| 241 | state.rstate.ba = new bool[0+1];
|
---|
| 242 | state.rstate.ra = new double[8+1];
|
---|
| 243 | state.rstate.stage = -1;
|
---|
| 244 |
|
---|
| 245 | //
|
---|
| 246 | // prepare internal structures
|
---|
| 247 | //
|
---|
| 248 | lmprepare(n, m, true, ref state);
|
---|
| 249 |
|
---|
| 250 | //
|
---|
| 251 | // initialize, check parameters
|
---|
| 252 | //
|
---|
| 253 | state.xupdated = false;
|
---|
| 254 | state.n = n;
|
---|
| 255 | state.m = m;
|
---|
| 256 | state.epsf = epsf;
|
---|
| 257 | state.epsx = epsx;
|
---|
| 258 | state.maxits = maxits;
|
---|
| 259 | state.flags = 0;
|
---|
| 260 | if( (double)(state.epsf)==(double)(0) & (double)(state.epsx)==(double)(0) & state.maxits==0 )
|
---|
| 261 | {
|
---|
| 262 | state.epsx = 1.0E-6;
|
---|
| 263 | }
|
---|
| 264 | state.usermode = lmmodefgj;
|
---|
| 265 | state.wrongparams = false;
|
---|
| 266 | if( n<1 | m<1 | (double)(epsf)<(double)(0) | (double)(epsx)<(double)(0) | maxits<0 )
|
---|
| 267 | {
|
---|
| 268 | state.wrongparams = true;
|
---|
| 269 | return;
|
---|
| 270 | }
|
---|
| 271 | for(i_=0; i_<=n-1;i_++)
|
---|
| 272 | {
|
---|
| 273 | state.x[i_] = x[i_];
|
---|
| 274 | }
|
---|
| 275 | }
|
---|
| 276 |
|
---|
| 277 |
|
---|
| 278 | /*************************************************************************
|
---|
| 279 | CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION
|
---|
| 280 |
|
---|
| 281 | Optimization using Jacobi matrix. Algorithm - classic Levenberg-Marquardt
|
---|
| 282 | method.
|
---|
| 283 |
|
---|
| 284 | Function F is represented as sum of squares:
|
---|
| 285 |
|
---|
| 286 | F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
|
---|
| 287 |
|
---|
| 288 | EXAMPLE
|
---|
| 289 |
|
---|
| 290 | See HTML-documentation.
|
---|
| 291 |
|
---|
| 292 | INPUT PARAMETERS:
|
---|
| 293 | N - dimension, N>1
|
---|
| 294 | M - number of functions f[i]
|
---|
| 295 | X - initial solution, array[0..N-1]
|
---|
| 296 | EpsF - stopping criterion. Algorithm stops if
|
---|
| 297 | |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
|
---|
| 298 | EpsX - stopping criterion. Algorithm stops if
|
---|
| 299 | |X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
|
---|
| 300 | MaxIts - stopping criterion. Algorithm stops after MaxIts iterations.
|
---|
| 301 | MaxIts=0 means no stopping criterion.
|
---|
| 302 |
|
---|
| 303 | OUTPUT PARAMETERS:
|
---|
| 304 | State - structure which stores algorithm state between subsequent
|
---|
| 305 | calls of MinLMIteration. Used for reverse communication.
|
---|
| 306 | This structure should be passed to MinLMIteration subroutine.
|
---|
| 307 |
|
---|
| 308 | See also MinLMIteration, MinLMResults.
|
---|
| 309 |
|
---|
| 310 | NOTE
|
---|
| 311 |
|
---|
| 312 | Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
|
---|
| 313 | stopping criterion selection (small EpsX).
|
---|
| 314 |
|
---|
| 315 | -- ALGLIB --
|
---|
| 316 | Copyright 30.03.2009 by Bochkanov Sergey
|
---|
| 317 | *************************************************************************/
|
---|
| 318 | public static void minlmfj(int n,
|
---|
| 319 | int m,
|
---|
| 320 | ref double[] x,
|
---|
| 321 | double epsf,
|
---|
| 322 | double epsx,
|
---|
| 323 | int maxits,
|
---|
| 324 | ref lmstate state)
|
---|
| 325 | {
|
---|
| 326 | int i_ = 0;
|
---|
| 327 |
|
---|
| 328 |
|
---|
| 329 | //
|
---|
| 330 | // Prepare RComm
|
---|
| 331 | //
|
---|
| 332 | state.rstate.ia = new int[3+1];
|
---|
| 333 | state.rstate.ba = new bool[0+1];
|
---|
| 334 | state.rstate.ra = new double[8+1];
|
---|
| 335 | state.rstate.stage = -1;
|
---|
| 336 |
|
---|
| 337 | //
|
---|
| 338 | // prepare internal structures
|
---|
| 339 | //
|
---|
| 340 | lmprepare(n, m, true, ref state);
|
---|
| 341 |
|
---|
| 342 | //
|
---|
| 343 | // initialize, check parameters
|
---|
| 344 | //
|
---|
| 345 | state.xupdated = false;
|
---|
| 346 | state.n = n;
|
---|
| 347 | state.m = m;
|
---|
| 348 | state.epsf = epsf;
|
---|
| 349 | state.epsx = epsx;
|
---|
| 350 | state.maxits = maxits;
|
---|
| 351 | state.flags = 0;
|
---|
| 352 | if( (double)(state.epsf)==(double)(0) & (double)(state.epsx)==(double)(0) & state.maxits==0 )
|
---|
| 353 | {
|
---|
| 354 | state.epsx = 1.0E-6;
|
---|
| 355 | }
|
---|
| 356 | state.usermode = lmmodefj;
|
---|
| 357 | state.wrongparams = false;
|
---|
| 358 | if( n<1 | m<1 | (double)(epsf)<(double)(0) | (double)(epsx)<(double)(0) | maxits<0 )
|
---|
| 359 | {
|
---|
| 360 | state.wrongparams = true;
|
---|
| 361 | return;
|
---|
| 362 | }
|
---|
| 363 | for(i_=0; i_<=n-1;i_++)
|
---|
| 364 | {
|
---|
| 365 | state.x[i_] = x[i_];
|
---|
| 366 | }
|
---|
| 367 | }
|
---|
| 368 |
|
---|
| 369 |
|
---|
| 370 | /*************************************************************************
|
---|
| 371 | One Levenberg-Marquardt iteration.
|
---|
| 372 |
|
---|
| 373 | Called after inialization of State structure with MinLMXXX subroutine.
|
---|
| 374 | See HTML docs for examples.
|
---|
| 375 |
|
---|
| 376 | Input parameters:
|
---|
| 377 | State - structure which stores algorithm state between subsequent
|
---|
| 378 | calls and which is used for reverse communication. Must be
|
---|
| 379 | initialized with MinLMXXX call first.
|
---|
| 380 |
|
---|
| 381 | If subroutine returned False, iterative algorithm has converged.
|
---|
| 382 |
|
---|
| 383 | If subroutine returned True, then:
|
---|
| 384 | * if State.NeedF=True, - function value F at State.X[0..N-1]
|
---|
| 385 | is required
|
---|
| 386 | * if State.NeedFG=True - function value F and gradient G
|
---|
| 387 | are required
|
---|
| 388 | * if State.NeedFiJ=True - function vector f[i] and Jacobi matrix J
|
---|
| 389 | are required
|
---|
| 390 | * if State.NeedFGH=True - function value F, gradient G and Hesian H
|
---|
| 391 | are required
|
---|
| 392 |
|
---|
| 393 | One and only one of this fields can be set at time.
|
---|
| 394 |
|
---|
| 395 | Results are stored:
|
---|
| 396 | * function value - in LMState.F
|
---|
| 397 | * gradient - in LMState.G[0..N-1]
|
---|
| 398 | * Jacobi matrix - in LMState.J[0..M-1,0..N-1]
|
---|
| 399 | * Hessian - in LMState.H[0..N-1,0..N-1]
|
---|
| 400 |
|
---|
| 401 | -- ALGLIB --
|
---|
| 402 | Copyright 10.03.2009 by Bochkanov Sergey
|
---|
| 403 | *************************************************************************/
|
---|
| 404 | public static bool minlmiteration(ref lmstate state)
|
---|
| 405 | {
|
---|
| 406 | bool result = new bool();
|
---|
| 407 | int n = 0;
|
---|
| 408 | int m = 0;
|
---|
| 409 | int i = 0;
|
---|
| 410 | double xnorm = 0;
|
---|
| 411 | double stepnorm = 0;
|
---|
| 412 | bool spd = new bool();
|
---|
| 413 | double fbase = 0;
|
---|
| 414 | double fnew = 0;
|
---|
| 415 | double lambda = 0;
|
---|
| 416 | double nu = 0;
|
---|
| 417 | double lambdaup = 0;
|
---|
| 418 | double lambdadown = 0;
|
---|
| 419 | int lbfgsflags = 0;
|
---|
| 420 | double v = 0;
|
---|
| 421 | int i_ = 0;
|
---|
| 422 |
|
---|
| 423 |
|
---|
| 424 | //
|
---|
| 425 | // Reverse communication preparations
|
---|
| 426 | // I know it looks ugly, but it works the same way
|
---|
| 427 | // anywhere from C++ to Python.
|
---|
| 428 | //
|
---|
| 429 | // This code initializes locals by:
|
---|
| 430 | // * random values determined during code
|
---|
| 431 | // generation - on first subroutine call
|
---|
| 432 | // * values from previous call - on subsequent calls
|
---|
| 433 | //
|
---|
| 434 | if( state.rstate.stage>=0 )
|
---|
| 435 | {
|
---|
| 436 | n = state.rstate.ia[0];
|
---|
| 437 | m = state.rstate.ia[1];
|
---|
| 438 | i = state.rstate.ia[2];
|
---|
| 439 | lbfgsflags = state.rstate.ia[3];
|
---|
| 440 | spd = state.rstate.ba[0];
|
---|
| 441 | xnorm = state.rstate.ra[0];
|
---|
| 442 | stepnorm = state.rstate.ra[1];
|
---|
| 443 | fbase = state.rstate.ra[2];
|
---|
| 444 | fnew = state.rstate.ra[3];
|
---|
| 445 | lambda = state.rstate.ra[4];
|
---|
| 446 | nu = state.rstate.ra[5];
|
---|
| 447 | lambdaup = state.rstate.ra[6];
|
---|
| 448 | lambdadown = state.rstate.ra[7];
|
---|
| 449 | v = state.rstate.ra[8];
|
---|
| 450 | }
|
---|
| 451 | else
|
---|
| 452 | {
|
---|
| 453 | n = -983;
|
---|
| 454 | m = -989;
|
---|
| 455 | i = -834;
|
---|
| 456 | lbfgsflags = 900;
|
---|
| 457 | spd = true;
|
---|
| 458 | xnorm = 364;
|
---|
| 459 | stepnorm = 214;
|
---|
| 460 | fbase = -338;
|
---|
| 461 | fnew = -686;
|
---|
| 462 | lambda = 912;
|
---|
| 463 | nu = 585;
|
---|
| 464 | lambdaup = 497;
|
---|
| 465 | lambdadown = -271;
|
---|
| 466 | v = -581;
|
---|
| 467 | }
|
---|
| 468 | if( state.rstate.stage==0 )
|
---|
| 469 | {
|
---|
| 470 | goto lbl_0;
|
---|
| 471 | }
|
---|
| 472 | if( state.rstate.stage==1 )
|
---|
| 473 | {
|
---|
| 474 | goto lbl_1;
|
---|
| 475 | }
|
---|
| 476 | if( state.rstate.stage==2 )
|
---|
| 477 | {
|
---|
| 478 | goto lbl_2;
|
---|
| 479 | }
|
---|
| 480 | if( state.rstate.stage==3 )
|
---|
| 481 | {
|
---|
| 482 | goto lbl_3;
|
---|
| 483 | }
|
---|
| 484 | if( state.rstate.stage==4 )
|
---|
| 485 | {
|
---|
| 486 | goto lbl_4;
|
---|
| 487 | }
|
---|
| 488 | if( state.rstate.stage==5 )
|
---|
| 489 | {
|
---|
| 490 | goto lbl_5;
|
---|
| 491 | }
|
---|
| 492 | if( state.rstate.stage==6 )
|
---|
| 493 | {
|
---|
| 494 | goto lbl_6;
|
---|
| 495 | }
|
---|
| 496 |
|
---|
| 497 | //
|
---|
| 498 | // Routine body
|
---|
| 499 | //
|
---|
| 500 | System.Diagnostics.Debug.Assert(state.usermode==lmmodefj | state.usermode==lmmodefgj | state.usermode==lmmodefgh, "LM: internal error");
|
---|
| 501 | if( state.wrongparams )
|
---|
| 502 | {
|
---|
| 503 | state.repterminationtype = -1;
|
---|
| 504 | result = false;
|
---|
| 505 | return result;
|
---|
| 506 | }
|
---|
| 507 |
|
---|
| 508 | //
|
---|
| 509 | // prepare params
|
---|
| 510 | //
|
---|
| 511 | n = state.n;
|
---|
| 512 | m = state.m;
|
---|
| 513 | lambdaup = 10;
|
---|
| 514 | lambdadown = 0.3;
|
---|
| 515 | nu = 2;
|
---|
| 516 | lbfgsflags = 0;
|
---|
| 517 |
|
---|
| 518 | //
|
---|
| 519 | // if we have F and G
|
---|
| 520 | //
|
---|
| 521 | if( ! ((state.usermode==lmmodefgj | state.usermode==lmmodefgh) & state.flags/lmflagnoprelbfgs%2==0) )
|
---|
| 522 | {
|
---|
| 523 | goto lbl_7;
|
---|
| 524 | }
|
---|
| 525 |
|
---|
| 526 | //
|
---|
| 527 | // First stage of the hybrid algorithm: LBFGS
|
---|
| 528 | //
|
---|
| 529 | lbfgs.minlbfgs(n, Math.Min(n, lmprelbfgsm), ref state.x, 0.0, 0.0, 0.0, Math.Max(5, n), 0, ref state.internalstate);
|
---|
| 530 | lbl_9:
|
---|
| 531 | if( ! lbfgs.minlbfgsiteration(ref state.internalstate) )
|
---|
| 532 | {
|
---|
| 533 | goto lbl_10;
|
---|
| 534 | }
|
---|
| 535 |
|
---|
| 536 | //
|
---|
| 537 | // RComm
|
---|
| 538 | //
|
---|
| 539 | for(i_=0; i_<=n-1;i_++)
|
---|
| 540 | {
|
---|
| 541 | state.x[i_] = state.internalstate.x[i_];
|
---|
| 542 | }
|
---|
| 543 | lmclearrequestfields(ref state);
|
---|
| 544 | state.needfg = true;
|
---|
| 545 | state.rstate.stage = 0;
|
---|
| 546 | goto lbl_rcomm;
|
---|
| 547 | lbl_0:
|
---|
| 548 | state.repnfunc = state.repnfunc+1;
|
---|
| 549 | state.repngrad = state.repngrad+1;
|
---|
| 550 |
|
---|
| 551 | //
|
---|
| 552 | // Call LBFGS
|
---|
| 553 | //
|
---|
| 554 | state.internalstate.f = state.f;
|
---|
| 555 | for(i_=0; i_<=n-1;i_++)
|
---|
| 556 | {
|
---|
| 557 | state.internalstate.g[i_] = state.g[i_];
|
---|
| 558 | }
|
---|
| 559 | goto lbl_9;
|
---|
| 560 | lbl_10:
|
---|
| 561 | lbfgs.minlbfgsresults(ref state.internalstate, ref state.x, ref state.internalrep);
|
---|
| 562 | lbl_7:
|
---|
| 563 |
|
---|
| 564 | //
|
---|
| 565 | // Second stage of the hybrid algorithm: LM
|
---|
| 566 | // Initialize quadratic model.
|
---|
| 567 | //
|
---|
| 568 | if( state.usermode!=lmmodefgh )
|
---|
| 569 | {
|
---|
| 570 | goto lbl_11;
|
---|
| 571 | }
|
---|
| 572 |
|
---|
| 573 | //
|
---|
| 574 | // RComm
|
---|
| 575 | //
|
---|
| 576 | lmclearrequestfields(ref state);
|
---|
| 577 | state.needfgh = true;
|
---|
| 578 | state.rstate.stage = 1;
|
---|
| 579 | goto lbl_rcomm;
|
---|
| 580 | lbl_1:
|
---|
| 581 | state.repnfunc = state.repnfunc+1;
|
---|
| 582 | state.repngrad = state.repngrad+1;
|
---|
| 583 | state.repnhess = state.repnhess+1;
|
---|
| 584 |
|
---|
| 585 | //
|
---|
| 586 | // generate raw quadratic model
|
---|
| 587 | //
|
---|
| 588 | for(i=0; i<=n-1; i++)
|
---|
| 589 | {
|
---|
| 590 | for(i_=0; i_<=n-1;i_++)
|
---|
| 591 | {
|
---|
| 592 | state.rawmodel[i,i_] = state.h[i,i_];
|
---|
| 593 | }
|
---|
| 594 | }
|
---|
| 595 | for(i_=0; i_<=n-1;i_++)
|
---|
| 596 | {
|
---|
| 597 | state.gbase[i_] = state.g[i_];
|
---|
| 598 | }
|
---|
| 599 | fbase = state.f;
|
---|
| 600 | lbl_11:
|
---|
| 601 | if( ! (state.usermode==lmmodefgj | state.usermode==lmmodefj) )
|
---|
| 602 | {
|
---|
| 603 | goto lbl_13;
|
---|
| 604 | }
|
---|
| 605 |
|
---|
| 606 | //
|
---|
| 607 | // RComm
|
---|
| 608 | //
|
---|
| 609 | lmclearrequestfields(ref state);
|
---|
| 610 | state.needfij = true;
|
---|
| 611 | state.rstate.stage = 2;
|
---|
| 612 | goto lbl_rcomm;
|
---|
| 613 | lbl_2:
|
---|
| 614 | state.repnfunc = state.repnfunc+1;
|
---|
| 615 | state.repnjac = state.repnjac+1;
|
---|
| 616 |
|
---|
| 617 | //
|
---|
| 618 | // generate raw quadratic model
|
---|
| 619 | //
|
---|
| 620 | blas.matrixmatrixmultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.j, 0, m-1, 0, n-1, false, 1.0, ref state.rawmodel, 0, n-1, 0, n-1, 0.0, ref state.work);
|
---|
| 621 | blas.matrixvectormultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.fi, 0, m-1, 1.0, ref state.gbase, 0, n-1, 0.0);
|
---|
| 622 | fbase = 0.0;
|
---|
| 623 | for(i_=0; i_<=m-1;i_++)
|
---|
| 624 | {
|
---|
| 625 | fbase += state.fi[i_]*state.fi[i_];
|
---|
| 626 | }
|
---|
| 627 | lbl_13:
|
---|
| 628 | lambda = 0.001;
|
---|
| 629 | lbl_15:
|
---|
| 630 | if( false )
|
---|
| 631 | {
|
---|
| 632 | goto lbl_16;
|
---|
| 633 | }
|
---|
| 634 |
|
---|
| 635 | //
|
---|
| 636 | // 1. Model = RawModel+lambda*I
|
---|
| 637 | // 2. Try to solve (RawModel+Lambda*I)*dx = -g.
|
---|
| 638 | // Increase lambda if left part is not positive definite.
|
---|
| 639 | //
|
---|
| 640 | for(i=0; i<=n-1; i++)
|
---|
| 641 | {
|
---|
| 642 | for(i_=0; i_<=n-1;i_++)
|
---|
| 643 | {
|
---|
| 644 | state.model[i,i_] = state.rawmodel[i,i_];
|
---|
| 645 | }
|
---|
| 646 | state.model[i,i] = state.model[i,i]+lambda;
|
---|
| 647 | }
|
---|
| 648 | spd = trfac.spdmatrixcholesky(ref state.model, n, true);
|
---|
| 649 | state.repncholesky = state.repncholesky+1;
|
---|
| 650 | if( !spd )
|
---|
| 651 | {
|
---|
| 652 | lambda = lambda*lambdaup*nu;
|
---|
| 653 | nu = nu*2;
|
---|
| 654 | goto lbl_15;
|
---|
| 655 | }
|
---|
| 656 | densesolver.spdmatrixcholeskysolve(ref state.model, n, true, ref state.gbase, ref state.solverinfo, ref state.solverrep, ref state.xdir);
|
---|
| 657 | if( state.solverinfo<0 )
|
---|
| 658 | {
|
---|
| 659 | lambda = lambda*lambdaup*nu;
|
---|
| 660 | nu = nu*2;
|
---|
| 661 | goto lbl_15;
|
---|
| 662 | }
|
---|
| 663 | for(i_=0; i_<=n-1;i_++)
|
---|
| 664 | {
|
---|
| 665 | state.xdir[i_] = -1*state.xdir[i_];
|
---|
| 666 | }
|
---|
| 667 |
|
---|
| 668 | //
|
---|
| 669 | // Candidate lambda found.
|
---|
| 670 | // 1. Save old w in WBase
|
---|
| 671 | // 1. Test some stopping criterions
|
---|
| 672 | // 2. If error(w+wdir)>error(w), increase lambda
|
---|
| 673 | //
|
---|
| 674 | for(i_=0; i_<=n-1;i_++)
|
---|
| 675 | {
|
---|
| 676 | state.xbase[i_] = state.x[i_];
|
---|
| 677 | }
|
---|
| 678 | for(i_=0; i_<=n-1;i_++)
|
---|
| 679 | {
|
---|
| 680 | state.x[i_] = state.x[i_] + state.xdir[i_];
|
---|
| 681 | }
|
---|
| 682 | xnorm = 0.0;
|
---|
| 683 | for(i_=0; i_<=n-1;i_++)
|
---|
| 684 | {
|
---|
| 685 | xnorm += state.xbase[i_]*state.xbase[i_];
|
---|
| 686 | }
|
---|
| 687 | stepnorm = 0.0;
|
---|
| 688 | for(i_=0; i_<=n-1;i_++)
|
---|
| 689 | {
|
---|
| 690 | stepnorm += state.xdir[i_]*state.xdir[i_];
|
---|
| 691 | }
|
---|
| 692 | xnorm = Math.Sqrt(xnorm);
|
---|
| 693 | stepnorm = Math.Sqrt(stepnorm);
|
---|
| 694 | if( (double)(stepnorm)<=(double)(state.epsx*(1+xnorm)) )
|
---|
| 695 | {
|
---|
| 696 |
|
---|
| 697 | //
|
---|
| 698 | // step size if small enough
|
---|
| 699 | //
|
---|
| 700 | state.repterminationtype = 2;
|
---|
| 701 | goto lbl_16;
|
---|
| 702 | }
|
---|
| 703 | lmclearrequestfields(ref state);
|
---|
| 704 | state.needf = true;
|
---|
| 705 | state.rstate.stage = 3;
|
---|
| 706 | goto lbl_rcomm;
|
---|
| 707 | lbl_3:
|
---|
| 708 | state.repnfunc = state.repnfunc+1;
|
---|
| 709 | fnew = state.f;
|
---|
| 710 | if( (double)(Math.Abs(fnew-fbase))<=(double)(state.epsf*Math.Max(1, Math.Max(Math.Abs(fbase), Math.Abs(fnew)))) )
|
---|
| 711 | {
|
---|
| 712 |
|
---|
| 713 | //
|
---|
| 714 | // function change is small enough
|
---|
| 715 | //
|
---|
| 716 | state.repterminationtype = 1;
|
---|
| 717 | goto lbl_16;
|
---|
| 718 | }
|
---|
| 719 | if( (double)(fnew)>(double)(fbase) )
|
---|
| 720 | {
|
---|
| 721 |
|
---|
| 722 | //
|
---|
| 723 | // restore state and continue out search for lambda
|
---|
| 724 | //
|
---|
| 725 | for(i_=0; i_<=n-1;i_++)
|
---|
| 726 | {
|
---|
| 727 | state.x[i_] = state.xbase[i_];
|
---|
| 728 | }
|
---|
| 729 | lambda = lambda*lambdaup*nu;
|
---|
| 730 | nu = nu*2;
|
---|
| 731 | goto lbl_15;
|
---|
| 732 | }
|
---|
| 733 | if( ! ((state.usermode==lmmodefgj | state.usermode==lmmodefgh) & state.flags/lmflagnointlbfgs%2==0) )
|
---|
| 734 | {
|
---|
| 735 | goto lbl_17;
|
---|
| 736 | }
|
---|
| 737 |
|
---|
| 738 | //
|
---|
| 739 | // Optimize using inv(cholesky(H)) as preconditioner
|
---|
| 740 | //
|
---|
| 741 | if( ! trinverse.rmatrixtrinverse(ref state.model, n, true, false) )
|
---|
| 742 | {
|
---|
| 743 | goto lbl_19;
|
---|
| 744 | }
|
---|
| 745 |
|
---|
| 746 | //
|
---|
| 747 | // if matrix can be inverted use it.
|
---|
| 748 | // just silently move to next iteration otherwise.
|
---|
| 749 | // (will be very, very rare, mostly for specially
|
---|
| 750 | // designed near-degenerate tasks)
|
---|
| 751 | //
|
---|
| 752 | for(i_=0; i_<=n-1;i_++)
|
---|
| 753 | {
|
---|
| 754 | state.xbase[i_] = state.x[i_];
|
---|
| 755 | }
|
---|
| 756 | for(i=0; i<=n-1; i++)
|
---|
| 757 | {
|
---|
| 758 | state.xprec[i] = 0;
|
---|
| 759 | }
|
---|
| 760 | lbfgs.minlbfgs(n, Math.Min(n, lmintlbfgsits), ref state.xprec, 0.0, 0.0, 0.0, lmintlbfgsits, lbfgsflags, ref state.internalstate);
|
---|
| 761 | lbl_21:
|
---|
| 762 | if( ! lbfgs.minlbfgsiteration(ref state.internalstate) )
|
---|
| 763 | {
|
---|
| 764 | goto lbl_22;
|
---|
| 765 | }
|
---|
| 766 |
|
---|
| 767 | //
|
---|
| 768 | // convert XPrec to unpreconditioned form, then call RComm.
|
---|
| 769 | //
|
---|
| 770 | for(i=0; i<=n-1; i++)
|
---|
| 771 | {
|
---|
| 772 | v = 0.0;
|
---|
| 773 | for(i_=i; i_<=n-1;i_++)
|
---|
| 774 | {
|
---|
| 775 | v += state.internalstate.x[i_]*state.model[i,i_];
|
---|
| 776 | }
|
---|
| 777 | state.x[i] = state.xbase[i]+v;
|
---|
| 778 | }
|
---|
| 779 | lmclearrequestfields(ref state);
|
---|
| 780 | state.needfg = true;
|
---|
| 781 | state.rstate.stage = 4;
|
---|
| 782 | goto lbl_rcomm;
|
---|
| 783 | lbl_4:
|
---|
| 784 | state.repnfunc = state.repnfunc+1;
|
---|
| 785 | state.repngrad = state.repngrad+1;
|
---|
| 786 |
|
---|
| 787 | //
|
---|
| 788 | // 1. pass State.F to State.InternalState.F
|
---|
| 789 | // 2. convert gradient back to preconditioned form
|
---|
| 790 | //
|
---|
| 791 | state.internalstate.f = state.f;
|
---|
| 792 | for(i=0; i<=n-1; i++)
|
---|
| 793 | {
|
---|
| 794 | state.internalstate.g[i] = 0;
|
---|
| 795 | }
|
---|
| 796 | for(i=0; i<=n-1; i++)
|
---|
| 797 | {
|
---|
| 798 | v = state.g[i];
|
---|
| 799 | for(i_=i; i_<=n-1;i_++)
|
---|
| 800 | {
|
---|
| 801 | state.internalstate.g[i_] = state.internalstate.g[i_] + v*state.model[i,i_];
|
---|
| 802 | }
|
---|
| 803 | }
|
---|
| 804 |
|
---|
| 805 | //
|
---|
| 806 | // next iteration
|
---|
| 807 | //
|
---|
| 808 | goto lbl_21;
|
---|
| 809 | lbl_22:
|
---|
| 810 |
|
---|
| 811 | //
|
---|
| 812 | // change LBFGS flags to NoRealloc.
|
---|
| 813 | // L-BFGS subroutine will use memory allocated from previous run.
|
---|
| 814 | // it is possible since all subsequent calls will be with same N/M.
|
---|
| 815 | //
|
---|
| 816 | lbfgsflags = lbfgsnorealloc;
|
---|
| 817 |
|
---|
| 818 | //
|
---|
| 819 | // back to unpreconditioned X
|
---|
| 820 | //
|
---|
| 821 | lbfgs.minlbfgsresults(ref state.internalstate, ref state.xprec, ref state.internalrep);
|
---|
| 822 | for(i=0; i<=n-1; i++)
|
---|
| 823 | {
|
---|
| 824 | v = 0.0;
|
---|
| 825 | for(i_=i; i_<=n-1;i_++)
|
---|
| 826 | {
|
---|
| 827 | v += state.xprec[i_]*state.model[i,i_];
|
---|
| 828 | }
|
---|
| 829 | state.x[i] = state.xbase[i]+v;
|
---|
| 830 | }
|
---|
| 831 | lbl_19:
|
---|
| 832 | lbl_17:
|
---|
| 833 |
|
---|
| 834 | //
|
---|
| 835 | // Accept new position.
|
---|
| 836 | // Calculate Hessian
|
---|
| 837 | //
|
---|
| 838 | if( state.usermode!=lmmodefgh )
|
---|
| 839 | {
|
---|
| 840 | goto lbl_23;
|
---|
| 841 | }
|
---|
| 842 |
|
---|
| 843 | //
|
---|
| 844 | // RComm
|
---|
| 845 | //
|
---|
| 846 | lmclearrequestfields(ref state);
|
---|
| 847 | state.needfgh = true;
|
---|
| 848 | state.rstate.stage = 5;
|
---|
| 849 | goto lbl_rcomm;
|
---|
| 850 | lbl_5:
|
---|
| 851 | state.repnfunc = state.repnfunc+1;
|
---|
| 852 | state.repngrad = state.repngrad+1;
|
---|
| 853 | state.repnhess = state.repnhess+1;
|
---|
| 854 |
|
---|
| 855 | //
|
---|
| 856 | // Update raw quadratic model
|
---|
| 857 | //
|
---|
| 858 | for(i=0; i<=n-1; i++)
|
---|
| 859 | {
|
---|
| 860 | for(i_=0; i_<=n-1;i_++)
|
---|
| 861 | {
|
---|
| 862 | state.rawmodel[i,i_] = state.h[i,i_];
|
---|
| 863 | }
|
---|
| 864 | }
|
---|
| 865 | for(i_=0; i_<=n-1;i_++)
|
---|
| 866 | {
|
---|
| 867 | state.gbase[i_] = state.g[i_];
|
---|
| 868 | }
|
---|
| 869 | fbase = state.f;
|
---|
| 870 | lbl_23:
|
---|
| 871 | if( ! (state.usermode==lmmodefgj | state.usermode==lmmodefj) )
|
---|
| 872 | {
|
---|
| 873 | goto lbl_25;
|
---|
| 874 | }
|
---|
| 875 |
|
---|
| 876 | //
|
---|
| 877 | // RComm
|
---|
| 878 | //
|
---|
| 879 | lmclearrequestfields(ref state);
|
---|
| 880 | state.needfij = true;
|
---|
| 881 | state.rstate.stage = 6;
|
---|
| 882 | goto lbl_rcomm;
|
---|
| 883 | lbl_6:
|
---|
| 884 | state.repnfunc = state.repnfunc+1;
|
---|
| 885 | state.repnjac = state.repnjac+1;
|
---|
| 886 |
|
---|
| 887 | //
|
---|
| 888 | // generate raw quadratic model
|
---|
| 889 | //
|
---|
| 890 | blas.matrixmatrixmultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.j, 0, m-1, 0, n-1, false, 1.0, ref state.rawmodel, 0, n-1, 0, n-1, 0.0, ref state.work);
|
---|
| 891 | blas.matrixvectormultiply(ref state.j, 0, m-1, 0, n-1, true, ref state.fi, 0, m-1, 1.0, ref state.gbase, 0, n-1, 0.0);
|
---|
| 892 | fbase = 0.0;
|
---|
| 893 | for(i_=0; i_<=m-1;i_++)
|
---|
| 894 | {
|
---|
| 895 | fbase += state.fi[i_]*state.fi[i_];
|
---|
| 896 | }
|
---|
| 897 | lbl_25:
|
---|
| 898 | state.repiterationscount = state.repiterationscount+1;
|
---|
| 899 | if( state.repiterationscount>=state.maxits & state.maxits>0 )
|
---|
| 900 | {
|
---|
| 901 | state.repterminationtype = 5;
|
---|
| 902 | goto lbl_16;
|
---|
| 903 | }
|
---|
| 904 |
|
---|
| 905 | //
|
---|
| 906 | // Update lambda
|
---|
| 907 | //
|
---|
| 908 | lambda = lambda*lambdadown;
|
---|
| 909 | nu = 2;
|
---|
| 910 | goto lbl_15;
|
---|
| 911 | lbl_16:
|
---|
| 912 | result = false;
|
---|
| 913 | return result;
|
---|
| 914 |
|
---|
| 915 | //
|
---|
| 916 | // Saving state
|
---|
| 917 | //
|
---|
| 918 | lbl_rcomm:
|
---|
| 919 | result = true;
|
---|
| 920 | state.rstate.ia[0] = n;
|
---|
| 921 | state.rstate.ia[1] = m;
|
---|
| 922 | state.rstate.ia[2] = i;
|
---|
| 923 | state.rstate.ia[3] = lbfgsflags;
|
---|
| 924 | state.rstate.ba[0] = spd;
|
---|
| 925 | state.rstate.ra[0] = xnorm;
|
---|
| 926 | state.rstate.ra[1] = stepnorm;
|
---|
| 927 | state.rstate.ra[2] = fbase;
|
---|
| 928 | state.rstate.ra[3] = fnew;
|
---|
| 929 | state.rstate.ra[4] = lambda;
|
---|
| 930 | state.rstate.ra[5] = nu;
|
---|
| 931 | state.rstate.ra[6] = lambdaup;
|
---|
| 932 | state.rstate.ra[7] = lambdadown;
|
---|
| 933 | state.rstate.ra[8] = v;
|
---|
| 934 | return result;
|
---|
| 935 | }
|
---|
| 936 |
|
---|
| 937 |
|
---|
| 938 | /*************************************************************************
|
---|
| 939 | Levenberg-Marquardt algorithm results
|
---|
| 940 |
|
---|
| 941 | Called after MinLMIteration returned False.
|
---|
| 942 |
|
---|
| 943 | Input parameters:
|
---|
| 944 | State - algorithm state (used by MinLMIteration).
|
---|
| 945 |
|
---|
| 946 | Output parameters:
|
---|
| 947 | X - array[0..N-1], solution
|
---|
| 948 | Rep - optimization report:
|
---|
| 949 | * Rep.TerminationType completetion code:
|
---|
| 950 | * -1 incorrect parameters were specified
|
---|
| 951 | * 1 relative function improvement is no more than
|
---|
| 952 | EpsF.
|
---|
| 953 | * 2 relative step is no more than EpsX.
|
---|
| 954 | * 5 MaxIts steps was taken
|
---|
| 955 | * Rep.IterationsCount contains iterations count
|
---|
| 956 | * Rep.NFunc - number of function calculations
|
---|
| 957 | * Rep.NJac - number of Jacobi matrix calculations
|
---|
| 958 | * Rep.NGrad - number of gradient calculations
|
---|
| 959 | * Rep.NHess - number of Hessian calculations
|
---|
| 960 | * Rep.NCholesky - number of Cholesky decomposition calculations
|
---|
| 961 |
|
---|
| 962 | -- ALGLIB --
|
---|
| 963 | Copyright 10.03.2009 by Bochkanov Sergey
|
---|
| 964 | *************************************************************************/
|
---|
| 965 | public static void minlmresults(ref lmstate state,
|
---|
| 966 | ref double[] x,
|
---|
| 967 | ref lmreport rep)
|
---|
| 968 | {
|
---|
| 969 | int i_ = 0;
|
---|
| 970 |
|
---|
| 971 | x = new double[state.n-1+1];
|
---|
| 972 | for(i_=0; i_<=state.n-1;i_++)
|
---|
| 973 | {
|
---|
| 974 | x[i_] = state.x[i_];
|
---|
| 975 | }
|
---|
| 976 | rep.iterationscount = state.repiterationscount;
|
---|
| 977 | rep.terminationtype = state.repterminationtype;
|
---|
| 978 | rep.nfunc = state.repnfunc;
|
---|
| 979 | rep.njac = state.repnjac;
|
---|
| 980 | rep.ngrad = state.repngrad;
|
---|
| 981 | rep.nhess = state.repnhess;
|
---|
| 982 | rep.ncholesky = state.repncholesky;
|
---|
| 983 | }
|
---|
| 984 |
|
---|
| 985 |
|
---|
| 986 | /*************************************************************************
|
---|
| 987 | Prepare internal structures (except for RComm).
|
---|
| 988 |
|
---|
| 989 | Note: M must be zero for FGH mode, non-zero for FJ/FGJ mode.
|
---|
| 990 | *************************************************************************/
|
---|
| 991 | private static void lmprepare(int n,
|
---|
| 992 | int m,
|
---|
| 993 | bool havegrad,
|
---|
| 994 | ref lmstate state)
|
---|
| 995 | {
|
---|
| 996 | state.repiterationscount = 0;
|
---|
| 997 | state.repterminationtype = 0;
|
---|
| 998 | state.repnfunc = 0;
|
---|
| 999 | state.repnjac = 0;
|
---|
| 1000 | state.repngrad = 0;
|
---|
| 1001 | state.repnhess = 0;
|
---|
| 1002 | state.repncholesky = 0;
|
---|
| 1003 | if( n<0 | m<0 )
|
---|
| 1004 | {
|
---|
| 1005 | return;
|
---|
| 1006 | }
|
---|
| 1007 | if( havegrad )
|
---|
| 1008 | {
|
---|
| 1009 | state.g = new double[n-1+1];
|
---|
| 1010 | }
|
---|
| 1011 | if( m!=0 )
|
---|
| 1012 | {
|
---|
| 1013 | state.j = new double[m-1+1, n-1+1];
|
---|
| 1014 | state.fi = new double[m-1+1];
|
---|
| 1015 | state.h = new double[0+1, 0+1];
|
---|
| 1016 | }
|
---|
| 1017 | else
|
---|
| 1018 | {
|
---|
| 1019 | state.j = new double[0+1, 0+1];
|
---|
| 1020 | state.fi = new double[0+1];
|
---|
| 1021 | state.h = new double[n-1+1, n-1+1];
|
---|
| 1022 | }
|
---|
| 1023 | state.x = new double[n-1+1];
|
---|
| 1024 | state.rawmodel = new double[n-1+1, n-1+1];
|
---|
| 1025 | state.model = new double[n-1+1, n-1+1];
|
---|
| 1026 | state.xbase = new double[n-1+1];
|
---|
| 1027 | state.xprec = new double[n-1+1];
|
---|
| 1028 | state.gbase = new double[n-1+1];
|
---|
| 1029 | state.xdir = new double[n-1+1];
|
---|
| 1030 | state.work = new double[Math.Max(n, m)+1];
|
---|
| 1031 | }
|
---|
| 1032 |
|
---|
| 1033 |
|
---|
| 1034 | /*************************************************************************
|
---|
| 1035 | Clears request fileds (to be sure that we don't forgot to clear something)
|
---|
| 1036 | *************************************************************************/
|
---|
| 1037 | private static void lmclearrequestfields(ref lmstate state)
|
---|
| 1038 | {
|
---|
| 1039 | state.needf = false;
|
---|
| 1040 | state.needfg = false;
|
---|
| 1041 | state.needfgh = false;
|
---|
| 1042 | state.needfij = false;
|
---|
| 1043 | }
|
---|
| 1044 | }
|
---|
| 1045 | }
|
---|