MATH-1648: Derivatives for "BicubicInterpolatingFunction".

commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunction.java

@@ -17,6 +17,8 @@
 package org.apache.commons.math4.legacy.analysis.interpolation;
 
 import java.util.Arrays;
+import java.util.function.DoubleBinaryOperator;
+import java.util.function.Function;
 
 import org.apache.commons.numbers.core.Sum;
 import org.apache.commons.math4.legacy.analysis.BivariateFunction;
@@ -92,6 +94,38 @@ public class BicubicInterpolatingFunction
         throws DimensionMismatchException,
                NoDataException,
                NonMonotonicSequenceException {
+        this(x, y, f, dFdX, dFdY, d2FdXdY, false);
+    }
+
+    /**
+     * @param x Sample values of the x-coordinate, in increasing order.
+     * @param y Sample values of the y-coordinate, in increasing order.
+     * @param f Values of the function on every grid point.
+     * @param dFdX Values of the partial derivative of function with respect
+     * to x on every grid point.
+     * @param dFdY Values of the partial derivative of function with respect
+     * to y on every grid point.
+     * @param d2FdXdY Values of the cross partial derivative of function on
+     * every grid point.
+     * @param initializeDerivatives Whether to initialize the internal data
+     * needed for calling any of the methods that compute the partial derivatives
+     * this function.
+     * @throws DimensionMismatchException if the various arrays do not contain
+     * the expected number of elements.
+     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
+     * not strictly increasing.
+     * @throws NoDataException if any of the arrays has zero length.
+     */
+    public BicubicInterpolatingFunction(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;
 
@@ -147,7 +181,10 @@ public class BicubicInterpolatingFunction
                     d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
                 };
 
-                splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta));
+                splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta),
+                                                    xR,
+                                                    yR,
+                                                    initializeDerivatives);
             }
         }
     }
@@ -181,6 +218,75 @@ public class BicubicInterpolatingFunction
             y > yval[yval.length - 1]);
     }
 
+    /**
+     * @return the first partial derivative respect to x.
+     * @throws NullPointerException if the internal data were not initialized
+     * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
+     *             double[][],double[][],double[][],boolean) constructor}).
+     */
+    public DoubleBinaryOperator partialDerivativeX() {
+        return partialDerivative(BicubicFunction::partialDerivativeX);
+    }
+
+    /**
+     * @return the first partial derivative respect to y.
+     * @throws NullPointerException if the internal data were not initialized
+     * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
+     *             double[][],double[][],double[][],boolean) constructor}).
+     */
+    public DoubleBinaryOperator partialDerivativeY() {
+        return partialDerivative(BicubicFunction::partialDerivativeY);
+    }
+
+    /**
+     * @return the second partial derivative respect to x.
+     * @throws NullPointerException if the internal data were not initialized
+     * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
+     *             double[][],double[][],double[][],boolean) constructor}).
+     */
+    public DoubleBinaryOperator partialDerivativeXX() {
+        return partialDerivative(BicubicFunction::partialDerivativeXX);
+    }
+
+    /**
+     * @return the second partial derivative respect to y.
+     * @throws NullPointerException if the internal data were not initialized
+     * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
+     *             double[][],double[][],double[][],boolean) constructor}).
+     */
+    public DoubleBinaryOperator partialDerivativeYY() {
+        return partialDerivative(BicubicFunction::partialDerivativeYY);
+    }
+
+    /**
+     * @return the second partial cross derivative.
+     * @throws NullPointerException if the internal data were not initialized
+     * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
+     *             double[][],double[][],double[][],boolean) constructor}).
+     */
+    public DoubleBinaryOperator partialDerivativeXY() {
+        return partialDerivative(BicubicFunction::partialDerivativeXY);
+    }
+
+    /**
+     * @param which derivative function to apply.
+     * @return the selected partial derivative.
+     * @throws NullPointerException if the internal data were not initialized
+     * (cf. {@link #BicubicInterpolatingFunction(double[],double[],double[][],
+     *             double[][],double[][],double[][],boolean) constructor}).
+     */
+    private DoubleBinaryOperator partialDerivative(Function<BicubicFunction, BivariateFunction> which) {
+        return (x, y) -> {
+            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 which.apply(splines[i][j]).value(xN, yN);
+        };
+    }
+
     /**
      * @param c Coordinate.
      * @param val Coordinate samples.
@@ -256,6 +362,7 @@ public class BicubicInterpolatingFunction
 
         return a;
     }
+
 }
 
 /**
@@ -266,13 +373,31 @@ class BicubicFunction implements BivariateFunction {
     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.
+     * @param xR x spacing.
+     * @param yR y spacing.
+     * @param initializeDerivatives Whether to initialize the internal data
+     * needed for calling any of the methods that compute the partial derivatives
+     * this function.
      */
