[2563] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 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 logit
|
---|
| 26 | {
|
---|
| 27 | public struct logitmodel
|
---|
| 28 | {
|
---|
| 29 | public double[] w;
|
---|
| 30 | };
|
---|
| 31 |
|
---|
| 32 |
|
---|
| 33 | public struct logitmcstate
|
---|
| 34 | {
|
---|
| 35 | public bool brackt;
|
---|
| 36 | public bool stage1;
|
---|
| 37 | public int infoc;
|
---|
| 38 | public double dg;
|
---|
| 39 | public double dgm;
|
---|
| 40 | public double dginit;
|
---|
| 41 | public double dgtest;
|
---|
| 42 | public double dgx;
|
---|
| 43 | public double dgxm;
|
---|
| 44 | public double dgy;
|
---|
| 45 | public double dgym;
|
---|
| 46 | public double finit;
|
---|
| 47 | public double ftest1;
|
---|
| 48 | public double fm;
|
---|
| 49 | public double fx;
|
---|
| 50 | public double fxm;
|
---|
| 51 | public double fy;
|
---|
| 52 | public double fym;
|
---|
| 53 | public double stx;
|
---|
| 54 | public double sty;
|
---|
| 55 | public double stmin;
|
---|
| 56 | public double stmax;
|
---|
| 57 | public double width;
|
---|
| 58 | public double width1;
|
---|
| 59 | public double xtrapf;
|
---|
| 60 | };
|
---|
| 61 |
|
---|
| 62 |
|
---|
| 63 | /*************************************************************************
|
---|
| 64 | MNLReport structure contains information about training process:
|
---|
| 65 | * NGrad - number of gradient calculations
|
---|
| 66 | * NHess - number of Hessian calculations
|
---|
| 67 | *************************************************************************/
|
---|
| 68 | public struct mnlreport
|
---|
| 69 | {
|
---|
| 70 | public int ngrad;
|
---|
| 71 | public int nhess;
|
---|
| 72 | };
|
---|
| 73 |
|
---|
| 74 |
|
---|
| 75 |
|
---|
| 76 |
|
---|
| 77 | public const double xtol = 100*AP.Math.MachineEpsilon;
|
---|
| 78 | public const double ftol = 0.0001;
|
---|
| 79 | public const double gtol = 0.3;
|
---|
| 80 | public const int maxfev = 20;
|
---|
| 81 | public const double stpmin = 1.0E-2;
|
---|
| 82 | public const double stpmax = 1.0E5;
|
---|
| 83 | public const int logitvnum = 6;
|
---|
| 84 |
|
---|
| 85 |
|
---|
| 86 | /*************************************************************************
|
---|
| 87 | This subroutine trains logit model.
|
---|
| 88 |
|
---|
| 89 | INPUT PARAMETERS:
|
---|
| 90 | XY - training set, array[0..NPoints-1,0..NVars]
|
---|
| 91 | First NVars columns store values of independent
|
---|
| 92 | variables, next column stores number of class (from 0
|
---|
| 93 | to NClasses-1) which dataset element belongs to. Fractional
|
---|
| 94 | values are rounded to nearest integer.
|
---|
| 95 | NPoints - training set size, NPoints>=1
|
---|
| 96 | NVars - number of independent variables, NVars>=1
|
---|
| 97 | NClasses - number of classes, NClasses>=2
|
---|
| 98 |
|
---|
| 99 | OUTPUT PARAMETERS:
|
---|
| 100 | Info - return code:
|
---|
| 101 | * -2, if there is a point with class number
|
---|
| 102 | outside of [0..NClasses-1].
|
---|
| 103 | * -1, if incorrect parameters was passed
|
---|
| 104 | (NPoints<NVars+2, NVars<1, NClasses<2).
|
---|
| 105 | * 1, if task has been solved
|
---|
| 106 | LM - model built
|
---|
| 107 | Rep - training report
|
---|
| 108 |
|
---|
| 109 | -- ALGLIB --
|
---|
| 110 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 111 | *************************************************************************/
|
---|
| 112 | public static void mnltrainh(ref double[,] xy,
|
---|
| 113 | int npoints,
|
---|
| 114 | int nvars,
|
---|
| 115 | int nclasses,
|
---|
| 116 | ref int info,
|
---|
| 117 | ref logitmodel lm,
|
---|
| 118 | ref mnlreport rep)
|
---|
| 119 | {
|
---|
| 120 | int i = 0;
|
---|
| 121 | int j = 0;
|
---|
| 122 | int k = 0;
|
---|
| 123 | int m = 0;
|
---|
| 124 | int n = 0;
|
---|
| 125 | int ssize = 0;
|
---|
| 126 | bool allsame = new bool();
|
---|
| 127 | int offs = 0;
|
---|
| 128 | double threshold = 0;
|
---|
| 129 | double wminstep = 0;
|
---|
| 130 | double decay = 0;
|
---|
| 131 | int wdim = 0;
|
---|
| 132 | int expoffs = 0;
|
---|
| 133 | double v = 0;
|
---|
| 134 | double s = 0;
|
---|
| 135 | mlpbase.multilayerperceptron network = new mlpbase.multilayerperceptron();
|
---|
| 136 | int nin = 0;
|
---|
| 137 | int nout = 0;
|
---|
| 138 | int wcount = 0;
|
---|
| 139 | double e = 0;
|
---|
| 140 | double[] g = new double[0];
|
---|
| 141 | double[,] h = new double[0,0];
|
---|
| 142 | bool spd = new bool();
|
---|
| 143 | int cvcnt = 0;
|
---|
| 144 | double[] x = new double[0];
|
---|
| 145 | double[] y = new double[0];
|
---|
| 146 | double[] wbase = new double[0];
|
---|
| 147 | double wstep = 0;
|
---|
| 148 | double[] wdir = new double[0];
|
---|
| 149 | double[] work = new double[0];
|
---|
| 150 | int mcstage = 0;
|
---|
| 151 | logitmcstate mcstate = new logitmcstate();
|
---|
| 152 | int mcinfo = 0;
|
---|
| 153 | int mcnfev = 0;
|
---|
| 154 | int i_ = 0;
|
---|
| 155 | int i1_ = 0;
|
---|
| 156 |
|
---|
| 157 | threshold = 1000*AP.Math.MachineEpsilon;
|
---|
| 158 | wminstep = 0.001;
|
---|
| 159 | decay = 0.001;
|
---|
| 160 |
|
---|
| 161 | //
|
---|
| 162 | // Test for inputs
|
---|
| 163 | //
|
---|
| 164 | if( npoints<nvars+2 | nvars<1 | nclasses<2 )
|
---|
| 165 | {
|
---|
| 166 | info = -1;
|
---|
| 167 | return;
|
---|
| 168 | }
|
---|
| 169 | for(i=0; i<=npoints-1; i++)
|
---|
| 170 | {
|
---|
| 171 | if( (int)Math.Round(xy[i,nvars])<0 | (int)Math.Round(xy[i,nvars])>=nclasses )
|
---|
| 172 | {
|
---|
| 173 | info = -2;
|
---|
| 174 | return;
|
---|
| 175 | }
|
---|
| 176 | }
|
---|
| 177 | info = 1;
|
---|
| 178 |
|
---|
| 179 | //
|
---|
| 180 | // Initialize data
|
---|
| 181 | //
|
---|
| 182 | rep.ngrad = 0;
|
---|
| 183 | rep.nhess = 0;
|
---|
| 184 |
|
---|
| 185 | //
|
---|
| 186 | // Allocate array
|
---|
| 187 | //
|
---|
| 188 | wdim = (nvars+1)*(nclasses-1);
|
---|
| 189 | offs = 5;
|
---|
| 190 | expoffs = offs+wdim;
|
---|
| 191 | ssize = 5+(nvars+1)*(nclasses-1)+nclasses;
|
---|
| 192 | lm.w = new double[ssize-1+1];
|
---|
| 193 | lm.w[0] = ssize;
|
---|
| 194 | lm.w[1] = logitvnum;
|
---|
| 195 | lm.w[2] = nvars;
|
---|
| 196 | lm.w[3] = nclasses;
|
---|
| 197 | lm.w[4] = offs;
|
---|
| 198 |
|
---|
| 199 | //
|
---|
| 200 | // Degenerate case: all outputs are equal
|
---|
| 201 | //
|
---|
| 202 | allsame = true;
|
---|
| 203 | for(i=1; i<=npoints-1; i++)
|
---|
| 204 | {
|
---|
| 205 | if( (int)Math.Round(xy[i,nvars])!=(int)Math.Round(xy[i-1,nvars]) )
|
---|
| 206 | {
|
---|
| 207 | allsame = false;
|
---|
| 208 | }
|
---|
| 209 | }
|
---|
| 210 | if( allsame )
|
---|
| 211 | {
|
---|
| 212 | for(i=0; i<=(nvars+1)*(nclasses-1)-1; i++)
|
---|
| 213 | {
|
---|
| 214 | lm.w[offs+i] = 0;
|
---|
| 215 | }
|
---|
| 216 | v = -(2*Math.Log(AP.Math.MinRealNumber));
|
---|
| 217 | k = (int)Math.Round(xy[0,nvars]);
|
---|
| 218 | if( k==nclasses-1 )
|
---|
| 219 | {
|
---|
| 220 | for(i=0; i<=nclasses-2; i++)
|
---|
| 221 | {
|
---|
| 222 | lm.w[offs+i*(nvars+1)+nvars] = -v;
|
---|
| 223 | }
|
---|
| 224 | }
|
---|
| 225 | else
|
---|
| 226 | {
|
---|
| 227 | for(i=0; i<=nclasses-2; i++)
|
---|
| 228 | {
|
---|
| 229 | if( i==k )
|
---|
| 230 | {
|
---|
| 231 | lm.w[offs+i*(nvars+1)+nvars] = +v;
|
---|
| 232 | }
|
---|
| 233 | else
|
---|
| 234 | {
|
---|
| 235 | lm.w[offs+i*(nvars+1)+nvars] = 0;
|
---|
| 236 | }
|
---|
| 237 | }
|
---|
| 238 | }
|
---|
| 239 | return;
|
---|
| 240 | }
|
---|
| 241 |
|
---|
| 242 | //
|
---|
| 243 | // General case.
|
---|
| 244 | // Prepare task and network. Allocate space.
|
---|
| 245 | //
|
---|
| 246 | mlpbase.mlpcreatec0(nvars, nclasses, ref network);
|
---|
| 247 | mlpbase.mlpinitpreprocessor(ref network, ref xy, npoints);
|
---|
| 248 | mlpbase.mlpproperties(ref network, ref nin, ref nout, ref wcount);
|
---|
| 249 | for(i=0; i<=wcount-1; i++)
|
---|
| 250 | {
|
---|
| 251 | network.weights[i] = (2*AP.Math.RandomReal()-1)/nvars;
|
---|
| 252 | }
|
---|
| 253 | g = new double[wcount-1+1];
|
---|
| 254 | h = new double[wcount-1+1, wcount-1+1];
|
---|
| 255 | wbase = new double[wcount-1+1];
|
---|
| 256 | wdir = new double[wcount-1+1];
|
---|
| 257 | work = new double[wcount-1+1];
|
---|
| 258 |
|
---|
| 259 | //
|
---|
| 260 | // First stage: optimize in gradient direction.
|
---|
| 261 | //
|
---|
| 262 | for(k=0; k<=wcount/3+10; k++)
|
---|
| 263 | {
|
---|
| 264 |
|
---|
| 265 | //
|
---|
| 266 | // Calculate gradient in starting point
|
---|
| 267 | //
|
---|
| 268 | mlpbase.mlpgradnbatch(ref network, ref xy, npoints, ref e, ref g);
|
---|
| 269 | v = 0.0;
|
---|
| 270 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 271 | {
|
---|
| 272 | v += network.weights[i_]*network.weights[i_];
|
---|
| 273 | }
|
---|
| 274 | e = e+0.5*decay*v;
|
---|
| 275 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 276 | {
|
---|
| 277 | g[i_] = g[i_] + decay*network.weights[i_];
|
---|
| 278 | }
|
---|
| 279 | rep.ngrad = rep.ngrad+1;
|
---|
| 280 |
|
---|
| 281 | //
|
---|
| 282 | // Setup optimization scheme
|
---|
| 283 | //
|
---|
| 284 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 285 | {
|
---|
| 286 | wdir[i_] = -g[i_];
|
---|
| 287 | }
|
---|
| 288 | v = 0.0;
|
---|
| 289 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 290 | {
|
---|
| 291 | v += wdir[i_]*wdir[i_];
|
---|
| 292 | }
|
---|
| 293 | wstep = Math.Sqrt(v);
|
---|
| 294 | v = 1/Math.Sqrt(v);
|
---|
| 295 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 296 | {
|
---|
| 297 | wdir[i_] = v*wdir[i_];
|
---|
| 298 | }
|
---|
| 299 | mcstage = 0;
|
---|
| 300 | mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
|
---|
| 301 | while( mcstage!=0 )
|
---|
| 302 | {
|
---|
| 303 | mlpbase.mlpgradnbatch(ref network, ref xy, npoints, ref e, ref g);
|
---|
| 304 | v = 0.0;
|
---|
| 305 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 306 | {
|
---|
| 307 | v += network.weights[i_]*network.weights[i_];
|
---|
| 308 | }
|
---|
| 309 | e = e+0.5*decay*v;
|
---|
| 310 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 311 | {
|
---|
| 312 | g[i_] = g[i_] + decay*network.weights[i_];
|
---|
| 313 | }
|
---|
| 314 | rep.ngrad = rep.ngrad+1;
|
---|
| 315 | mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
|
---|
| 316 | }
|
---|
| 317 | }
|
---|
| 318 |
|
---|
| 319 | //
|
---|
| 320 | // Second stage: use Hessian when we are close to the minimum
|
---|
| 321 | //
|
---|
| 322 | while( true )
|
---|
| 323 | {
|
---|
| 324 |
|
---|
| 325 | //
|
---|
| 326 | // Calculate and update E/G/H
|
---|
| 327 | //
|
---|
| 328 | mlpbase.mlphessiannbatch(ref network, ref xy, npoints, ref e, ref g, ref h);
|
---|
| 329 | v = 0.0;
|
---|
| 330 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 331 | {
|
---|
| 332 | v += network.weights[i_]*network.weights[i_];
|
---|
| 333 | }
|
---|
| 334 | e = e+0.5*decay*v;
|
---|
| 335 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 336 | {
|
---|
| 337 | g[i_] = g[i_] + decay*network.weights[i_];
|
---|
| 338 | }
|
---|
| 339 | for(k=0; k<=wcount-1; k++)
|
---|
| 340 | {
|
---|
| 341 | h[k,k] = h[k,k]+decay;
|
---|
| 342 | }
|
---|
| 343 | rep.nhess = rep.nhess+1;
|
---|
| 344 |
|
---|
| 345 | //
|
---|
| 346 | // Select step direction
|
---|
| 347 | // NOTE: it is important to use lower-triangle Cholesky
|
---|
| 348 | // factorization since it is much faster than higher-triangle version.
|
---|
| 349 | //
|
---|
| 350 | spd = cholesky.spdmatrixcholesky(ref h, wcount, false);
|
---|
| 351 | if( spd )
|
---|
| 352 | {
|
---|
| 353 | spd = spdsolve.spdmatrixcholeskysolve(ref h, g, wcount, false, ref wdir);
|
---|
| 354 | }
|
---|
| 355 | if( spd )
|
---|
| 356 | {
|
---|
| 357 |
|
---|
| 358 | //
|
---|
| 359 | // H is positive definite.
|
---|
| 360 | // Step in Newton direction.
|
---|
| 361 | //
|
---|
| 362 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 363 | {
|
---|
| 364 | wdir[i_] = -1*wdir[i_];
|
---|
| 365 | }
|
---|
| 366 | spd = true;
|
---|
| 367 | }
|
---|
| 368 | else
|
---|
| 369 | {
|
---|
| 370 |
|
---|
| 371 | //
|
---|
| 372 | // H is indefinite.
|
---|
| 373 | // Step in gradient direction.
|
---|
| 374 | //
|
---|
| 375 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 376 | {
|
---|
| 377 | wdir[i_] = -g[i_];
|
---|
| 378 | }
|
---|
| 379 | spd = false;
|
---|
| 380 | }
|
---|
| 381 |
|
---|
| 382 | //
|
---|
| 383 | // Optimize in WDir direction
|
---|
| 384 | //
|
---|
| 385 | v = 0.0;
|
---|
| 386 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 387 | {
|
---|
| 388 | v += wdir[i_]*wdir[i_];
|
---|
| 389 | }
|
---|
| 390 | wstep = Math.Sqrt(v);
|
---|
| 391 | v = 1/Math.Sqrt(v);
|
---|
| 392 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 393 | {
|
---|
| 394 | wdir[i_] = v*wdir[i_];
|
---|
| 395 | }
|
---|
| 396 | mcstage = 0;
|
---|
| 397 | mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
|
---|
| 398 | while( mcstage!=0 )
|
---|
| 399 | {
|
---|
| 400 | mlpbase.mlpgradnbatch(ref network, ref xy, npoints, ref e, ref g);
|
---|
| 401 | v = 0.0;
|
---|
| 402 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 403 | {
|
---|
| 404 | v += network.weights[i_]*network.weights[i_];
|
---|
| 405 | }
|
---|
| 406 | e = e+0.5*decay*v;
|
---|
| 407 | for(i_=0; i_<=wcount-1;i_++)
|
---|
| 408 | {
|
---|
| 409 | g[i_] = g[i_] + decay*network.weights[i_];
|
---|
| 410 | }
|
---|
| 411 | rep.ngrad = rep.ngrad+1;
|
---|
| 412 | mnlmcsrch(wcount, ref network.weights, ref e, ref g, ref wdir, ref wstep, ref mcinfo, ref mcnfev, ref work, ref mcstate, ref mcstage);
|
---|
| 413 | }
|
---|
| 414 | if( spd & (mcinfo==2 | mcinfo==4 | mcinfo==6) )
|
---|
| 415 | {
|
---|
| 416 | break;
|
---|
| 417 | }
|
---|
| 418 | }
|
---|
| 419 |
|
---|
| 420 | //
|
---|
| 421 | // Convert from NN format to MNL format
|
---|
| 422 | //
|
---|
| 423 | i1_ = (0) - (offs);
|
---|
| 424 | for(i_=offs; i_<=offs+wcount-1;i_++)
|
---|
| 425 | {
|
---|
| 426 | lm.w[i_] = network.weights[i_+i1_];
|
---|
| 427 | }
|
---|
| 428 | for(k=0; k<=nvars-1; k++)
|
---|
| 429 | {
|
---|
| 430 | for(i=0; i<=nclasses-2; i++)
|
---|
| 431 | {
|
---|
| 432 | s = network.columnsigmas[k];
|
---|
| 433 | if( (double)(s)==(double)(0) )
|
---|
| 434 | {
|
---|
| 435 | s = 1;
|
---|
| 436 | }
|
---|
| 437 | j = offs+(nvars+1)*i;
|
---|
| 438 | v = lm.w[j+k];
|
---|
| 439 | lm.w[j+k] = v/s;
|
---|
| 440 | lm.w[j+nvars] = lm.w[j+nvars]+v*network.columnmeans[k]/s;
|
---|
| 441 | }
|
---|
| 442 | }
|
---|
| 443 | for(k=0; k<=nclasses-2; k++)
|
---|
| 444 | {
|
---|
| 445 | lm.w[offs+(nvars+1)*k+nvars] = -lm.w[offs+(nvars+1)*k+nvars];
|
---|
| 446 | }
|
---|
| 447 | }
|
---|
| 448 |
|
---|
| 449 |
|
---|
| 450 | /*************************************************************************
|
---|
| 451 | Procesing
|
---|
| 452 |
|
---|
| 453 | INPUT PARAMETERS:
|
---|
| 454 | LM - logit model, passed by non-constant reference
|
---|
| 455 | (some fields of structure are used as temporaries
|
---|
| 456 | when calculating model output).
|
---|
| 457 | X - input vector, array[0..NVars-1].
|
---|
| 458 |
|
---|
| 459 | OUTPUT PARAMETERS:
|
---|
| 460 | Y - result, array[0..NClasses-1]
|
---|
| 461 | Vector of posterior probabilities for classification task.
|
---|
| 462 | Subroutine does not allocate memory for this vector, it is
|
---|
| 463 | responsibility of a caller to allocate it. Array must be
|
---|
| 464 | at least [0..NClasses-1].
|
---|
| 465 |
|
---|
| 466 | -- ALGLIB --
|
---|
| 467 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 468 | *************************************************************************/
|
---|
| 469 | public static void mnlprocess(ref logitmodel lm,
|
---|
| 470 | ref double[] x,
|
---|
| 471 | ref double[] y)
|
---|
| 472 | {
|
---|
| 473 | int nvars = 0;
|
---|
| 474 | int nclasses = 0;
|
---|
| 475 | int offs = 0;
|
---|
| 476 | int i = 0;
|
---|
| 477 | int i1 = 0;
|
---|
| 478 | double v = 0;
|
---|
| 479 | double mx = 0;
|
---|
| 480 | double s = 0;
|
---|
| 481 |
|
---|
| 482 | System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLProcess: unexpected model version");
|
---|
| 483 | nvars = (int)Math.Round(lm.w[2]);
|
---|
| 484 | nclasses = (int)Math.Round(lm.w[3]);
|
---|
| 485 | offs = (int)Math.Round(lm.w[4]);
|
---|
| 486 | mnliexp(ref lm.w, ref x);
|
---|
| 487 | s = 0;
|
---|
| 488 | i1 = offs+(nvars+1)*(nclasses-1);
|
---|
| 489 | for(i=i1; i<=i1+nclasses-1; i++)
|
---|
| 490 | {
|
---|
| 491 | s = s+lm.w[i];
|
---|
| 492 | }
|
---|
| 493 | for(i=0; i<=nclasses-1; i++)
|
---|
| 494 | {
|
---|
| 495 | y[i] = lm.w[i1+i]/s;
|
---|
| 496 | }
|
---|
| 497 | }
|
---|
| 498 |
|
---|
| 499 |
|
---|
| 500 | /*************************************************************************
|
---|
| 501 | Unpacks coefficients of logit model. Logit model have form:
|
---|
| 502 |
|
---|
| 503 | P(class=i) = S(i) / (S(0) + S(1) + ... +S(M-1))
|
---|
| 504 | S(i) = Exp(A[i,0]*X[0] + ... + A[i,N-1]*X[N-1] + A[i,N]), when i<M-1
|
---|
| 505 | S(M-1) = 1
|
---|
| 506 |
|
---|
| 507 | INPUT PARAMETERS:
|
---|
| 508 | LM - logit model in ALGLIB format
|
---|
| 509 |
|
---|
| 510 | OUTPUT PARAMETERS:
|
---|
| 511 | V - coefficients, array[0..NClasses-2,0..NVars]
|
---|
| 512 | NVars - number of independent variables
|
---|
| 513 | NClasses - number of classes
|
---|
| 514 |
|
---|
| 515 | -- ALGLIB --
|
---|
| 516 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 517 | *************************************************************************/
|
---|
| 518 | public static void mnlunpack(ref logitmodel lm,
|
---|
| 519 | ref double[,] a,
|
---|
| 520 | ref int nvars,
|
---|
| 521 | ref int nclasses)
|
---|
| 522 | {
|
---|
| 523 | int offs = 0;
|
---|
| 524 | int i = 0;
|
---|
| 525 | int i_ = 0;
|
---|
| 526 | int i1_ = 0;
|
---|
| 527 |
|
---|
| 528 | System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLUnpack: unexpected model version");
|
---|
| 529 | nvars = (int)Math.Round(lm.w[2]);
|
---|
| 530 | nclasses = (int)Math.Round(lm.w[3]);
|
---|
| 531 | offs = (int)Math.Round(lm.w[4]);
|
---|
| 532 | a = new double[nclasses-2+1, nvars+1];
|
---|
| 533 | for(i=0; i<=nclasses-2; i++)
|
---|
| 534 | {
|
---|
| 535 | i1_ = (offs+i*(nvars+1)) - (0);
|
---|
| 536 | for(i_=0; i_<=nvars;i_++)
|
---|
| 537 | {
|
---|
| 538 | a[i,i_] = lm.w[i_+i1_];
|
---|
| 539 | }
|
---|
| 540 | }
|
---|
| 541 | }
|
---|
| 542 |
|
---|
| 543 |
|
---|
| 544 | /*************************************************************************
|
---|
| 545 | "Packs" coefficients and creates logit model in ALGLIB format (MNLUnpack
|
---|
| 546 | reversed).
|
---|
| 547 |
|
---|
| 548 | INPUT PARAMETERS:
|
---|
| 549 | A - model (see MNLUnpack)
|
---|
| 550 | NVars - number of independent variables
|
---|
| 551 | NClasses - number of classes
|
---|
| 552 |
|
---|
| 553 | OUTPUT PARAMETERS:
|
---|
| 554 | LM - logit model.
|
---|
| 555 |
|
---|
| 556 | -- ALGLIB --
|
---|
| 557 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 558 | *************************************************************************/
|
---|
| 559 | public static void mnlpack(ref double[,] a,
|
---|
| 560 | int nvars,
|
---|
| 561 | int nclasses,
|
---|
| 562 | ref logitmodel lm)
|
---|
| 563 | {
|
---|
| 564 | int offs = 0;
|
---|
| 565 | int i = 0;
|
---|
| 566 | int wdim = 0;
|
---|
| 567 | int ssize = 0;
|
---|
| 568 | int i_ = 0;
|
---|
| 569 | int i1_ = 0;
|
---|
| 570 |
|
---|
| 571 | wdim = (nvars+1)*(nclasses-1);
|
---|
| 572 | offs = 5;
|
---|
| 573 | ssize = 5+(nvars+1)*(nclasses-1)+nclasses;
|
---|
| 574 | lm.w = new double[ssize-1+1];
|
---|
| 575 | lm.w[0] = ssize;
|
---|
| 576 | lm.w[1] = logitvnum;
|
---|
| 577 | lm.w[2] = nvars;
|
---|
| 578 | lm.w[3] = nclasses;
|
---|
| 579 | lm.w[4] = offs;
|
---|
| 580 | for(i=0; i<=nclasses-2; i++)
|
---|
| 581 | {
|
---|
| 582 | i1_ = (0) - (offs+i*(nvars+1));
|
---|
| 583 | for(i_=offs+i*(nvars+1); i_<=offs+i*(nvars+1)+nvars;i_++)
|
---|
| 584 | {
|
---|
| 585 | lm.w[i_] = a[i,i_+i1_];
|
---|
| 586 | }
|
---|
| 587 | }
|
---|
| 588 | }
|
---|
| 589 |
|
---|
| 590 |
|
---|
| 591 | /*************************************************************************
|
---|
| 592 | Copying of LogitModel strucure
|
---|
| 593 |
|
---|
| 594 | INPUT PARAMETERS:
|
---|
| 595 | LM1 - original
|
---|
| 596 |
|
---|
| 597 | OUTPUT PARAMETERS:
|
---|
| 598 | LM2 - copy
|
---|
| 599 |
|
---|
| 600 | -- ALGLIB --
|
---|
| 601 | Copyright 15.03.2009 by Bochkanov Sergey
|
---|
| 602 | *************************************************************************/
|
---|
| 603 | public static void mnlcopy(ref logitmodel lm1,
|
---|
| 604 | ref logitmodel lm2)
|
---|
| 605 | {
|
---|
| 606 | int k = 0;
|
---|
| 607 | int i_ = 0;
|
---|
| 608 |
|
---|
| 609 | k = (int)Math.Round(lm1.w[0]);
|
---|
| 610 | lm2.w = new double[k-1+1];
|
---|
| 611 | for(i_=0; i_<=k-1;i_++)
|
---|
| 612 | {
|
---|
| 613 | lm2.w[i_] = lm1.w[i_];
|
---|
| 614 | }
|
---|
| 615 | }
|
---|
| 616 |
|
---|
| 617 |
|
---|
| 618 | /*************************************************************************
|
---|
| 619 | Serialization of LogitModel strucure
|
---|
| 620 |
|
---|
| 621 | INPUT PARAMETERS:
|
---|
| 622 | LM - original
|
---|
| 623 |
|
---|
| 624 | OUTPUT PARAMETERS:
|
---|
| 625 | RA - array of real numbers which stores model,
|
---|
| 626 | array[0..RLen-1]
|
---|
| 627 | RLen - RA lenght
|
---|
| 628 |
|
---|
| 629 | -- ALGLIB --
|
---|
| 630 | Copyright 15.03.2009 by Bochkanov Sergey
|
---|
| 631 | *************************************************************************/
|
---|
| 632 | public static void mnlserialize(ref logitmodel lm,
|
---|
| 633 | ref double[] ra,
|
---|
| 634 | ref int rlen)
|
---|
| 635 | {
|
---|
| 636 | int i_ = 0;
|
---|
| 637 | int i1_ = 0;
|
---|
| 638 |
|
---|
| 639 | rlen = (int)Math.Round(lm.w[0])+1;
|
---|
| 640 | ra = new double[rlen-1+1];
|
---|
| 641 | ra[0] = logitvnum;
|
---|
| 642 | i1_ = (0) - (1);
|
---|
| 643 | for(i_=1; i_<=rlen-1;i_++)
|
---|
| 644 | {
|
---|
| 645 | ra[i_] = lm.w[i_+i1_];
|
---|
| 646 | }
|
---|
| 647 | }
|
---|
| 648 |
|
---|
| 649 |
|
---|
| 650 | /*************************************************************************
|
---|
| 651 | Unserialization of LogitModel strucure
|
---|
| 652 |
|
---|
| 653 | INPUT PARAMETERS:
|
---|
| 654 | RA - real array which stores model
|
---|
| 655 |
|
---|
| 656 | OUTPUT PARAMETERS:
|
---|
| 657 | LM - restored model
|
---|
| 658 |
|
---|
| 659 | -- ALGLIB --
|
---|
| 660 | Copyright 15.03.2009 by Bochkanov Sergey
|
---|
| 661 | *************************************************************************/
|
---|
| 662 | public static void mnlunserialize(ref double[] ra,
|
---|
| 663 | ref logitmodel lm)
|
---|
| 664 | {
|
---|
| 665 | int i_ = 0;
|
---|
| 666 | int i1_ = 0;
|
---|
| 667 |
|
---|
| 668 | System.Diagnostics.Debug.Assert((int)Math.Round(ra[0])==logitvnum, "MNLUnserialize: incorrect array!");
|
---|
| 669 | lm.w = new double[(int)Math.Round(ra[1])-1+1];
|
---|
| 670 | i1_ = (1) - (0);
|
---|
| 671 | for(i_=0; i_<=(int)Math.Round(ra[1])-1;i_++)
|
---|
| 672 | {
|
---|
| 673 | lm.w[i_] = ra[i_+i1_];
|
---|
| 674 | }
|
---|
| 675 | }
|
---|
| 676 |
|
---|
| 677 |
|
---|
| 678 | /*************************************************************************
|
---|
| 679 | Average cross-entropy (in bits per element) on the test set
|
---|
| 680 |
|
---|
| 681 | INPUT PARAMETERS:
|
---|
| 682 | LM - logit model
|
---|
| 683 | XY - test set
|
---|
| 684 | NPoints - test set size
|
---|
| 685 |
|
---|
| 686 | RESULT:
|
---|
| 687 | CrossEntropy/(NPoints*ln(2)).
|
---|
| 688 |
|
---|
| 689 | -- ALGLIB --
|
---|
| 690 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 691 | *************************************************************************/
|
---|
| 692 | public static double mnlavgce(ref logitmodel lm,
|
---|
| 693 | ref double[,] xy,
|
---|
| 694 | int npoints)
|
---|
| 695 | {
|
---|
| 696 | double result = 0;
|
---|
| 697 | int nvars = 0;
|
---|
| 698 | int nclasses = 0;
|
---|
| 699 | int i = 0;
|
---|
| 700 | int j = 0;
|
---|
| 701 | double[] workx = new double[0];
|
---|
| 702 | double[] worky = new double[0];
|
---|
| 703 | int nmax = 0;
|
---|
| 704 | int i_ = 0;
|
---|
| 705 |
|
---|
| 706 | System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLClsError: unexpected model version");
|
---|
| 707 | nvars = (int)Math.Round(lm.w[2]);
|
---|
| 708 | nclasses = (int)Math.Round(lm.w[3]);
|
---|
| 709 | workx = new double[nvars-1+1];
|
---|
| 710 | worky = new double[nclasses-1+1];
|
---|
| 711 | result = 0;
|
---|
| 712 | for(i=0; i<=npoints-1; i++)
|
---|
| 713 | {
|
---|
| 714 | System.Diagnostics.Debug.Assert((int)Math.Round(xy[i,nvars])>=0 & (int)Math.Round(xy[i,nvars])<nclasses, "MNLAvgCE: incorrect class number!");
|
---|
| 715 |
|
---|
| 716 | //
|
---|
| 717 | // Process
|
---|
| 718 | //
|
---|
| 719 | for(i_=0; i_<=nvars-1;i_++)
|
---|
| 720 | {
|
---|
| 721 | workx[i_] = xy[i,i_];
|
---|
| 722 | }
|
---|
| 723 | mnlprocess(ref lm, ref workx, ref worky);
|
---|
| 724 | if( (double)(worky[(int)Math.Round(xy[i,nvars])])>(double)(0) )
|
---|
| 725 | {
|
---|
| 726 | result = result-Math.Log(worky[(int)Math.Round(xy[i,nvars])]);
|
---|
| 727 | }
|
---|
| 728 | else
|
---|
| 729 | {
|
---|
| 730 | result = result-Math.Log(AP.Math.MinRealNumber);
|
---|
| 731 | }
|
---|
| 732 | }
|
---|
| 733 | result = result/(npoints*Math.Log(2));
|
---|
| 734 | return result;
|
---|
| 735 | }
|
---|
| 736 |
|
---|
| 737 |
|
---|
| 738 | /*************************************************************************
|
---|
| 739 | Relative classification error on the test set
|
---|
| 740 |
|
---|
| 741 | INPUT PARAMETERS:
|
---|
| 742 | LM - logit model
|
---|
| 743 | XY - test set
|
---|
| 744 | NPoints - test set size
|
---|
| 745 |
|
---|
| 746 | RESULT:
|
---|
| 747 | percent of incorrectly classified cases.
|
---|
| 748 |
|
---|
| 749 | -- ALGLIB --
|
---|
| 750 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 751 | *************************************************************************/
|
---|
| 752 | public static double mnlrelclserror(ref logitmodel lm,
|
---|
| 753 | ref double[,] xy,
|
---|
| 754 | int npoints)
|
---|
| 755 | {
|
---|
| 756 | double result = 0;
|
---|
| 757 |
|
---|
| 758 | result = (double)(mnlclserror(ref lm, ref xy, npoints))/(double)(npoints);
|
---|
| 759 | return result;
|
---|
| 760 | }
|
---|
| 761 |
|
---|
| 762 |
|
---|
| 763 | /*************************************************************************
|
---|
| 764 | RMS error on the test set
|
---|
| 765 |
|
---|
| 766 | INPUT PARAMETERS:
|
---|
| 767 | LM - logit model
|
---|
| 768 | XY - test set
|
---|
| 769 | NPoints - test set size
|
---|
| 770 |
|
---|
| 771 | RESULT:
|
---|
| 772 | root mean square error (error when estimating posterior probabilities).
|
---|
| 773 |
|
---|
| 774 | -- ALGLIB --
|
---|
| 775 | Copyright 30.08.2008 by Bochkanov Sergey
|
---|
| 776 | *************************************************************************/
|
---|
| 777 | public static double mnlrmserror(ref logitmodel lm,
|
---|
| 778 | ref double[,] xy,
|
---|
| 779 | int npoints)
|
---|
| 780 | {
|
---|
| 781 | double result = 0;
|
---|
| 782 | double relcls = 0;
|
---|
| 783 | double avgce = 0;
|
---|
| 784 | double rms = 0;
|
---|
| 785 | double avg = 0;
|
---|
| 786 | double avgrel = 0;
|
---|
| 787 |
|
---|
| 788 | System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNLRMSError: Incorrect MNL version!");
|
---|
| 789 | mnlallerrors(ref lm, ref xy, npoints, ref relcls, ref avgce, ref rms, ref avg, ref avgrel);
|
---|
| 790 | result = rms;
|
---|
| 791 | return result;
|
---|
| 792 | }
|
---|
| 793 |
|
---|
| 794 |
|
---|
| 795 | /*************************************************************************
|
---|
| 796 | Average error on the test set
|
---|
| 797 |
|
---|
| 798 | INPUT PARAMETERS:
|
---|
| 799 | LM - logit model
|
---|
| 800 | XY - test set
|
---|
| 801 | NPoints - test set size
|
---|
| 802 |
|
---|
| 803 | RESULT:
|
---|
| 804 | average error (error when estimating posterior probabilities).
|
---|
| 805 |
|
---|
| 806 | -- ALGLIB --
|
---|
| 807 | Copyright 30.08.2008 by Bochkanov Sergey
|
---|
| 808 | *************************************************************************/
|
---|
| 809 | public static double mnlavgerror(ref logitmodel lm,
|
---|
| 810 | ref double[,] xy,
|
---|
| 811 | int npoints)
|
---|
| 812 | {
|
---|
| 813 | double result = 0;
|
---|
| 814 | double relcls = 0;
|
---|
| 815 | double avgce = 0;
|
---|
| 816 | double rms = 0;
|
---|
| 817 | double avg = 0;
|
---|
| 818 | double avgrel = 0;
|
---|
| 819 |
|
---|
| 820 | System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNLRMSError: Incorrect MNL version!");
|
---|
| 821 | mnlallerrors(ref lm, ref xy, npoints, ref relcls, ref avgce, ref rms, ref avg, ref avgrel);
|
---|
| 822 | result = avg;
|
---|
| 823 | return result;
|
---|
| 824 | }
|
---|
| 825 |
|
---|
| 826 |
|
---|
| 827 | /*************************************************************************
|
---|
| 828 | Average relative error on the test set
|
---|
| 829 |
|
---|
| 830 | INPUT PARAMETERS:
|
---|
| 831 | LM - logit model
|
---|
| 832 | XY - test set
|
---|
| 833 | NPoints - test set size
|
---|
| 834 |
|
---|
| 835 | RESULT:
|
---|
| 836 | average relative error (error when estimating posterior probabilities).
|
---|
| 837 |
|
---|
| 838 | -- ALGLIB --
|
---|
| 839 | Copyright 30.08.2008 by Bochkanov Sergey
|
---|
| 840 | *************************************************************************/
|
---|
| 841 | public static double mnlavgrelerror(ref logitmodel lm,
|
---|
| 842 | ref double[,] xy,
|
---|
| 843 | int ssize)
|
---|
| 844 | {
|
---|
| 845 | double result = 0;
|
---|
| 846 | double relcls = 0;
|
---|
| 847 | double avgce = 0;
|
---|
| 848 | double rms = 0;
|
---|
| 849 | double avg = 0;
|
---|
| 850 | double avgrel = 0;
|
---|
| 851 |
|
---|
| 852 | System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNLRMSError: Incorrect MNL version!");
|
---|
| 853 | mnlallerrors(ref lm, ref xy, ssize, ref relcls, ref avgce, ref rms, ref avg, ref avgrel);
|
---|
| 854 | result = avgrel;
|
---|
| 855 | return result;
|
---|
| 856 | }
|
---|
| 857 |
|
---|
| 858 |
|
---|
| 859 | /*************************************************************************
|
---|
| 860 | Classification error on test set = MNLRelClsError*NPoints
|
---|
| 861 |
|
---|
| 862 | -- ALGLIB --
|
---|
| 863 | Copyright 10.09.2008 by Bochkanov Sergey
|
---|
| 864 | *************************************************************************/
|
---|
| 865 | public static int mnlclserror(ref logitmodel lm,
|
---|
| 866 | ref double[,] xy,
|
---|
| 867 | int npoints)
|
---|
| 868 | {
|
---|
| 869 | int result = 0;
|
---|
| 870 | int nvars = 0;
|
---|
| 871 | int nclasses = 0;
|
---|
| 872 | int i = 0;
|
---|
| 873 | int j = 0;
|
---|
| 874 | double[] workx = new double[0];
|
---|
| 875 | double[] worky = new double[0];
|
---|
| 876 | int nmax = 0;
|
---|
| 877 | int i_ = 0;
|
---|
| 878 |
|
---|
| 879 | System.Diagnostics.Debug.Assert((double)(lm.w[1])==(double)(logitvnum), "MNLClsError: unexpected model version");
|
---|
| 880 | nvars = (int)Math.Round(lm.w[2]);
|
---|
| 881 | nclasses = (int)Math.Round(lm.w[3]);
|
---|
| 882 | workx = new double[nvars-1+1];
|
---|
| 883 | worky = new double[nclasses-1+1];
|
---|
| 884 | result = 0;
|
---|
| 885 | for(i=0; i<=npoints-1; i++)
|
---|
| 886 | {
|
---|
| 887 |
|
---|
| 888 | //
|
---|
| 889 | // Process
|
---|
| 890 | //
|
---|
| 891 | for(i_=0; i_<=nvars-1;i_++)
|
---|
| 892 | {
|
---|
| 893 | workx[i_] = xy[i,i_];
|
---|
| 894 | }
|
---|
| 895 | mnlprocess(ref lm, ref workx, ref worky);
|
---|
| 896 |
|
---|
| 897 | //
|
---|
| 898 | // Logit version of the answer
|
---|
| 899 | //
|
---|
| 900 | nmax = 0;
|
---|
| 901 | for(j=0; j<=nclasses-1; j++)
|
---|
| 902 | {
|
---|
| 903 | if( (double)(worky[j])>(double)(worky[nmax]) )
|
---|
| 904 | {
|
---|
| 905 | nmax = j;
|
---|
| 906 | }
|
---|
| 907 | }
|
---|
| 908 |
|
---|
| 909 | //
|
---|
| 910 | // compare
|
---|
| 911 | //
|
---|
| 912 | if( nmax!=(int)Math.Round(xy[i,nvars]) )
|
---|
| 913 | {
|
---|
| 914 | result = result+1;
|
---|
| 915 | }
|
---|
| 916 | }
|
---|
| 917 | return result;
|
---|
| 918 | }
|
---|
| 919 |
|
---|
| 920 |
|
---|
| 921 | /*************************************************************************
|
---|
| 922 | Internal subroutine. Places exponents of the anti-overflow shifted
|
---|
| 923 | internal linear outputs into the service part of the W array.
|
---|
| 924 | *************************************************************************/
|
---|
| 925 | private static void mnliexp(ref double[] w,
|
---|
| 926 | ref double[] x)
|
---|
| 927 | {
|
---|
| 928 | int nvars = 0;
|
---|
| 929 | int nclasses = 0;
|
---|
| 930 | int offs = 0;
|
---|
| 931 | int i = 0;
|
---|
| 932 | int i1 = 0;
|
---|
| 933 | double v = 0;
|
---|
| 934 | double mx = 0;
|
---|
| 935 | int i_ = 0;
|
---|
| 936 | int i1_ = 0;
|
---|
| 937 |
|
---|
| 938 | System.Diagnostics.Debug.Assert((double)(w[1])==(double)(logitvnum), "LOGIT: unexpected model version");
|
---|
| 939 | nvars = (int)Math.Round(w[2]);
|
---|
| 940 | nclasses = (int)Math.Round(w[3]);
|
---|
| 941 | offs = (int)Math.Round(w[4]);
|
---|
| 942 | i1 = offs+(nvars+1)*(nclasses-1);
|
---|
| 943 | for(i=0; i<=nclasses-2; i++)
|
---|
| 944 | {
|
---|
| 945 | i1_ = (0)-(offs+i*(nvars+1));
|
---|
| 946 | v = 0.0;
|
---|
| 947 | for(i_=offs+i*(nvars+1); i_<=offs+i*(nvars+1)+nvars-1;i_++)
|
---|
| 948 | {
|
---|
| 949 | v += w[i_]*x[i_+i1_];
|
---|
| 950 | }
|
---|
| 951 | w[i1+i] = v+w[offs+i*(nvars+1)+nvars];
|
---|
| 952 | }
|
---|
| 953 | w[i1+nclasses-1] = 0;
|
---|
| 954 | mx = 0;
|
---|
| 955 | for(i=i1; i<=i1+nclasses-1; i++)
|
---|
| 956 | {
|
---|
| 957 | mx = Math.Max(mx, w[i]);
|
---|
| 958 | }
|
---|
| 959 | for(i=i1; i<=i1+nclasses-1; i++)
|
---|
| 960 | {
|
---|
| 961 | w[i] = Math.Exp(w[i]-mx);
|
---|
| 962 | }
|
---|
| 963 | }
|
---|
| 964 |
|
---|
| 965 |
|
---|
| 966 | /*************************************************************************
|
---|
| 967 | Calculation of all types of errors
|
---|
| 968 |
|
---|
| 969 | -- ALGLIB --
|
---|
| 970 | Copyright 30.08.2008 by Bochkanov Sergey
|
---|
| 971 | *************************************************************************/
|
---|
| 972 | private static void mnlallerrors(ref logitmodel lm,
|
---|
| 973 | ref double[,] xy,
|
---|
| 974 | int npoints,
|
---|
| 975 | ref double relcls,
|
---|
| 976 | ref double avgce,
|
---|
| 977 | ref double rms,
|
---|
| 978 | ref double avg,
|
---|
| 979 | ref double avgrel)
|
---|
| 980 | {
|
---|
| 981 | int nvars = 0;
|
---|
| 982 | int nclasses = 0;
|
---|
| 983 | int i = 0;
|
---|
| 984 | double[] buf = new double[0];
|
---|
| 985 | double[] workx = new double[0];
|
---|
| 986 | double[] y = new double[0];
|
---|
| 987 | double[] dy = new double[0];
|
---|
| 988 | int i_ = 0;
|
---|
| 989 |
|
---|
| 990 | System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==logitvnum, "MNL unit: Incorrect MNL version!");
|
---|
| 991 | nvars = (int)Math.Round(lm.w[2]);
|
---|
| 992 | nclasses = (int)Math.Round(lm.w[3]);
|
---|
| 993 | workx = new double[nvars-1+1];
|
---|
| 994 | y = new double[nclasses-1+1];
|
---|
| 995 | dy = new double[0+1];
|
---|
| 996 | bdss.dserrallocate(nclasses, ref buf);
|
---|
| 997 | for(i=0; i<=npoints-1; i++)
|
---|
| 998 | {
|
---|
| 999 | for(i_=0; i_<=nvars-1;i_++)
|
---|
| 1000 | {
|
---|
| 1001 | workx[i_] = xy[i,i_];
|
---|
| 1002 | }
|
---|
| 1003 | mnlprocess(ref lm, ref workx, ref y);
|
---|
| 1004 | dy[0] = xy[i,nvars];
|
---|
| 1005 | bdss.dserraccumulate(ref buf, ref y, ref dy);
|
---|
| 1006 | }
|
---|
| 1007 | bdss.dserrfinish(ref buf);
|
---|
| 1008 | relcls = buf[0];
|
---|
| 1009 | avgce = buf[1];
|
---|
| 1010 | rms = buf[2];
|
---|
| 1011 | avg = buf[3];
|
---|
| 1012 | avgrel = buf[4];
|
---|
| 1013 | }
|
---|
| 1014 |
|
---|
| 1015 |
|
---|
| 1016 | /*************************************************************************
|
---|
| 1017 | THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES A SUFFICIENT
|
---|
| 1018 | DECREASE CONDITION AND A CURVATURE CONDITION.
|
---|
| 1019 |
|
---|
| 1020 | AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF UNCERTAINTY WITH
|
---|
| 1021 | ENDPOINTS STX AND STY. THE INTERVAL OF UNCERTAINTY IS INITIALLY CHOSEN
|
---|
| 1022 | SO THAT IT CONTAINS A MINIMIZER OF THE MODIFIED FUNCTION
|
---|
| 1023 |
|
---|
| 1024 | F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
|
---|
| 1025 |
|
---|
| 1026 | IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE
|
---|
| 1027 | FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, THEN THE INTERVAL OF
|
---|
| 1028 | UNCERTAINTY IS CHOSEN SO THAT IT CONTAINS A MINIMIZER OF F(X+STP*S).
|
---|
| 1029 |
|
---|
| 1030 | THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES THE SUFFICIENT
|
---|
| 1031 | DECREASE CONDITION
|
---|
| 1032 |
|
---|
| 1033 | F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
|
---|
| 1034 |
|
---|
| 1035 | AND THE CURVATURE CONDITION
|
---|
| 1036 |
|
---|
| 1037 | ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
|
---|
| 1038 |
|
---|
| 1039 | IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION IS BOUNDED
|
---|
| 1040 | BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES BOTH CONDITIONS.
|
---|
| 1041 | IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH CONDITIONS, THEN THE
|
---|
| 1042 | ALGORITHM USUALLY STOPS WHEN ROUNDING ERRORS PREVENT FURTHER PROGRESS.
|
---|
| 1043 | IN THIS CASE STP ONLY SATISFIES THE SUFFICIENT DECREASE CONDITION.
|
---|
| 1044 |
|
---|
| 1045 | PARAMETERS DESCRIPRION
|
---|
| 1046 |
|
---|
| 1047 | N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER OF VARIABLES.
|
---|
| 1048 |
|
---|
| 1049 | X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE BASE POINT FOR
|
---|
| 1050 | THE LINE SEARCH. ON OUTPUT IT CONTAINS X+STP*S.
|
---|
| 1051 |
|
---|
| 1052 | F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F AT X. ON OUTPUT
|
---|
| 1053 | IT CONTAINS THE VALUE OF F AT X + STP*S.
|
---|
| 1054 |
|
---|
| 1055 | G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE GRADIENT OF F AT X.
|
---|
| 1056 | ON OUTPUT IT CONTAINS THE GRADIENT OF F AT X + STP*S.
|
---|
| 1057 |
|
---|
| 1058 | S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE SEARCH DIRECTION.
|
---|
| 1059 |
|
---|
| 1060 | STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN INITIAL ESTIMATE
|
---|
| 1061 | OF A SATISFACTORY STEP. ON OUTPUT STP CONTAINS THE FINAL ESTIMATE.
|
---|
| 1062 |
|
---|
| 1063 | FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION OCCURS WHEN THE
|
---|
| 1064 | SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
|
---|
| 1065 | SATISFIED.
|
---|
| 1066 |
|
---|
| 1067 | XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS WHEN THE RELATIVE
|
---|
| 1068 | WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL.
|
---|
| 1069 |
|
---|
| 1070 | STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH SPECIFY LOWER AND
|
---|
| 1071 | UPPER BOUNDS FOR THE STEP.
|
---|
| 1072 |
|
---|
| 1073 | MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION OCCURS WHEN THE
|
---|
| 1074 | NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV BY THE END OF AN ITERATION.
|
---|
| 1075 |
|
---|
| 1076 | INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
|
---|
| 1077 | INFO = 0 IMPROPER INPUT PARAMETERS.
|
---|
| 1078 |
|
---|
| 1079 | INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE
|
---|
| 1080 | DIRECTIONAL DERIVATIVE CONDITION HOLD.
|
---|
| 1081 |
|
---|
| 1082 | INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
|
---|
| 1083 | IS AT MOST XTOL.
|
---|
| 1084 |
|
---|
| 1085 | INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
|
---|
| 1086 |
|
---|
| 1087 | INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN.
|
---|
| 1088 |
|
---|
| 1089 | INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX.
|
---|
| 1090 |
|
---|
| 1091 | INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
|
---|
| 1092 | THERE MAY NOT BE A STEP WHICH SATISFIES THE
|
---|
| 1093 | SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
|
---|
| 1094 | TOLERANCES MAY BE TOO SMALL.
|
---|
| 1095 |
|
---|
| 1096 | NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF CALLS TO FCN.
|
---|
| 1097 |
|
---|
| 1098 | WA IS A WORK ARRAY OF LENGTH N.
|
---|
| 1099 |
|
---|
| 1100 | ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
|
---|
| 1101 | JORGE J. MORE', DAVID J. THUENTE
|
---|
| 1102 | *************************************************************************/
|
---|
| 1103 | private static void mnlmcsrch(int n,
|
---|
| 1104 | ref double[] x,
|
---|
| 1105 | ref double f,
|
---|
| 1106 | ref double[] g,
|
---|
| 1107 | ref double[] s,
|
---|
| 1108 | ref double stp,
|
---|
| 1109 | ref int info,
|
---|
| 1110 | ref int nfev,
|
---|
| 1111 | ref double[] wa,
|
---|
| 1112 | ref logitmcstate state,
|
---|
| 1113 | ref int stage)
|
---|
| 1114 | {
|
---|
| 1115 | double v = 0;
|
---|
| 1116 | double p5 = 0;
|
---|
| 1117 | double p66 = 0;
|
---|
| 1118 | double zero = 0;
|
---|
| 1119 | int i_ = 0;
|
---|
| 1120 |
|
---|
| 1121 |
|
---|
| 1122 | //
|
---|
| 1123 | // init
|
---|
| 1124 | //
|
---|
| 1125 | p5 = 0.5;
|
---|
| 1126 | p66 = 0.66;
|
---|
| 1127 | state.xtrapf = 4.0;
|
---|
| 1128 | zero = 0;
|
---|
| 1129 |
|
---|
| 1130 | //
|
---|
| 1131 | // Main cycle
|
---|
| 1132 | //
|
---|
| 1133 | while( true )
|
---|
| 1134 | {
|
---|
| 1135 | if( stage==0 )
|
---|
| 1136 | {
|
---|
| 1137 |
|
---|
| 1138 | //
|
---|
| 1139 | // NEXT
|
---|
| 1140 | //
|
---|
| 1141 | stage = 2;
|
---|
| 1142 | continue;
|
---|
| 1143 | }
|
---|
| 1144 | if( stage==2 )
|
---|
| 1145 | {
|
---|
| 1146 | state.infoc = 1;
|
---|
| 1147 | info = 0;
|
---|
| 1148 |
|
---|
| 1149 | //
|
---|
| 1150 | // CHECK THE INPUT PARAMETERS FOR ERRORS.
|
---|
| 1151 | //
|
---|
| 1152 | 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 )
|
---|
| 1153 | {
|
---|
| 1154 | stage = 0;
|
---|
| 1155 | return;
|
---|
| 1156 | }
|
---|
| 1157 |
|
---|
| 1158 | //
|
---|
| 1159 | // COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
|
---|
| 1160 | // AND CHECK THAT S IS A DESCENT DIRECTION.
|
---|
| 1161 | //
|
---|
| 1162 | v = 0.0;
|
---|
| 1163 | for(i_=0; i_<=n-1;i_++)
|
---|
| 1164 | {
|
---|
| 1165 | v += g[i_]*s[i_];
|
---|
| 1166 | }
|
---|
| 1167 | state.dginit = v;
|
---|
| 1168 | if( (double)(state.dginit)>=(double)(0) )
|
---|
| 1169 | {
|
---|
| 1170 | stage = 0;
|
---|
| 1171 | return;
|
---|
| 1172 | }
|
---|
| 1173 |
|
---|
| 1174 | //
|
---|
| 1175 | // INITIALIZE LOCAL VARIABLES.
|
---|
| 1176 | //
|
---|
| 1177 | state.brackt = false;
|
---|
| 1178 | state.stage1 = true;
|
---|
| 1179 | nfev = 0;
|
---|
| 1180 | state.finit = f;
|
---|
| 1181 | state.dgtest = ftol*state.dginit;
|
---|
| 1182 | state.width = stpmax-stpmin;
|
---|
| 1183 | state.width1 = state.width/p5;
|
---|
| 1184 | for(i_=0; i_<=n-1;i_++)
|
---|
| 1185 | {
|
---|
| 1186 | wa[i_] = x[i_];
|
---|
| 1187 | }
|
---|
| 1188 |
|
---|
| 1189 | //
|
---|
| 1190 | // THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
|
---|
| 1191 | // FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
|
---|
| 1192 | // THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
|
---|
| 1193 | // FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
|
---|
| 1194 | // THE INTERVAL OF UNCERTAINTY.
|
---|
| 1195 | // THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
|
---|
| 1196 | // FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
|
---|
| 1197 | //
|
---|
| 1198 | state.stx = 0;
|
---|
| 1199 | state.fx = state.finit;
|
---|
| 1200 | state.dgx = state.dginit;
|
---|
| 1201 | state.sty = 0;
|
---|
| 1202 | state.fy = state.finit;
|
---|
| 1203 | state.dgy = state.dginit;
|
---|
| 1204 |
|
---|
| 1205 | //
|
---|
| 1206 | // NEXT
|
---|
| 1207 | //
|
---|
| 1208 | stage = 3;
|
---|
| 1209 | continue;
|
---|
| 1210 | }
|
---|
| 1211 | if( stage==3 )
|
---|
| 1212 | {
|
---|
| 1213 |
|
---|
| 1214 | //
|
---|
| 1215 | // START OF ITERATION.
|
---|
| 1216 | //
|
---|
| 1217 | // SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
|
---|
| 1218 | // TO THE PRESENT INTERVAL OF UNCERTAINTY.
|
---|
| 1219 | //
|
---|
| 1220 | if( state.brackt )
|
---|
| 1221 | {
|
---|
| 1222 | if( (double)(state.stx)<(double)(state.sty) )
|
---|
| 1223 | {
|
---|
| 1224 | state.stmin = state.stx;
|
---|
| 1225 | state.stmax = state.sty;
|
---|
| 1226 | }
|
---|
| 1227 | else
|
---|
| 1228 | {
|
---|
| 1229 | state.stmin = state.sty;
|
---|
| 1230 | state.stmax = state.stx;
|
---|
| 1231 | }
|
---|
| 1232 | }
|
---|
| 1233 | else
|
---|
| 1234 | {
|
---|
| 1235 | state.stmin = state.stx;
|
---|
| 1236 | state.stmax = stp+state.xtrapf*(stp-state.stx);
|
---|
| 1237 | }
|
---|
| 1238 |
|
---|
| 1239 | //
|
---|
| 1240 | // FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
|
---|
| 1241 | //
|
---|
| 1242 | if( (double)(stp)>(double)(stpmax) )
|
---|
| 1243 | {
|
---|
| 1244 | stp = stpmax;
|
---|
| 1245 | }
|
---|
| 1246 | if( (double)(stp)<(double)(stpmin) )
|
---|
| 1247 | {
|
---|
| 1248 | stp = stpmin;
|
---|
| 1249 | }
|
---|
| 1250 |
|
---|
| 1251 | //
|
---|
| 1252 | // IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
|
---|
| 1253 | // STP BE THE LOWEST POINT OBTAINED SO FAR.
|
---|
| 1254 | //
|
---|
| 1255 | 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) )
|
---|
| 1256 | {
|
---|
| 1257 | stp = state.stx;
|
---|
| 1258 | }
|
---|
| 1259 |
|
---|
| 1260 | //
|
---|
| 1261 | // EVALUATE THE FUNCTION AND GRADIENT AT STP
|
---|
| 1262 | // AND COMPUTE THE DIRECTIONAL DERIVATIVE.
|
---|
| 1263 | //
|
---|
| 1264 | for(i_=0; i_<=n-1;i_++)
|
---|
| 1265 | {
|
---|
| 1266 | x[i_] = wa[i_];
|
---|
| 1267 | }
|
---|
| 1268 | for(i_=0; i_<=n-1;i_++)
|
---|
| 1269 | {
|
---|
| 1270 | x[i_] = x[i_] + stp*s[i_];
|
---|
| 1271 | }
|
---|
| 1272 |
|
---|
| 1273 | //
|
---|
| 1274 | // NEXT
|
---|
| 1275 | //
|
---|
| 1276 | stage = 4;
|
---|
| 1277 | return;
|
---|
| 1278 | }
|
---|
| 1279 | if( stage==4 )
|
---|
| 1280 | {
|
---|
| 1281 | info = 0;
|
---|
| 1282 | nfev = nfev+1;
|
---|
| 1283 | v = 0.0;
|
---|
| 1284 | for(i_=0; i_<=n-1;i_++)
|
---|
| 1285 | {
|
---|
| 1286 | v += g[i_]*s[i_];
|
---|
| 1287 | }
|
---|
| 1288 | state.dg = v;
|
---|
| 1289 | state.ftest1 = state.finit+stp*state.dgtest;
|
---|
| 1290 |
|
---|
| 1291 | //
|
---|
| 1292 | // TEST FOR CONVERGENCE.
|
---|
| 1293 | //
|
---|
| 1294 | if( state.brackt & ((double)(stp)<=(double)(state.stmin) | (double)(stp)>=(double)(state.stmax)) | state.infoc==0 )
|
---|
| 1295 | {
|
---|
| 1296 | info = 6;
|
---|
| 1297 | }
|
---|
| 1298 | if( (double)(stp)==(double)(stpmax) & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)<=(double)(state.dgtest) )
|
---|
| 1299 | {
|
---|
| 1300 | info = 5;
|
---|
| 1301 | }
|
---|
| 1302 | if( (double)(stp)==(double)(stpmin) & ((double)(f)>(double)(state.ftest1) | (double)(state.dg)>=(double)(state.dgtest)) )
|
---|
| 1303 | {
|
---|
| 1304 | info = 4;
|
---|
| 1305 | }
|
---|
| 1306 | if( nfev>=maxfev )
|
---|
| 1307 | {
|
---|
| 1308 | info = 3;
|
---|
| 1309 | }
|
---|
| 1310 | if( state.brackt & (double)(state.stmax-state.stmin)<=(double)(xtol*state.stmax) )
|
---|
| 1311 | {
|
---|
| 1312 | info = 2;
|
---|
| 1313 | }
|
---|
| 1314 | if( (double)(f)<=(double)(state.ftest1) & (double)(Math.Abs(state.dg))<=(double)(-(gtol*state.dginit)) )
|
---|
| 1315 | {
|
---|
| 1316 | info = 1;
|
---|
| 1317 | }
|
---|
| 1318 |
|
---|
| 1319 | //
|
---|
| 1320 | // CHECK FOR TERMINATION.
|
---|
| 1321 | //
|
---|
| 1322 | if( info!=0 )
|
---|
| 1323 | {
|
---|
| 1324 | stage = 0;
|
---|
| 1325 | return;
|
---|
| 1326 | }
|
---|
| 1327 |
|
---|
| 1328 | //
|
---|
| 1329 | // IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
|
---|
| 1330 | // FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
|
---|
| 1331 | //
|
---|
| 1332 | if( state.stage1 & (double)(f)<=(double)(state.ftest1) & (double)(state.dg)>=(double)(Math.Min(ftol, gtol)*state.dginit) )
|
---|
| 1333 | {
|
---|
| 1334 | state.stage1 = false;
|
---|
| 1335 | }
|
---|
| 1336 |
|
---|
| 1337 | //
|
---|
| 1338 | // A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
|
---|
| 1339 | // WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
|
---|
| 1340 | // FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
|
---|
| 1341 | // DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
|
---|
| 1342 | // OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
|
---|
| 1343 | //
|
---|
| 1344 | if( state.stage1 & (double)(f)<=(double)(state.fx) & (double)(f)>(double)(state.ftest1) )
|
---|
| 1345 | {
|
---|
| 1346 |
|
---|
| 1347 | //
|
---|
| 1348 | // DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
|
---|
| 1349 | //
|
---|
| 1350 | state.fm = f-stp*state.dgtest;
|
---|
| 1351 | state.fxm = state.fx-state.stx*state.dgtest;
|
---|
| 1352 | state.fym = state.fy-state.sty*state.dgtest;
|
---|
| 1353 | state.dgm = state.dg-state.dgtest;
|
---|
| 1354 | state.dgxm = state.dgx-state.dgtest;
|
---|
| 1355 | state.dgym = state.dgy-state.dgtest;
|
---|
| 1356 |
|
---|
| 1357 | //
|
---|
| 1358 | // CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
|
---|
| 1359 | // AND TO COMPUTE THE NEW STEP.
|
---|
| 1360 | //
|
---|
| 1361 | mnlmcstep(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);
|
---|
| 1362 |
|
---|
| 1363 | //
|
---|
| 1364 | // RESET THE FUNCTION AND GRADIENT VALUES FOR F.
|
---|
| 1365 | //
|
---|
| 1366 | state.fx = state.fxm+state.stx*state.dgtest;
|
---|
| 1367 | state.fy = state.fym+state.sty*state.dgtest;
|
---|
| 1368 | state.dgx = state.dgxm+state.dgtest;
|
---|
| 1369 | state.dgy = state.dgym+state.dgtest;
|
---|
| 1370 | }
|
---|
| 1371 | else
|
---|
| 1372 | {
|
---|
| 1373 |
|
---|
| 1374 | //
|
---|
| 1375 | // CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
|
---|
| 1376 | // AND TO COMPUTE THE NEW STEP.
|
---|
| 1377 | //
|
---|
| 1378 | mnlmcstep(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);
|
---|
| 1379 | }
|
---|
| 1380 |
|
---|
| 1381 | //
|
---|
| 1382 | // FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
|
---|
| 1383 | // INTERVAL OF UNCERTAINTY.
|
---|
| 1384 | //
|
---|
| 1385 | if( state.brackt )
|
---|
| 1386 | {
|
---|
| 1387 | if( (double)(Math.Abs(state.sty-state.stx))>=(double)(p66*state.width1) )
|
---|
| 1388 | {
|
---|
| 1389 | stp = state.stx+p5*(state.sty-state.stx);
|
---|
| 1390 | }
|
---|
| 1391 | state.width1 = state.width;
|
---|
| 1392 | state.width = Math.Abs(state.sty-state.stx);
|
---|
| 1393 | }
|
---|
| 1394 |
|
---|
| 1395 | //
|
---|
| 1396 | // NEXT.
|
---|
| 1397 | //
|
---|
| 1398 | stage = 3;
|
---|
| 1399 | continue;
|
---|
| 1400 | }
|
---|
| 1401 | }
|
---|
| 1402 | }
|
---|
| 1403 |
|
---|
| 1404 |
|
---|
| 1405 | private static void mnlmcstep(ref double stx,
|
---|
| 1406 | ref double fx,
|
---|
| 1407 | ref double dx,
|
---|
| 1408 | ref double sty,
|
---|
| 1409 | ref double fy,
|
---|
| 1410 | ref double dy,
|
---|
| 1411 | ref double stp,
|
---|
| 1412 | double fp,
|
---|
| 1413 | double dp,
|
---|
| 1414 | ref bool brackt,
|
---|
| 1415 | double stmin,
|
---|
| 1416 | double stmax,
|
---|
| 1417 | ref int info)
|
---|
| 1418 | {
|
---|
| 1419 | bool bound = new bool();
|
---|
| 1420 | double gamma = 0;
|
---|
| 1421 | double p = 0;
|
---|
| 1422 | double q = 0;
|
---|
| 1423 | double r = 0;
|
---|
| 1424 | double s = 0;
|
---|
| 1425 | double sgnd = 0;
|
---|
| 1426 | double stpc = 0;
|
---|
| 1427 | double stpf = 0;
|
---|
| 1428 | double stpq = 0;
|
---|
| 1429 | double theta = 0;
|
---|
| 1430 |
|
---|
| 1431 | info = 0;
|
---|
| 1432 |
|
---|
| 1433 | //
|
---|
| 1434 | // CHECK THE INPUT PARAMETERS FOR ERRORS.
|
---|
| 1435 | //
|
---|
| 1436 | 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) )
|
---|
| 1437 | {
|
---|
| 1438 | return;
|
---|
| 1439 | }
|
---|
| 1440 |
|
---|
| 1441 | //
|
---|
| 1442 | // DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
|
---|
| 1443 | //
|
---|
| 1444 | sgnd = dp*(dx/Math.Abs(dx));
|
---|
| 1445 |
|
---|
| 1446 | //
|
---|
| 1447 | // FIRST CASE. A HIGHER FUNCTION VALUE.
|
---|
| 1448 | // THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
|
---|
| 1449 | // TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
|
---|
| 1450 | // ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
|
---|
| 1451 | //
|
---|
| 1452 | if( (double)(fp)>(double)(fx) )
|
---|
| 1453 | {
|
---|
| 1454 | info = 1;
|
---|
| 1455 | bound = true;
|
---|
| 1456 | theta = 3*(fx-fp)/(stp-stx)+dx+dp;
|
---|
| 1457 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
|
---|
| 1458 | gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dx/s*(dp/s));
|
---|
| 1459 | if( (double)(stp)<(double)(stx) )
|
---|
| 1460 | {
|
---|
| 1461 | gamma = -gamma;
|
---|
| 1462 | }
|
---|
| 1463 | p = gamma-dx+theta;
|
---|
| 1464 | q = gamma-dx+gamma+dp;
|
---|
| 1465 | r = p/q;
|
---|
| 1466 | stpc = stx+r*(stp-stx);
|
---|
| 1467 | stpq = stx+dx/((fx-fp)/(stp-stx)+dx)/2*(stp-stx);
|
---|
| 1468 | if( (double)(Math.Abs(stpc-stx))<(double)(Math.Abs(stpq-stx)) )
|
---|
| 1469 | {
|
---|
| 1470 | stpf = stpc;
|
---|
| 1471 | }
|
---|
| 1472 | else
|
---|
| 1473 | {
|
---|
| 1474 | stpf = stpc+(stpq-stpc)/2;
|
---|
| 1475 | }
|
---|
| 1476 | brackt = true;
|
---|
| 1477 | }
|
---|
| 1478 | else
|
---|
| 1479 | {
|
---|
| 1480 | if( (double)(sgnd)<(double)(0) )
|
---|
| 1481 | {
|
---|
| 1482 |
|
---|
| 1483 | //
|
---|
| 1484 | // SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
|
---|
| 1485 | // OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
|
---|
| 1486 | // STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
|
---|
| 1487 | // THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
|
---|
| 1488 | //
|
---|
| 1489 | info = 2;
|
---|
| 1490 | bound = false;
|
---|
| 1491 | theta = 3*(fx-fp)/(stp-stx)+dx+dp;
|
---|
| 1492 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
|
---|
| 1493 | gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dx/s*(dp/s));
|
---|
| 1494 | if( (double)(stp)>(double)(stx) )
|
---|
| 1495 | {
|
---|
| 1496 | gamma = -gamma;
|
---|
| 1497 | }
|
---|
| 1498 | p = gamma-dp+theta;
|
---|
| 1499 | q = gamma-dp+gamma+dx;
|
---|
| 1500 | r = p/q;
|
---|
| 1501 | stpc = stp+r*(stx-stp);
|
---|
| 1502 | stpq = stp+dp/(dp-dx)*(stx-stp);
|
---|
| 1503 | if( (double)(Math.Abs(stpc-stp))>(double)(Math.Abs(stpq-stp)) )
|
---|
| 1504 | {
|
---|
| 1505 | stpf = stpc;
|
---|
| 1506 | }
|
---|
| 1507 | else
|
---|
| 1508 | {
|
---|
| 1509 | stpf = stpq;
|
---|
| 1510 | }
|
---|
| 1511 | brackt = true;
|
---|
| 1512 | }
|
---|
| 1513 | else
|
---|
| 1514 | {
|
---|
| 1515 | if( (double)(Math.Abs(dp))<(double)(Math.Abs(dx)) )
|
---|
| 1516 | {
|
---|
| 1517 |
|
---|
| 1518 | //
|
---|
| 1519 | // THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
|
---|
| 1520 | // SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
|
---|
| 1521 | // THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
|
---|
| 1522 | // IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
|
---|
| 1523 | // IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
|
---|
| 1524 | // EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
|
---|
| 1525 | // COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
|
---|
| 1526 | // CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
|
---|
| 1527 | //
|
---|
| 1528 | info = 3;
|
---|
| 1529 | bound = true;
|
---|
| 1530 | theta = 3*(fx-fp)/(stp-stx)+dx+dp;
|
---|
| 1531 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dx), Math.Abs(dp)));
|
---|
| 1532 |
|
---|
| 1533 | //
|
---|
| 1534 | // THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
|
---|
| 1535 | // TO INFINITY IN THE DIRECTION OF THE STEP.
|
---|
| 1536 | //
|
---|
| 1537 | gamma = s*Math.Sqrt(Math.Max(0, AP.Math.Sqr(theta/s)-dx/s*(dp/s)));
|
---|
| 1538 | if( (double)(stp)>(double)(stx) )
|
---|
| 1539 | {
|
---|
| 1540 | gamma = -gamma;
|
---|
| 1541 | }
|
---|
| 1542 | p = gamma-dp+theta;
|
---|
| 1543 | q = gamma+(dx-dp)+gamma;
|
---|
| 1544 | r = p/q;
|
---|
| 1545 | if( (double)(r)<(double)(0) & (double)(gamma)!=(double)(0) )
|
---|
| 1546 | {
|
---|
| 1547 | stpc = stp+r*(stx-stp);
|
---|
| 1548 | }
|
---|
| 1549 | else
|
---|
| 1550 | {
|
---|
| 1551 | if( (double)(stp)>(double)(stx) )
|
---|
| 1552 | {
|
---|
| 1553 | stpc = stmax;
|
---|
| 1554 | }
|
---|
| 1555 | else
|
---|
| 1556 | {
|
---|
| 1557 | stpc = stmin;
|
---|
| 1558 | }
|
---|
| 1559 | }
|
---|
| 1560 | stpq = stp+dp/(dp-dx)*(stx-stp);
|
---|
| 1561 | if( brackt )
|
---|
| 1562 | {
|
---|
| 1563 | if( (double)(Math.Abs(stp-stpc))<(double)(Math.Abs(stp-stpq)) )
|
---|
| 1564 | {
|
---|
| 1565 | stpf = stpc;
|
---|
| 1566 | }
|
---|
| 1567 | else
|
---|
| 1568 | {
|
---|
| 1569 | stpf = stpq;
|
---|
| 1570 | }
|
---|
| 1571 | }
|
---|
| 1572 | else
|
---|
| 1573 | {
|
---|
| 1574 | if( (double)(Math.Abs(stp-stpc))>(double)(Math.Abs(stp-stpq)) )
|
---|
| 1575 | {
|
---|
| 1576 | stpf = stpc;
|
---|
| 1577 | }
|
---|
| 1578 | else
|
---|
| 1579 | {
|
---|
| 1580 | stpf = stpq;
|
---|
| 1581 | }
|
---|
| 1582 | }
|
---|
| 1583 | }
|
---|
| 1584 | else
|
---|
| 1585 | {
|
---|
| 1586 |
|
---|
| 1587 | //
|
---|
| 1588 | // FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
|
---|
| 1589 | // SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
|
---|
| 1590 | // NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
|
---|
| 1591 | // IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
|
---|
| 1592 | //
|
---|
| 1593 | info = 4;
|
---|
| 1594 | bound = false;
|
---|
| 1595 | if( brackt )
|
---|
| 1596 | {
|
---|
| 1597 | theta = 3*(fp-fy)/(sty-stp)+dy+dp;
|
---|
| 1598 | s = Math.Max(Math.Abs(theta), Math.Max(Math.Abs(dy), Math.Abs(dp)));
|
---|
| 1599 | gamma = s*Math.Sqrt(AP.Math.Sqr(theta/s)-dy/s*(dp/s));
|
---|
| 1600 | if( (double)(stp)>(double)(sty) )
|
---|
| 1601 | {
|
---|
| 1602 | gamma = -gamma;
|
---|
| 1603 | }
|
---|
| 1604 | p = gamma-dp+theta;
|
---|
| 1605 | q = gamma-dp+gamma+dy;
|
---|
| 1606 | r = p/q;
|
---|
| 1607 | stpc = stp+r*(sty-stp);
|
---|
| 1608 | stpf = stpc;
|
---|
| 1609 | }
|
---|
| 1610 | else
|
---|
| 1611 | {
|
---|
| 1612 | if( (double)(stp)>(double)(stx) )
|
---|
| 1613 | {
|
---|
| 1614 | stpf = stmax;
|
---|
| 1615 | }
|
---|
| 1616 | else
|
---|
| 1617 | {
|
---|
| 1618 | stpf = stmin;
|
---|
| 1619 | }
|
---|
| 1620 | }
|
---|
| 1621 | }
|
---|
| 1622 | }
|
---|
| 1623 | }
|
---|
| 1624 |
|
---|
| 1625 | //
|
---|
| 1626 | // UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
|
---|
| 1627 | // DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
|
---|
| 1628 | //
|
---|
| 1629 | if( (double)(fp)>(double)(fx) )
|
---|
| 1630 | {
|
---|
| 1631 | sty = stp;
|
---|
| 1632 | fy = fp;
|
---|
| 1633 | dy = dp;
|
---|
| 1634 | }
|
---|
| 1635 | else
|
---|
| 1636 | {
|
---|
| 1637 | if( (double)(sgnd)<(double)(0.0) )
|
---|
| 1638 | {
|
---|
| 1639 | sty = stx;
|
---|
| 1640 | fy = fx;
|
---|
| 1641 | dy = dx;
|
---|
| 1642 | }
|
---|
| 1643 | stx = stp;
|
---|
| 1644 | fx = fp;
|
---|
| 1645 | dx = dp;
|
---|
| 1646 | }
|
---|
| 1647 |
|
---|
| 1648 | //
|
---|
| 1649 | // COMPUTE THE NEW STEP AND SAFEGUARD IT.
|
---|
| 1650 | //
|
---|
| 1651 | stpf = Math.Min(stmax, stpf);
|
---|
| 1652 | stpf = Math.Max(stmin, stpf);
|
---|
| 1653 | stp = stpf;
|
---|
| 1654 | if( brackt & bound )
|
---|
| 1655 | {
|
---|
| 1656 | if( (double)(sty)>(double)(stx) )
|
---|
| 1657 | {
|
---|
| 1658 | stp = Math.Min(stx+0.66*(sty-stx), stp);
|
---|
| 1659 | }
|
---|
| 1660 | else
|
---|
| 1661 | {
|
---|
| 1662 | stp = Math.Max(stx+0.66*(sty-stx), stp);
|
---|
| 1663 | }
|
---|
| 1664 | }
|
---|
| 1665 | }
|
---|
| 1666 | }
|
---|
| 1667 | }
|
---|