source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/interpolation/PolyCubicSpline.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: 15.6 KB
Line 
1/**********************************************************
2*
3* PolyCubicSpline.java
4*
5* Class for performing an interpolation on the tabulated
6* function y = f(x1,x2, x3 .... xn) using a natural cubic splines
7* Assumes second derivatives at end points = 0 (natural spines)
8*
9* WRITTEN BY: Dr Michael Thomas Flanagan
10*
11* DATE: 15 February 2006,
12* UPDATES: 9 June 2007, 27 July 2007, 4 December 2007, 21 September 2008, 12 October 2009, 31 October 2009
13*
14* DOCUMENTATION:
15* See Michael Thomas Flanagan's Java library on-line web page:
16* http://www.ee.ucl.ac.uk/~mflanaga/java/PolyCubicSpline.html
17* http://www.ee.ucl.ac.uk/~mflanaga/java/
18*
19* Copyright (c) 2006 - 2009 Michael Thomas Flanagan
20*
21* PERMISSION TO COPY:
22*
23* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
24* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
25* and associated documentation or publications.
26*
27* Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice, this list of conditions
28* and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
29*
30* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
31* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission from the Michael Thomas Flanagan:
32*
33* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
34* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
35* or its derivatives.
36*
37***************************************************************************************/
38
39
40package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
41
42import java.lang.reflect.Array;
43
44import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
45
46public class PolyCubicSpline{
47
48 private int nDimensions = 0; // number of the dimensions of the tabulated points array, y=f(x1,x2,x3 . . xn), i.e. n
49 private Object fOfX = null; // tabulated values of y = f(x1,x2,x3 . . fn)
50 // as a multidimensional array of double [x1 length][x2 length] ... [xn length]
51 private Object internalDeriv = null; // Object to store second derivatives
52 // as a multidimensional array of Objects the innermost three layers being double arrays
53 private Object xArrays = null; // The variable arrays x1, x2, x3 . . . xn
54 // packed as an Object as a multidimensional array of double [][]
55 // where xArrays[0] = array of x1 values, xArrays[1] = array of x2 values etc
56 private Object method = null; // interpolation method
57 private double[][] xArray = null; // The variable arrays x1, x2, x3 . . . xn
58 private double[] csArray = null; // array for final cubic spline interpolation
59 private PolyCubicSpline[] pcs = null; // array of PolyCubicSplines for use with recursive step
60 private int dimOne = 0; // xArray dimension in a recursive step
61
62 private double yValue = 0.0D; // returned interpolated value
63 private double[] xMin = null; // minimum values of the x arrays
64 private double[] xMax = null; // maximum values of the x arrays
65 private boolean calculationDone = false; // = true when derivatives calculated
66 private boolean averageIdenticalAbscissae = false; // if true: the the ordinate values for identical abscissae are averaged
67 // If false: the abscissae values are separated by 0.001 of the total abscissae range;
68 private static double potentialRoundingError = 5e-15; // potential rounding error used in checking wheter a value lies within the interpolation bounds (static value)
69 private static boolean roundingCheck = true; // = true: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit (static value)
70
71
72 // Constructor
73 public PolyCubicSpline(Object xArrays, Object fOfX){
74
75 this.fOfX = Fmath.copyObject(fOfX);
76 this.xArrays = Fmath.copyObject(xArrays);
77
78 // Calculate fOfX array dimension number
79 Object internalArrays = Fmath.copyObject(fOfX);
80 this.nDimensions = 1;
81 while(!((internalArrays = Array.get(internalArrays, 0)) instanceof Double))this.nDimensions++;
82
83 // Repack xArrays as 2 dimensional array if entered a single dimensioned array for a simple cubic spline
84 if(this.xArrays instanceof double[] && this.nDimensions == 1){
85 double[][] xArraysTemp = new double[1][];
86 xArraysTemp[0] = (double[])this.xArrays;
87 this.xArrays = (Object)xArraysTemp;
88 }
89 else{
90 if(!(this.xArrays instanceof double[][]))throw new IllegalArgumentException("xArrays should be a two dimensional array of doubles");
91 }
92
93 // x -arrays and their limits
94 this.xArray = (double[][])this.xArrays;
95 this.limits();
96
97 // Select interpolation method
98 switch(this.nDimensions){
99 case 0: throw new IllegalArgumentException("data array must have at least one dimension");
100 case 1: // If fOfX is one dimensional perform simple cubic spline
101 CubicSpline cs = new CubicSpline(this.xArray[0], (double[])this.fOfX);
102 if(this.averageIdenticalAbscissae)cs.averageIdenticalAbscissae();
103 this.internalDeriv = (Object)(cs.getDeriv());
104 this.method = (Object)cs;
105 this.calculationDone = true;
106 break;
107 case 2: // If fOfX is two dimensional perform bicubic spline
108 BiCubicSpline bcs = new BiCubicSpline(this.xArray[0], this.xArray[1], (double[][])this.fOfX);
109 if(this.averageIdenticalAbscissae)bcs.averageIdenticalAbscissae();
110 this.internalDeriv = (Object)(bcs.getDeriv());
111 this.method = (Object)bcs;
112 this.calculationDone = true;
113 break;
114 case 3: // If fOfX is three dimensional perform tricubic spline
115 TriCubicSpline tcs = new TriCubicSpline(xArray[0], xArray[1], xArray[2], (double[][][])this.fOfX);
116 if(this.averageIdenticalAbscissae)tcs.averageIdenticalAbscissae();
117 this.internalDeriv = (Object)(tcs.getDeriv());
118 this.method = (Object)tcs;
119 this.calculationDone = true;
120 break;
121 case 4: // If fOfX is four dimensional perform quadricubic spline
122 QuadriCubicSpline qcs = new QuadriCubicSpline(xArray[0], xArray[1], xArray[2], xArray[3], (double[][][][])this.fOfX);
123 if(this.averageIdenticalAbscissae)qcs.averageIdenticalAbscissae();
124 this.internalDeriv = (Object)(qcs.getDeriv());
125 this.method = (Object)qcs;
126 this.calculationDone = true;
127 break;
128 default: // If fOfX is greater than four dimensional, recursively call PolyCubicSpline
129 // with, as arguments, the n1 fOfX sub-arrays, each of (number of dimensions - 1) dimensions,
130 // where n1 is the number of x1 variables.
131 Object obj = fOfX;
132 this.dimOne = Array.getLength(obj);
133 this.csArray = new double [dimOne];
134 double[][] newXarrays= new double[this.nDimensions-1][];
135 for(int i=0; i<this.nDimensions-1; i++){
136 newXarrays[i] = xArray[i+1];
137 }
138
139 Object[] objDeriv = new Object[dimOne];
140 if(calculationDone)objDeriv = (Object[])this.internalDeriv;
141 this.pcs = new PolyCubicSpline[dimOne];
142 for(int i=0; i<dimOne; i++){
143 Object objT = (Object)Array.get(obj, i);
144 this.pcs[i] = new PolyCubicSpline(newXarrays, objT);
145 if(this.averageIdenticalAbscissae)pcs[i].averageIdenticalAbscissae();
146 if(this.calculationDone)pcs[i].setDeriv(objDeriv[i]);
147 if(!this.calculationDone)objDeriv[i] = pcs[i].getDeriv();
148 }
149 this.internalDeriv = (Object)objDeriv;
150 this.calculationDone = true;
151 }
152
153 }
154
155 // Reset rounding error check option
156 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
157 // This method causes this check to be ignored and an exception to be thrown if any point lies outside the interpolation bounds
158 public static void noRoundingErrorCheck(){
159 PolyCubicSpline.roundingCheck = false;
160 QuadriCubicSpline.noRoundingErrorCheck();
161 TriCubicSpline.noRoundingErrorCheck();
162 BiCubicSpline.noRoundingErrorCheck();
163 CubicSpline.noRoundingErrorCheck();
164 }
165
166 // Reset potential rounding error value
167 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
168 // The default value for the potential rounding error is 5e-15*times the 10^exponent of the value outside the bounds
169 // This method allows the 5e-15 to be reset
170 public static void potentialRoundingError(double potentialRoundingError){
171 PolyCubicSpline.potentialRoundingError = potentialRoundingError;
172 QuadriCubicSpline.potentialRoundingError(potentialRoundingError);
173 TriCubicSpline.potentialRoundingError(potentialRoundingError);
174 BiCubicSpline.potentialRoundingError(potentialRoundingError);
175 CubicSpline.potentialRoundingError(potentialRoundingError);
176 }
177
178 // Limits to x arrays
179 private void limits(){
180 this.xMin = new double[this.nDimensions];
181 this.xMax = new double[this.nDimensions];
182 for(int i=0; i<this.nDimensions; i++){
183 this.xMin[i] = Fmath.minimum(xArray[i]);
184 this.xMax[i] = Fmath.maximum(xArray[i]);
185 }
186 }
187
188 // Get minimum limits
189 public double[] getXmin(){
190 return this.xMin;
191 }
192
193 // Get maximum limits
194 public double[] getXmax(){
195 return this.xMax;
196 }
197
198 // Get number of dimensions
199 public int getNumberOfDimensions(){
200 return this.nDimensions;
201 }
202
203 // Get limits to x
204 public double[] getLimits(){
205 double[] limits = new double[2*this.nDimensions];
206 int j = 0;
207 for(int i=0; i<this.nDimensions; i++){
208 limits[j] = this.xMin[i];
209 j++;
210 limits[j] = this.xMax[i];
211 j++;
212 }
213 return limits;
214 }
215
216 // Display limits to x
217 public void displayLimits(){
218 System.out.println(" ");
219 for(int i=0; i<this.nDimensions; i++){
220 System.out.println("The limits to the x array " + i + " are " + xMin[i] + " and " + xMax[i]);
221 }
222 System.out.println(" ");
223 }
224
225 // Reset the default handing of identical abscissae with different ordinates
226 // from the default option of separating the two relevant abscissae by 0.001 of the range
227 // to avraging the relevant ordinates
228 public void averageIdenticalAbscissae(){
229 this.averageIdenticalAbscissae = true;
230 }
231
232 // Interpolation method
233 public double interpolate(double[] unknownCoord){
234
235 int nUnknown = unknownCoord.length;
236 if(nUnknown!=this.nDimensions)throw new IllegalArgumentException("Number of unknown value coordinates, " + nUnknown + ", does not equal the number of tabulated data dimensions, " + this.nDimensions);
237
238 switch(this.nDimensions){
239 case 0: throw new IllegalArgumentException("data array must have at least one dimension");
240 case 1: // If fOfX is one dimensional perform simple cubic spline
241 this.yValue = ((CubicSpline)(this.method)).interpolate(unknownCoord[0]);
242 break;
243 case 2: // If fOfX is two dimensional perform bicubic spline
244 this.yValue = ((BiCubicSpline)(this.method)).interpolate(unknownCoord[0], unknownCoord[1]);
245 break;
246 case 3: // If fOfX is three dimensional perform tricubic spline
247 this.yValue = ((TriCubicSpline)(this.method)).interpolate(unknownCoord[0], unknownCoord[1], unknownCoord[2]);
248 break;
249 case 4: // If fOfX is four dimensional perform quadricubic spline
250 this.yValue = ((QuadriCubicSpline)(this.method)).interpolate(unknownCoord[0], unknownCoord[1], unknownCoord[2], unknownCoord[2]);
251 break;
252 default: // If fOfX is greater than four dimensional, recursively call PolyCubicSpline
253 // with, as arguments, the n1 fOfX sub-arrays, each of (number of dimensions - 1) dimensions,
254 // where n1 is the number of x1 variables.
255 double[] newCoord = new double[this.nDimensions-1];
256 for(int i=0; i<this.nDimensions-1; i++){
257 newCoord[i] = unknownCoord[i+1];
258 }
259 for(int i=0; i<this.dimOne; i++){
260 csArray[i] = pcs[i].interpolate(newCoord);
261 }
262
263 // Perform simple cubic spline on the array of above returned interpolates
264 CubicSpline ncs = new CubicSpline(this.xArray[0], this.csArray);
265 this.yValue = ncs.interpolate(unknownCoord[0]);
266 }
267
268 return this.yValue;
269 }
270
271 // Set derivatives (internal array)
272 public void setDeriv(Object internalDeriv){
273 this.internalDeriv = internalDeriv;
274 }
275
276 // Get derivatives (internal array)
277 public Object getDeriv(){
278 return this.internalDeriv;
279 }
280
281
282}
283
Note: See TracBrowser for help on using the repository browser.