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 |
|
---|
38 | package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
|
---|
39 |
|
---|
40 | import java.util.ArrayList;
|
---|
41 |
|
---|
42 | import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
|
---|
43 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
|
---|
44 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
|
---|
45 |
|
---|
46 | public 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 |
|
---|