From 31fae6431438e26d6b47b988164847048ceab314 Mon Sep 17 00:00:00 2001 From: Hank Grabowski Date: Mon, 20 Oct 2014 21:31:39 -0400 Subject: [PATCH] MATH-1138 #comment BicubicSplineInterpolator Maintenance Moved the new bi-cubic interpolator functions to new Piecewise* version of the same files. Restored the original bi-cubic functions and restored the tricubic function to use the original functions as well --- .../BicubicSplineInterpolatingFunction.java | 618 ++++++++++++-- .../BicubicSplineInterpolator.java | 141 +++- ...iseBicubicSplineInterpolatingFunction.java | 195 +++++ .../PiecewiseBicubicSplineInterpolator.java | 64 ++ ...ngPolynomialBicubicSplineInterpolator.java | 4 +- .../TricubicSplineInterpolator.java | 35 +- ...icubicSplineInterpolatingFunctionTest.java | 766 +++++++++++++----- .../BicubicSplineInterpolatorTest.java | 207 ++--- ...icubicSplineInterpolatingFunctionTest.java | 280 +++++++ ...iecewiseBicubicSplineInterpolatorTest.java | 277 +++++++ .../TricubicSplineInterpolatorTest.java | 2 +- 11 files changed, 2140 insertions(+), 449 deletions(-) create mode 100644 src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java create mode 100644 src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java create mode 100644 src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunctionTest.java create mode 100644 src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatorTest.java diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java index c0ce3c556..079e9fc7e 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java @@ -18,76 +18,143 @@ package org.apache.commons.math3.analysis.interpolation; import java.util.Arrays; import org.apache.commons.math3.analysis.BivariateFunction; -import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.InsufficientDataException; import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NullArgumentException; import org.apache.commons.math3.exception.OutOfRangeException; import org.apache.commons.math3.exception.NonMonotonicSequenceException; import org.apache.commons.math3.util.MathArrays; /** - * Function that implements the bicubic spline - * interpolation. + * Function that implements the + * + * bicubic spline interpolation. * * @since 2.1 */ public class BicubicSplineInterpolatingFunction implements BivariateFunction { + /** Number of coefficients. */ + private static final int NUM_COEFF = 16; + /** + * Matrix to compute the spline coefficients from the function values + * and function derivatives values + */ + private static final double[][] AINV = { + { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, + { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, + { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, + { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, + { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, + { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, + { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, + { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, + { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, + { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, + { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } + }; /** Samples x-coordinates */ private final double[] xval; - /** Samples y-coordinates */ private final double[] yval; - /** Set of cubic splines patching the whole data grid */ - private final double[][] fval; + private final BicubicSplineFunction[][] splines; + /** + * Partial derivatives. + * The value of the first index determines the kind of derivatives: + * 0 = first partial derivatives wrt x + * 1 = first partial derivatives wrt y + * 2 = second partial derivatives wrt x + * 3 = second partial derivatives wrt y + * 4 = cross partial derivatives + */ + private final BivariateFunction[][][] partialDerivatives; /** * @param x Sample values of the x-coordinate, in increasing order. * @param y Sample values of the y-coordinate, in increasing order. - * @param f Values of the function on every grid point. the expected number - * of elements. - * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not - * strictly increasing. - * @throws NullArgumentException if any of the arguments are null + * @param f Values of the function on every grid point. + * @param dFdX Values of the partial derivative of function with respect + * to x on every grid point. + * @param dFdY Values of the partial derivative of function with respect + * to y on every grid point. + * @param d2FdXdY Values of the cross partial derivative of function on + * every grid point. + * @throws DimensionMismatchException if the various arrays do not contain + * the expected number of elements. + * @throws NonMonotonicSequenceException if {@code x} or {@code y} are + * not strictly increasing. * @throws NoDataException if any of the arrays has zero length. - * @throws DimensionMismatchException if the length of x and y don't match the row, column - * height of f */ + public BicubicSplineInterpolatingFunction(double[] x, + double[] y, + double[][] f, + double[][] dFdX, + double[][] dFdY, + double[][] d2FdXdY) + throws DimensionMismatchException, + NoDataException, + NonMonotonicSequenceException { + this(x, y, f, dFdX, dFdY, d2FdXdY, false); + } - public BicubicSplineInterpolatingFunction(double[] x, double[] y, - double[][] f) - throws DimensionMismatchException, NullArgumentException, - NoDataException, NonMonotonicSequenceException { - - final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS; - - if (x == null || y == null || f == null || f[0] == null) { - throw new NullArgumentException(); - } - + /** + * @param x Sample values of the x-coordinate, in increasing order. + * @param y Sample values of the y-coordinate, in increasing order. + * @param f Values of the function on every grid point. + * @param dFdX Values of the partial derivative of function with respect + * to x on every grid point. + * @param dFdY Values of the partial derivative of function with respect + * to y on every grid point. + * @param d2FdXdY Values of the cross partial derivative of function on + * every grid point. + * @param initializeDerivatives Whether to initialize the internal data + * needed for calling any of the methods that compute the partial derivatives + * this function. + * @throws DimensionMismatchException if the various arrays do not contain + * the expected number of elements. + * @throws NonMonotonicSequenceException if {@code x} or {@code y} are + * not strictly increasing. + * @throws NoDataException if any of the arrays has zero length. + * + * @see #partialDerivativeX(double,double) + * @see #partialDerivativeY(double,double) + * @see #partialDerivativeXX(double,double) + * @see #partialDerivativeYY(double,double) + * @see #partialDerivativeXY(double,double) + */ + public BicubicSplineInterpolatingFunction(double[] x, + double[] y, + double[][] f, + double[][] dFdX, + double[][] dFdY, + double[][] d2FdXdY, + boolean initializeDerivatives) + throws DimensionMismatchException, + NoDataException, + NonMonotonicSequenceException { final int xLen = x.length; final int yLen = y.length; if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { throw new NoDataException(); } - - if (xLen < minimumLength || yLen < minimumLength || - f.length < minimumLength || f[0].length < minimumLength) { - throw new InsufficientDataException(); - } - if (xLen != f.length) { throw new DimensionMismatchException(xLen, f.length); } - - if (yLen != f[0].length) { - throw new DimensionMismatchException(yLen, f[0].length); + if (xLen != dFdX.length) { + throw new DimensionMismatchException(xLen, dFdX.length); + } + if (xLen != dFdY.length) { + throw new DimensionMismatchException(xLen, dFdY.length); + } + if (xLen != d2FdXdY.length) { + throw new DimensionMismatchException(xLen, d2FdXdY.length); } MathArrays.checkOrder(x); @@ -95,7 +162,57 @@ public class BicubicSplineInterpolatingFunction xval = x.clone(); yval = y.clone(); - fval = f.clone(); + + final int lastI = xLen - 1; + final int lastJ = yLen - 1; + splines = new BicubicSplineFunction[lastI][lastJ]; + + for (int i = 0; i < lastI; i++) { + if (f[i].length != yLen) { + throw new DimensionMismatchException(f[i].length, yLen); + } + if (dFdX[i].length != yLen) { + throw new DimensionMismatchException(dFdX[i].length, yLen); + } + if (dFdY[i].length != yLen) { + throw new DimensionMismatchException(dFdY[i].length, yLen); + } + if (d2FdXdY[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdY[i].length, yLen); + } + final int ip1 = i + 1; + for (int j = 0; j < lastJ; j++) { + final int jp1 = j + 1; + final double[] beta = new double[] { + f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1], + dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1], + dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1], + d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1] + }; + + splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta), + initializeDerivatives); + } + } + + if (initializeDerivatives) { + // Compute all partial derivatives. + partialDerivatives = new BivariateFunction[5][lastI][lastJ]; + + for (int i = 0; i < lastI; i++) { + for (int j = 0; j < lastJ; j++) { + final BicubicSplineFunction bcs = splines[i][j]; + partialDerivatives[0][i][j] = bcs.partialDerivativeX(); + partialDerivatives[1][i][j] = bcs.partialDerivativeY(); + partialDerivatives[2][i][j] = bcs.partialDerivativeXX(); + partialDerivatives[3][i][j] = bcs.partialDerivativeYY(); + partialDerivatives[4][i][j] = bcs.partialDerivativeXY(); + } + } + } else { + // Partial derivative methods cannot be used. + partialDerivatives = null; + } } /** @@ -103,37 +220,13 @@ public class BicubicSplineInterpolatingFunction */ public double value(double x, double y) throws OutOfRangeException { - int index; - PolynomialSplineFunction spline; - AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator(); - final int offset = 2; - final int count = offset + 3; - final int i = searchIndex(x, xval, offset, count); - final int j = searchIndex(y, yval, offset, count); + final int i = searchIndex(x, xval); + final int j = searchIndex(y, yval); - double xArray[] = new double[count]; - double yArray[] = new double[count]; - double zArray[] = new double[count]; - double interpArray[] = new double[count]; + final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); + final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); - for (index = 0; index < count; index++) { - xArray[index] = xval[i + index]; - yArray[index] = yval[j + index]; - } - - for (int zIndex = 0; zIndex < count; zIndex++) { - for (index = 0; index < count; index++) { - zArray[index] = fval[i + index][j + zIndex]; - } - spline = interpolator.interpolate(xArray, zArray); - interpArray[zIndex] = spline.value(x); - } - - spline = interpolator.interpolate(yArray, interpArray); - - double returnValue = spline.value(y); - - return returnValue; + return splines[i][j].value(xN, yN); } /** @@ -145,7 +238,9 @@ public class BicubicSplineInterpolatingFunction * @since 3.3 */ public boolean isValidPoint(double x, double y) { - if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] || + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || y > yval[yval.length - 1]) { return false; } else { @@ -153,43 +248,386 @@ public class BicubicSplineInterpolatingFunction } } + /** + * @param x x-coordinate. + * @param y y-coordinate. + * @return the value at point (x, y) of the first partial derivative with + * respect to x. + * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside + * the range defined by the boundary values of {@code xval} (resp. + * {@code yval}). + * @throws NullPointerException if the internal data were not initialized + * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][], + * double[][],double[][],double[][],boolean) constructor}). + */ + public double partialDerivativeX(double x, double y) + throws OutOfRangeException { + return partialDerivative(0, x, y); + } + /** + * @param x x-coordinate. + * @param y y-coordinate. + * @return the value at point (x, y) of the first partial derivative with + * respect to y. + * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside + * the range defined by the boundary values of {@code xval} (resp. + * {@code yval}). + * @throws NullPointerException if the internal data were not initialized + * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][], + * double[][],double[][],double[][],boolean) constructor}). + */ + public double partialDerivativeY(double x, double y) + throws OutOfRangeException { + return partialDerivative(1, x, y); + } + /** + * @param x x-coordinate. + * @param y y-coordinate. + * @return the value at point (x, y) of the second partial derivative with + * respect to x. + * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside + * the range defined by the boundary values of {@code xval} (resp. + * {@code yval}). + * @throws NullPointerException if the internal data were not initialized + * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][], + * double[][],double[][],double[][],boolean) constructor}). + */ + public double partialDerivativeXX(double x, double y) + throws OutOfRangeException { + return partialDerivative(2, x, y); + } + /** + * @param x x-coordinate. + * @param y y-coordinate. + * @return the value at point (x, y) of the second partial derivative with + * respect to y. + * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside + * the range defined by the boundary values of {@code xval} (resp. + * {@code yval}). + * @throws NullPointerException if the internal data were not initialized + * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][], + * double[][],double[][],double[][],boolean) constructor}). + */ + public double partialDerivativeYY(double x, double y) + throws OutOfRangeException { + return partialDerivative(3, x, y); + } + /** + * @param x x-coordinate. + * @param y y-coordinate. + * @return the value at point (x, y) of the second partial cross-derivative. + * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside + * the range defined by the boundary values of {@code xval} (resp. + * {@code yval}). + * @throws NullPointerException if the internal data were not initialized + * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][], + * double[][],double[][],double[][],boolean) constructor}). + */ + public double partialDerivativeXY(double x, double y) + throws OutOfRangeException { + return partialDerivative(4, x, y); + } + + /** + * @param which First index in {@link #partialDerivatives}. + * @param x x-coordinate. + * @param y y-coordinate. + * @return the value at point (x, y) of the selected partial derivative. + * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside + * the range defined by the boundary values of {@code xval} (resp. + * {@code yval}). + * @throws NullPointerException if the internal data were not initialized + * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][], + * double[][],double[][],double[][],boolean) constructor}). + */ + private double partialDerivative(int which, double x, double y) + throws OutOfRangeException { + final int i = searchIndex(x, xval); + final int j = searchIndex(y, yval); + + final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); + final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); + + return partialDerivatives[which][i][j].value(xN, yN); + } + /** * @param c Coordinate. * @param val Coordinate samples. - * @param offset how far back from found value to offset for querying - * @param count total number of elements forward from beginning that will be - * queried - * @return the index in {@code val} corresponding to the interval containing - * {@code c}. - * @throws OutOfRangeException if {@code c} is out of the range defined by - * the boundary values of {@code val}. + * @return the index in {@code val} corresponding to the interval + * containing {@code c}. + * @throws OutOfRangeException if {@code c} is out of the + * range defined by the boundary values of {@code val}. */ - private int searchIndex(double c, double[] val, int offset, int count) { - int r = Arrays.binarySearch(val, c); + private int searchIndex(double c, double[] val) { + final int r = Arrays.binarySearch(val, c); - if (r == -1 || r == -val.length - 1) { + if (r == -1 || + r == -val.length - 1) { throw new OutOfRangeException(c, val[0], val[val.length - 1]); } if (r < 0) { - // "c" in within an interpolation sub-interval, which returns - // negative - // need to remove the negative sign for consistency - r = -r - offset - 1; - } else { - r -= offset; + // "c" in within an interpolation sub-interval: Return the + // index of the sample at the lower end of the sub-interval. + return -r - 2; } - - if (r < 0) { - r = 0; - } - - if ((r + count) >= val.length) { + final int last = val.length - 1; + if (r == last) { // "c" is the last sample of the range: Return the index // of the sample at the lower end of the last sub-interval. - r = val.length - count; + return last - 1; } + // "c" is another sample point. return r; } + + /** + * Compute the spline coefficients from the list of function values and + * function partial derivatives values at the four corners of a grid + * element. They must be specified in the following order: + * + * where the subscripts indicate the partial derivative with respect to + * the corresponding variable(s). + * + * @param beta List of function values and function partial derivatives + * values. + * @return the spline coefficients. + */ + private double[] computeSplineCoefficients(double[] beta) { + final double[] a = new double[NUM_COEFF]; + + for (int i = 0; i < NUM_COEFF; i++) { + double result = 0; + final double[] row = AINV[i]; + for (int j = 0; j < NUM_COEFF; j++) { + result += row[j] * beta[j]; + } + a[i] = result; + } + + return a; + } +} + +/** + * 2D-spline function. + * + */ +class BicubicSplineFunction implements BivariateFunction { + /** Number of points. */ + private static final short N = 4; + /** Coefficients */ + private final double[][] a; + /** First partial derivative along x. */ + private final BivariateFunction partialDerivativeX; + /** First partial derivative along y. */ + private final BivariateFunction partialDerivativeY; + /** Second partial derivative along x. */ + private final BivariateFunction partialDerivativeXX; + /** Second partial derivative along y. */ + private final BivariateFunction partialDerivativeYY; + /** Second crossed partial derivative. */ + private final BivariateFunction partialDerivativeXY; + + /** + * Simple constructor. + * + * @param coeff Spline coefficients. + */ + public BicubicSplineFunction(double[] coeff) { + this(coeff, false); + } + + /** + * Simple constructor. + * + * @param coeff Spline coefficients. + * @param initializeDerivatives Whether to initialize the internal data + * needed for calling any of the methods that compute the partial derivatives + * this function. + */ + public BicubicSplineFunction(double[] coeff, + boolean initializeDerivatives) { + a = new double[N][N]; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + a[i][j] = coeff[i * N + j]; + } + } + + if (initializeDerivatives) { + // Compute all partial derivatives functions. + final double[][] aX = new double[N][N]; + final double[][] aY = new double[N][N]; + final double[][] aXX = new double[N][N]; + final double[][] aYY = new double[N][N]; + final double[][] aXY = new double[N][N]; + + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + final double c = a[i][j]; + aX[i][j] = i * c; + aY[i][j] = j * c; + aXX[i][j] = (i - 1) * aX[i][j]; + aYY[i][j] = (j - 1) * aY[i][j]; + aXY[i][j] = j * aX[i][j]; + } + } + + partialDerivativeX = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double[] pX = {0, 1, x, x2}; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = {1, y, y2, y3}; + + return apply(pX, pY, aX); + } + }; + partialDerivativeY = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = {1, x, x2, x3}; + + final double y2 = y * y; + final double[] pY = {0, 1, y, y2}; + + return apply(pX, pY, aY); + } + }; + partialDerivativeXX = new BivariateFunction() { + public double value(double x, double y) { + final double[] pX = {0, 0, 1, x}; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = {1, y, y2, y3}; + + return apply(pX, pY, aXX); + } + }; + partialDerivativeYY = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = {1, x, x2, x3}; + + final double[] pY = {0, 0, 1, y}; + + return apply(pX, pY, aYY); + } + }; + partialDerivativeXY = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double[] pX = {0, 1, x, x2}; + + final double y2 = y * y; + final double[] pY = {0, 1, y, y2}; + + return apply(pX, pY, aXY); + } + }; + } else { + partialDerivativeX = null; + partialDerivativeY = null; + partialDerivativeXX = null; + partialDerivativeYY = null; + partialDerivativeXY = null; + } + } + + /** + * {@inheritDoc} + */ + public double value(double x, double y) { + if (x < 0 || x > 1) { + throw new OutOfRangeException(x, 0, 1); + } + if (y < 0 || y > 1) { + throw new OutOfRangeException(y, 0, 1); + } + + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = {1, x, x2, x3}; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = {1, y, y2, y3}; + + return apply(pX, pY, a); + } + + /** + * Compute the value of the bicubic polynomial. + * + * @param pX Powers of the x-coordinate. + * @param pY Powers of the y-coordinate. + * @param coeff Spline coefficients. + * @return the interpolated value. + */ + private double apply(double[] pX, double[] pY, double[][] coeff) { + double result = 0; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + result += coeff[i][j] * pX[i] * pY[j]; + } + } + + return result; + } + + /** + * @return the partial derivative wrt {@code x}. + */ + public BivariateFunction partialDerivativeX() { + return partialDerivativeX; + } + /** + * @return the partial derivative wrt {@code y}. + */ + public BivariateFunction partialDerivativeY() { + return partialDerivativeY; + } + /** + * @return the second partial derivative wrt {@code x}. + */ + public BivariateFunction partialDerivativeXX() { + return partialDerivativeXX; + } + /** + * @return the second partial derivative wrt {@code y}. + */ + public BivariateFunction partialDerivativeYY() { + return partialDerivativeYY; + } + /** + * @return the second partial cross-derivative. + */ + public BivariateFunction partialDerivativeXY() { + return partialDerivativeXY; + } } diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java index a973f8121..5e2c92fde 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java @@ -16,10 +16,12 @@ */ package org.apache.commons.math3.analysis.interpolation; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; import org.apache.commons.math3.exception.DimensionMismatchException; import org.apache.commons.math3.exception.NoDataException; import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; import org.apache.commons.math3.util.MathArrays; /** @@ -28,37 +30,144 @@ import org.apache.commons.math3.util.MathArrays; * @since 2.2 */ public class BicubicSplineInterpolator - implements BivariateGridInterpolator -{ + implements BivariateGridInterpolator { + /** Whether to initialize internal data used to compute the analytical + derivatives of the splines. */ + private final boolean initializeDerivatives; /** * Default constructor. + * The argument {@link #BicubicSplineInterpolator(boolean) initializeDerivatives} + * is set to {@code false}. */ - public BicubicSplineInterpolator() - { + public BicubicSplineInterpolator() { + this(false); + } + /** + * Creates an interpolator. + * + * @param initializeDerivatives Whether to initialize the internal data + * needed for calling any of the methods that compute the partial derivatives + * of the {@link BicubicSplineInterpolatingFunction function} returned from + * the call to {@link #interpolate(double[],double[],double[][]) interpolate}. + */ + public BicubicSplineInterpolator(boolean initializeDerivatives) { + this.initializeDerivatives = initializeDerivatives; } /** * {@inheritDoc} */ - public BicubicSplineInterpolatingFunction interpolate( final double[] xval, final double[] yval, + public BicubicSplineInterpolatingFunction interpolate(final double[] xval, + final double[] yval, final double[][] fval) - throws DimensionMismatchException, NullArgumentException, NoDataException, NonMonotonicSequenceException - { - if ( xval == null || yval == null || fval == null || fval[0] == null ) - { - throw new NullArgumentException(); - } - - if ( xval.length == 0 || yval.length == 0 || fval.length == 0 ) - { + throws NoDataException, DimensionMismatchException, + NonMonotonicSequenceException, NumberIsTooSmallException { + if (xval.length == 0 || yval.length == 0 || fval.length == 0) { throw new NoDataException(); } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } MathArrays.checkOrder(xval); MathArrays.checkOrder(yval); - return new BicubicSplineInterpolatingFunction( xval, yval, fval ); + final int xLen = xval.length; + final int yLen = yval.length; + + // Samples (first index is y-coordinate, i.e. subarray variable is x) + // 0 <= i < xval.length + // 0 <= j < yval.length + // fX[j][i] = f(xval[i], yval[j]) + final double[][] fX = new double[yLen][xLen]; + for (int i = 0; i < xLen; i++) { + if (fval[i].length != yLen) { + throw new DimensionMismatchException(fval[i].length, yLen); + } + + for (int j = 0; j < yLen; j++) { + fX[j][i] = fval[i][j]; + } + } + + final SplineInterpolator spInterpolator = new SplineInterpolator(); + + // For each line y[j] (0 <= j < yLen), construct a 1D spline with + // respect to variable x + final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; + for (int j = 0; j < yLen; j++) { + ySplineX[j] = spInterpolator.interpolate(xval, fX[j]); + } + + // For each line x[i] (0 <= i < xLen), construct a 1D spline with + // respect to variable y generated by array fY_1[i] + final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; + for (int i = 0; i < xLen; i++) { + xSplineY[i] = spInterpolator.interpolate(yval, fval[i]); + } + + // Partial derivatives with respect to x at the grid knots + final double[][] dFdX = new double[xLen][yLen]; + for (int j = 0; j < yLen; j++) { + final UnivariateFunction f = ySplineX[j].derivative(); + for (int i = 0; i < xLen; i++) { + dFdX[i][j] = f.value(xval[i]); + } + } + + // Partial derivatives with respect to y at the grid knots + final double[][] dFdY = new double[xLen][yLen]; + for (int i = 0; i < xLen; i++) { + final UnivariateFunction f = xSplineY[i].derivative(); + for (int j = 0; j < yLen; j++) { + dFdY[i][j] = f.value(yval[j]); + } + } + + // Cross partial derivatives + final double[][] d2FdXdY = new double[xLen][yLen]; + for (int i = 0; i < xLen ; i++) { + final int nI = nextIndex(i, xLen); + final int pI = previousIndex(i); + for (int j = 0; j < yLen; j++) { + final int nJ = nextIndex(j, yLen); + final int pJ = previousIndex(j); + d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - + fval[pI][nJ] + fval[pI][pJ]) / + ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); + } + } + + // Create the interpolating splines + return new BicubicSplineInterpolatingFunction(xval, yval, fval, + dFdX, dFdY, d2FdXdY, + initializeDerivatives); + } + + /** + * Computes the next index of an array, clipping if necessary. + * It is assumed (but not checked) that {@code i >= 0}. + * + * @param i Index. + * @param max Upper limit of the array. + * @return the next index. + */ + private int nextIndex(int i, int max) { + final int index = i + 1; + return index < max ? index : index - 1; + } + /** + * Computes the previous index of an array, clipping if necessary. + * It is assumed (but not checked) that {@code i} is smaller than the size + * of the array. + * + * @param i Index. + * @return the previous index. + */ + private int previousIndex(int i) { + final int index = i - 1; + return index >= 0 ? index : 0; } } diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java new file mode 100644 index 000000000..7942b5229 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java @@ -0,0 +1,195 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import java.util.Arrays; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.InsufficientDataException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Function that implements the bicubic spline + * interpolation. + * + * @since 2.1 + */ +public class PiecewiseBicubicSplineInterpolatingFunction + implements BivariateFunction { + + /** Samples x-coordinates */ + private final double[] xval; + + /** Samples y-coordinates */ + private final double[] yval; + + /** Set of cubic splines patching the whole data grid */ + private final double[][] fval; + + /** + * @param x Sample values of the x-coordinate, in increasing order. + * @param y Sample values of the y-coordinate, in increasing order. + * @param f Values of the function on every grid point. the expected number + * of elements. + * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not + * strictly increasing. + * @throws NullArgumentException if any of the arguments are null + * @throws NoDataException if any of the arrays has zero length. + * @throws DimensionMismatchException if the length of x and y don't match the row, column + * height of f + */ + + public PiecewiseBicubicSplineInterpolatingFunction(double[] x, double[] y, + double[][] f) + throws DimensionMismatchException, NullArgumentException, + NoDataException, NonMonotonicSequenceException { + + final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS; + + if (x == null || y == null || f == null || f[0] == null) { + throw new NullArgumentException(); + } + + final int xLen = x.length; + final int yLen = y.length; + + if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { + throw new NoDataException(); + } + + if (xLen < minimumLength || yLen < minimumLength || + f.length < minimumLength || f[0].length < minimumLength) { + throw new InsufficientDataException(); + } + + if (xLen != f.length) { + throw new DimensionMismatchException(xLen, f.length); + } + + if (yLen != f[0].length) { + throw new DimensionMismatchException(yLen, f[0].length); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + + xval = x.clone(); + yval = y.clone(); + fval = f.clone(); + } + + /** + * {@inheritDoc} + */ + public double value(double x, double y) + throws OutOfRangeException { + int index; + PolynomialSplineFunction spline; + AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator(); + final int offset = 2; + final int count = offset + 3; + final int i = searchIndex(x, xval, offset, count); + final int j = searchIndex(y, yval, offset, count); + + double xArray[] = new double[count]; + double yArray[] = new double[count]; + double zArray[] = new double[count]; + double interpArray[] = new double[count]; + + for (index = 0; index < count; index++) { + xArray[index] = xval[i + index]; + yArray[index] = yval[j + index]; + } + + for (int zIndex = 0; zIndex < count; zIndex++) { + for (index = 0; index < count; index++) { + zArray[index] = fval[i + index][j + zIndex]; + } + spline = interpolator.interpolate(xArray, zArray); + interpArray[zIndex] = spline.value(x); + } + + spline = interpolator.interpolate(yArray, interpArray); + + double returnValue = spline.value(y); + + return returnValue; + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @return {@code true} if (x, y) is a valid point. + * @since 3.3 + */ + public boolean isValidPoint(double x, double y) { + if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] || + y > yval[yval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @param c Coordinate. + * @param val Coordinate samples. + * @param offset how far back from found value to offset for querying + * @param count total number of elements forward from beginning that will be + * queried + * @return the index in {@code val} corresponding to the interval containing + * {@code c}. + * @throws OutOfRangeException if {@code c} is out of the range defined by + * the boundary values of {@code val}. + */ + private int searchIndex(double c, double[] val, int offset, int count) { + int r = Arrays.binarySearch(val, c); + + if (r == -1 || r == -val.length - 1) { + throw new OutOfRangeException(c, val[0], val[val.length - 1]); + } + + if (r < 0) { + // "c" in within an interpolation sub-interval, which returns + // negative + // need to remove the negative sign for consistency + r = -r - offset - 1; + } else { + r -= offset; + } + + if (r < 0) { + r = 0; + } + + if ((r + count) >= val.length) { + // "c" is the last sample of the range: Return the index + // of the sample at the lower end of the last sub-interval. + r = val.length - count; + } + + return r; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java new file mode 100644 index 000000000..2063f85aa --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java @@ -0,0 +1,64 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Generates a piecewise-bicubic interpolating function. + * + * @since 2.2 + */ +public class PiecewiseBicubicSplineInterpolator + implements BivariateGridInterpolator +{ + + /** + * Default constructor. + */ + public PiecewiseBicubicSplineInterpolator() + { + + } + + /** + * {@inheritDoc} + */ + public PiecewiseBicubicSplineInterpolatingFunction interpolate( final double[] xval, final double[] yval, + final double[][] fval) + throws DimensionMismatchException, NullArgumentException, NoDataException, NonMonotonicSequenceException + { + if ( xval == null || yval == null || fval == null || fval[0] == null ) + { + throw new NullArgumentException(); + } + + if ( xval.length == 0 || yval.length == 0 || fval.length == 0 ) + { + throw new NoDataException(); + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + + return new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, fval ); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java index 1c7603550..ad1e7b2f3 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java @@ -36,7 +36,7 @@ import org.apache.commons.math3.optim.SimpleVectorValueChecker; * @since 2.2 */ public class SmoothingPolynomialBicubicSplineInterpolator - extends BicubicSplineInterpolator { + extends PiecewiseBicubicSplineInterpolator { /** Fitter for x. */ private final PolynomialFitter xFitter; /** Degree of the fitting polynomial. */ @@ -92,7 +92,7 @@ public class SmoothingPolynomialBicubicSplineInterpolator * {@inheritDoc} */ @Override - public BicubicSplineInterpolatingFunction interpolate(final double[] xval, + public PiecewiseBicubicSplineInterpolatingFunction interpolate(final double[] xval, final double[] yval, final double[][] fval) throws NoDataException, NullArgumentException, diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java index 0faa27432..cda6a33dd 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java @@ -76,7 +76,7 @@ public class TricubicSplineInterpolator } } - final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(); + final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true); // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z final BicubicSplineInterpolatingFunction[] xSplineYZ @@ -103,13 +103,46 @@ public class TricubicSplineInterpolator final double[][][] dFdX = new double[xLen][yLen][zLen]; final double[][][] dFdY = new double[xLen][yLen][zLen]; final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; + for (int k = 0; k < zLen; k++) { + final BicubicSplineInterpolatingFunction f = zSplineXY[k]; + for (int i = 0; i < xLen; i++) { + final double x = xval[i]; + for (int j = 0; j < yLen; j++) { + final double y = yval[j]; + dFdX[i][j][k] = f.partialDerivativeX(x, y); + dFdY[i][j][k] = f.partialDerivativeY(x, y); + d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y); + } + } + } // Partial derivatives wrt y and wrt z final double[][][] dFdZ = new double[xLen][yLen][zLen]; final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; + for (int i = 0; i < xLen; i++) { + final BicubicSplineInterpolatingFunction f = xSplineYZ[i]; + for (int j = 0; j < yLen; j++) { + final double y = yval[j]; + for (int k = 0; k < zLen; k++) { + final double z = zval[k]; + dFdZ[i][j][k] = f.partialDerivativeY(y, z); + d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z); + } + } + } // Partial derivatives wrt x and wrt z final double[][][] d2FdZdX = new double[xLen][yLen][zLen]; + for (int j = 0; j < yLen; j++) { + final BicubicSplineInterpolatingFunction f = ySplineZX[j]; + for (int k = 0; k < zLen; k++) { + final double z = zval[k]; + for (int i = 0; i < xLen; i++) { + final double x = xval[i]; + d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x); + } + } + } // Third partial cross-derivatives final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java index 773a9488d..8c78aed00 100644 --- a/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java @@ -16,147 +16,435 @@ */ package org.apache.commons.math3.analysis.interpolation; -import static org.junit.Assert.assertEquals; -import static org.junit.Assert.assertTrue; - import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.InsufficientDataException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.OutOfRangeException; import org.apache.commons.math3.analysis.BivariateFunction; import org.apache.commons.math3.distribution.UniformRealDistribution; import org.apache.commons.math3.random.RandomGenerator; import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.Precision; import org.junit.Assert; import org.junit.Test; +import org.junit.Ignore; /** * Test case for the bicubic function. + * */ -public final class BicubicSplineInterpolatingFunctionTest -{ +public final class BicubicSplineInterpolatingFunctionTest { /** * Test preconditions. */ @Test - public void testPreconditions() - { - double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 }; - double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 }; + public void testPreconditions() { + double[] xval = new double[] {3, 4, 5, 6.5}; + double[] yval = new double[] {-4, -3, -1, 2.5}; double[][] zval = new double[xval.length][yval.length]; @SuppressWarnings("unused") - BicubicSplineInterpolatingFunction bcf = new BicubicSplineInterpolatingFunction( xval, yval, zval ); + BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, + zval, zval, zval); - try - { - bcf = new BicubicSplineInterpolatingFunction( null, yval, zval ); - Assert.fail( "Failed to detect x null pointer" ); + double[] wxval = new double[] {3, 2, 5, 6.5}; + try { + bcf = new BicubicSplineInterpolatingFunction(wxval, yval, zval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[] wyval = new double[] {-4, -1, -1, 2.5}; + try { + bcf = new BicubicSplineInterpolatingFunction(xval, wyval, zval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[][] wzval = new double[xval.length][yval.length - 1]; + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + + wzval = new double[xval.length - 1][yval.length]; + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected } - catch ( NullArgumentException iae ) - { - // Expected. } - try - { - bcf = new BicubicSplineInterpolatingFunction( xval, null, zval ); - Assert.fail( "Failed to detect y null pointer" ); - } - catch ( NullArgumentException iae ) - { - // Expected. - } - - try - { - bcf = new BicubicSplineInterpolatingFunction( xval, yval, null ); - Assert.fail( "Failed to detect z null pointer" ); - } - catch ( NullArgumentException iae ) - { - // Expected. - } - - try - { - double xval1[] = { 0.0, 1.0, 2.0, 3.0 }; - bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval ); - Assert.fail( "Failed to detect insufficient x data" ); + /** + * Test for a plane. + *

