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