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