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

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 <a
+ * Function that implements the
- * href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline
+ * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
- * interpolation</a>.
+ * bicubic spline interpolation</a>.
  *
  * @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
+     * @param f Values of the function on every grid point.
-     *        of elements.
+     * @param dFdX Values of the partial derivative of function with respect
-     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not
+     * to x on every grid point.
-     *         strictly increasing.
+     * @param dFdY Values of the partial derivative of function with respect
-     * @throws NullArgumentException if any of the arguments are null
+     * 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)
+     * @param x Sample values of the x-coordinate, in increasing order.
-        throws DimensionMismatchException, NullArgumentException,
+     * @param y Sample values of the y-coordinate, in increasing order.
-        NoDataException, NonMonotonicSequenceException {
+     * @param f Values of the function on every grid point.
-
+     * @param dFdX Values of the partial derivative of function with respect
-        final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS;
+     * to x on every grid point.
-
+     * @param dFdY Values of the partial derivative of function with respect
-        if (x == null || y == null || f == null || f[0] == null) {
+     * to y on every grid point.
-            throw new NullArgumentException();
+     * @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 (xLen != dFdX.length) {
-        if (yLen != f[0].length) {
+            throw new DimensionMismatchException(xLen, dFdX.length);
-            throw new DimensionMismatchException(yLen, f[0].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;
+        final int i = searchIndex(x, xval);
-        PolynomialSplineFunction spline;
+        final int j = searchIndex(y, yval);
-        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];
+        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
-        double yArray[] = new double[count];
+        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
-        double zArray[] = new double[count];
-        double interpArray[] = new double[count];
 
-        for (index = 0; index < count; index++) {
+        return splines[i][j].value(xN, yN);
-            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;
     }
 
     /**
@@ -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
+     * @return the index in {@code val} corresponding to the interval
-     * @param count total number of elements forward from beginning that will be
+     * containing {@code c}.
-     *        queried
+     * @throws OutOfRangeException if {@code c} is out of the
-     * @return the index in {@code val} corresponding to the interval containing
+     * range defined by the boundary values of {@code val}.
-     *         {@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) {
+    private int searchIndex(double c, double[] val) {
-        int r = Arrays.binarySearch(val, c);
+        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
+            // "c" in within an interpolation sub-interval: Return the
-            // negative
+            // index of the sample at the lower end of the sub-interval.
-            // need to remove the negative sign for consistency
+            return -r - 2;
-            r = -r - offset - 1;
-        } else {
-            r -= offset;
         }
-
+        final int last = val.length - 1;
-        if (r < 0) {
+        if (r == last) {
-            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 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:
+     * <ul>
+     *  <li>f(0,0)</li>
+     *  <li>f(1,0)</li>
+     *  <li>f(0,1)</li>
+     *  <li>f(1,1)</li>
+     *  <li>f<sub>x</sub>(0,0)</li>
+     *  <li>f<sub>x</sub>(1,0)</li>
+     *  <li>f<sub>x</sub>(0,1)</li>
+     *  <li>f<sub>x</sub>(1,1)</li>
+     *  <li>f<sub>y</sub>(0,0)</li>
+     *  <li>f<sub>y</sub>(1,0)</li>
+     *  <li>f<sub>y</sub>(0,1)</li>
+     *  <li>f<sub>y</sub>(1,1)</li>
+     *  <li>f<sub>xy</sub>(0,0)</li>
+     *  <li>f<sub>xy</sub>(1,0)</li>
+     *  <li>f<sub>xy</sub>(0,1)</li>
+     *  <li>f<sub>xy</sub>(1,1)</li>
+     * </ul>
+     * 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;
+    }
 }

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
+        throws NoDataException, DimensionMismatchException,
-    {
+               NonMonotonicSequenceException, NumberIsTooSmallException {
-        if ( xval == null || yval == null || fval == null || fval[0] == null )
+        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
-        {
-            throw new NullArgumentException();
-        }
-
-        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;
     }
 }

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 <a
+ * href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline
+ * interpolation</a>.
+ *
+ * @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;
+    }
+}

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 );
+    }
+}

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,

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];

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.MathIllegalArgumentException;
-import org.apache.commons.math3.exception.NonMonotonicSequenceException;
+import org.apache.commons.math3.exception.OutOfRangeException;
-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;
+import org.junit.Ignore;
 
 /**
  * Test case for the bicubic function.
+ * 
  */
-public final class BicubicSplineInterpolatingFunctionTest
+public final class BicubicSplineInterpolatingFunctionTest {
-{
     /**
      * Test preconditions.
      */
     @Test
-    public void testPreconditions()
+    public void testPreconditions() {
-    {
+        double[] xval = new double[] {3, 4, 5, 6.5};
-        double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
+        double[] yval = new double[] {-4, -3, -1, 2.5};
-        double[] yval = new double[] { -4, -3, -1, 2.5, 3.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
+        double[] wxval = new double[] {3, 2, 5, 6.5};
-        {
+        try {
-            bcf = new BicubicSplineInterpolatingFunction( null, yval, zval );
+            bcf = new BicubicSplineInterpolatingFunction(wxval, yval, zval, zval, zval, zval);
-            Assert.fail( "Failed to detect x null pointer" );
+            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
+    /**
-        {
+     * Test for a plane.
-            bcf = new BicubicSplineInterpolatingFunction( xval, null, zval );
+     * <p>
-            Assert.fail( "Failed to detect y null pointer" );
+     * z = 2 x - 3 y + 5
-        }
+     */
-        catch ( NullArgumentException iae )
+    @Ignore@Test
-        {
+    public void testPlane() {
-            // Expected.
+        double[] xval = new double[] {3, 4, 5, 6.5};
-    }
+        double[] yval = new double[] {-4, -3, -1, 2, 2.5};
-
+        // Function values
-        try
+        BivariateFunction f = new BivariateFunction() {
-        {
+                public double value(double x, double y) {
-            bcf = new BicubicSplineInterpolatingFunction( xval, yval, null );
+                    return 2 * x - 3 * y + 5;
-            Assert.fail( "Failed to detect z null pointer" );
+                }
-        }
+            };
-        catch ( NullArgumentException iae )
+        double[][] zval = new double[xval.length][yval.length];
-        {
+        for (int i = 0; i < xval.length; i++) {
-            // Expected.
+            for (int j = 0; j < yval.length; j++) {
-    }
+                zval[i][j] = f.value(xval[i], yval[j]);
-
+            }
-        try
+        }
-        {
+        // Partial derivatives with respect to x
-            double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
+        double[][] dZdX = new double[xval.length][yval.length];
-            bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
+        for (int i = 0; i < xval.length; i++) {
-            Assert.fail( "Failed to detect insufficient x data" );
+            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
+        BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
-        {
+                                                                           dZdX, dZdY, dZdXdY);
-            double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
+        double x, y;
-            bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
+        double expected, result;
-            Assert.fail( "Failed to detect insufficient y data" );
-                }
-        catch ( InsufficientDataException iae )
-        {
-            // Expected.
-                }
 
-        try
+        x = 4;
-        {
+        y = -3;
-            double zval1[][] = new double[4][4];
+        expected = f.value(x, y);
-            bcf = new BicubicSplineInterpolatingFunction( xval, yval, zval1 );
+        result = bcf.value(x, y);
-            Assert.fail( "Failed to detect insufficient z data" );
+        Assert.assertEquals("On sample point",
-                }
+                            expected, result, 1e-15);
-        catch ( InsufficientDataException iae )
-        {
-            // Expected.
-                }
 
-        try
+        x = 4.5;
-        {
+        y = -1.5;
-            double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+        expected = f.value(x, y);
-            bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
+        result = bcf.value(x, y);
-            Assert.fail( "Failed to detect data set array with different sizes." );
+        Assert.assertEquals("Half-way between sample points (middle of the patch)",
-                }
+                            expected, result, 0.3);
-        catch ( DimensionMismatchException iae )
+
-        {
+        x = 3.5;
-            // Expected.
+        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
+    /**
-        {
+     * Test for a paraboloid.
-            double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+     * <p>
-            bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
+     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
-            Assert.fail( "Failed to detect data set array with different sizes." );
+     */
-        }
+    @Ignore@Test
-        catch ( DimensionMismatchException iae )
+    public void testParaboloid() {
-        {
+        double[] xval = new double[] {3, 4, 5, 6.5};
-            // Expected.
+        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.
+        double[][] zval = new double[xval.length][yval.length];
-        try
+        for (int i = 0; i < xval.length; i++) {
-        {
+            for (int j = 0; j < yval.length; j++) {
-            double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
+                zval[i][j] = f.value(xval[i], yval[j]);
-            bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
+            }
-            Assert.fail( "Failed to detect unsorted x arguments." );
+        }
+        // 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.
+        BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
-        try
+                                                                           dZdX, dZdY, dZdXdY);
-        {
+        double x, y;
-            double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
+        double expected, result;
-            bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
+        
-            Assert.fail( "Failed to detect unsorted y arguments." );
+        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}.
+     * <p>
+     * f(x, y) = &Sigma;<sub>i</sub>&Sigma;<sub>j</sub> (i+1) (j+2) x<sup>i</sup> y<sup>j</sup>
+     */
+    @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.
+     * <p>
+     * f(x, y) = 5
+     *           - 3 x + 2 y
+     *           - x y + 2 x<sup>2</sup> - 3 y<sup>2</sup>
+     *           + 4 x<sup>2</sup> y - x y<sup>2</sup> - 3 x<sup>3</sup> + y<sup>3</sup>
+     */
+    @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()
+    public void testInterpolation1() {
-    {
+        final int sz = 21;
-        final int numberOfElements = 10;
+        double[] xval = new double[sz];
-        final double minimumX = -10;
+        double[] yval = new double[sz];
-        final double maximumX = 10;
+        // Coordinate values
-        final double minimumY = -10;
+        final double delta = 1d / (sz - 1);
-        final double maximumY = 10;
+        for (int i = 0; i < sz; i++) {
-        final int numberOfSamples = 100;
+            xval[i] = -1 + 15 * i * delta;
-        final double interpolationTolerance = 7e-15;
+            yval[i] = -20 + 30 * i * delta;
-        final double maxTolerance = 6e-14;
+        }
 
         // Function values
-        BivariateFunction f = new BivariateFunction()
+        BivariateFunction f = new BivariateFunction() {
-        {
+                public double value(double x, double y) {
-            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,
+        final BivariateFunction bcf
-                           interpolationTolerance, maxTolerance );
+            = 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 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
      */
     @Test
-    public void testInterpolationParabaloid()
+    public void testInterpolation2() {
-    {
+        final int sz = 21;
-        final int numberOfElements = 10;
+        double[] xval = new double[sz];
-        final double minimumX = -10;
+        double[] yval = new double[sz];
-        final double maximumX = 10;
+        // Coordinate values
-        final double minimumY = -10;
+        final double delta = 1d / (sz - 1);
-        final double maximumY = 10;
+        for (int i = 0; i < sz; i++) {
-        final int numberOfSamples = 100;
+            xval[i] = -1 + 15 * i * delta;
-        final double interpolationTolerance = 2e-14;
+            yval[i] = -20 + 30 * i * delta;
-        final double maxTolerance = 6e-14;
+        }
 
         // Function values
-        BivariateFunction f = new BivariateFunction()
+        BivariateFunction f = new BivariateFunction() {
-        {
+                public double value(double x, double y) {
-            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];
-        testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f,
+        for (int i = 0; i < xval.length; i++) {
-                           interpolationTolerance, maxTolerance );
+            for (int j = 0; j < yval.length; j++) {
+                zval[i][j] = f.value(xval[i], yval[j]);
+            }
         }
-
+        // Partial derivatives with respect to x
-    private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY,
+        double[][] dZdX = new double[xval.length][yval.length];
-                                    int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance,
+        BivariateFunction dfdX = new BivariateFunction() {
-                                    double maxTolerance )
+                public double value(double x, double y) {
-    {
+                    return 4 * (x + y);
-        double expected;
+                }
-        double actual;
+            };
-        double currentX;
+        for (int i = 0; i < xval.length; i++) {
-        double currentY;
+            for (int j = 0; j < yval.length; j++) {
-        final double deltaX = ( maximumX - minimumX ) / ( (double) numberOfElements );
+                dZdX[i][j] = dfdX.value(xval[i], yval[j]);
-        final double deltaY = ( maximumY - minimumY ) / ( (double) numberOfElements );
+            }
-        double xValues[] = new double[numberOfElements];
+        }
-        double yValues[] = new double[numberOfElements];
+        // Partial derivatives with respect to y
-        double zValues[][] = new double[numberOfElements][numberOfElements];
+        double[][] dZdY = new double[xval.length][yval.length];
-
+        BivariateFunction dfdY = new BivariateFunction() {
-        for ( int i = 0; i < numberOfElements; i++ )
+                public double value(double x, double y) {
-        {
+                    return 4 * x - 6 * y;
-            xValues[i] = minimumX + deltaX * (double) i;
+                }
-            for ( int j = 0; j < numberOfElements; j++ )
+            };
-            {
+        for (int i = 0; i < xval.length; i++) {
-                yValues[j] = minimumY + deltaY * (double) j;
+            for (int j = 0; j < yval.length; j++) {
-                zValues[i][j] = f.value( xValues[i], yValues[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++ )
+        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
-        {
+        final UniformRealDistribution distX
-            currentX = xValues[i];
+            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
-            for ( int j = 0; j < numberOfElements; j++ )
+        final UniformRealDistribution distY
-            {
+            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
-                currentY = yValues[j];
+
-                expected = f.value( currentX, currentY );
+        final double tol = 224;
-                actual = interpolation.value( currentX, currentY );
+        for (int i = 0; i < sz; i++) {
-                assertTrue( Precision.equals( expected, actual ) );
+            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) {}
     }
 }

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()
+    public void testPreconditions() {
-    {
+        double[] xval = new double[] {3, 4, 5, 6.5};
-        double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
+        double[] yval = new double[] {-4, -3, -1, 2.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 BicubicSplineInterpolator();
         
-        try
+        @SuppressWarnings("unused")
-        {
+        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
-            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[] wxval = new double[] {3, 2, 5, 6.5};
-        {
+        try {
-            double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
+            p = interpolator.interpolate(wxval, yval, zval);
-            interpolator.interpolate( xval, yval1, zval );
+            Assert.fail("an exception should have been thrown");
-            Assert.fail( "Failed to detect data set array with different sizes." );
+        } catch (MathIllegalArgumentException e) {
-        }
+            // Expected
-        catch ( DimensionMismatchException iae )
-        {
-            // Expected.
         }
 
-        // X values not sorted.
+        double[] wyval = new double[] {-4, -3, -1, -1};
-        try
+        try {
-        {
+            p = interpolator.interpolate(xval, wyval, zval);
-            double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
+            Assert.fail("an exception should have been thrown");
-            interpolator.interpolate( xval1, yval, zval );
+        } catch (MathIllegalArgumentException e) {
-            Assert.fail( "Failed to detect unsorted x arguments." );
+            // Expected
-        }
-        catch ( NonMonotonicSequenceException iae )
-        {
-            // Expected.
         }
 
-        // Y values not sorted.
+        double[][] wzval = new double[xval.length][yval.length + 1];
-        try
+        try {
-        {
+            p = interpolator.interpolate(xval, yval, wzval);
-            double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
+            Assert.fail("an exception should have been thrown");
-            interpolator.interpolate( xval, yval1, zval );
+        } catch (DimensionMismatchException e) {
-            Assert.fail( "Failed to detect unsorted y arguments." );
+            // Expected
         }
-        catch ( NonMonotonicSequenceException iae )
+        wzval = new double[xval.length - 1][yval.length];
-        {
+        try {
-            // Expected.
+            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()
+        BivariateFunction f = new BivariateFunction() {
-        {
+                public double value(double x, double y) {
-            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 i = 0; i < xval.length; i++) {
-        {
+            for (int j = 0; j < yval.length; j++) {
-            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 distX
-        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );
+            = 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;
+        final double tol = 6;
-        for ( int i = 0; i < numSamples; i++ )
+        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 x<sup>2</sup> - 3 y<sup>2</sup> + 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()
+        BivariateFunction f = new BivariateFunction() {
-        {
+                public double value(double x, double y) {
-            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 i = 0; i < xval.length; i++) {
-        {
+            for (int j = 0; j < yval.length; j++) {
-            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 distX
-        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );
+            = 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;
+        final double tol = 251;
-        for ( int i = 0; i < numSamples; i++ )
+        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);

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.
+     * <p>
+     * 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.
+     * <p>
+     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 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 );
+    }
+}

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.
+     * <p>
+     * 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.
+     * <p>
+     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 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();
+        }
+    }
+}

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};