From 8e6dee595832951e4d74e795a76e5aff236d4181 Mon Sep 17 00:00:00 2001 From: amoscatelli Date: Fri, 22 Jul 2022 15:21:49 +0200 Subject: [PATCH] MATH-1648: Derivatives for "BicubicInterpolatingFunction". --- .../BicubicInterpolatingFunction.java | 239 +++++++++++++++++- .../interpolation/BicubicInterpolator.java | 29 ++- .../BicubicInterpolatingFunctionTest.java | 228 +++++++++++++++++ 3 files changed, 493 insertions(+), 3 deletions(-) diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunction.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunction.java index fda6349c6..2543f0b8f 100644 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunction.java +++ b/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 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; + } + } diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolator.java index 839d11181..0b3cedb23 100644 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolator.java +++ b/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) { diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunctionTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunctionTest.java index 173ff54c5..6a17b44d9 100644 --- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/BicubicInterpolatingFunctionTest.java +++ b/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}. + *

+ * f(x, y) = ΣiΣj (i+1) (j+2) xi yj + */ + @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. + *

+ * f(x, y) = 5 + * - 3 x + 2 y + * - x y + 2 x2 - 3 y2 + * + 4 x2 y - x y2 - 3 x3 + y3 + */ + @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.