source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/interpolation/CubicInterpolation.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.7 KB
Line 
1/**********************************************************
2*
3* CubicInterpolation.java
4*
5* Class for performing an interpolation on the tabulated
6* function y = f(x1) using a cubic interploation procedure
7**
8* WRITTEN BY: Dr Michael Thomas Flanagan
9*
10* DATE: 15-16 January 2011
11* UPDATE:
12*
13* DOCUMENTATION:
14* See Michael Thomas Flanagan's Java library on-line web page:
15* http://www.ee.ucl.ac.uk/~mflanaga/java/CubicInterpolation.html
16* http://www.ee.ucl.ac.uk/~mflanaga/java/
17*
18* Copyright (c) 2011 Michael Thomas Flanagan
19*
20* PERMISSION TO COPY:
21*
22* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
23* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
24* and associated documentation or publications.
25*
26* 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
27* and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
28*
29* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
30* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission from the Michael Thomas Flanagan:
31*
32* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
33* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
34* or its derivatives.
35*
36***************************************************************************************/
37
38package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
39
40import java.util.ArrayList;
41
42import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
43import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
44import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
45
46public class CubicInterpolation{
47
48 private int nPoints = 0; // no. of x tabulated points
49 private double[] x = null; // x in tabulated function f(x1,x2)
50 private double[] y = null; // y=f(x) tabulated function
51 private double[] dydx = null; // dy/dx
52 private boolean derivCalculated = false; // = true when the derivatives have been calculated or entered
53 private CubicSpline cs = null; // cubic spline used in calculating the derivatives
54 private double incrX = 0; // x increment used in calculating the derivatives
55
56 double[][] coeff = null; // cubic coefficients
57
58 private double xx = Double.NaN; // value of x at which an interpolated y value is required
59
60 // Weights used in calculating the grid square coefficients
61 private double[][] weights = {{1.0,0.0,0.0,0.0},{0.0,0.0,1.0,0.0},{-3.0,3.0,-2.0,-1.0},{2.0,-2.0,1.0,1.0}};
62
63 private int[] xIndices = null; // x data indices before ordering
64
65 private double xMin = 0.0; // minimum value of x
66 private double xMax = 0.0; // maximum value of x
67
68 private double interpolatedValue = Double.NaN; // interpolated value of y
69 private double interpolatedDydx = Double.NaN; // interpolated value of dydx
70
71 private boolean numerDiffFlag = true; // = true: if numerical differentiation performed h1 and h2 calculated using delta
72 // = false: if numerical differentiation performed h1 and h2 calculated only provided data points
73 private static double delta = 1e-3; // fractional step factor used in calculating the derivatives
74 private static double potentialRoundingError = 5e-15; // potential rounding error used in checking wheter a value lies within the interpolation bounds (static value)
75 private static boolean roundingCheck = false; // = true: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit (static value)
76
77 // Constructor without derivatives
78 public CubicInterpolation(double[] x, double[] y, int numerDiffOption){
79 // set numerical differencing option
80 if(numerDiffOption==0){
81 this.numerDiffFlag = false;
82 }
83 else{
84 if(numerDiffOption==1){
85 this.numerDiffFlag = true;
86 }
87 else{
88 throw new IllegalArgumentException("The numerical differencing option, " + numerDiffOption + ", must be 0 or 1");
89 }
90 }
91
92 // initialize the data
93 this.initialize(Conv.copy(x), Conv.copy(y));
94
95 // calculate the derivatives
96 this.calcDeriv();
97
98 // calculate grid coefficients for all grid squares
99 this.gridCoefficients();
100 }
101
102 // Constructor with derivatives
103 public CubicInterpolation(double[] x, double[] y, double[] dydx){
104
105 // initialize the data
106 this.initialize(Conv.copy(x), Conv.copy(y), Conv.copy(dydx));
107
108 // calculate grid coefficients for all grid squares
109 this.gridCoefficients();
110 }
111
112 // Initialize the data
113 private void initialize(double[] x, double[] y){
114 this.initialize(x, y, null, false);
115 }
116
117 private void initialize(double[] x, double[] y, double[] dydx){
118 this.initialize(x, y, dydx, true);
119 }
120
121 private void initialize(double[] x, double[] y, double[] dydx, boolean flag){
122 int nl = 3;
123 if(flag)nl = 2;
124 int nPoints=x.length;
125 if(nPoints!=y.length)throw new IllegalArgumentException("Arrays x and y-row are of different length " + nPoints + " " + y.length);
126 if(nPoints<nl)throw new IllegalArgumentException("The data matrix must have a minimum size of " + nl + " X " + nl);
127
128 // order data
129 ArrayMaths am = new ArrayMaths(x);
130 am = am.sort();
131 this.xIndices = am.originalIndices();
132 x = am.array();
133 double[] hold = new double[nPoints];
134
135 for(int i=0; i<nPoints; i++){
136 hold[i] = y[this.xIndices[i]];
137 }
138 for(int i=0; i<nPoints; i++){
139 y[i] = hold[i];
140 }
141
142 if(flag){
143 for(int i=0; i<nPoints; i++){
144 hold[i] = dydx[this.xIndices[i]];
145 }
146 for(int i=0; i<nPoints; i++){
147 dydx[i] = hold[i];
148 }
149 }
150
151 // check for identical x values
152 for(int i=1; i<nPoints; i++){
153 if(x[i]==x[i-1]){
154 System.out.println("x["+this.xIndices[i]+"] and x["+this.xIndices[i+1]+"] are identical, " + x[i]);
155 System.out.println("The y values have been averaged and one point has been deleted");
156 y[i-1] = (y[i-1] + y[i])/2.0;
157 for(int j=i; j<nPoints-1; j++){
158 x[j]=x[j+1];
159 y[j]=y[j+1];
160 this.xIndices[j]=this.xIndices[j+1];
161 }
162 if(flag){
163 dydx[i-1] = (dydx[i-1] + dydx[i])/2.0;
164 for(int j=i; j<nPoints-1; j++){
165 dydx[j]=dydx[j+1];
166 }
167 }
168 nPoints--;
169 }
170 }
171
172 // assign variables
173 this.nPoints = nPoints;
174 this.x = new double[this.nPoints];
175 this.y = new double[this.nPoints];
176 this.dydx = new double[this.nPoints];
177
178 for(int i=0; i<this.nPoints; i++){
179 this.x[i]=x[i];
180 this.y[i]=y[i];
181 }
182 if(flag){
183 for(int j=0; j<this.nPoints; j++){
184 this.dydx[j]=dydx[j];
185 }
186 this.derivCalculated = true;
187 }
188
189 // limits
190 this.xMin = Fmath.minimum(this.x);
191 this.xMax = Fmath.maximum(this.x);
192
193 if(!flag && this.numerDiffFlag){
194 // numerical difference increments
195 double range = this.xMax - this.xMin;
196 double averageSeparation = range/this.nPoints;
197 double minSep = this.x[1] - this.x[0];
198 double minimumSeparation = minSep;
199 for(int i=2; i<this.nPoints; i++){
200 minSep = this.x[i] - this.x[i-1];
201 if(minSep<minimumSeparation)minimumSeparation = minSep;
202 }
203
204 this.incrX = range*CubicInterpolation.delta;
205 double defaultIncr = minimumSeparation;
206 if(minimumSeparation<averageSeparation/10.0)defaultIncr = averageSeparation/10.0;
207 if(this.incrX>averageSeparation)this.incrX = defaultIncr;
208 }
209 }
210
211 // Calculate the derivatives
212 private void calcDeriv(){
213
214 if(this.numerDiffFlag){
215 // Numerical differentiation using delta and interpolation
216 this.cs = new CubicSpline(this.x, this.y);
217 double[] xjp1 = new double[this.nPoints];
218 double[] xjm1 = new double[this.nPoints];
219 for(int i=0; i<this.nPoints; i++){
220 xjp1[i] = this.x[i] + this.incrX;
221 if(xjp1[i]>this.x[this.nPoints-1])xjp1[i] = this.x[this.nPoints-1];
222 xjm1[i] = this.x[i] - this.incrX;
223 if(xjm1[i]<this.x[0])xjm1[i] = this.x[0];
224 }
225
226 for(int i=0; i<this.nPoints; i++){
227 this.dydx[i] = (cs.interpolate(xjp1[i]) - cs.interpolate(xjm1[i]))/(xjp1[i] - xjm1[i]);
228 }
229 }
230 else{
231 // Numerical differentiation using provided data points
232 int iip =0;
233 int iim =0;
234 for(int i=0; i<this.nPoints; i++){
235 iip = i+1;
236 if(iip>=this.nPoints)iip = this.nPoints-1;
237 iim = i-1;
238 if(iim<0)iim = 0;
239 this.dydx[i] = (this.y[iip] - this.y[iim])/(this.x[iip] - this.x[iim]);
240 }
241 }
242 this.derivCalculated = true;
243 }
244
245 // Grid coefficients
246 private void gridCoefficients(){
247
248 double[] xt = new double[4];
249 this.coeff = new double[this.nPoints][4];
250 double d1 = 0.0;
251 for(int i=0; i<this.nPoints-1; i++){
252 d1 = this.x[i+1] - this.x[i];
253 xt[0] = this.y[i];
254 xt[1] = this.y[i+1];
255 xt[2] = this.dydx[i]*d1;
256 xt[3] = this.dydx[i+1]*d1;
257
258 double xh = 0.0;
259 for(int k=0; k<4; k++){
260 for(int kk=0; kk<4; kk++){
261 xh += this.weights[k][kk]*xt[kk];
262 }
263 this.coeff[i][k] = xh;
264 xh = 0.0;
265 }
266 }
267 }
268
269 // Returns an interpolated value of y for a value of x
270 // from a tabulated function y=f(x,x2)
271 public double interpolate(double xx){
272 // check that xx and xx2 are within the limits
273 if(xx<x[0]){
274 if(xx>=x[0]-CubicInterpolation.potentialRoundingError){
275 xx=this.x[0];
276 }
277 else{
278 throw new IllegalArgumentException(xx + " is outside the limits, " + x[0] + " - " + x[this.nPoints-1]);
279 }
280 }
281
282 if(xx>this.x[this.nPoints-1]){
283 if(xx<=this.x[this.nPoints-1]+CubicInterpolation.potentialRoundingError){
284 xx=this.x[this.nPoints-1];
285 }
286 else{
287 throw new IllegalArgumentException(xx + " is outside the limits, " + this.x[0] + " - " + this.x[this.nPoints-1]);
288 }
289 }
290
291 // assign variables
292 this.xx = xx;
293
294 // Find grid surrounding the interpolation point
295 int gridn = 0;
296 int counter = 1;
297 boolean test = true;
298 while(test){
299 if(xx<this.x[counter]){
300 gridn = counter - 1;
301 test = false;
302 }
303 else{
304 counter++;
305 if(counter>=this.nPoints){
306 gridn = this.nPoints-2;
307 test = false;
308 }
309 }
310 }
311
312 // interpolation
313 double xNormalised = (xx - x[gridn])/(x[gridn+1] - x[gridn]);
314 this.interpolatedValue = 0.0;
315 for(int i=0; i<4; i++){
316 this.interpolatedValue += this.coeff[gridn][i]*Math.pow(xNormalised, i);
317 }
318 this.interpolatedDydx = 0.0;
319 for(int i=1; i<4; i++){
320 this.interpolatedDydx += i*this.coeff[gridn][i]*Math.pow(xNormalised, i-1);
321 }
322
323 return this.interpolatedValue;
324 }
325
326 // Return last interpolated value and the interpolated gradients
327 public double[] getInterpolatedValues(){
328 double[] ret = new double[3];
329 ret[0] = this.interpolatedValue;
330 ret[1] = this.interpolatedDydx;
331 ret[2] = this.xx;
332 return ret;
333 }
334
335 // Return grid point values of dydx
336 public double[] getGridDydx(){
337 double[] ret = new double[this.nPoints];
338 for(int i=0; i<this.nPoints; i++){
339 ret[this.xIndices[i]] = this.dydx[i];
340 }
341 return ret;
342 }
343
344 // Reset the numerical differentiation incremental factor delta
345 public static void resetDelta(double delta){
346 CubicInterpolation.delta = delta;
347 }
348
349 // Reset rounding error check option
350 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
351 // This method causes this check to be ignored and an exception to be thrown if any point lies outside the interpolation bounds
352 public static void noRoundingErrorCheck(){
353 CubicInterpolation.roundingCheck = false;
354 CubicInterpolation.potentialRoundingError = 0.0;
355 }
356
357 // Reset potential rounding error value
358 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
359 // The default value for the potential rounding error is 5e-15*times the 10^exponent of the value outside the bounds
360 // This method allows the 5e-15 to be reset
361 public static void potentialRoundingError(double potentialRoundingError){
362 CubicInterpolation.potentialRoundingError = potentialRoundingError;
363 }
364
365 // Get minimum limit
366 public double getXmin(){
367 return this.xMin;
368 }
369
370 // Get maximum limit
371 public double getXmax(){
372 return this.xMax;
373 }
374
375 // Get limits to x
376 public double[] getLimits(){
377 double[] limits = {this.xMin, this.xMax};
378 return limits;
379 }
380
381 // Display limits to x
382 public void displayLimits(){
383 System.out.println(" ");
384 System.out.println("The limits to the x array are " + this.xMin + " and " + this.xMax);
385 System.out.println(" ");
386 }
387
388}
389
Note: See TracBrowser for help on using the repository browser.