-    BicubicFunction(double[] coeff) {
+    BicubicFunction(double[] coeff,
+                    double xR,
+                    double yR,
+                    boolean initializeDerivatives) {
         a = new double[N][N];
         for (int j = 0; j < N; j++) {
             final double[] aJ = a[j];
@@ -280,6 +405,80 @@ class BicubicFunction implements BivariateFunction {
                 aJ[i] = 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 = (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) / xR;
+            };
+            partialDerivativeY = (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) / yR;
+            };
+            partialDerivativeXX = (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) / (xR * xR);
+            };
+            partialDerivativeYY = (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) / (yR * yR);
+            };
+            partialDerivativeXY = (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) / (xR * yR);
+            };
+        } else {
+            partialDerivativeX = null;
+            partialDerivativeY = null;
+            partialDerivativeXX = null;
+            partialDerivativeYY = null;
+            partialDerivativeXY = null;
+        }
     }
 
     /**
@@ -322,4 +521,40 @@ class BicubicFunction implements BivariateFunction {
 
         return result;
     }
+
+    /**
+     * @return the partial derivative wrt {@code x}.
+     */
+    BivariateFunction partialDerivativeX() {
+        return partialDerivativeX;
+    }
+
+    /**
+     * @return the partial derivative wrt {@code y}.
+     */
+    BivariateFunction partialDerivativeY() {
+        return partialDerivativeY;
+    }
+
+    /**
+     * @return the second partial derivative wrt {@code x}.
+     */
+    BivariateFunction partialDerivativeXX() {
+        return partialDerivativeXX;
+    }
+
+    /**
+     * @return the second partial derivative wrt {@code y}.
+     */
+    BivariateFunction partialDerivativeYY() {
+        return partialDerivativeYY;
+    }
+
+    /**
+     * @return the second partial cross-derivative.
+     */
+    BivariateFunction partialDerivativeXY() {
+        return partialDerivativeXY;
+    }
+
 }

commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolator.java

@@ -41,6 +41,32 @@ import org.apache.commons.math4.legacy.core.MathArrays;
  */
 public class BicubicInterpolator
     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 #BicubicInterpolator(boolean) initializeDerivatives}