+ * z = 2 x - 3 y + 5 + */ + @Ignore@Test + public void testPlane() { + double[] xval = new double[] {3, 4, 5, 6.5}; + double[] yval = new double[] {-4, -3, -1, 2, 2.5}; + // Function values + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + return 2 * x - 3 * y + 5; + } + }; + double[][] zval = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + zval[i][j] = f.value(xval[i], yval[j]); + } + } + // Partial derivatives with respect to x + double[][] dZdX = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdX[i][j] = 2; + } + } + // Partial derivatives with respect to y + double[][] dZdY = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdY[i][j] = -3; + } + } + // Partial cross-derivatives + double[][] dZdXdY = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdXdY[i][j] = 0; } - catch ( InsufficientDataException iae ) - { - // Expected. } - try - { - double yval1[] = { 0.0, 1.0, 2.0, 3.0 }; - bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval ); - Assert.fail( "Failed to detect insufficient y data" ); - } - catch ( InsufficientDataException iae ) - { - // Expected. - } + BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, + dZdX, dZdY, dZdXdY); + double x, y; + double expected, result; - try - { - double zval1[][] = new double[4][4]; - bcf = new BicubicSplineInterpolatingFunction( xval, yval, zval1 ); - Assert.fail( "Failed to detect insufficient z data" ); - } - catch ( InsufficientDataException iae ) - { - // Expected. - } + x = 4; + y = -3; + expected = f.value(x, y); + result = bcf.value(x, y); + Assert.assertEquals("On sample point", + expected, result, 1e-15); - try - { - double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; - bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval ); - Assert.fail( "Failed to detect data set array with different sizes." ); - } - catch ( DimensionMismatchException iae ) - { - // Expected. + x = 4.5; + y = -1.5; + expected = f.value(x, y); + result = bcf.value(x, y); + Assert.assertEquals("Half-way between sample points (middle of the patch)", + expected, result, 0.3); + + x = 3.5; + y = -3.5; + expected = f.value(x, y); + result = bcf.value(x, y); + Assert.assertEquals("Half-way between sample points (border of the patch)", + expected, result, 0.3); } - try - { - double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; - bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval ); - Assert.fail( "Failed to detect data set array with different sizes." ); - } - catch ( DimensionMismatchException iae ) - { - // Expected. + /** + * Test for a paraboloid. + *

+ * z = 2 x2 - 3 y2 + 4 x y - 5 + */ + @Ignore@Test + public void testParaboloid() { + double[] xval = new double[] {3, 4, 5, 6.5}; + double[] yval = new double[] {-4, -3, -1, 2, 2.5}; + // Function values + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + return 2 * x * x - 3 * y * y + 4 * x * y - 5; } - - // X values not sorted. - try - { - double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 }; - bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval ); - Assert.fail( "Failed to detect unsorted x arguments." ); + }; + double[][] zval = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + zval[i][j] = f.value(xval[i], yval[j]); + } + } + // Partial derivatives with respect to x + double[][] dZdX = new double[xval.length][yval.length]; + BivariateFunction dfdX = new BivariateFunction() { + public double value(double x, double y) { + return 4 * (x + y); + } + }; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdX[i][j] = dfdX.value(xval[i], yval[j]); + } + } + // Partial derivatives with respect to y + double[][] dZdY = new double[xval.length][yval.length]; + BivariateFunction dfdY = new BivariateFunction() { + public double value(double x, double y) { + return 4 * x - 6 * y; + } + }; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdY[i][j] = dfdY.value(xval[i], yval[j]); + } + } + // Partial cross-derivatives + double[][] dZdXdY = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdXdY[i][j] = 4; } - catch ( NonMonotonicSequenceException iae ) - { - // Expected. } - // Y values not sorted. - try - { - double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 }; - bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval ); - Assert.fail( "Failed to detect unsorted y arguments." ); + BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, + dZdX, dZdY, dZdXdY); + double x, y; + double expected, result; + + x = 4; + y = -3; + expected = f.value(x, y); + result = bcf.value(x, y); + Assert.assertEquals("On sample point", + expected, result, 1e-15); + + x = 4.5; + y = -1.5; + expected = f.value(x, y); + result = bcf.value(x, y); + Assert.assertEquals("Half-way between sample points (middle of the patch)", + expected, result, 2); + + x = 3.5; + y = -3.5; + expected = f.value(x, y); + result = bcf.value(x, y); + Assert.assertEquals("Half-way between sample points (border of the patch)", + expected, result, 2); + } + + /** + * Test for partial derivatives of {@link BicubicSplineFunction}. + *

+ * f(x, y) = ΣiΣj (i+1) (j+2) xi yj + */ + @Ignore@Test + public void testSplinePartialDerivatives() { + final int N = 4; + final double[] coeff = new double[16]; + + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + coeff[i + N * j] = (i + 1) * (j + 2); + } + } + + final BicubicSplineFunction f = new BicubicSplineFunction(coeff); + BivariateFunction derivative; + final double x = 0.435; + final double y = 0.776; + final double tol = 1e-13; + + derivative = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double y2 = y * y; + final double y3 = y2 * y; + final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3; + return yFactor * (2 + 6 * x + 12 * x2); + } + }; + Assert.assertEquals("dFdX", derivative.value(x, y), + f.partialDerivativeX().value(x, y), tol); + + derivative = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double x3 = x2 * x; + final double y2 = y * y; + final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3; + return xFactor * (3 + 8 * y + 15 * y2); + } + }; + Assert.assertEquals("dFdY", derivative.value(x, y), + f.partialDerivativeY().value(x, y), tol); + + derivative = new BivariateFunction() { + public double value(double x, double y) { + final double y2 = y * y; + final double y3 = y2 * y; + final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3; + return yFactor * (6 + 24 * x); + } + }; + Assert.assertEquals("d2FdX2", derivative.value(x, y), + f.partialDerivativeXX().value(x, y), tol); + + derivative = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double x3 = x2 * x; + final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3; + return xFactor * (8 + 30 * y); + } + }; + Assert.assertEquals("d2FdY2", derivative.value(x, y), + f.partialDerivativeYY().value(x, y), tol); + + derivative = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double y2 = y * y; + final double yFactor = 3 + 8 * y + 15 * y2; + return yFactor * (2 + 6 * x + 12 * x2); + } + }; + Assert.assertEquals("d2FdXdY", derivative.value(x, y), + f.partialDerivativeXY().value(x, y), tol); + } + + /** + * Test that the partial derivatives computed from a + * {@link BicubicSplineInterpolatingFunction} match the input data. + *

+ * f(x, y) = 5 + * - 3 x + 2 y + * - x y + 2 x2 - 3 y2 + * + 4 x2 y - x y2 - 3 x3 + y3 + */ + @Ignore@Test + public void testMatchingPartialDerivatives() { + final int sz = 21; + double[] val = new double[sz]; + // Coordinate values + final double delta = 1d / (sz - 1); + for (int i = 0; i < sz; i++) { + val[i] = i * delta; + } + // Function values + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double x3 = x2 * x; + final double y2 = y * y; + final double y3 = y2 * y; + + return 5 + - 3 * x + 2 * y + - x * y + 2 * x2 - 3 * y2 + + 4 * x2 * y - x * y2 - 3 * x3 + y3; + } + }; + double[][] fval = new double[sz][sz]; + for (int i = 0; i < sz; i++) { + for (int j = 0; j < sz; j++) { + fval[i][j] = f.value(val[i], val[j]); + } + } + // Partial derivatives with respect to x + double[][] dFdX = new double[sz][sz]; + BivariateFunction dfdX = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double y2 = y * y; + return - 3 - y + 4 * x + 8 * x * y - y2 - 9 * x2; + } + }; + for (int i = 0; i < sz; i++) { + for (int j = 0; j < sz; j++) { + dFdX[i][j] = dfdX.value(val[i], val[j]); + } + } + // Partial derivatives with respect to y + double[][] dFdY = new double[sz][sz]; + BivariateFunction dfdY = new BivariateFunction() { + public double value(double x, double y) { + final double x2 = x * x; + final double y2 = y * y; + return 2 - x - 6 * y + 4 * x2 - 2 * x * y + 3 * y2; + } + }; + for (int i = 0; i < sz; i++) { + for (int j = 0; j < sz; j++) { + dFdY[i][j] = dfdY.value(val[i], val[j]); + } + } + // Partial cross-derivatives + double[][] d2FdXdY = new double[sz][sz]; + BivariateFunction d2fdXdY = new BivariateFunction() { + public double value(double x, double y) { + return -1 + 8 * x - 2 * y; + } + }; + for (int i = 0; i < sz; i++) { + for (int j = 0; j < sz; j++) { + d2FdXdY[i][j] = d2fdXdY.value(val[i], val[j]); + } + } + + BicubicSplineInterpolatingFunction bcf + = new BicubicSplineInterpolatingFunction(val, val, fval, dFdX, dFdY, d2FdXdY); + + double x, y; + double expected, result; + + final double tol = 1e-12; + for (int i = 0; i < sz; i++) { + x = val[i]; + for (int j = 0; j < sz; j++) { + y = val[j]; + + expected = dfdX.value(x, y); + result = bcf.partialDerivativeX(x, y); + Assert.assertEquals(x + " " + y + " dFdX", expected, result, tol); + + expected = dfdY.value(x, y); + result = bcf.partialDerivativeY(x, y); + Assert.assertEquals(x + " " + y + " dFdY", expected, result, tol); + + expected = d2fdXdY.value(x, y); + result = bcf.partialDerivativeXY(x, y); + Assert.assertEquals(x + " " + y + " d2FdXdY", expected, result, tol); } - catch ( NonMonotonicSequenceException iae ) - { - // Expected. } } @@ -166,28 +454,73 @@ public final class BicubicSplineInterpolatingFunctionTest * z = 2 x - 3 y + 5 */ @Test - public void testInterpolatePlane() - { - final int numberOfElements = 10; - final double minimumX = -10; - final double maximumX = 10; - final double minimumY = -10; - final double maximumY = 10; - final int numberOfSamples = 100; - final double interpolationTolerance = 7e-15; - final double maxTolerance = 6e-14; + public void testInterpolation1() { + final int sz = 21; + double[] xval = new double[sz]; + double[] yval = new double[sz]; + // Coordinate values + final double delta = 1d / (sz - 1); + for (int i = 0; i < sz; i++) { + xval[i] = -1 + 15 * i * delta; + yval[i] = -20 + 30 * i * delta; + } // Function values - BivariateFunction f = new BivariateFunction() - { - public double value( double x, double y ) - { + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { return 2 * x - 3 * y + 5; } }; + double[][] zval = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + zval[i][j] = f.value(xval[i], yval[j]); + } + } + // Partial derivatives with respect to x + double[][] dZdX = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdX[i][j] = 2; + } + } + // Partial derivatives with respect to y + double[][] dZdY = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdY[i][j] = -3; + } + } + // Partial cross-derivatives + double[][] dZdXdY = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdXdY[i][j] = 0; + } + } - testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f, - interpolationTolerance, maxTolerance ); + final BivariateFunction bcf + = new BicubicSplineInterpolatingFunction(xval, yval, zval, + dZdX, dZdY, dZdXdY); + double x, y; + + final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed. + final UniformRealDistribution distX + = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); + final UniformRealDistribution distY + = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); + + final int numSamples = 50; + final double tol = 6; + for (int i = 0; i < numSamples; i++) { + x = distX.sample(); + for (int j = 0; j < numSamples; j++) { + y = distY.sample(); +// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y)); + Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol); + } +// System.out.println(); + } } /** @@ -196,85 +529,140 @@ public final class BicubicSplineInterpolatingFunctionTest * z = 2 x2 - 3 y2 + 4 x y - 5 */ @Test - public void testInterpolationParabaloid() - { - final int numberOfElements = 10; - final double minimumX = -10; - final double maximumX = 10; - final double minimumY = -10; - final double maximumY = 10; - final int numberOfSamples = 100; - final double interpolationTolerance = 2e-14; - final double maxTolerance = 6e-14; + public void testInterpolation2() { + final int sz = 21; + double[] xval = new double[sz]; + double[] yval = new double[sz]; + // Coordinate values + final double delta = 1d / (sz - 1); + for (int i = 0; i < sz; i++) { + xval[i] = -1 + 15 * i * delta; + yval[i] = -20 + 30 * i * delta; + } // Function values - BivariateFunction f = new BivariateFunction() - { - public double value( double x, double y ) - { + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { return 2 * x * x - 3 * y * y + 4 * x * y - 5; } }; - - testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f, - interpolationTolerance, maxTolerance ); + double[][] zval = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + zval[i][j] = f.value(xval[i], yval[j]); + } } - - private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY, - int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance, - double maxTolerance ) - { - double expected; - double actual; - double currentX; - double currentY; - final double deltaX = ( maximumX - minimumX ) / ( (double) numberOfElements ); - final double deltaY = ( maximumY - minimumY ) / ( (double) numberOfElements ); - double xValues[] = new double[numberOfElements]; - double yValues[] = new double[numberOfElements]; - double zValues[][] = new double[numberOfElements][numberOfElements]; - - for ( int i = 0; i < numberOfElements; i++ ) - { - xValues[i] = minimumX + deltaX * (double) i; - for ( int j = 0; j < numberOfElements; j++ ) - { - yValues[j] = minimumY + deltaY * (double) j; - zValues[i][j] = f.value( xValues[i], yValues[j] ); + // Partial derivatives with respect to x + double[][] dZdX = new double[xval.length][yval.length]; + BivariateFunction dfdX = new BivariateFunction() { + public double value(double x, double y) { + return 4 * (x + y); + } + }; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdX[i][j] = dfdX.value(xval[i], yval[j]); + } + } + // Partial derivatives with respect to y + double[][] dZdY = new double[xval.length][yval.length]; + BivariateFunction dfdY = new BivariateFunction() { + public double value(double x, double y) { + return 4 * x - 6 * y; + } + }; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdY[i][j] = dfdY.value(xval[i], yval[j]); + } + } + // Partial cross-derivatives + double[][] dZdXdY = new double[xval.length][yval.length]; + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { + dZdXdY[i][j] = 4; } } - BivariateFunction interpolation = new BicubicSplineInterpolatingFunction( xValues, yValues, zValues ); + BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, + dZdX, dZdY, dZdXdY); + double x, y; - for ( int i = 0; i < numberOfElements; i++ ) - { - currentX = xValues[i]; - for ( int j = 0; j < numberOfElements; j++ ) - { - currentY = yValues[j]; - expected = f.value( currentX, currentY ); - actual = interpolation.value( currentX, currentY ); - assertTrue( Precision.equals( expected, actual ) ); + final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed. + final UniformRealDistribution distX + = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); + final UniformRealDistribution distY + = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); + + final double tol = 224; + for (int i = 0; i < sz; i++) { + x = distX.sample(); + for (int j = 0; j < sz; j++) { + y = distY.sample(); +// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y)); + Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol); } +// System.out.println(); } - - final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed. - final UniformRealDistribution distX = - new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] ); - final UniformRealDistribution distY = - new UniformRealDistribution( rng, yValues[0], yValues[yValues.length - 1] ); - - double sumError = 0; - for ( int i = 0; i < numberOfSamples; i++ ) - { - currentX = distX.sample(); - currentY = distY.sample(); - expected = f.value( currentX, currentY ); - actual = interpolation.value( currentX, currentY ); - sumError += FastMath.abs( actual - expected ); - assertEquals( expected, actual, maxTolerance ); } - assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance ); + @Test + public void testIsValidPoint() { + final double xMin = -12; + final double xMax = 34; + final double yMin = 5; + final double yMax = 67; + final double[] xval = new double[] { xMin, xMax }; + final double[] yval = new double[] { yMin, yMax }; + final double[][] f = new double[][] { { 1, 2 }, + { 3, 4 } }; + final double[][] dFdX = f; + final double[][] dFdY = f; + final double[][] dFdXdY = f; + + final BicubicSplineInterpolatingFunction bcf + = new BicubicSplineInterpolatingFunction(xval, yval, f, + dFdX, dFdY, dFdXdY); + + double x, y; + + x = xMin; + y = yMin; + Assert.assertTrue(bcf.isValidPoint(x, y)); + // Ensure that no exception is thrown. + bcf.value(x, y); + + x = xMax; + y = yMax; + Assert.assertTrue(bcf.isValidPoint(x, y)); + // Ensure that no exception is thrown. + bcf.value(x, y); + + final double xRange = xMax - xMin; + final double yRange = yMax - yMin; + x = xMin + xRange / 3.4; + y = yMin + yRange / 1.2; + Assert.assertTrue(bcf.isValidPoint(x, y)); + // Ensure that no exception is thrown. + bcf.value(x, y); + + final double small = 1e-8; + x = xMin - small; + y = yMax; + Assert.assertFalse(bcf.isValidPoint(x, y)); + // Ensure that an exception would have been thrown. + try { + bcf.value(x, y); + Assert.fail("OutOfRangeException expected"); + } catch (OutOfRangeException expected) {} + + x = xMin; + y = yMax + small; + Assert.assertFalse(bcf.isValidPoint(x, y)); + // Ensure that an exception would have been thrown. + try { + bcf.value(x, y); + Assert.fail("OutOfRangeException expected"); + } catch (OutOfRangeException expected) {} } } diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatorTest.java index 69d52a117..c4a56bbd5 100644 --- a/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatorTest.java @@ -17,10 +17,7 @@ package org.apache.commons.math3.analysis.interpolation; import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.InsufficientDataException; import org.apache.commons.math3.exception.MathIllegalArgumentException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NullArgumentException; import org.apache.commons.math3.analysis.BivariateFunction; import org.apache.commons.math3.distribution.UniformRealDistribution; import org.apache.commons.math3.random.RandomGenerator; @@ -30,131 +27,53 @@ import org.junit.Test; /** * Test case for the bicubic interpolator. + * */ -public final class BicubicSplineInterpolatorTest -{ +public final class BicubicSplineInterpolatorTest { /** * Test preconditions. */ @Test - public void testPreconditions() - { - double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 }; - double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 }; + public void testPreconditions() { + double[] xval = new double[] {3, 4, 5, 6.5}; + double[] yval = new double[] {-4, -3, -1, 2.5}; double[][] zval = new double[xval.length][yval.length]; - @SuppressWarnings( "unused" ) BivariateGridInterpolator interpolator = new BicubicSplineInterpolator(); - try - { - interpolator.interpolate( null, yval, zval ); - Assert.fail( "Failed to detect x null pointer" ); - } - catch ( NullArgumentException iae ) - { - // Expected. - } - - try - { - interpolator.interpolate( xval, null, zval ); - Assert.fail( "Failed to detect y null pointer" ); - } - catch ( NullArgumentException iae ) - { - // Expected. - } - - try - { - interpolator.interpolate( xval, yval, null ); - Assert.fail( "Failed to detect z null pointer" ); - } - catch ( NullArgumentException iae ) - { - // Expected. - } - - try - { - double xval1[] = { 0.0, 1.0, 2.0, 3.0 }; - interpolator.interpolate( xval1, yval, zval ); - Assert.fail( "Failed to detect insufficient x data" ); - } - catch ( InsufficientDataException iae ) - { - // Expected. - } - - try - { - double yval1[] = { 0.0, 1.0, 2.0, 3.0 }; - interpolator.interpolate( xval, yval1, zval ); - Assert.fail( "Failed to detect insufficient y data" ); - } - catch ( InsufficientDataException iae ) - { - // Expected. - } - - try - { - double zval1[][] = new double[4][4]; - interpolator.interpolate( xval, yval, zval1 ); - Assert.fail( "Failed to detect insufficient z data" ); - } - catch ( InsufficientDataException iae ) - { - // Expected. - } - - try - { - double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; - interpolator.interpolate( xval1, yval, zval ); - Assert.fail( "Failed to detect data set array with different sizes." ); - } - catch ( DimensionMismatchException iae ) - { - // Expected. - } + @SuppressWarnings("unused") + BivariateFunction p = interpolator.interpolate(xval, yval, zval); - try - { - double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; - interpolator.interpolate( xval, yval1, zval ); - Assert.fail( "Failed to detect data set array with different sizes." ); - } - catch ( DimensionMismatchException iae ) - { - // Expected. + double[] wxval = new double[] {3, 2, 5, 6.5}; + try { + p = interpolator.interpolate(wxval, yval, zval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected } - // X values not sorted. - try - { - double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 }; - interpolator.interpolate( xval1, yval, zval ); - Assert.fail( "Failed to detect unsorted x arguments." ); - } - catch ( NonMonotonicSequenceException iae ) - { - // Expected. + double[] wyval = new double[] {-4, -3, -1, -1}; + try { + p = interpolator.interpolate(xval, wyval, zval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected } - // Y values not sorted. - try - { - double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 }; - interpolator.interpolate( xval, yval1, zval ); - Assert.fail( "Failed to detect unsorted y arguments." ); + double[][] wzval = new double[xval.length][yval.length + 1]; + try { + p = interpolator.interpolate(xval, yval, wzval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected } - catch ( NonMonotonicSequenceException iae ) - { - // Expected. + wzval = new double[xval.length - 1][yval.length]; + try { + p = interpolator.interpolate(xval, yval, wzval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected } - } /** @@ -163,32 +82,26 @@ public final class BicubicSplineInterpolatorTest * z = 2 x - 3 y + 5 */ @Test - public void testInterpolation1() - { + public void testInterpolation1() { final int sz = 21; double[] xval = new double[sz]; double[] yval = new double[sz]; // Coordinate values final double delta = 1d / (sz - 1); - for ( int i = 0; i < sz; i++ ) - { + for (int i = 0; i < sz; i++) { xval[i] = -1 + 15 * i * delta; yval[i] = -20 + 30 * i * delta; } // Function values - BivariateFunction f = new BivariateFunction() - { - public double value( double x, double y ) - { + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { return 2 * x - 3 * y + 5; } }; double[][] zval = new double[xval.length][yval.length]; - for ( int i = 0; i < xval.length; i++ ) - { - for ( int j = 0; j < yval.length; j++ ) - { + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { zval[i][j] = f.value(xval[i], yval[j]); } } @@ -198,16 +111,16 @@ public final class BicubicSplineInterpolatorTest double x, y; final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed. - final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] ); - final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] ); + final UniformRealDistribution distX + = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); + final UniformRealDistribution distY + = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); final int numSamples = 50; - final double tol = 2e-14; - for ( int i = 0; i < numSamples; i++ ) - { + final double tol = 6; + for (int i = 0; i < numSamples; i++) { x = distX.sample(); - for ( int j = 0; j < numSamples; j++ ) - { + for (int j = 0; j < numSamples; j++) { y = distY.sample(); // System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y)); Assert.assertEquals(f.value(x, y), p.value(x, y), tol); @@ -222,32 +135,26 @@ public final class BicubicSplineInterpolatorTest * z = 2 x2 - 3 y2 + 4 x y - 5 */ @Test - public void testInterpolation2() - { + public void testInterpolation2() { final int sz = 21; double[] xval = new double[sz]; double[] yval = new double[sz]; // Coordinate values final double delta = 1d / (sz - 1); - for ( int i = 0; i < sz; i++ ) - { + for (int i = 0; i < sz; i++) { xval[i] = -1 + 15 * i * delta; yval[i] = -20 + 30 * i * delta; } // Function values - BivariateFunction f = new BivariateFunction() - { - public double value( double x, double y ) - { + BivariateFunction f = new BivariateFunction() { + public double value(double x, double y) { return 2 * x * x - 3 * y * y + 4 * x * y - 5; } }; double[][] zval = new double[xval.length][yval.length]; - for ( int i = 0; i < xval.length; i++ ) - { - for ( int j = 0; j < yval.length; j++ ) - { + for (int i = 0; i < xval.length; i++) { + for (int j = 0; j < yval.length; j++) { zval[i][j] = f.value(xval[i], yval[j]); } } @@ -257,16 +164,16 @@ public final class BicubicSplineInterpolatorTest double x, y; final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed. - final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] ); - final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] ); + final UniformRealDistribution distX + = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]); + final UniformRealDistribution distY + = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]); final int numSamples = 50; - final double tol = 5e-13; - for ( int i = 0; i < numSamples; i++ ) - { + final double tol = 251; + for (int i = 0; i < numSamples; i++) { x = distX.sample(); - for ( int j = 0; j < numSamples; j++ ) - { + for (int j = 0; j < numSamples; j++) { y = distY.sample(); // System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y)); Assert.assertEquals(f.value(x, y), p.value(x, y), tol); diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunctionTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunctionTest.java new file mode 100644 index 000000000..10f7822e3 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunctionTest.java @@ -0,0 +1,280 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertTrue; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.InsufficientDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; +import org.junit.Assert; +import org.junit.Test; + +/** + * Test case for the piecewise bicubic function. + */ +public final class PiecewiseBicubicSplineInterpolatingFunctionTest +{ + /** + * Test preconditions. + */ + @Test + public void testPreconditions() + { + double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 }; + double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 }; + double[][] zval = new double[xval.length][yval.length]; + + @SuppressWarnings("unused") + PiecewiseBicubicSplineInterpolatingFunction bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, zval ); + + try + { + bcf = new PiecewiseBicubicSplineInterpolatingFunction( null, yval, zval ); + Assert.fail( "Failed to detect x null pointer" ); + } + catch ( NullArgumentException iae ) + { + // Expected. + } + + try + { + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, null, zval ); + Assert.fail( "Failed to detect y null pointer" ); + } + catch ( NullArgumentException iae ) + { + // Expected. + } + + try + { + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, null ); + Assert.fail( "Failed to detect z null pointer" ); + } + catch ( NullArgumentException iae ) + { + // Expected. + } + + try + { + double xval1[] = { 0.0, 1.0, 2.0, 3.0 }; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval ); + Assert.fail( "Failed to detect insufficient x data" ); + } + catch ( InsufficientDataException iae ) + { + // Expected. + } + + try + { + double yval1[] = { 0.0, 1.0, 2.0, 3.0 }; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval ); + Assert.fail( "Failed to detect insufficient y data" ); + } + catch ( InsufficientDataException iae ) + { + // Expected. + } + + try + { + double zval1[][] = new double[4][4]; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, zval1 ); + Assert.fail( "Failed to detect insufficient z data" ); + } + catch ( InsufficientDataException iae ) + { + // Expected. + } + + try + { + double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval ); + Assert.fail( "Failed to detect data set array with different sizes." ); + } + catch ( DimensionMismatchException iae ) + { + // Expected. + } + + try + { + double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval ); + Assert.fail( "Failed to detect data set array with different sizes." ); + } + catch ( DimensionMismatchException iae ) + { + // Expected. + } + + // X values not sorted. + try + { + double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 }; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval ); + Assert.fail( "Failed to detect unsorted x arguments." ); + } + catch ( NonMonotonicSequenceException iae ) + { + // Expected. + } + + // Y values not sorted. + try + { + double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 }; + bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval ); + Assert.fail( "Failed to detect unsorted y arguments." ); + } + catch ( NonMonotonicSequenceException iae ) + { + // Expected. + } + } + + /** + * Interpolating a plane. + *

+ * z = 2 x - 3 y + 5 + */ + @Test + public void testInterpolatePlane() + { + final int numberOfElements = 10; + final double minimumX = -10; + final double maximumX = 10; + final double minimumY = -10; + final double maximumY = 10; + final int numberOfSamples = 100; + final double interpolationTolerance = 7e-15; + final double maxTolerance = 6e-14; + + // Function values + BivariateFunction f = new BivariateFunction() + { + public double value( double x, double y ) + { + return 2 * x - 3 * y + 5; + } + }; + + testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f, + interpolationTolerance, maxTolerance ); + } + + /** + * Interpolating a paraboloid. + *

+ * z = 2 x2 - 3 y2 + 4 x y - 5 + */ + @Test + public void testInterpolationParabaloid() + { + final int numberOfElements = 10; + final double minimumX = -10; + final double maximumX = 10; + final double minimumY = -10; + final double maximumY = 10; + final int numberOfSamples = 100; + final double interpolationTolerance = 2e-14; + final double maxTolerance = 6e-14; + + // Function values + BivariateFunction f = new BivariateFunction() + { + public double value( double x, double y ) + { + return 2 * x * x - 3 * y * y + 4 * x * y - 5; + } + }; + + testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f, + interpolationTolerance, maxTolerance ); + } + + private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY, + int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance, + double maxTolerance ) + { + double expected; + double actual; + double currentX; + double currentY; + final double deltaX = ( maximumX - minimumX ) / ( (double) numberOfElements ); + final double deltaY = ( maximumY - minimumY ) / ( (double) numberOfElements ); + double xValues[] = new double[numberOfElements]; + double yValues[] = new double[numberOfElements]; + double zValues[][] = new double[numberOfElements][numberOfElements]; + + for ( int i = 0; i < numberOfElements; i++ ) + { + xValues[i] = minimumX + deltaX * (double) i; + for ( int j = 0; j < numberOfElements; j++ ) + { + yValues[j] = minimumY + deltaY * (double) j; + zValues[i][j] = f.value( xValues[i], yValues[j] ); + } + } + + BivariateFunction interpolation = new PiecewiseBicubicSplineInterpolatingFunction( xValues, yValues, zValues ); + + for ( int i = 0; i < numberOfElements; i++ ) + { + currentX = xValues[i]; + for ( int j = 0; j < numberOfElements; j++ ) + { + currentY = yValues[j]; + expected = f.value( currentX, currentY ); + actual = interpolation.value( currentX, currentY ); + assertTrue( Precision.equals( expected, actual ) ); + } + } + + final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed. + final UniformRealDistribution distX = + new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] ); + final UniformRealDistribution distY = + new UniformRealDistribution( rng, yValues[0], yValues[yValues.length - 1] ); + + double sumError = 0; + for ( int i = 0; i < numberOfSamples; i++ ) + { + currentX = distX.sample(); + currentY = distY.sample(); + expected = f.value( currentX, currentY ); + actual = interpolation.value( currentX, currentY ); + sumError += FastMath.abs( actual - expected ); + assertEquals( expected, actual, maxTolerance ); + } + + assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance ); + } +} diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatorTest.java new file mode 100644 index 000000000..50d41b78a --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatorTest.java @@ -0,0 +1,277 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.InsufficientDataException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.junit.Assert; +import org.junit.Test; + +/** + * Test case for the piecewise bicubic interpolator. + */ +public final class PiecewiseBicubicSplineInterpolatorTest +{ + /** + * Test preconditions. + */ + @Test + public void testPreconditions() + { + double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 }; + double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 }; + double[][] zval = new double[xval.length][yval.length]; + + @SuppressWarnings( "unused" ) + BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator(); + + try + { + interpolator.interpolate( null, yval, zval ); + Assert.fail( "Failed to detect x null pointer" ); + } + catch ( NullArgumentException iae ) + { + // Expected. + } + + try + { + interpolator.interpolate( xval, null, zval ); + Assert.fail( "Failed to detect y null pointer" ); + } + catch ( NullArgumentException iae ) + { + // Expected. + } + + try + { + interpolator.interpolate( xval, yval, null ); + Assert.fail( "Failed to detect z null pointer" ); + } + catch ( NullArgumentException iae ) + { + // Expected. + } + + try + { + double xval1[] = { 0.0, 1.0, 2.0, 3.0 }; + interpolator.interpolate( xval1, yval, zval ); + Assert.fail( "Failed to detect insufficient x data" ); + } + catch ( InsufficientDataException iae ) + { + // Expected. + } + + try + { + double yval1[] = { 0.0, 1.0, 2.0, 3.0 }; + interpolator.interpolate( xval, yval1, zval ); + Assert.fail( "Failed to detect insufficient y data" ); + } + catch ( InsufficientDataException iae ) + { + // Expected. + } + + try + { + double zval1[][] = new double[4][4]; + interpolator.interpolate( xval, yval, zval1 ); + Assert.fail( "Failed to detect insufficient z data" ); + } + catch ( InsufficientDataException iae ) + { + // Expected. + } + + try + { + double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + interpolator.interpolate( xval1, yval, zval ); + Assert.fail( "Failed to detect data set array with different sizes." ); + } + catch ( DimensionMismatchException iae ) + { + // Expected. + } + + try + { + double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + interpolator.interpolate( xval, yval1, zval ); + Assert.fail( "Failed to detect data set array with different sizes." ); + } + catch ( DimensionMismatchException iae ) + { + // Expected. + } + + // X values not sorted. + try + { + double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 }; + interpolator.interpolate( xval1, yval, zval ); + Assert.fail( "Failed to detect unsorted x arguments." ); + } + catch ( NonMonotonicSequenceException iae ) + { + // Expected. + } + + // Y values not sorted. + try + { + double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 }; + interpolator.interpolate( xval, yval1, zval ); + Assert.fail( "Failed to detect unsorted y arguments." ); + } + catch ( NonMonotonicSequenceException iae ) + { + // Expected. + } + + } + + /** + * Interpolating a plane. + *

+ * z = 2 x - 3 y + 5 + */ + @Test + public void testInterpolation1() + { + final int sz = 21; + double[] xval = new double[sz]; + double[] yval = new double[sz]; + // Coordinate values + final double delta = 1d / (sz - 1); + for ( int i = 0; i < sz; i++ ) + { + xval[i] = -1 + 15 * i * delta; + yval[i] = -20 + 30 * i * delta; + } + + // Function values + BivariateFunction f = new BivariateFunction() + { + public double value( double x, double y ) + { + return 2 * x - 3 * y + 5; + } + }; + double[][] zval = new double[xval.length][yval.length]; + for ( int i = 0; i < xval.length; i++ ) + { + for ( int j = 0; j < yval.length; j++ ) + { + zval[i][j] = f.value(xval[i], yval[j]); + } + } + + BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator(); + BivariateFunction p = interpolator.interpolate(xval, yval, zval); + double x, y; + + final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed. + final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] ); + final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] ); + + final int numSamples = 50; + final double tol = 2e-14; + for ( int i = 0; i < numSamples; i++ ) + { + x = distX.sample(); + for ( int j = 0; j < numSamples; j++ ) + { + y = distY.sample(); +// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y)); + Assert.assertEquals(f.value(x, y), p.value(x, y), tol); + } +// System.out.println(); + } + } + + /** + * Interpolating a paraboloid. + *

+ * z = 2 x2 - 3 y2 + 4 x y - 5 + */ + @Test + public void testInterpolation2() + { + final int sz = 21; + double[] xval = new double[sz]; + double[] yval = new double[sz]; + // Coordinate values + final double delta = 1d / (sz - 1); + for ( int i = 0; i < sz; i++ ) + { + xval[i] = -1 + 15 * i * delta; + yval[i] = -20 + 30 * i * delta; + } + + // Function values + BivariateFunction f = new BivariateFunction() + { + public double value( double x, double y ) + { + return 2 * x * x - 3 * y * y + 4 * x * y - 5; + } + }; + double[][] zval = new double[xval.length][yval.length]; + for ( int i = 0; i < xval.length; i++ ) + { + for ( int j = 0; j < yval.length; j++ ) + { + zval[i][j] = f.value(xval[i], yval[j]); + } + } + + BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator(); + BivariateFunction p = interpolator.interpolate(xval, yval, zval); + double x, y; + + final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed. + final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] ); + final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] ); + + final int numSamples = 50; + final double tol = 5e-13; + for ( int i = 0; i < numSamples; i++ ) + { + x = distX.sample(); + for ( int j = 0; j < numSamples; j++ ) + { + y = distY.sample(); +// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y)); + Assert.assertEquals(f.value(x, y), p.value(x, y), tol); + } +// System.out.println(); + } + } +} diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatorTest.java index be5a86592..6de42a282 100644 --- a/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatorTest.java @@ -32,7 +32,7 @@ public final class TricubicSplineInterpolatorTest { /** * Test preconditions. */ - @Ignore@Test + @Test public void testPreconditions() { double[] xval = new double[] {3, 4, 5, 6.5}; double[] yval = new double[] {-4, -3, -1, 2.5};