This is an automated email from the ASF dual-hosted git repository.
erans pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-math.git
commit 8e6dee595832951e4d74e795a76e5aff236d4181
Author: amoscatelli <alessandro.moscate...@live.com>
AuthorDate: Fri Jul 22 15:21:49 2022 +0200
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<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;
+ }
+
}
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}.
+ * <p>
+ * f(x, y) = Σ<sub>i</sub>Σ<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.
