source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/math/Minimisation.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 40.8 KB
Line 
1/*
2* Class Minimisation
3*
4* Contains methods for finding the values of the
5* function parameters that minimise that function
6* using the Nelder and Mead Simplex method.
7*
8* The function needed by the minimisation method
9* is supplied by though the interface, MinimisationFunction or Minimization
10*
11* WRITTEN BY: Dr Michael Thomas Flanagan
12*
13* DATE: April 2003
14* MODIFIED: 29 December 2005, 18 February 2006, 28 December 2007, 10/12 May 2008,
15* July 2008, 12 January 2010, 23 September 2010
16*
17* DOCUMENTATION:
18* See Michael Thomas Flanagan's Java library on-line web page:
19* http://www.ee.ucl.ac.uk/~mflanaga/java/Minimisation.html
20* http://www.ee.ucl.ac.uk/~mflanaga/java/
21*
22* Copyright (c) 2003 - 2008
23*
24* PERMISSION TO COPY:
25* Permission to use, copy and modify this software and its documentation for
26* NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement
27* to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.
28*
29* Dr Michael Thomas Flanagan makes no representations about the suitability
30* or fitness of the software for any or for a particular purpose.
31* Michael Thomas Flanagan shall not be liable for any damages suffered
32* as a result of using, modifying or distributing this software or its derivatives.
33*
34***************************************************************************************/
35
36package agents.anac.y2015.agentBuyogV2.flanagan.math;
37
38import java.util.*;
39
40import agents.anac.y2015.agentBuyogV2.flanagan.io.FileOutput;
41import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
42
43
44// Minimisation class
45public class Minimisation{
46
47 protected boolean iseOption = true; // = true if minimis... spelling used
48 // = false if minimiz... spelling used
49 protected int nParam=0; // number of unknown parameters to be estimated
50 protected double[]paramValue = null; // function parameter values (returned at function minimum)
51 protected String[] paraName = null; // names of parameters, eg, c[0], c[1], c[2] . . .
52 protected double functValue = 0.0D; // current value of the function to be minimised
53 protected double lastFunctValNoCnstrnt=0.0D;// Last function value with no constraint penalty
54 protected double minimum = 0.0D; // value of the function to be minimised at the minimum
55 protected int prec = 4; // number of places to which double variables are truncated on output to text files
56 protected int field = 13; // field width on output to text files
57 protected boolean convStatus = false; // Status of minimisation on exiting minimisation method
58 // = true - convergence criterion was met
59 // = false - convergence criterion not met - current estimates returned
60 protected boolean suppressNoConvergenceMessage = false; // if true - suppress the print out of a message saying that convergence was not achieved.
61 protected int scaleOpt=0; // if = 0; no scaling of initial estimates
62 // if = 1; initial simplex estimates scaled to unity
63 // if = 2; initial estimates scaled by user provided values in scale[]
64 // (default = 0)
65 protected double[] scale = null; // values to scale initial estimate (see scaleOpt above)
66 protected boolean penalty = false; // true if single parameter penalty function is included
67 protected boolean sumPenalty = false; // true if multiple parameter penalty function is included
68 protected int nConstraints = 0; // number of single parameter constraints
69 protected int nSumConstraints = 0; // number of multiple parameter constraints
70 protected int maxConstraintIndex = -1; // maximum index of constrained parameter/s
71 protected ArrayList<Object> penalties = new ArrayList<Object>(); // method index,
72 // number of single parameter constraints,
73 // then repeated for each constraint:
74 // penalty parameter index,
75 // below or above constraint flag,
76 // constraint boundary value
77 protected ArrayList<Object> sumPenalties = new ArrayList<Object>(); // constraint method index,
78 // number of multiple parameter constraints,
79 // then repeated for each constraint:
80 // number of parameters in summation
81 // penalty parameter indices,
82 // summation signs
83 // below or above constraint flag,
84 // constraint boundary value
85 protected int[] penaltyCheck = null; // = -1 values below the single constraint boundary not allowed
86 // = +1 values above the single constraint boundary not allowed
87 protected int[] sumPenaltyCheck = null; // = -1 values below the multiple constraint boundary not allowed
88 // = +1 values above the multiple constraint boundary not allowed
89 protected double penaltyWeight = 1.0e30; // weight for the penalty functions
90 protected int[] penaltyParam = null; // indices of paramaters subject to single parameter constraint
91 protected int[][] sumPenaltyParam = null; // indices of paramaters subject to multiple parameter constraint
92 protected double[][] sumPlusOrMinus = null; // value before each parameter in multiple parameter summation
93 protected int[] sumPenaltyNumber = null; // number of paramaters in each multiple parameter constraint
94
95 protected double[] constraints = null; // single parameter constraint values
96 protected double constraintTolerance = 1e-4; // tolerance in constraining parameter/s to a fixed value
97 protected double[] sumConstraints = null; // multiple parameter constraint values
98 protected int constraintMethod = 0; // constraint method number
99 // =0: cliff to the power two (only method at present)
100 protected int nMax = 3000; // Nelder and Mead simplex maximum number of iterations
101 protected int nIter = 0; // Nelder and Mead simplex number of iterations performed
102 protected int konvge = 3; // Nelder and Mead simplex number of restarts allowed
103 protected int kRestart = 0; // Nelder and Mead simplex number of restarts taken
104 protected double fTol = 1e-13; // Nelder and Mead simplex convergence tolerance
105 protected double rCoeff = 1.0D; // Nelder and Mead simplex reflection coefficient
106 protected double eCoeff = 2.0D; // Nelder and Mead simplex extension coefficient
107 protected double cCoeff = 0.5D; // Nelder and Mead simplex contraction coefficient
108 protected double[] startH = null; // Nelder and Mead simplex initial estimates
109 protected double[] step = null; // Nelder and Mead simplex step values
110 protected double dStep = 0.5D; // Nelder and Mead simplex default step value
111 protected int minTest = 0; // Nelder and Mead minimum test
112 // = 0; tests simplex sd < fTol
113 // allows options for further tests to be added later
114 protected double simplexSd = 0.0D; // simplex standard deviation
115
116
117 //Constructor
118 public Minimisation(){
119 this.iseOption = true;
120 }
121
122 // Suppress the print out of a message saying that convergence was not achieved.
123 public void suppressNoConvergenceMessage(){
124 this.suppressNoConvergenceMessage = true;
125 }
126
127 // Suppress the print out of a message saying that convergence was not achieved.
128 public void supressNoConvergenceMessage(){
129 this.suppressNoConvergenceMessage = true;
130 }
131
132 // Nelder and Mead Simplex minimisation
133 public void nelderMead(MinimisationFunction gg, double[] start, double[] step, double fTol, int nMax){
134 Object g = (Object)gg;
135 this.nelderMead(g, start, step, fTol, nMax);
136 }
137
138 // Nelder and Mead Simplex minimisation
139 public void nelderMead(MinimizationFunction gg, double[] start, double[] step, double fTol, int nMax){
140 Object g = (Object)gg;
141 this.nelderMead(g, start, step, fTol, nMax);
142 }
143
144 // Nelder and Mead Simplex minimisation
145 public void nelderMead(Object g, double[] start, double[] step, double fTol, int nMax){
146
147 boolean testContract=false; // test whether a simplex contraction has been performed
148 int np = start.length; // number of unknown parameters;
149 if(this.maxConstraintIndex>=np)throw new IllegalArgumentException("You have entered more constrained parameters ("+this.maxConstraintIndex+") than minimisation parameters (" + np + ")");
150 this.nParam = np;
151 this.convStatus = true;
152 int nnp = np+1; // Number of simplex apices
153 this.lastFunctValNoCnstrnt=0.0D;
154
155 if(this.scaleOpt<2)this.scale = new double[np];
156 if(scaleOpt==2 && scale.length!=start.length)throw new IllegalArgumentException("scale array and initial estimate array are of different lengths");
157 if(step.length!=start.length)throw new IllegalArgumentException("step array length " + step.length + " and initial estimate array length " + start.length + " are of different");
158
159 // check for zero step sizes
160 for(int i=0; i<np; i++){
161 if(step[i]==0.0D){
162 if(start[i]!=0.0){
163 step[i] = start[i]*0.1;
164 }
165 else{
166 step[i] = 1.0;
167 System.out.println("As no step size has been entered for an itial estimate of zero an initial step size of unity has been used");
168 System.out.println("You are advised to repeat the minimization using one of the methods allowing the setting of an appropriate non-zero initial step size");
169 }
170 }
171 }
172
173 // set up arrays
174 this.paramValue = new double[np];
175 this.startH = new double[np];
176 this.step = new double[np];
177 double[]pmin = new double[np]; //Nelder and Mead Pmin
178
179 double[][] pp = new double[nnp][nnp]; //Nelder and Mead P
180 double[] yy = new double[nnp]; //Nelder and Mead y
181 double[] pbar = new double[nnp]; //Nelder and Mead P with bar superscript
182 double[] pstar = new double[nnp]; //Nelder and Mead P*
183 double[] p2star = new double[nnp]; //Nelder and Mead P**
184
185 // Set any single parameter constraint parameters
186 if(this.penalty){
187 Integer itemp = (Integer)this.penalties.get(1);
188 this.nConstraints = itemp.intValue();
189 this.penaltyParam = new int[this.nConstraints];
190 this.penaltyCheck = new int[this.nConstraints];
191 this.constraints = new double[this.nConstraints];
192 Double dtemp = null;
193 int j=2;
194 for(int i=0;i<this.nConstraints;i++){
195 itemp = (Integer)this.penalties.get(j);
196 this.penaltyParam[i] = itemp.intValue();
197 j++;
198 itemp = (Integer)this.penalties.get(j);
199 this.penaltyCheck[i] = itemp.intValue();
200 j++;
201 dtemp = (Double)this.penalties.get(j);
202 this.constraints[i] = dtemp.doubleValue();
203 j++;
204 }
205 }
206
207 // Set any multiple parameter constraint parameters
208 if(this.sumPenalty){
209 Integer itemp = (Integer)this.sumPenalties.get(1);
210 this.nSumConstraints = itemp.intValue();
211 this.sumPenaltyParam = new int[this.nSumConstraints][];
212 this.sumPlusOrMinus = new double[this.nSumConstraints][];
213 this.sumPenaltyCheck = new int[this.nSumConstraints];
214 this.sumPenaltyNumber = new int[this.nSumConstraints];
215 this.sumConstraints = new double[this.nSumConstraints];
216 int[] itempArray = null;
217 double[] dtempArray = null;
218 Double dtemp = null;
219 int j=2;
220 for(int i=0;i<this.nSumConstraints;i++){
221 itemp = (Integer)this.sumPenalties.get(j);
222 this.sumPenaltyNumber[i] = itemp.intValue();
223 j++;
224 itempArray = (int[])this.sumPenalties.get(j);
225 this.sumPenaltyParam[i] = itempArray;
226 j++;
227 dtempArray = (double[])this.sumPenalties.get(j);
228 this.sumPlusOrMinus[i] = dtempArray;
229 j++;
230 itemp = (Integer)this.sumPenalties.get(j);
231 this.sumPenaltyCheck[i] = itemp.intValue();
232 j++;
233 dtemp = (Double)this.sumPenalties.get(j);
234 this.sumConstraints[i] = dtemp.doubleValue();
235 j++;
236 }
237 }
238
239 // Store unscaled start values
240 for(int i=0; i<np; i++)this.startH[i]=start[i];
241
242 // scale initial estimates and step sizes
243 if(this.scaleOpt>0){
244 boolean testzero=false;
245 for(int i=0; i<np; i++)if(start[i]==0.0D)testzero=true;
246 if(testzero){
247 System.out.println("Neler and Mead Simplex: a start value of zero precludes scaling");
248 System.out.println("Regression performed without scaling");
249 this.scaleOpt=0;
250 }
251 }
252 switch(this.scaleOpt){
253 case 0: for(int i=0; i<np; i++)scale[i]=1.0D;
254 break;
255 case 1: for(int i=0; i<np; i++){
256 scale[i]=1.0/start[i];
257 step[i]=step[i]/start[i];
258 start[i]=1.0D;
259 }
260 break;
261 case 2: for(int i=0; i<np; i++){
262 step[i]*=scale[i];
263 start[i]*= scale[i];
264 }
265 break;
266 }
267
268 // set class member values
269 this.fTol=fTol;
270 this.nMax=nMax;
271 this.nIter=0;
272 for(int i=0; i<np; i++){
273 this.step[i]=step[i];
274 this.scale[i]=scale[i];
275 }
276
277 // initial simplex
278 double sho=0.0D;
279 for (int i=0; i<np; ++i){
280 sho=start[i];
281 pstar[i]=sho;
282 p2star[i]=sho;
283 pmin[i]=sho;
284 }
285
286 int jcount=this.konvge; // count of number of restarts still available
287
288 for (int i=0; i<np; ++i){
289 pp[i][nnp-1]=start[i];
290 }
291 yy[nnp-1]=this.functionValue(g, start);
292 for (int j=0; j<np; ++j){
293 start[j]=start[j]+step[j];
294
295 for (int i=0; i<np; ++i)pp[i][j]=start[i];
296 yy[j]=this.functionValue(g, start);
297 start[j]=start[j]-step[j];
298 }
299
300 // loop over allowed iterations
301 double ynewlo=0.0D; // current value lowest y
302 double ystar = 0.0D; // Nelder and Mead y*
303 double y2star = 0.0D; // Nelder and Mead y**
304 double ylo = 0.0D; // Nelder and Mead y(low)
305 double fMin; // function value at minimum
306 // variables used in calculating the variance of the simplex at a putative minimum
307 double curMin = 00D, sumnm = 0.0D, summnm = 0.0D, zn = 0.0D;
308 int ilo=0; // index of low apex
309 int ihi=0; // index of high apex
310 int ln=0; // counter for a check on low and high apices
311 boolean test = true; // test becomes false on reaching minimum
312
313 while(test){
314 // Determine h
315 ylo=yy[0];
316 ynewlo=ylo;
317 ilo=0;
318 ihi=0;
319 for (int i=1; i<nnp; ++i){
320 if (yy[i]<ylo){
321 ylo=yy[i];
322 ilo=i;
323 }
324 if (yy[i]>ynewlo){
325 ynewlo=yy[i];
326 ihi=i;
327 }
328 }
329 // Calculate pbar
330 for (int i=0; i<np; ++i){
331 zn=0.0D;
332 for (int j=0; j<nnp; ++j){
333 zn += pp[i][j];
334 }
335 zn -= pp[i][ihi];
336 pbar[i] = zn/np;
337 }
338
339 // Calculate p=(1+alpha).pbar-alpha.ph {Reflection}
340 for (int i=0; i<np; ++i)pstar[i]=(1.0 + this.rCoeff)*pbar[i]-this.rCoeff*pp[i][ihi];
341
342 // Calculate y*
343 ystar=this.functionValue(g, pstar);
344
345 ++this.nIter;
346
347 // check for y*<yi
348 if(ystar < ylo){
349 // Form p**=(1+gamma).p*-gamma.pbar {Extension}
350 for (int i=0; i<np; ++i)p2star[i]=pstar[i]*(1.0D + this.eCoeff)-this.eCoeff*pbar[i];
351 // Calculate y**
352 y2star=this.functionValue(g, p2star);
353 ++this.nIter;
354 if(y2star < ylo){
355 // Replace ph by p**
356 for (int i=0; i<np; ++i)pp[i][ihi] = p2star[i];
357 yy[ihi] = y2star;
358 }
359 else{
360 //Replace ph by p*
361 for (int i=0; i<np; ++i)pp[i][ihi]=pstar[i];
362 yy[ihi]=ystar;
363 }
364 }
365 else{
366 // Check y*>yi, i!=h
367 ln=0;
368 for (int i=0; i<nnp; ++i)if (i!=ihi && ystar > yy[i]) ++ln;
369 if (ln==np ){
370 // y*>= all yi; Check if y*>yh
371 if(ystar<=yy[ihi]){
372 // Replace ph by p*
373 for (int i=0; i<np; ++i)pp[i][ihi]=pstar[i];
374 yy[ihi]=ystar;
375 }
376 // Calculate p** =beta.ph+(1-beta)pbar {Contraction}
377 for (int i=0; i<np; ++i)p2star[i]=this.cCoeff*pp[i][ihi] + (1.0 - this.cCoeff)*pbar[i];
378 // Calculate y**
379 y2star=this.functionValue(g, p2star);
380 ++this.nIter;
381 // Check if y**>yh
382 if(y2star>yy[ihi]){
383 //Replace all pi by (pi+pl)/2
384
385 for (int j=0; j<nnp; ++j){
386 for (int i=0; i<np; ++i){
387 pp[i][j]=0.5*(pp[i][j] + pp[i][ilo]);
388 pmin[i]=pp[i][j];
389 }
390 yy[j]=this.functionValue(g, pmin);
391 }
392 this.nIter += nnp;
393 }
394 else{
395 // Replace ph by p**
396 for (int i=0; i<np; ++i)pp[i][ihi] = p2star[i];
397 yy[ihi] = y2star;
398 }
399 }
400 else{
401 // replace ph by p*
402 for (int i=0; i<np; ++i)pp[i][ihi]=pstar[i];
403 yy[ihi]=ystar;
404 }
405 }
406
407 // test for convergence
408 // calculte sd of simplex and minimum point
409 sumnm=0.0;
410 ynewlo=yy[0];
411 ilo=0;
412 for (int i=0; i<nnp; ++i){
413 sumnm += yy[i];
414 if(ynewlo>yy[i]){
415 ynewlo=yy[i];
416 ilo=i;
417 }
418 }
419 sumnm /= (double)(nnp);
420 summnm=0.0;
421 for (int i=0; i<nnp; ++i){
422 zn=yy[i]-sumnm;
423 summnm += zn*zn;
424 }
425 curMin=Math.sqrt(summnm/np);
426
427 // test simplex sd
428 switch(this.minTest){
429 case 0:
430 if(curMin<fTol)test=false;
431 break;
432 }
433 this.minimum=ynewlo;
434 if(!test){
435 // store parameter values
436 for (int i=0; i<np; ++i)pmin[i]=pp[i][ilo];
437 yy[nnp-1]=ynewlo;
438 // store simplex sd
439 this.simplexSd = curMin;
440 // test for restart
441 --jcount;
442 if(jcount>0){
443 test=true;
444 for (int j=0; j<np; ++j){
445 pmin[j]=pmin[j]+step[j];
446 for (int i=0; i<np; ++i)pp[i][j]=pmin[i];
447 yy[j]=this.functionValue(g, pmin);
448 pmin[j]=pmin[j]-step[j];
449 }
450 }
451 }
452
453 if(test && this.nIter>this.nMax){
454 if(!this.suppressNoConvergenceMessage){
455 System.out.println("Maximum iteration number reached, in Minimisation.simplex(...)");
456 System.out.println("without the convergence criterion being satisfied");
457 System.out.println("Current parameter estimates and function value returned");
458 }
459 this.convStatus = false;
460 // store current estimates
461 for (int i=0; i<np; ++i)pmin[i]=pp[i][ilo];
462 yy[nnp-1]=ynewlo;
463 test=false;
464 }
465 }
466
467 for (int i=0; i<np; ++i){
468 pmin[i] = pp[i][ilo];
469 paramValue[i] = pmin[i]/this.scale[i];
470 }
471 this.minimum=ynewlo;
472 this.kRestart=this.konvge-jcount;
473
474 }
475
476 // Nelder and Mead simplex
477 // Default maximum iterations
478 public void nelderMead(MinimisationFunction g, double[] start, double[] step, double fTol){
479 int nMaxx = this.nMax;
480 this.nelderMead(g, start, step, fTol, nMaxx);
481 }
482
483 public void nelderMead(MinimizationFunction g, double[] start, double[] step, double fTol){
484 int nMaxx = this.nMax;
485 this.nelderMead(g, start, step, fTol, nMaxx);
486 }
487
488 // Nelder and Mead simplex
489 // Default tolerance
490 public void nelderMead(MinimisationFunction g, double[] start, double[] step, int nMax){
491 double fToll = this.fTol;
492 this.nelderMead(g, start, step, fToll, nMax);
493 }
494
495 public void nelderMead(MinimizationFunction g, double[] start, double[] step, int nMax){
496 double fToll = this.fTol;
497 this.nelderMead(g, start, step, fToll, nMax);
498 }
499
500 // Nelder and Mead simplex
501 // Default tolerance
502 // Default maximum iterations
503 public void nelderMead(MinimisationFunction g, double[] start, double[] step){
504 double fToll = this.fTol;
505 int nMaxx = this.nMax;
506 this.nelderMead(g, start, step, fToll, nMaxx);
507 }
508
509 public void nelderMead(MinimizationFunction g, double[] start, double[] step){
510 double fToll = this.fTol;
511 int nMaxx = this.nMax;
512 this.nelderMead(g, start, step, fToll, nMaxx);
513 }
514
515 // Nelder and Mead simplex
516 // Default step option - all step[i] = dStep
517 public void nelderMead(MinimisationFunction g, double[] start, double fTol, int nMax){
518 int n=start.length;
519 double[] stepp = new double[n];
520 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
521 this.nelderMead(g, start, stepp, fTol, nMax);
522 }
523
524 public void nelderMead(MinimizationFunction g, double[] start, double fTol, int nMax){
525 int n=start.length;
526 double[] stepp = new double[n];
527 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
528 this.nelderMead(g, start, stepp, fTol, nMax);
529 }
530
531 // Nelder and Mead simplex
532 // Default maximum iterations
533 // Default step option - all step[i] = dStep
534 public void nelderMead(MinimisationFunction g, double[] start, double fTol){
535 int n=start.length;
536 int nMaxx = this.nMax;
537 double[] stepp = new double[n];
538 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
539 this.nelderMead(g, start, stepp, fTol, nMaxx);
540 }
541
542 public void nelderMead(MinimizationFunction g, double[] start, double fTol){
543 int n=start.length;
544 int nMaxx = this.nMax;
545 double[] stepp = new double[n];
546 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
547 this.nelderMead(g, start, stepp, fTol, nMaxx);
548 }
549
550 // Nelder and Mead simplex
551 // Default tolerance
552 // Default step option - all step[i] = dStep
553 public void nelderMead(MinimisationFunction g, double[] start, int nMax){
554 int n=start.length;
555 double fToll = this.fTol;
556 double[] stepp = new double[n];
557 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
558 this.nelderMead(g, start, stepp, fToll, nMax);
559 }
560
561 public void nelderMead(MinimizationFunction g, double[] start, int nMax){
562 int n=start.length;
563 double fToll = this.fTol;
564 double[] stepp = new double[n];
565 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
566 this.nelderMead(g, start, stepp, fToll, nMax);
567 }
568
569 // Nelder and Mead simplex
570 // Default tolerance
571 // Default maximum iterations
572 // Default step option - all step[i] = dStep
573 public void nelderMead(MinimisationFunction g, double[] start){
574 int n=start.length;
575 int nMaxx = this.nMax;
576 double fToll = this.fTol;
577 double[] stepp = new double[n];
578 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
579 this.nelderMead(g, start, stepp, fToll, nMaxx);
580 }
581
582 public void nelderMead(MinimizationFunction g, double[] start){
583 int n=start.length;
584 int nMaxx = this.nMax;
585 double fToll = this.fTol;
586 double[] stepp = new double[n];
587 for(int i=0; i<n;i++)stepp[i]=this.dStep*start[i];
588 this.nelderMead(g, start, stepp, fToll, nMaxx);
589 }
590
591
592 // Calculate the function value for minimisation
593 protected double functionValue(Object g, double[] x){
594 if(this.iseOption){
595 return functionValue((MinimisationFunction) g, x);
596 }
597 else{
598 return functionValue((MinimizationFunction) g, x);
599 }
600 }
601
602
603 // Calculate the function value for minimisation
604 protected double functionValue(MinimisationFunction g, double[] x){
605
606 double[] param = new double[this.nParam];
607 // rescale
608 for(int i=0; i<this.nParam; i++)param[i]=x[i]/scale[i];
609
610 boolean test = this.functionValueCommon(x, param);
611
612 if(test){
613 this.functValue = g.function(param);
614 this.lastFunctValNoCnstrnt = this.functValue;
615 }
616 return this.functValue;
617 }
618
619 // Calculate the function value for minimisation
620 protected double functionValue(MinimizationFunction g, double[] x){
621
622 double[] param = new double[this.nParam];
623 // rescale
624 for(int i=0; i<this.nParam; i++)param[i]=x[i]/scale[i];
625
626 boolean test = this.functionValueCommon(x, param);
627
628 if(test){
629 this.functValue = g.function(param);
630 this.lastFunctValNoCnstrnt = this.functValue;
631 }
632 return this.functValue;
633 }
634
635 // Common method for the functionValue(..) methods
636 protected boolean functionValueCommon(double[] x, double[] param){
637
638 // single parameter penalty functions
639 double tempFunctVal = this.lastFunctValNoCnstrnt;
640 boolean test=true;
641 if(this.penalty){
642 int k=0;
643 for(int i=0; i<this.nConstraints; i++){
644 k = this.penaltyParam[i];
645 switch(penaltyCheck[i]){
646 case -1: if(param[k]<constraints[i]){
647 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(constraints[i]-param[k]);
648 test=false;
649 }
650 break;
651 case 0: if(param[k]<constraints[i]*(1.0-this.constraintTolerance)){
652 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(constraints[i]*(1.0-this.constraintTolerance)-param[k]);
653 test=false;
654 }
655 if(param[k]>constraints[i]*(1.0+this.constraintTolerance)){
656 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(param[k]-constraints[i]*(1.0+this.constraintTolerance));
657 test=false;
658 }
659 break;
660 case 1: if(param[k]>constraints[i]){
661 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(param[k]-constraints[i]);
662 test=false;
663 }
664 break;
665 }
666 }
667 }
668
669 // multiple parameter penalty functions
670 if(this.sumPenalty){
671 int kk = 0;
672 double pSign = 0;
673 for(int i=0; i<this.nSumConstraints; i++){
674 double sumPenaltySum = 0.0D;
675 for(int j=0; j<this.sumPenaltyNumber[i]; j++){
676 kk = this.sumPenaltyParam[i][j];
677 pSign = this.sumPlusOrMinus[i][j];
678 sumPenaltySum += param[kk]*pSign;
679 }
680 switch(this.sumPenaltyCheck[i]){
681 case -1: if(sumPenaltySum<sumConstraints[i]){
682 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(sumConstraints[i]-sumPenaltySum);
683 test=false;
684 }
685 break;
686 case 0: if(sumPenaltySum<sumConstraints[i]*(1.0-this.constraintTolerance)){
687 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(sumConstraints[i]*(1.0-this.constraintTolerance)-sumPenaltySum);
688 test=false;
689 }
690 if(sumPenaltySum>sumConstraints[i]*(1.0+this.constraintTolerance)){
691 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(sumPenaltySum-sumConstraints[i]*(1.0+this.constraintTolerance));
692 test=false;
693 }
694 break;
695 case 1: if(sumPenaltySum>sumConstraints[i]){
696 this.functValue = tempFunctVal + this.penaltyWeight*Fmath.square(sumPenaltySum-sumConstraints[i]);
697 test=false;
698 }
699 break;
700 }
701 }
702 }
703 return test;
704 }
705
706
707 // add a single parameter constraint boundary for the minimisation
708 public void addConstraint(int paramIndex, int conDir, double constraint){
709 this.penalty=true;
710 // First element reserved for method number if other methods than 'cliff' are added later
711 if(this.penalties.isEmpty())this.penalties.add(new Integer(this.constraintMethod));
712
713 // add constraint
714 if(penalties.size()==1){
715 this.penalties.add(new Integer(1));
716 }
717 else{
718 int nPC = ((Integer)this.penalties.get(1)).intValue();
719 nPC++;
720 this.penalties.set(1, new Integer(nPC));
721 }
722 this.penalties.add(new Integer(paramIndex));
723 this.penalties.add(new Integer(conDir));
724 this.penalties.add(new Double(constraint));
725 if(paramIndex>this.maxConstraintIndex)this.maxConstraintIndex = paramIndex;
726
727 }
728
729 // add a multiple parameter constraint boundary for the non-linear regression
730 public void addConstraint(int[] paramIndices, int[] plusOrMinus, int conDir, double constraint){
731 ArrayMaths am = new ArrayMaths(plusOrMinus);
732 double[] dpom = am.getArray_as_double();
733 addConstraint(paramIndices, dpom, conDir, constraint);
734 }
735
736 // add a multiple parameter constraint boundary for the minimisation
737 public void addConstraint(int[] paramIndices, double[] plusOrMinus, int conDir, double constraint){
738 int nCon = paramIndices.length;
739 int nPorM = plusOrMinus.length;
740 if(nCon!=nPorM)throw new IllegalArgumentException("num of parameters, " + nCon + ", does not equal number of parameter signs, " + nPorM);
741 this.sumPenalty=true;
742 // First element reserved for method number if other methods than 'cliff' are added later
743 if(this.sumPenalties.isEmpty())this.sumPenalties.add(new Integer(this.constraintMethod));
744
745 // add constraint
746 if(sumPenalties.size()==1){
747 this.sumPenalties.add(new Integer(1));
748 }
749 else{
750 int nPC = ((Integer)this.sumPenalties.get(1)).intValue();
751 nPC++;
752 this.sumPenalties.set(1, new Integer(nPC));
753 }
754 this.sumPenalties.add(new Integer(nCon));
755 this.sumPenalties.add(paramIndices);
756 this.sumPenalties.add(plusOrMinus);
757 this.sumPenalties.add(new Integer(conDir));
758 this.sumPenalties.add(new Double(constraint));
759 ArrayMaths am = new ArrayMaths(paramIndices);
760 int maxI = am.getMaximum_as_int();
761 if(maxI>this.maxConstraintIndex)this.maxConstraintIndex = maxI;
762 }
763
764 // Set constraint method
765 public void setConstraintMethod(int conMeth){
766 this.constraintMethod = conMeth;
767 if(!this.penalties.isEmpty())this.penalties.set(0, new Integer(this.constraintMethod));
768 }
769
770 // remove all constraint boundaries for the minimisation
771 public void removeConstraints(){
772
773 // check if single parameter constraints already set
774 if(!this.penalties.isEmpty()){
775 int m=this.penalties.size();
776
777 // remove single parameter constraints
778 for(int i=m-1; i>=0; i--){
779 this.penalties.remove(i);
780 }
781 }
782 this.penalty = false;
783 this.nConstraints = 0;
784
785 // check if mutiple parameter constraints already set
786 if(!this.sumPenalties.isEmpty()){
787 int m=this.sumPenalties.size();
788
789 // remove multiple parameter constraints
790 for(int i=m-1; i>=0; i--){
791 this.sumPenalties.remove(i);
792 }
793 }
794 this.sumPenalty = false;
795 this.nSumConstraints = 0;
796 this.maxConstraintIndex = -1;
797 }
798
799 // Reset the tolerance used in a fixed value constraint
800 public void setConstraintTolerance(double tolerance){
801 this.constraintTolerance = tolerance;
802 }
803
804 // Print the results of the minimisation
805 // File name provided
806 // prec = truncation precision
807 public void print(String filename, int prec){
808 this.prec = prec;
809 this.print(filename);
810 }
811
812 // Print the results of the minimisation
813 // No file name provided
814 // prec = truncation precision
815 public void print(int prec){
816 this.prec = prec;
817 String filename="MinimisationOutput.txt";
818 this.print(filename);
819 }
820
821 // Print the results of the minimisation
822 // File name provided
823 // prec = truncation precision
824 public void print(String filename){
825
826 if(filename.indexOf('.')==-1)filename = filename+".txt";
827 FileOutput fout = new FileOutput(filename, 'n');
828 fout.dateAndTimeln(filename);
829 fout.println(" ");
830 fout.println("Simplex minimisation, using the method of Nelder and Mead,");
831 fout.println("of the function y = f(c[0], c[1], c[2] . . .)");
832 this.paraName = new String[this.nParam];
833 for(int i=0;i<this.nParam;i++)this.paraName[i]="c["+i+"]";
834
835 fout.println();
836 if(!this.convStatus){
837 fout.println("Convergence criterion was not satisfied");
838 fout.println("The following results are the current estimates on exiting the minimisation method");
839 fout.println();
840 }
841
842 fout.println("Value of parameters at the minimum");
843 fout.println(" ");
844
845 fout.printtab(" ", this.field);
846 fout.printtab("Value at", this.field);
847 fout.printtab("Initial", this.field);
848 fout.println("Initial");
849
850 fout.printtab(" ", this.field);
851 fout.printtab("mimium", this.field);
852 fout.printtab("estimate", this.field);
853 fout.println("step");
854
855 for(int i=0; i<this.nParam; i++){
856 fout.printtab(this.paraName[i], this.field);
857 fout.printtab(Fmath.truncate(paramValue[i],this.prec), this.field);
858 fout.printtab(Fmath.truncate(this.startH[i],this.prec), this.field);
859 fout.println(Fmath.truncate(this.step[i],this.prec));
860 }
861 fout.println();
862
863 fout.println(" ");
864
865 fout.printtab("Number of paramaters");
866 fout.println(this.nParam);
867 fout.printtab("Number of iterations taken");
868 fout.println(this.nIter);
869 fout.printtab("Maximum number of iterations allowed");
870 fout.println(this.nMax);
871 fout.printtab("Number of restarts taken");
872 fout.println(this.kRestart);
873 fout.printtab("Maximum number of restarts allowed");
874 fout.println(this.konvge);
875 fout.printtab("Standard deviation of the simplex at the minimum");
876 fout.println(Fmath.truncate(this.simplexSd, this.prec));
877 fout.printtab("Convergence tolerance");
878 fout.println(this.fTol);
879 switch(minTest){
880 case 0: if(this.convStatus){
881 fout.println("simplex sd < the tolerance");
882 }
883 else{
884 fout.println("NOTE!!! simplex sd > the tolerance");
885 }
886 break;
887 }
888 fout.println();
889 fout.println("End of file");
890 fout.close();
891 }
892
893 // Print the results of the minimisation
894 // No file name provided
895 public void print(){
896 String filename="MinimisationOutput.txt";
897 this.print(filename);
898 }
899
900 // Get the minimisation status
901 // true if convergence was achieved
902 // false if convergence not achieved before maximum number of iterations
903 // current values then returned
904 public boolean getConvStatus(){
905 return this.convStatus;
906 }
907
908 // Reset scaling factors (scaleOpt 0 and 1, see below for scaleOpt 2)
909 public void setScale(int n){
910 if(n<0 || n>1)throw new IllegalArgumentException("The argument must be 0 (no scaling) 1(initial estimates all scaled to unity) or the array of scaling factors");
911 this.scaleOpt=n;
912 }
913
914 // Reset scaling factors (scaleOpt 2, see above for scaleOpt 0 and 1)
915 public void setScale(double[] sc){
916 this.scale=sc;
917 this.scaleOpt=2;
918 }
919
920 // Get scaling factors
921 public double[] getScale(){
922 return this.scale;
923 }
924
925 // Reset the minimisation convergence test option
926 public void setMinTest(int n){
927 if(n<0 || n>1)throw new IllegalArgumentException("minTest must be 0 or 1");
928 this.minTest=n;
929 }
930
931 // Get the minimisation convergence test option
932 public int getMinTest(){
933 return this.minTest;
934 }
935
936 // Get the simplex sd at the minimum
937 public double getSimplexSd(){
938 return this.simplexSd;
939 }
940
941 // Get the values of the parameters at the minimum
942 public double[] getParamValues(){
943 return this.paramValue;
944 }
945
946 // Get the function value at minimum
947 public double getMinimum(){
948 return this.minimum;
949 }
950
951 // Get the number of iterations in Nelder and Mead
952 public int getNiter(){
953 return this.nIter;
954 }
955
956
957 // Set the maximum number of iterations allowed in Nelder and Mead
958 public void setNmax(int nmax){
959 this.nMax = nmax;
960 }
961
962 // Get the maximum number of iterations allowed in Nelder and Mead
963 public int getNmax(){
964 return this.nMax;
965 }
966
967 // Get the number of restarts in Nelder and Mead
968 public int getNrestarts(){
969 return this.kRestart;
970 }
971
972 // Set the maximum number of restarts allowed in Nelder and Mead
973 public void setNrestartsMax(int nrs){
974 this.konvge = nrs;
975 }
976
977 // Get the maximum number of restarts allowed in Nelder amd Mead
978 public int getNrestartsMax(){
979 return this.konvge;
980 }
981
982 // Reset the Nelder and Mead reflection coefficient [alpha]
983 public void setNMreflect(double refl){
984 this.rCoeff = refl;
985 }
986
987 // Get the Nelder and Mead reflection coefficient [alpha]
988 public double getNMreflect(){
989 return this.rCoeff;
990 }
991
992 // Reset the Nelder and Mead extension coefficient [beta]
993 public void setNMextend(double ext){
994 this.eCoeff = ext;
995 }
996 // Get the Nelder and Mead extension coefficient [beta]
997 public double getNMextend(){
998 return this.eCoeff;
999 }
1000
1001 // Reset the Nelder and Mead contraction coefficient [gamma]
1002 public void setNMcontract(double con){
1003 this.cCoeff = con;
1004 }
1005
1006 // Get the Nelder and Mead contraction coefficient [gamma]
1007 public double getNMcontract(){
1008 return cCoeff;
1009 }
1010
1011 // Set the minimisation tolerance
1012 public void setTolerance(double tol){
1013 this.fTol = tol;
1014 }
1015
1016
1017 // Get the minimisation tolerance
1018 public double getTolerance(){
1019 return this.fTol;
1020 }
1021
1022}
Note: See TracBrowser for help on using the repository browser.