source: src/main/java/agents/anac/y2019/harddealer/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 24.0 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package agents.anac.y2019.harddealer.math3.analysis.interpolation;
18
19import java.util.Arrays;
20import agents.anac.y2019.harddealer.math3.analysis.BivariateFunction;
21import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
22import agents.anac.y2019.harddealer.math3.exception.NoDataException;
23import agents.anac.y2019.harddealer.math3.exception.OutOfRangeException;
24import agents.anac.y2019.harddealer.math3.exception.NonMonotonicSequenceException;
25import agents.anac.y2019.harddealer.math3.util.MathArrays;
26
27/**
28 * Function that implements the
29 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
30 * bicubic spline interpolation</a>. Due to numerical accuracy issues this should not
31 * be used.
32 *
33 * @since 2.1
34 * @deprecated as of 3.4 replaced by
35 * {@link agents.anac.y2019.harddealer.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction}
36 */
37@Deprecated
38public class BicubicSplineInterpolatingFunction
39 implements BivariateFunction {
40 /** Number of coefficients. */
41 private static final int NUM_COEFF = 16;
42 /**
43 * Matrix to compute the spline coefficients from the function values
44 * and function derivatives values
45 */
46 private static final double[][] AINV = {
47 { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
48 { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
49 { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
50 { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
51 { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
52 { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
53 { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
54 { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
55 { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
56 { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
57 { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
58 { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
59 { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
60 { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
61 { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
62 { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
63 };
64
65 /** Samples x-coordinates */
66 private final double[] xval;
67 /** Samples y-coordinates */
68 private final double[] yval;
69 /** Set of cubic splines patching the whole data grid */
70 private final BicubicSplineFunction[][] splines;
71 /**
72 * Partial derivatives.
73 * The value of the first index determines the kind of derivatives:
74 * 0 = first partial derivatives wrt x
75 * 1 = first partial derivatives wrt y
76 * 2 = second partial derivatives wrt x
77 * 3 = second partial derivatives wrt y
78 * 4 = cross partial derivatives
79 */
80 private final BivariateFunction[][][] partialDerivatives;
81
82 /**
83 * @param x Sample values of the x-coordinate, in increasing order.
84 * @param y Sample values of the y-coordinate, in increasing order.
85 * @param f Values of the function on every grid point.
86 * @param dFdX Values of the partial derivative of function with respect
87 * to x on every grid point.
88 * @param dFdY Values of the partial derivative of function with respect
89 * to y on every grid point.
90 * @param d2FdXdY Values of the cross partial derivative of function on
91 * every grid point.
92 * @throws DimensionMismatchException if the various arrays do not contain
93 * the expected number of elements.
94 * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
95 * not strictly increasing.
96 * @throws NoDataException if any of the arrays has zero length.
97 */
98 public BicubicSplineInterpolatingFunction(double[] x,
99 double[] y,
100 double[][] f,
101 double[][] dFdX,
102 double[][] dFdY,
103 double[][] d2FdXdY)
104 throws DimensionMismatchException,
105 NoDataException,
106 NonMonotonicSequenceException {
107 this(x, y, f, dFdX, dFdY, d2FdXdY, false);
108 }
109
110 /**
111 * @param x Sample values of the x-coordinate, in increasing order.
112 * @param y Sample values of the y-coordinate, in increasing order.
113 * @param f Values of the function on every grid point.
114 * @param dFdX Values of the partial derivative of function with respect
115 * to x on every grid point.
116 * @param dFdY Values of the partial derivative of function with respect
117 * to y on every grid point.
118 * @param d2FdXdY Values of the cross partial derivative of function on
119 * every grid point.
120 * @param initializeDerivatives Whether to initialize the internal data
121 * needed for calling any of the methods that compute the partial derivatives
122 * this function.
123 * @throws DimensionMismatchException if the various arrays do not contain
124 * the expected number of elements.
125 * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
126 * not strictly increasing.
127 * @throws NoDataException if any of the arrays has zero length.
128 *
129 * @see #partialDerivativeX(double,double)
130 * @see #partialDerivativeY(double,double)
131 * @see #partialDerivativeXX(double,double)
132 * @see #partialDerivativeYY(double,double)
133 * @see #partialDerivativeXY(double,double)
134 */
135 public BicubicSplineInterpolatingFunction(double[] x,
136 double[] y,
137 double[][] f,
138 double[][] dFdX,
139 double[][] dFdY,
140 double[][] d2FdXdY,
141 boolean initializeDerivatives)
142 throws DimensionMismatchException,
143 NoDataException,
144 NonMonotonicSequenceException {
145 final int xLen = x.length;
146 final int yLen = y.length;
147
148 if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
149 throw new NoDataException();
150 }
151 if (xLen != f.length) {
152 throw new DimensionMismatchException(xLen, f.length);
153 }
154 if (xLen != dFdX.length) {
155 throw new DimensionMismatchException(xLen, dFdX.length);
156 }
157 if (xLen != dFdY.length) {
158 throw new DimensionMismatchException(xLen, dFdY.length);
159 }
160 if (xLen != d2FdXdY.length) {
161 throw new DimensionMismatchException(xLen, d2FdXdY.length);
162 }
163
164 MathArrays.checkOrder(x);
165 MathArrays.checkOrder(y);
166
167 xval = x.clone();
168 yval = y.clone();
169
170 final int lastI = xLen - 1;
171 final int lastJ = yLen - 1;
172 splines = new BicubicSplineFunction[lastI][lastJ];
173
174 for (int i = 0; i < lastI; i++) {
175 if (f[i].length != yLen) {
176 throw new DimensionMismatchException(f[i].length, yLen);
177 }
178 if (dFdX[i].length != yLen) {
179 throw new DimensionMismatchException(dFdX[i].length, yLen);
180 }
181 if (dFdY[i].length != yLen) {
182 throw new DimensionMismatchException(dFdY[i].length, yLen);
183 }
184 if (d2FdXdY[i].length != yLen) {
185 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
186 }
187 final int ip1 = i + 1;
188 for (int j = 0; j < lastJ; j++) {
189 final int jp1 = j + 1;
190 final double[] beta = new double[] {
191 f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
192 dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
193 dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
194 d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
195 };
196
197 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta),
198 initializeDerivatives);
199 }
200 }
201
202 if (initializeDerivatives) {
203 // Compute all partial derivatives.
204 partialDerivatives = new BivariateFunction[5][lastI][lastJ];
205
206 for (int i = 0; i < lastI; i++) {
207 for (int j = 0; j < lastJ; j++) {
208 final BicubicSplineFunction bcs = splines[i][j];
209 partialDerivatives[0][i][j] = bcs.partialDerivativeX();
210 partialDerivatives[1][i][j] = bcs.partialDerivativeY();
211 partialDerivatives[2][i][j] = bcs.partialDerivativeXX();
212 partialDerivatives[3][i][j] = bcs.partialDerivativeYY();
213 partialDerivatives[4][i][j] = bcs.partialDerivativeXY();
214 }
215 }
216 } else {
217 // Partial derivative methods cannot be used.
218 partialDerivatives = null;
219 }
220 }
221
222 /**
223 * {@inheritDoc}
224 */
225 public double value(double x, double y)
226 throws OutOfRangeException {
227 final int i = searchIndex(x, xval);
228 final int j = searchIndex(y, yval);
229
230 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
231 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
232
233 return splines[i][j].value(xN, yN);
234 }
235
236 /**
237 * Indicates whether a point is within the interpolation range.
238 *
239 * @param x First coordinate.
240 * @param y Second coordinate.
241 * @return {@code true} if (x, y) is a valid point.
242 * @since 3.3
243 */
244 public boolean isValidPoint(double x, double y) {
245 if (x < xval[0] ||
246 x > xval[xval.length - 1] ||
247 y < yval[0] ||
248 y > yval[yval.length - 1]) {
249 return false;
250 } else {
251 return true;
252 }
253 }
254
255 /**
256 * @param x x-coordinate.
257 * @param y y-coordinate.
258 * @return the value at point (x, y) of the first partial derivative with
259 * respect to x.
260 * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
261 * the range defined by the boundary values of {@code xval} (resp.
262 * {@code yval}).
263 * @throws NullPointerException if the internal data were not initialized
264 * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
265 * double[][],double[][],double[][],boolean) constructor}).
266 */
267 public double partialDerivativeX(double x, double y)
268 throws OutOfRangeException {
269 return partialDerivative(0, x, y);
270 }
271 /**
272 * @param x x-coordinate.
273 * @param y y-coordinate.
274 * @return the value at point (x, y) of the first partial derivative with
275 * respect to y.
276 * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
277 * the range defined by the boundary values of {@code xval} (resp.
278 * {@code yval}).
279 * @throws NullPointerException if the internal data were not initialized
280 * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
281 * double[][],double[][],double[][],boolean) constructor}).
282 */
283 public double partialDerivativeY(double x, double y)
284 throws OutOfRangeException {
285 return partialDerivative(1, x, y);
286 }
287 /**
288 * @param x x-coordinate.
289 * @param y y-coordinate.
290 * @return the value at point (x, y) of the second partial derivative with
291 * respect to x.
292 * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
293 * the range defined by the boundary values of {@code xval} (resp.
294 * {@code yval}).
295 * @throws NullPointerException if the internal data were not initialized
296 * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
297 * double[][],double[][],double[][],boolean) constructor}).
298 */
299 public double partialDerivativeXX(double x, double y)
300 throws OutOfRangeException {
301 return partialDerivative(2, x, y);
302 }
303 /**
304 * @param x x-coordinate.
305 * @param y y-coordinate.
306 * @return the value at point (x, y) of the second partial derivative with
307 * respect to y.
308 * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
309 * the range defined by the boundary values of {@code xval} (resp.
310 * {@code yval}).
311 * @throws NullPointerException if the internal data were not initialized
312 * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
313 * double[][],double[][],double[][],boolean) constructor}).
314 */
315 public double partialDerivativeYY(double x, double y)
316 throws OutOfRangeException {
317 return partialDerivative(3, x, y);
318 }
319 /**
320 * @param x x-coordinate.
321 * @param y y-coordinate.
322 * @return the value at point (x, y) of the second partial cross-derivative.
323 * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
324 * the range defined by the boundary values of {@code xval} (resp.
325 * {@code yval}).
326 * @throws NullPointerException if the internal data were not initialized
327 * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
328 * double[][],double[][],double[][],boolean) constructor}).
329 */
330 public double partialDerivativeXY(double x, double y)
331 throws OutOfRangeException {
332 return partialDerivative(4, x, y);
333 }
334
335 /**
336 * @param which First index in {@link #partialDerivatives}.
337 * @param x x-coordinate.
338 * @param y y-coordinate.
339 * @return the value at point (x, y) of the selected partial derivative.
340 * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
341 * the range defined by the boundary values of {@code xval} (resp.
342 * {@code yval}).
343 * @throws NullPointerException if the internal data were not initialized
344 * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
345 * double[][],double[][],double[][],boolean) constructor}).
346 */
347 private double partialDerivative(int which, double x, double y)
348 throws OutOfRangeException {
349 final int i = searchIndex(x, xval);
350 final int j = searchIndex(y, yval);
351
352 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
353 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
354
355 return partialDerivatives[which][i][j].value(xN, yN);
356 }
357
358 /**
359 * @param c Coordinate.
360 * @param val Coordinate samples.
361 * @return the index in {@code val} corresponding to the interval
362 * containing {@code c}.
363 * @throws OutOfRangeException if {@code c} is out of the
364 * range defined by the boundary values of {@code val}.
365 */
366 private int searchIndex(double c, double[] val) {
367 final int r = Arrays.binarySearch(val, c);
368
369 if (r == -1 ||
370 r == -val.length - 1) {
371 throw new OutOfRangeException(c, val[0], val[val.length - 1]);
372 }
373
374 if (r < 0) {
375 // "c" in within an interpolation sub-interval: Return the
376 // index of the sample at the lower end of the sub-interval.
377 return -r - 2;
378 }
379 final int last = val.length - 1;
380 if (r == last) {
381 // "c" is the last sample of the range: Return the index
382 // of the sample at the lower end of the last sub-interval.
383 return last - 1;
384 }
385
386 // "c" is another sample point.
387 return r;
388 }
389
390 /**
391 * Compute the spline coefficients from the list of function values and
392 * function partial derivatives values at the four corners of a grid
393 * element. They must be specified in the following order:
394 * <ul>
395 * <li>f(0,0)</li>
396 * <li>f(1,0)</li>
397 * <li>f(0,1)</li>
398 * <li>f(1,1)</li>
399 * <li>f<sub>x</sub>(0,0)</li>
400 * <li>f<sub>x</sub>(1,0)</li>
401 * <li>f<sub>x</sub>(0,1)</li>
402 * <li>f<sub>x</sub>(1,1)</li>
403 * <li>f<sub>y</sub>(0,0)</li>
404 * <li>f<sub>y</sub>(1,0)</li>
405 * <li>f<sub>y</sub>(0,1)</li>
406 * <li>f<sub>y</sub>(1,1)</li>
407 * <li>f<sub>xy</sub>(0,0)</li>
408 * <li>f<sub>xy</sub>(1,0)</li>
409 * <li>f<sub>xy</sub>(0,1)</li>
410 * <li>f<sub>xy</sub>(1,1)</li>
411 * </ul>
412 * where the subscripts indicate the partial derivative with respect to
413 * the corresponding variable(s).
414 *
415 * @param beta List of function values and function partial derivatives
416 * values.
417 * @return the spline coefficients.
418 */
419 private double[] computeSplineCoefficients(double[] beta) {
420 final double[] a = new double[NUM_COEFF];
421
422 for (int i = 0; i < NUM_COEFF; i++) {
423 double result = 0;
424 final double[] row = AINV[i];
425 for (int j = 0; j < NUM_COEFF; j++) {
426 result += row[j] * beta[j];
427 }
428 a[i] = result;
429 }
430
431 return a;
432 }
433}
434
435/**
436 * 2D-spline function.
437 *
438 */
439class BicubicSplineFunction implements BivariateFunction {
440 /** Number of points. */
441 private static final short N = 4;
442 /** Coefficients */
443 private final double[][] a;
444 /** First partial derivative along x. */
445 private final BivariateFunction partialDerivativeX;
446 /** First partial derivative along y. */
447 private final BivariateFunction partialDerivativeY;
448 /** Second partial derivative along x. */
449 private final BivariateFunction partialDerivativeXX;
450 /** Second partial derivative along y. */
451 private final BivariateFunction partialDerivativeYY;
452 /** Second crossed partial derivative. */
453 private final BivariateFunction partialDerivativeXY;
454
455 /**
456 * Simple constructor.
457 *
458 * @param coeff Spline coefficients.
459 */
460 BicubicSplineFunction(double[] coeff) {
461 this(coeff, false);
462 }
463
464 /**
465 * Simple constructor.
466 *
467 * @param coeff Spline coefficients.
468 * @param initializeDerivatives Whether to initialize the internal data
469 * needed for calling any of the methods that compute the partial derivatives
470 * this function.
471 */
472 BicubicSplineFunction(double[] coeff, boolean initializeDerivatives) {
473 a = new double[N][N];
474 for (int i = 0; i < N; i++) {
475 for (int j = 0; j < N; j++) {
476 a[i][j] = coeff[i * N + j];
477 }
478 }
479
480 if (initializeDerivatives) {
481 // Compute all partial derivatives functions.
482 final double[][] aX = new double[N][N];
483 final double[][] aY = new double[N][N];
484 final double[][] aXX = new double[N][N];
485 final double[][] aYY = new double[N][N];
486 final double[][] aXY = new double[N][N];
487
488 for (int i = 0; i < N; i++) {
489 for (int j = 0; j < N; j++) {
490 final double c = a[i][j];
491 aX[i][j] = i * c;
492 aY[i][j] = j * c;
493 aXX[i][j] = (i - 1) * aX[i][j];
494 aYY[i][j] = (j - 1) * aY[i][j];
495 aXY[i][j] = j * aX[i][j];
496 }
497 }
498
499 partialDerivativeX = new BivariateFunction() {
500 /** {@inheritDoc} */
501 public double value(double x, double y) {
502 final double x2 = x * x;
503 final double[] pX = {0, 1, x, x2};
504
505 final double y2 = y * y;
506 final double y3 = y2 * y;
507 final double[] pY = {1, y, y2, y3};
508
509 return apply(pX, pY, aX);
510 }
511 };
512 partialDerivativeY = new BivariateFunction() {
513 /** {@inheritDoc} */
514 public double value(double x, double y) {
515 final double x2 = x * x;
516 final double x3 = x2 * x;
517 final double[] pX = {1, x, x2, x3};
518
519 final double y2 = y * y;
520 final double[] pY = {0, 1, y, y2};
521
522 return apply(pX, pY, aY);
523 }
524 };
525 partialDerivativeXX = new BivariateFunction() {
526 /** {@inheritDoc} */
527 public double value(double x, double y) {
528 final double[] pX = {0, 0, 1, x};
529
530 final double y2 = y * y;
531 final double y3 = y2 * y;
532 final double[] pY = {1, y, y2, y3};
533
534 return apply(pX, pY, aXX);
535 }
536 };
537 partialDerivativeYY = new BivariateFunction() {
538 /** {@inheritDoc} */
539 public double value(double x, double y) {
540 final double x2 = x * x;
541 final double x3 = x2 * x;
542 final double[] pX = {1, x, x2, x3};
543
544 final double[] pY = {0, 0, 1, y};
545
546 return apply(pX, pY, aYY);
547 }
548 };
549 partialDerivativeXY = new BivariateFunction() {
550 /** {@inheritDoc} */
551 public double value(double x, double y) {
552 final double x2 = x * x;
553 final double[] pX = {0, 1, x, x2};
554
555 final double y2 = y * y;
556 final double[] pY = {0, 1, y, y2};
557
558 return apply(pX, pY, aXY);
559 }
560 };
561 } else {
562 partialDerivativeX = null;
563 partialDerivativeY = null;
564 partialDerivativeXX = null;
565 partialDerivativeYY = null;
566 partialDerivativeXY = null;
567 }
568 }
569
570 /**
571 * {@inheritDoc}
572 */
573 public double value(double x, double y) {
574 if (x < 0 || x > 1) {
575 throw new OutOfRangeException(x, 0, 1);
576 }
577 if (y < 0 || y > 1) {
578 throw new OutOfRangeException(y, 0, 1);
579 }
580
581 final double x2 = x * x;
582 final double x3 = x2 * x;
583 final double[] pX = {1, x, x2, x3};
584
585 final double y2 = y * y;
586 final double y3 = y2 * y;
587 final double[] pY = {1, y, y2, y3};
588
589 return apply(pX, pY, a);
590 }
591
592 /**
593 * Compute the value of the bicubic polynomial.
594 *
595 * @param pX Powers of the x-coordinate.
596 * @param pY Powers of the y-coordinate.
597 * @param coeff Spline coefficients.
598 * @return the interpolated value.
599 */
600 private double apply(double[] pX, double[] pY, double[][] coeff) {
601 double result = 0;
602 for (int i = 0; i < N; i++) {
603 for (int j = 0; j < N; j++) {
604 result += coeff[i][j] * pX[i] * pY[j];
605 }
606 }
607
608 return result;
609 }
610
611 /**
612 * @return the partial derivative wrt {@code x}.
613 */
614 public BivariateFunction partialDerivativeX() {
615 return partialDerivativeX;
616 }
617 /**
618 * @return the partial derivative wrt {@code y}.
619 */
620 public BivariateFunction partialDerivativeY() {
621 return partialDerivativeY;
622 }
623 /**
624 * @return the second partial derivative wrt {@code x}.
625 */
626 public BivariateFunction partialDerivativeXX() {
627 return partialDerivativeXX;
628 }
629 /**
630 * @return the second partial derivative wrt {@code y}.
631 */
632 public BivariateFunction partialDerivativeYY() {
633 return partialDerivativeYY;
634 }
635 /**
636 * @return the second partial cross-derivative.
637 */
638 public BivariateFunction partialDerivativeXY() {
639 return partialDerivativeXY;
640 }
641}
Note: See TracBrowser for help on using the repository browser.