+     * is set to {@code false}.
+     */
+    public BicubicInterpolator() {
+        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 BicubicInterpolatingFunction function} returned from
+     * the call to {@link #interpolate(double[],double[],double[][]) interpolate}.
+     */
+    public BicubicInterpolator(boolean initializeDerivatives) {
+        this.initializeDerivatives = initializeDerivatives;
+    }
     /**
      * {@inheritDoc}
      */
@@ -96,7 +122,8 @@ public class BicubicInterpolator
 
         // Create the interpolating function.
         return new BicubicInterpolatingFunction(xval, yval, fval,
-                                                dFdX, dFdY, d2FdXdY) {
+                                                dFdX, dFdY, d2FdXdY,
+                                                initializeDerivatives) {
             /** {@inheritDoc} */
             @Override
             public boolean isValidPoint(double x, double y) {

commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunctionTest.java

@@ -16,6 +16,7 @@
  */
 package org.apache.commons.math4.legacy.analysis.interpolation;
 
+import java.util.function.DoubleBinaryOperator;
 import org.apache.commons.math4.legacy.analysis.BivariateFunction;
 import org.apache.commons.statistics.distribution.ContinuousDistribution;
 import org.apache.commons.statistics.distribution.UniformContinuousDistribution;
@@ -286,6 +287,233 @@ public final class BicubicInterpolatingFunctionTest {
                           maxTolerance,
                           false);
     }
+    
+    /**
+     * Test for partial derivatives of {@link BicubicFunction}.
+     * <p>
+     * f(x, y) = &Sigma;<sub>i</sub>&Sigma;<sub>j</sub> (i+1) (j+2) x<sup>i</sup> y<sup>j</sup>
+     */
+    @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 BicubicFunction f = new BicubicFunction(coeff, 1, 1, true);
+        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 BicubicInterpolatingFunction} 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>
+     */
+    @Test
+    public void testMatchingPartialDerivatives() {
+        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] = i * delta;
+            yval[i] = i * delta / 3;
+        }
+        // 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(xval[i], yval[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(xval[i], yval[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(xval[i], yval[j]);
+            }
+        }
+        // Second partial derivatives with respect to x
+        double[][] d2Fd2X = new double[sz][sz];
+        BivariateFunction d2fd2X = new BivariateFunction() {
+                public double value(double x, double y) {
+                    return 4 + 8 * y - 18 * x;
+                }
+            };
+        for (int i = 0; i < sz; i++) {
+            for (int j = 0; j < sz; j++) {
+                d2Fd2X[i][j] = d2fd2X.value(xval[i], yval[j]);
+            }
+        }
+        // Second partial derivatives with respect to y
+        double[][] d2Fd2Y = new double[sz][sz];
+        BivariateFunction d2fd2Y = new BivariateFunction() {
+                public double value(double x, double y) {
+                    return - 6 - 2 * x + 6 * y;
+                }
+            };
+        for (int i = 0; i < sz; i++) {
+            for (int j = 0; j < sz; j++) {
+                d2Fd2Y[i][j] = d2fd2Y.value(xval[i], yval[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(xval[i], yval[j]);
+            }
+        }
+
+        BicubicInterpolatingFunction bcf
+            = new BicubicInterpolatingFunction(xval, yval, fval, dFdX, dFdY, d2FdXdY, true);
+        DoubleBinaryOperator partialDerivativeX = bcf.partialDerivativeX();
+        DoubleBinaryOperator partialDerivativeY = bcf.partialDerivativeY();
+        DoubleBinaryOperator partialDerivativeXX = bcf.partialDerivativeXX();
+        DoubleBinaryOperator partialDerivativeYY = bcf.partialDerivativeYY();
+        DoubleBinaryOperator partialDerivativeXY = bcf.partialDerivativeXY();
+
+        double x;
+        double y;
+        double expected;
+        double result;
+
+        final double tol = 1e-10;
+        for (int i = 0; i < sz; i++) {
+            x = xval[i];
+            for (int j = 0; j < sz; j++) {
+                y = yval[j];
+
+                expected = dfdX.value(x, y);
+                result = partialDerivativeX.applyAsDouble(x, y);
+                Assert.assertEquals(x + " " + y + " dFdX", expected, result, tol);
+
+                expected = dfdY.value(x, y);
+                result = partialDerivativeY.applyAsDouble(x, y);
+                Assert.assertEquals(x + " " + y + " dFdY", expected, result, tol);
+
+                expected = d2fd2X.value(x, y);
+                result = partialDerivativeXX.applyAsDouble(x, y);
+                Assert.assertEquals(x + " " + y + " d2Fd2X", expected, result, tol);
+
+                expected = d2fd2Y.value(x, y);
+                result = partialDerivativeYY.applyAsDouble(x, y);
+                Assert.assertEquals(x + " " + y + " d2Fd2Y", expected, result, tol);
+
+                expected = d2fdXdY.value(x, y);
+                result = partialDerivativeXY.applyAsDouble(x, y);
+                Assert.assertEquals(x + " " + y + " d2FdXdY", expected, result, tol);
+            }
+        }
+    }
 
     /**
      * @param minimumX Lower bound of interpolation range along the x-coordinate.