source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/interpolation/CubicSpline.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: 26.7 KB
Line 
1/**********************************************************
2*
3* Class CubicSpline
4*
5* Class for performing an interpolation using a cubic spline
6* setTabulatedArrays and interpolate adapted, with modification to
7* an object-oriented approach, from Numerical Recipes in C (http://www.nr.com/)
8*
9*
10* WRITTEN BY: Dr Michael Thomas Flanagan
11*
12* DATE: May 2002
13* UPDATE: 29 April 2005, 17 February 2006, 21 September 2006, 4 December 2007
14* 24 March 2008 (Thanks to Peter Neuhaus, Florida Institute for Human and Machine Cognition)
15* 21 September 2008
16* 14 January 2009 - point deletion and check for 3 points reordered (Thanks to Jan Sacha, Vrije Universiteit Amsterdam)
17* 31 October 2009, 11-14 January 2010
18* 31 August 2010 - overriding natural spline error, introduced in earlier update, corrected (Thanks to Dagmara Oszkiewicz, University of Helsinki)
19* 2 February 2011
20*
21* DOCUMENTATION:
22* See Michael Thomas Flanagan's Java library on-line web page:
23* http://www.ee.ucl.ac.uk/~mflanaga/java/CubicSpline.html
24* http://www.ee.ucl.ac.uk/~mflanaga/java/
25*
26* Copyright (c) 2002 - 2010 Michael Thomas Flanagan
27*
28* PERMISSION TO COPY:
29*
30* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
31* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
32* and associated documentation or publications.
33*
34* Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice,
35* this list of conditions and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
36*
37* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
38* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission
39* from the Michael Thomas Flanagan:
40*
41* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
42* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
43* or its derivatives.
44*
45***************************************************************************************/
46
47
48package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
49
50import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
51import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
52
53public class CubicSpline{
54
55 private int nPoints = 0; // no. of tabulated points
56 private int nPointsOriginal = 0; // no. of tabulated points after any deletions of identical points
57 private double[] y = null; // y=f(x) tabulated function
58 private double[] x = null; // x in tabulated function f(x)
59 private double yy = Double.NaN; // interpolated value of y
60 private double dydx = Double.NaN; // interpolated value of the first derivative, dy/dx
61 private int[]newAndOldIndices; // record of indices on ordering x into ascending order
62 private double xMin = Double.NaN; // minimum x value
63 private double xMax = Double.NaN; // maximum x value
64 private double range = Double.NaN; // xMax - xMin
65 private double[] d2ydx2 = null; // second derivatives of y
66 private double yp1 = Double.NaN; // second derivative at point one
67 // default value = NaN (natural spline)
68 private double ypn = Double.NaN; // second derivative at point n
69 // default value = NaN (natural spline)
70 private boolean derivCalculated = false; // = true when the derivatives have been calculated
71
72 private boolean checkPoints = false; // = true when points checked for identical values
73 private static boolean supress = false; // if true: warning messages supressed
74
75 private static boolean averageIdenticalAbscissae = false; // if true: the the ordinate values for identical abscissae are averaged
76 // if false: the abscissae values are separated by 0.001 of the total abscissae range;
77 private static double potentialRoundingError = 5e-15; // potential rounding error used in checking wheter a value lies within the interpolation bounds
78 private static boolean roundingCheck = true; // = true: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
79
80
81 // Constructors
82 // Constructor with data arrays initialised to arrays x and y
83 public CubicSpline(double[] x, double[] y){
84 this.nPoints=x.length;
85 this.nPointsOriginal = this.nPoints;
86 if(this.nPoints!=y.length)throw new IllegalArgumentException("Arrays x and y are of different length "+ this.nPoints + " " + y.length);
87 if(this.nPoints<3)throw new IllegalArgumentException("A minimum of three data points is needed");
88 this.x = new double[nPoints];
89 this.y = new double[nPoints];
90 this.d2ydx2 = new double[nPoints];
91 for(int i=0; i<this.nPoints; i++){
92 this.x[i]=x[i];
93 this.y[i]=y[i];
94 }
95 this.orderPoints();
96 this.checkForIdenticalPoints();
97 this.calcDeriv();
98
99 }
100
101 // Constructor with data arrays initialised to zero
102 // Primarily for use by BiCubicSpline
103 public CubicSpline(int nPoints){
104 this.nPoints=nPoints;
105 this.nPointsOriginal = this.nPoints;
106 if(this.nPoints<3)throw new IllegalArgumentException("A minimum of three data points is needed");
107 this.x = new double[nPoints];
108 this.y = new double[nPoints];
109 this.d2ydx2 = new double[nPoints];
110 }
111
112 // METHODS
113 // Reset rounding error check option
114 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
115 // This method causes this check to be ignored and an exception to be thrown if any poit lies outside the interpolation bounds
116 public static void noRoundingErrorCheck(){
117 CubicSpline.roundingCheck = false;
118 }
119
120 // Reset potential rounding error value
121 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
122 // The default value for the potential rounding error is 5e-15*times the 10^exponent of the value outside the bounds
123 // This method allows the 5e-15 to be reset
124 public static void potentialRoundingError(double potentialRoundingError){
125 CubicSpline.potentialRoundingError = potentialRoundingError;
126 }
127
128 // Resets the x y data arrays - primarily for use in BiCubicSpline
129 public void resetData(double[] x, double[] y){
130 this.nPoints = this.nPointsOriginal;
131 if(x.length!=y.length)throw new IllegalArgumentException("Arrays x and y are of different length");
132 if(this.nPoints!=x.length)throw new IllegalArgumentException("Original array length not matched by new array length");
133
134 for(int i=0; i<this.nPoints; i++){
135 this.x[i]=x[i];
136 this.y[i]=y[i];
137 }
138 this.orderPoints();
139 this.checkForIdenticalPoints();
140 }
141
142 // Reset the default handing of identical abscissae with different ordinates
143 // from the default option of separating the two relevant abscissae by 0.001 of the range
144 // to averaging the relevant ordinates
145 public static void averageIdenticalAbscissae(){
146 CubicSpline.averageIdenticalAbscissae = true;
147 }
148
149 // Sort points into an ascending abscissa order
150 public void orderPoints(){
151 double[] dummy = new double[nPoints];
152 this.newAndOldIndices = new int[nPoints];
153 // Sort x into ascending order storing indices changes
154 Fmath.selectionSort(this.x, dummy, this.newAndOldIndices);
155 // Sort x into ascending order and make y match the new order storing both new x and new y
156 Fmath.selectionSort(this.x, this.y, this.x, this.y);
157
158 // Minimum and maximum values and range
159 this.xMin = Fmath.minimum(this.x);
160 this.xMax = Fmath.maximum(this.x);
161 range = xMax - xMin;
162 }
163
164 // get the maximum value
165 public double getXmax(){
166 return this.xMax;
167 }
168
169 // get the minimum value
170 public double getXmin(){
171 return this.xMin;
172 }
173
174 // get the limits of x
175 public double[] getLimits(){
176 double[] limits = {this.xMin, this.xMax};
177 return limits;
178 }
179
180 // print to screen the limis of x
181 public void displayLimits(){
182 System.out.println("\nThe limits of the abscissae (x-values) are " + this.xMin + " and " + this.xMax +"\n");
183 }
184
185
186 // Checks for and removes all but one of identical points
187 // Checks and appropriately handles identical abscissae with differing ordinates
188 public void checkForIdenticalPoints(){
189 int nP = this.nPoints;
190 boolean test1 = true;
191 int ii = 0;
192 while(test1){
193 boolean test2 = true;
194 int jj = ii+1;
195 while(test2){
196 if(this.x[ii]==this.x[jj]){
197 if(this.y[ii]==this.y[jj]){
198 if(!CubicSpline.supress){
199 System.out.print("CubicSpline: Two identical points, " + this.x[ii] + ", " + this.y[ii]);
200 System.out.println(", in data array at indices " + this.newAndOldIndices[ii] + " and " + this.newAndOldIndices[jj] + ", latter point removed");
201 }
202 double[] xx = new double[this.nPoints-1];
203 double[] yy = new double[this.nPoints-1];
204 int[] naoi = new int[this.nPoints-1];
205 for(int i=0; i<jj; i++){
206 xx[i] = this.x[i];
207 yy[i] = this.y[i];
208 naoi[i] = this.newAndOldIndices[i];
209 }
210 for(int i=jj; i<this.nPoints-1; i++){
211 xx[i] = this.x[i+1];
212 yy[i] = this.y[i+1];
213 naoi[i] = this.newAndOldIndices[i+1];
214 }
215 this.nPoints--;
216 this.x = Conv.copy(xx);
217 this.y = Conv.copy(yy);
218 this.newAndOldIndices = Conv.copy(naoi);
219 }
220 else{
221 if(CubicSpline.averageIdenticalAbscissae==true){
222 if(!CubicSpline.supress){
223 System.out.print("CubicSpline: Two identical points on the absicca (x-axis) with different ordinate (y-axis) values, " + x[ii] + ": " + y[ii] + ", " + y[jj]);
224 System.out.println(", average of the ordinates taken");
225 }
226 this.y[ii] = (this.y[ii] + this.y[jj])/2.0D;
227 double[] xx = new double[this.nPoints-1];
228 double[] yy = new double[this.nPoints-1];
229 int[] naoi = new int[this.nPoints-1];
230 for(int i=0; i<jj; i++){
231 xx[i] = this.x[i];
232 yy[i] = this.y[i];
233 naoi[i] = this.newAndOldIndices[i];
234 }
235 for(int i=jj; i<this.nPoints-1; i++){
236 xx[i] = this.x[i+1];
237 yy[i] = this.y[i+1];
238 naoi[i] = this.newAndOldIndices[i+1];
239 }
240 this.nPoints--;
241 this.x = Conv.copy(xx);
242 this.y = Conv.copy(yy);
243 this.newAndOldIndices = Conv.copy(naoi);
244 }
245 else{
246 double sepn = range*0.0005D;
247 if(!CubicSpline.supress)System.out.print("CubicSpline: Two identical points on the absicca (x-axis) with different ordinate (y-axis) values, " + x[ii] + ": " + y[ii] + ", " + y[jj]);
248 boolean check = false;
249 if(ii==0){
250 if(x[2]-x[1]<=sepn)sepn = (x[2]-x[1])/2.0D;
251 if(this.y[0]>this.y[1]){
252 if(this.y[1]>this.y[2]){
253 check = stay(ii, jj, sepn);
254 }
255 else{
256 check = swap(ii, jj, sepn);
257 }
258 }
259 else{
260 if(this.y[2]<=this.y[1]){
261 check = swap(ii, jj, sepn);
262 }
263 else{
264 check = stay(ii, jj, sepn);
265 }
266 }
267 }
268 if(jj==this.nPoints-1){
269 if(x[nP-2]-x[nP-3]<=sepn)sepn = (x[nP-2]-x[nP-3])/2.0D;
270 if(this.y[ii]<=this.y[jj]){
271 if(this.y[ii-1]<=this.y[ii]){
272 check = stay(ii, jj, sepn);
273 }
274 else{
275 check = swap(ii, jj, sepn);
276 }
277 }
278 else{
279 if(this.y[ii-1]<=this.y[ii]){
280 check = swap(ii, jj, sepn);
281 }
282 else{
283 check = stay(ii, jj, sepn);
284 }
285 }
286 }
287 if(ii!=0 && jj!=this.nPoints-1){
288 if(x[ii]-x[ii-1]<=sepn)sepn = (x[ii]-x[ii-1])/2;
289 if(x[jj+1]-x[jj]<=sepn)sepn = (x[jj+1]-x[jj])/2;
290 if(this.y[ii]>this.y[ii-1]){
291 if(this.y[jj]>this.y[ii]){
292 if(this.y[jj]>this.y[jj+1]){
293 if(this.y[ii-1]<=this.y[jj+1]){
294 check = stay(ii, jj, sepn);
295 }
296 else{
297 check = swap(ii, jj, sepn);
298 }
299 }
300 else{
301 check = stay(ii, jj, sepn);
302 }
303 }
304 else{
305 if(this.y[jj+1]>this.y[jj]){
306 if(this.y[jj+1]>this.y[ii-1] && this.y[jj+1]>this.y[ii-1]){
307 check = stay(ii, jj, sepn);
308 }
309 }
310 else{
311 check = swap(ii, jj, sepn);
312 }
313 }
314 }
315 else{
316 if(this.y[jj]>this.y[ii]){
317 if(this.y[jj+1]>this.y[jj]){
318 check = stay(ii, jj, sepn);
319 }
320 }
321 else{
322 if(this.y[jj+1]>this.y[ii-1]){
323 check = stay(ii, jj, sepn);
324 }
325 else{
326 check = swap(ii, jj, sepn);
327 }
328 }
329 }
330 }
331
332 if(check==false){
333 check = stay(ii, jj, sepn);
334 }
335 if(!CubicSpline.supress)System.out.println(", the two abscissae have been separated by a distance " + sepn);
336 jj++;
337 }
338 }
339 if((this.nPoints-1)==ii)test2 = false;
340 }
341 else{
342 jj++;
343 }
344 if(jj>=this.nPoints)test2 = false;
345 }
346 ii++;
347 if(ii>=this.nPoints-1)test1 = false;
348 }
349 if(this.nPoints<3)throw new IllegalArgumentException("Removal of duplicate points has reduced the number of points to less than the required minimum of three data points");
350
351 this.checkPoints = true;
352 }
353
354 // Swap method for checkForIdenticalPoints procedure
355 private boolean swap(int ii, int jj, double sepn){
356 this.x[ii] += sepn;
357 this.x[jj] -= sepn;
358 double hold = this.x[ii];
359 this.x[ii] = this.x[jj];
360 this.x[jj] = hold;
361 hold = this.y[ii];
362 this.y[ii] = this.y[jj];
363 this.y[jj] = hold;
364 return true;
365 }
366
367 // Stay method for checkForIdenticalPoints procedure
368 private boolean stay(int ii, int jj, double sepn){
369 this.x[ii] -= sepn;
370 this.x[jj] += sepn;
371 return true;
372 }
373
374 // Supress warning messages in the identifiaction of duplicate points
375 public static void supress(){
376 CubicSpline.supress = true;
377 }
378
379 // Unsupress warning messages in the identifiaction of duplicate points
380 public static void unsupress(){
381 CubicSpline.supress = false;
382 }
383
384 // Returns a new CubicSpline setting array lengths to n and all array values to zero with natural spline default
385 // Primarily for use in BiCubicSpline
386 public static CubicSpline zero(int n){
387 if(n<3)throw new IllegalArgumentException("A minimum of three data points is needed");
388 CubicSpline aa = new CubicSpline(n);
389 return aa;
390 }
391
392 // Create a one dimensional array of cubic spline objects of length n each of array length m
393 // Primarily for use in BiCubicSpline
394 public static CubicSpline[] oneDarray(int n, int m){
395 if(m<3)throw new IllegalArgumentException("A minimum of three data points is needed");
396 CubicSpline[] a =new CubicSpline[n];
397 for(int i=0; i<n; i++){
398 a[i]=CubicSpline.zero(m);
399 }
400 return a;
401 }
402
403 // Enters the first derivatives of the cubic spline at
404 // the first and last point of the tabulated data
405 // Overrides a natural spline
406 public void setDerivLimits(double yp1, double ypn){
407 this.yp1=yp1;
408 this.ypn=ypn;
409 this.calcDeriv();
410 }
411
412 // Resets a natural spline
413 // Use above - this kept for backward compatibility
414 public void setDerivLimits(){
415 this.yp1=Double.NaN;
416 this.ypn=Double.NaN;
417 this.calcDeriv();
418 }
419
420 // Enters the first derivatives of the cubic spline at
421 // the first and last point of the tabulated data
422 // Overrides a natural spline
423 // Use setDerivLimits(double yp1, double ypn) - this kept for backward compatibility
424 public void setDeriv(double yp1, double ypn){
425 this.yp1=yp1;
426 this.ypn=ypn;
427 this.calcDeriv();
428 }
429
430 // Returns the internal array of second derivatives
431 public double[] getDeriv(){
432 if(!this.derivCalculated)this.calcDeriv();
433 return this.d2ydx2;
434 }
435
436 // Sets the internal array of second derivatives
437 // Used primarily with BiCubicSpline
438 public void setDeriv(double[] deriv){
439 this.d2ydx2 = deriv;
440 this.derivCalculated = true;
441 }
442
443 // Calculates the second derivatives of the tabulated function
444 // for use by the cubic spline interpolation method (.interpolate)
445 // This method follows the procedure in Numerical Methods C language procedure for calculating second derivatives
446 public void calcDeriv(){
447 double p=0.0D,qn=0.0D,sig=0.0D,un=0.0D;
448 double[] u = new double[nPoints];
449
450 if (Double.isNaN(this.yp1)){
451 d2ydx2[0]=u[0]=0.0;
452 }
453 else{
454 this.d2ydx2[0] = -0.5;
455 u[0]=(3.0/(this.x[1]-this.x[0]))*((this.y[1]-this.y[0])/(this.x[1]-this.x[0])-this.yp1);
456 }
457
458 for(int i=1;i<=this.nPoints-2;i++){
459 sig=(this.x[i]-this.x[i-1])/(this.x[i+1]-this.x[i-1]);
460 p=sig*this.d2ydx2[i-1]+2.0;
461 this.d2ydx2[i]=(sig-1.0)/p;
462 u[i]=(this.y[i+1]-this.y[i])/(this.x[i+1]-this.x[i]) - (this.y[i]-this.y[i-1])/(this.x[i]-this.x[i-1]);
463 u[i]=(6.0*u[i]/(this.x[i+1]-this.x[i-1])-sig*u[i-1])/p;
464 }
465
466 if (Double.isNaN(this.ypn)){
467 qn=un=0.0;
468 }
469 else{
470 qn=0.5;
471 un=(3.0/(this.x[nPoints-1]-this.x[this.nPoints-2]))*(this.ypn-(this.y[this.nPoints-1]-this.y[this.nPoints-2])/(this.x[this.nPoints-1]-x[this.nPoints-2]));
472 }
473
474 this.d2ydx2[this.nPoints-1]=(un-qn*u[this.nPoints-2])/(qn*this.d2ydx2[this.nPoints-2]+1.0);
475 for(int k=this.nPoints-2;k>=0;k--){
476 this.d2ydx2[k]=this.d2ydx2[k]*this.d2ydx2[k+1]+u[k];
477 }
478 this.derivCalculated = true;
479 }
480
481 // INTERPOLATE
482 // Returns an interpolated value of y for a value of x from a tabulated function y=f(x)
483 // after the data has been entered via a constructor.
484 // The derivatives are calculated, bt calcDeriv(), on the first call to this method ands are
485 // then stored for use on all subsequent calls
486 public double interpolate(double xx){
487
488 // Check for violation of interpolation bounds
489 if (xx<this.x[0]){
490 // if violation is less than potntial rounding error - amend to lie with bounds
491 if(CubicSpline.roundingCheck && Math.abs(x[0]-xx)<=Math.pow(10, Math.floor(Math.log10(Math.abs(this.x[0]))))*CubicSpline.potentialRoundingError){
492 xx = x[0];
493 }
494 else{
495 throw new IllegalArgumentException("x ("+xx+") is outside the range of data points ("+x[0]+" to "+x[this.nPoints-1] + ")");
496 }
497 }
498 if (xx>this.x[this.nPoints-1]){
499 if(CubicSpline.roundingCheck && Math.abs(xx-this.x[this.nPoints-1])<=Math.pow(10, Math.floor(Math.log10(Math.abs(this.x[this.nPoints-1]))))*CubicSpline.potentialRoundingError){
500 xx = this.x[this.nPoints-1];
501 }
502 else{
503 throw new IllegalArgumentException("x ("+xx+") is outside the range of data points ("+x[0]+" to "+x[this.nPoints-1] + ")");
504 }
505 }
506
507 double h=0.0D,b=0.0D,a=0.0D, yy=0.0D;
508 int k=0;
509 int klo=0;
510 int khi=this.nPoints-1;
511 while (khi-klo > 1){
512 k=(khi+klo) >> 1;
513 if(this.x[k] > xx){
514 khi=k;
515 }
516 else{
517 klo=k;
518 }
519 }
520 h=this.x[khi]-this.x[klo];
521
522 if (h == 0.0){
523 throw new IllegalArgumentException("Two values of x are identical: point "+klo+ " ("+this.x[klo]+") and point "+khi+ " ("+this.x[khi]+")" );
524 }
525 else{
526 a=(this.x[khi]-xx)/h;
527 b=(xx-this.x[klo])/h;
528 this.yy = a*this.y[klo]+b*this.y[khi]+((a*a*a-a)*this.d2ydx2[klo]+(b*b*b-b)*this.d2ydx2[khi])*(h*h)/6.0;
529 this.dydx = (y[khi] - y[klo])/h - ((3*a*a - 1.0)*this.d2ydx2[klo] - (3*b*b - 1.0)*this.d2ydx2[khi])*h/6.0;
530 }
531 return this.yy;
532 }
533
534
535 // Returns an interpolated value of y and of the first derivative dy/dx for a value of x from a tabulated function y=f(x)
536 // after the data has been entered via a constructor.
537 public double[] interpolate_for_y_and_dydx(double xx){
538 this.interpolate(xx);
539 double[] ret = {this.yy, this.dydx};
540 return ret;
541 }
542
543
544 // Returns an interpolated value of y for a value of x (xx) from a tabulated function y=f(x)
545 // after the derivatives (deriv) have been calculated independently of calcDeriv().
546 public static double interpolate(double xx, double[] x, double[] y, double[] deriv){
547
548 if(((x.length != y.length) || (x.length != deriv.length)) || (y.length != deriv.length)){
549 throw new IllegalArgumentException("array lengths are not all equal");
550 }
551 int n = x.length;
552 double h=0.0D, b=0.0D, a=0.0D, yy = 0.0D;
553
554 int k=0;
555 int klo=0;
556 int khi=n-1;
557 while(khi-klo > 1){
558 k=(khi+klo) >> 1;
559 if(x[k] > xx){
560 khi=k;
561 }
562 else{
563 klo=k;
564 }
565 }
566 h=x[khi]-x[klo];
567
568 if (h == 0.0){
569 throw new IllegalArgumentException("Two values of x are identical");
570 }
571 else{
572 a=(x[khi]-xx)/h;
573 b=(xx-x[klo])/h;
574 yy=a*y[klo]+b*y[khi]+((a*a*a-a)*deriv[klo]+(b*b*b-b)*deriv[khi])*(h*h)/6.0;
575 }
576 return yy;
577 }
578
579}
Note: See TracBrowser for help on using the repository browser.