Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-366 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@937080 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski 2010-04-22 22:09:21 +00:00
parent d9f545ea68
commit dafd0976ad
11 changed files with 2089 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java Normal file
View File

@ -0,0 +1,40 @@
/*
* 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.math.analysis;
import org.apache.commons.math.FunctionEvaluationException;
/**
* An interface representing a trivariate real function.
*
* @since 2.2
* @version $Revision$ $Date$
*/
public interface TrivariateRealFunction {
/**
* Compute the value for the function.
*
* @param x x-coordinate for which the function value should be computed.
* @param y y-coordinate for which the function value should be computed.
* @param z z-coordinate for which the function value should be computed.
* @return the value.
* @throws FunctionEvaluationException if the function evaluation fails.
*/
public double value(double x, double y, double z)
throws FunctionEvaluationException;
}

422
src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java
View File

@ -18,6 +18,7 @@ package org.apache.commons.math.analysis.interpolation;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.DimensionMismatchException;
import org.apache.commons.math.analysis.BivariateRealFunction;
@ -58,19 +59,29 @@ public class BicubicSplineInterpolatingFunction
private final double[] xval;
/** Samples y-coordinates */
private final double[] yval;
/** Set of cubic splines pacthing the whole data grid */
/** Set of cubic splines patching the whole data grid */
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 BivariateRealFunction[][][] partialDerivatives = null;
/**
* @param x Sample values of the x-coordinate, in increasing order
* @param y Sample values of the y-coordinate, in increasing order
* @param z Values of the function on every grid point
* @param dZdX Values of the partial derivative of function with respect
* to x on every grid point
* @param dZdY Values of the partial derivative of function with respect
* to y on every grid point
* @param dZdXdY Values of the cross partial derivative of function on
* every grid point
* @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.
* @throws DimensionMismatchException if the various arrays do not contain
* the expected number of elements.
* @throws IllegalArgumentException if {@code x} or {@code y} are not strictly
@ -78,28 +89,28 @@ public class BicubicSplineInterpolatingFunction
*/
public BicubicSplineInterpolatingFunction(double[] x,
double[] y,
double[][] z,
double[][] dZdX,
double[][] dZdY,
double[][] dZdXdY)
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY)
throws DimensionMismatchException {
final int xLen = x.length;
final int yLen = y.length;
if (xLen == 0 || yLen == 0 || z.length == 0 || z[0].length == 0) {
if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
throw MathRuntimeException.createIllegalArgumentException("no data");
}
if (xLen != z.length) {
throw new DimensionMismatchException(xLen, z.length);
if (xLen != f.length) {
throw new DimensionMismatchException(xLen, f.length);
}
if (xLen != dZdX.length) {
throw new DimensionMismatchException(xLen, dZdX.length);
if (xLen != dFdX.length) {
throw new DimensionMismatchException(xLen, dFdX.length);
}
if (xLen != dZdY.length) {
throw new DimensionMismatchException(xLen, dZdY.length);
if (xLen != dFdY.length) {
throw new DimensionMismatchException(xLen, dFdY.length);
}
if (xLen != dZdXdY.length) {
throw new DimensionMismatchException(xLen, dZdXdY.length);
if (xLen != d2FdXdY.length) {
throw new DimensionMismatchException(xLen, d2FdXdY.length);
}
MathUtils.checkOrder(x, 1, true);
@ -113,26 +124,26 @@ public class BicubicSplineInterpolatingFunction
splines = new BicubicSplineFunction[lastI][lastJ];
for (int i = 0; i < lastI; i++) {
if (z[i].length != yLen) {
throw new DimensionMismatchException(z[i].length, yLen);
if (f[i].length != yLen) {
throw new DimensionMismatchException(f[i].length, yLen);
}
if (dZdX[i].length != yLen) {
throw new DimensionMismatchException(dZdX[i].length, yLen);
if (dFdX[i].length != yLen) {
throw new DimensionMismatchException(dFdX[i].length, yLen);
}
if (dZdY[i].length != yLen) {
throw new DimensionMismatchException(dZdY[i].length, yLen);
if (dFdY[i].length != yLen) {
throw new DimensionMismatchException(dFdY[i].length, yLen);
}
if (dZdXdY[i].length != yLen) {
throw new DimensionMismatchException(dZdXdY[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[] {
z[i][j], z[ip1][j], z[i][jp1], z[ip1][jp1],
dZdX[i][j], dZdX[ip1][j], dZdX[i][jp1], dZdX[ip1][jp1],
dZdY[i][j], dZdY[ip1][j], dZdY[i][jp1], dZdY[ip1][jp1],
dZdXdY[i][j], dZdXdY[ip1][j], dZdXdY[i][jp1], dZdXdY[ip1][jp1]
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));
@ -162,18 +173,119 @@ public class BicubicSplineInterpolatingFunction
}
/**
* @param c coordinate
* @param val coordinate samples
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the first partial derivative with
* respect to x.
*/
public double partialDerivativeX(double x, double y) {
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.
*/
public double partialDerivativeY(double x, double y) {
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.
*/
public double partialDerivativeXX(double x, double y) {
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.
*/
public double partialDerivativeYY(double x, double y) {
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.
*/
public double partialDerivativeXY(double x, double y) {
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.
*/
private double partialDerivative(int which, double x, double y) {
if (partialDerivatives == null) {
computePartialDerivatives();
}
final int i = searchIndex(x, xval);
if (i == -1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
x, xval[0], xval[xval.length - 1]);
}
final int j = searchIndex(y, yval);
if (j == -1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
y, yval[0], yval[yval.length - 1]);
}
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
double result = Double.NaN;
try { // Workaround to avoid carrying around "try" blocks for an
// exception that will never happen
result = partialDerivatives[which][i][j].value(xN, yN);
} catch (FunctionEvaluationException e) {
// Will never happen
}
return result;
}
/**
* Compute all partial derivatives.
*/
private void computePartialDerivatives() {
final int lastI = xval.length - 1;
final int lastJ = yval.length - 1;
partialDerivatives = new BivariateRealFunction[5][lastI][lastJ];
for (int i = 0; i < lastI; i++) {
for (int j = 0; j < lastJ; j++) {
final BicubicSplineFunction f = splines[i][j];
partialDerivatives[0][i][j] = f.partialDerivativeX();
partialDerivatives[1][i][j] = f.partialDerivativeY();
partialDerivatives[2][i][j] = f.partialDerivativeXX();
partialDerivatives[3][i][j] = f.partialDerivativeYY();
partialDerivatives[4][i][j] = f.partialDerivativeXY();
}
}
}
/**
* @param c Coordinate.
* @param val Coordinate samples.
* @return the index in {@code val} corresponding to the interval
* containing {@code c}, or {@code -1} if {@code c} is out of the
* range defined by the end values of {@code val}
* range defined by the end values of {@code val}.
*/
private int searchIndex(double c, double[] val) {
if (c < val[0]) {
return -1;
}
int max = val.length;
final int max = val.length;
for (int i = 1; i < max; i++) {
if (c <= val[i]) {
return i - 1;
@ -192,22 +304,25 @@ public class BicubicSplineInterpolatingFunction
* <li>f(1,0)</li>
* <li>f(0,1)</li>
* <li>f(1,1)</li>
* <li>fx(0,0)</li>
* <li>fx(1,0)</li>
* <li>fx(0,1)</li>
* <li>fx(1,1)</li>
* <li>fy(0,0)</li>
* <li>fy(1,0)</li>
* <li>fy(0,1)</li>
* <li>fy(1,1)</li>
* <li>fxy(0,0)</li>
* <li>fxy(1,0)</li>
* <li>fxy(0,1)</li>
* <li>fxy(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
* values.
* @return the spline coefficients.
*/
private double[] computeSplineCoefficients(double[] beta) {
final double[] a = new double[16];
@ -224,6 +339,7 @@ public class BicubicSplineInterpolatingFunction
return a;
}
}
/**
* 2D-spline function.
*
@ -231,41 +347,29 @@ public class BicubicSplineInterpolatingFunction
*/
class BicubicSplineFunction
implements BivariateRealFunction {
//CHECKSTYLE: stop MultipleVariableDeclarations
private static final short N = 4;
/** Coefficients */
private final double
a00, a01, a02, a03,
a10, a11, a12, a13,
a20, a21, a22, a23,
a30, a31, a32, a33;
//CHECKSTYLE: resume MultipleVariableDeclarations
private final double[][] a = new double[N][N];
/** Partial derivatives */
BivariateRealFunction partialDerivativeX = null;
BivariateRealFunction partialDerivativeY = null;
BivariateRealFunction partialDerivativeXX = null;
BivariateRealFunction partialDerivativeYY = null;
BivariateRealFunction partialDerivativeXY = null;
/**
* @param a Spline coefficients
*/
public BicubicSplineFunction(double[] a) {
this.a00 = a[0];
this.a10 = a[1];
this.a20 = a[2];
this.a30 = a[3];
this.a01 = a[4];
this.a11 = a[5];
this.a21 = a[6];
this.a31 = a[7];
this.a02 = a[8];
this.a12 = a[9];
this.a22 = a[10];
this.a32 = a[11];
this.a03 = a[12];
this.a13 = a[13];
this.a23 = a[14];
this.a33 = a[15];
public BicubicSplineFunction(double[] aV) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
a[i][j] = aV[i + N * j];
}
}
}
/**
* @param x x-coordinate of the interpolation point
* @param y y-coordinate of the interpolation point
* @return the interpolated value.
* {@inheritDoc}
*/
public double value(double x, double y) {
if (x < 0 || x > 1) {
@ -279,12 +383,162 @@ class BicubicSplineFunction
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 a00 + a01 * y + a02 * y2 + a03 * y3 +
a10 * x + a11 * x * y + a12 * x * y2 + a13 * x * y3 +
a20 * x2 + a21 * x2 * y + a22 * x2 * y2 + a23 * x2 * y3 +
a30 * x3 + a31 * x3 * y + a32 * x3 * y2 + a33 * x3 * 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 BivariateRealFunction partialDerivativeX() {
if (partialDerivativeX == null) {
computePartialDerivatives();
}
return partialDerivativeX;
}
/**
* @return the partial derivative wrt {@code y}.
*/
public BivariateRealFunction partialDerivativeY() {
if (partialDerivativeY == null) {
computePartialDerivatives();
}
return partialDerivativeY;
}
/**
* @return the second partial derivative wrt {@code x}.
*/
public BivariateRealFunction partialDerivativeXX() {
if (partialDerivativeXX == null) {
computePartialDerivatives();
}
return partialDerivativeXX;
}
/**
* @return the second partial derivative wrt {@code y}.
*/
public BivariateRealFunction partialDerivativeYY() {
if (partialDerivativeYY == null) {
computePartialDerivatives();
}
return partialDerivativeYY;
}
/**
* @return the second partial cross-derivative.
*/
public BivariateRealFunction partialDerivativeXY() {
if (partialDerivativeXY == null) {
computePartialDerivatives();
}
return partialDerivativeXY;
}
/**
* Compute all partial derivatives functions.
*/
private void computePartialDerivatives() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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);
}
};
}
}

2
src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java
View File

@ -40,7 +40,7 @@ public class SmoothingPolynomialBicubicSplineInterpolator
private final PolynomialFitter yFitter;
/**
* Default constructor. The degree of the fitting polynomials are set to 3.
* Default constructor. The degree of the fitting polynomials is set to 3.
*/
public SmoothingPolynomialBicubicSplineInterpolator() {
this(3);

487
src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java Normal file
View File

@ -0,0 +1,487 @@
/*
* 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.math.analysis.interpolation;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.DimensionMismatchException;
import org.apache.commons.math.analysis.TrivariateRealFunction;
/**
* Function that implements the
* <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
* tricubic spline interpolation</a>, as proposed in
* <quote>
* Tricubic interpolation in three dimensions<br/>
* F. Lekien and J. Marsden<br/>
* <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
* </quote>
*
* @version $Revision$ $Date$
* @since 2.2
*/
public class TricubicSplineInterpolatingFunction
implements TrivariateRealFunction {
/**
* 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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 },
{-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,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 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 },
{ 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 },
{ -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 },
{ 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 },
{ -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 },
{ 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 },
{ -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 },
{ 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 },
{ -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 },
{ 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 },
{ -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 },
{ 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 },
{ -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 },
{ 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 }
};
/** Samples x-coordinates */
private final double[] xval;
/** Samples y-coordinates */
private final double[] yval;
/** Samples z-coordinates */
private final double[] zval;
/** Set of cubic splines pacthing the whole data grid */
private final TricubicSplineFunction[][][] splines;
/**
* @param x Sample values of the x-coordinate, in increasing order.
* @param y Sample values of the y-coordinate, in increasing order.
* @param z 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 d2FdXdZ Values of the cross partial derivative of function on
* every grid point.
* @param d2FdYdZ Values of the cross partial derivative of function on
* every grid point.
* @param d3FdXdYdZ 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 IllegalArgumentException if {@code x}, {@code y} or {@code z}
* are not strictly increasing.
*/
public TricubicSplineInterpolatingFunction(double[] x,
double[] y,
double[] z,
double[][][] f,
double[][][] dFdX,
double[][][] dFdY,
double[][][] dFdZ,
double[][][] d2FdXdY,
double[][][] d2FdXdZ,
double[][][] d2FdYdZ,
double[][][] d3FdXdYdZ)
throws DimensionMismatchException {
final int xLen = x.length;
final int yLen = y.length;
final int zLen = z.length;
if (xLen == 0 || yLen == 0 || z.length == 0
|| f.length == 0 || f[0].length == 0) {
throw MathRuntimeException.createIllegalArgumentException("no data");
}
if (xLen != f.length) {
throw new DimensionMismatchException(xLen, f.length);
}
if (xLen != dFdX.length) {
throw new DimensionMismatchException(xLen, dFdX.length);
}
if (xLen != dFdY.length) {
throw new DimensionMismatchException(xLen, dFdY.length);
}
if (xLen != dFdZ.length) {
throw new DimensionMismatchException(xLen, dFdZ.length);
}
if (xLen != d2FdXdY.length) {
throw new DimensionMismatchException(xLen, d2FdXdY.length);
}
if (xLen != d2FdXdZ.length) {
throw new DimensionMismatchException(xLen, d2FdXdZ.length);
}
if (xLen != d2FdYdZ.length) {
throw new DimensionMismatchException(xLen, d2FdYdZ.length);
}
if (xLen != d3FdXdYdZ.length) {
throw new DimensionMismatchException(xLen, d3FdXdYdZ.length);
}
MathUtils.checkOrder(x, 1, true);
MathUtils.checkOrder(y, 1, true);
MathUtils.checkOrder(z, 1, true);
xval = x.clone();
yval = y.clone();
zval = z.clone();
final int lastI = xLen - 1;
final int lastJ = yLen - 1;
final int lastK = zLen - 1;
splines = new TricubicSplineFunction[lastI][lastJ][lastK];
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 (dFdZ[i].length != yLen) {
throw new DimensionMismatchException(dFdZ[i].length, yLen);
}
if (d2FdXdY[i].length != yLen) {
throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
}
if (d2FdXdZ[i].length != yLen) {
throw new DimensionMismatchException(d2FdXdZ[i].length, yLen);
}
if (d2FdYdZ[i].length != yLen) {
throw new DimensionMismatchException(d2FdYdZ[i].length, yLen);
}
if (d3FdXdYdZ[i].length != yLen) {
throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen);
}
final int ip1 = i + 1;
for (int j = 0; j < lastJ; j++) {
if (f[i][j].length != zLen) {
throw new DimensionMismatchException(f[i][j].length, zLen);
}
if (dFdX[i][j].length != zLen) {
throw new DimensionMismatchException(dFdX[i][j].length, zLen);
}
if (dFdY[i][j].length != zLen) {
throw new DimensionMismatchException(dFdY[i][j].length, zLen);
}
if (dFdZ[i][j].length != zLen) {
throw new DimensionMismatchException(dFdZ[i][j].length, zLen);
}
if (d2FdXdY[i][j].length != zLen) {
throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen);
}
if (d2FdXdZ[i][j].length != zLen) {
throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen);
}
if (d2FdYdZ[i][j].length != zLen) {
throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen);
}
if (d3FdXdYdZ[i][j].length != zLen) {
throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen);
}
final int jp1 = j + 1;
for (int k = 0; k < lastK; k++) {
final int kp1 = k + 1;
final double[] beta = new double[] {
f[i][j][k], f[ip1][j][k],
f[i][jp1][k], f[ip1][jp1][k],
f[i][j][kp1], f[ip1][j][kp1],
f[i][jp1][kp1], f[ip1][jp1][kp1],
dFdX[i][j][k], dFdX[ip1][j][k],
dFdX[i][jp1][k], dFdX[ip1][jp1][k],
dFdX[i][j][kp1], dFdX[ip1][j][kp1],
dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1],
dFdY[i][j][k], dFdY[ip1][j][k],
dFdY[i][jp1][k], dFdY[ip1][jp1][k],
dFdY[i][j][kp1], dFdY[ip1][j][kp1],
dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1],
dFdZ[i][j][k], dFdZ[ip1][j][k],
dFdZ[i][jp1][k], dFdZ[ip1][jp1][k],
dFdZ[i][j][kp1], dFdZ[ip1][j][kp1],
dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1],
d2FdXdY[i][j][k], d2FdXdY[ip1][j][k],
d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k],
d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1],
d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1],
d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k],
d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k],
d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1],
d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1],
d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k],
d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k],
d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1],
d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1],
d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k],
d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k],
d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1],
d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1],
};
splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta));
}
}
}
}
/**
* {@inheritDoc}
*/
public double value(double x, double y, double z) {
final int i = searchIndex(x, xval);
if (i == -1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
x, xval[0], xval[xval.length - 1]);
}
final int j = searchIndex(y, yval);
if (j == -1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
y, yval[0], yval[yval.length - 1]);
}
final int k = searchIndex(z, zval);
if (k == -1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
z, zval[0], zval[zval.length - 1]);
}
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);
return splines[i][j][k].value(xN, yN, zN);
}
/**
* @param c Coordinate.
* @param val Coordinate samples.
* @return the index in {@code val} corresponding to the interval
* containing {@code c}, or {@code -1} if {@code c} is out of the
* range defined by the end values of {@code val}.
*/
private int searchIndex(double c, double[] val) {
if (c < val[0]) {
return -1;
}
final int max = val.length;
for (int i = 1; i < max; i++) {
if (c <= val[i]) {
return i - 1;
}
}
return -1;
}
/**
* 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,0)</li>
* <li>f(1,0,0)</li>
* <li>f(0,1,0)</li>
* <li>f(1,1,0)</li>
* <li>f(0,0,1)</li>
* <li>f(1,0,1)</li>
* <li>f(0,1,1)</li>
* <li>f(1,1,1)</li>
*
* <li>f<sub>x</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>x</sub>(1,1,1)</li>
*
* <li>f<sub>y</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>y</sub>(1,1,1)</li>
*
* <li>f<sub>z</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>z</sub>(1,1,1)</li>
*
* <li>f<sub>xy</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>xy</sub>(1,1,1)</li>
*
* <li>f<sub>xz</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>xz</sub>(1,1,1)</li>
*
* <li>f<sub>yz</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>yz</sub>(1,1,1)</li>
*
* <li>f<sub>xyz</sub>(0,0,0)</li>
* <li>... <em>(same order as above)</em></li>
* <li>f<sub>xyz</sub>(1,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 int sz = 64;
final double[] a = new double[sz];
for (int i = 0; i < sz; i++) {
double result = 0;
final double[] row = AINV[i];
for (int j = 0; j < sz; j++) {
result += row[j] * beta[j];
}
a[i] = result;
}
return a;
}
}
/**
* 3D-spline function.
*
* @version $Revision$ $Date$
*/
class TricubicSplineFunction
implements TrivariateRealFunction {
private static final short N = 4;
private static final short N2 = N * N;
/** Coefficients */
private final double[][][] a = new double[N][N][N];
/**
* @param aV List of spline coefficients.
*/
public TricubicSplineFunction(double[] aV) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
for (int k = 0; k < N; k++) {
a[i][j][k] = aV[i + N * j + N2 * k];
}
}
}
}
/**
* @param x x-coordinate of the interpolation point.
* @param y y-coordinate of the interpolation point.
* @param z z-coordinate of the interpolation point.
* @return the interpolated value.
*/
public double value(double x, double y, double z) {
if (x < 0 || x > 1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
x, 0, 1);
}
if (y < 0 || y > 1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
y, 0, 1);
}
if (z < 0 || z > 1) {
throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
z, 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 };
final double z2 = z * z;
final double z3 = z2 * z;
final double[] pZ = { 1, z, z2, z3 };
double result = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
for (int k = 0; k < N; k++) {
result += a[i][j][k] * pX[i] * pY[j] * pZ[k];
}
}
}
return result;
}
}

200
src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java Normal file
View File

@ -0,0 +1,200 @@
/*
* 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.math.analysis.interpolation;
import org.apache.commons.math.DimensionMismatchException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
/**
* Generates a tricubic interpolating function.
*
* @version $Revision$ $Date$
* @since 2.2
*/
public class TricubicSplineInterpolator
implements TrivariateRealGridInterpolator {
/**
* {@inheritDoc}
*/
public TricubicSplineInterpolatingFunction interpolate(final double[] xval,
final double[] yval,
final double[] zval,
final double[][][] fval)
throws MathException, IllegalArgumentException {
if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
throw MathRuntimeException.createIllegalArgumentException("no data");
}
if (xval.length != fval.length) {
throw new DimensionMismatchException(xval.length, fval.length);
}
MathUtils.checkOrder(xval, 1, true);
MathUtils.checkOrder(yval, 1, true);
MathUtils.checkOrder(zval, 1, true);
final int xLen = xval.length;
final int yLen = yval.length;
final int zLen = zval.length;
// Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
// fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
// fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
final double[][][] fvalXY = new double[zLen][xLen][yLen];
final double[][][] fvalZX = new double[yLen][zLen][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++) {
if (fval[i][j].length != zLen) {
throw new DimensionMismatchException(fval[i][j].length, zLen);
}
for (int k = 0; k < zLen; k++) {
final double v = fval[i][j][k];
fvalXY[k][i][j] = v;
fvalZX[j][k][i] = v;
}
}
}
final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();
// For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
final BicubicSplineInterpolatingFunction[] xSplineYZ
= new BicubicSplineInterpolatingFunction[xLen];
for (int i = 0; i < xLen; i++) {
xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]);
}
// For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
final BicubicSplineInterpolatingFunction[] ySplineZX
= new BicubicSplineInterpolatingFunction[yLen];
for (int j = 0; j < yLen; j++) {
ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]);
}
// For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
final BicubicSplineInterpolatingFunction[] zSplineXY
= new BicubicSplineInterpolatingFunction[zLen];
for (int k = 0; k < zLen; k++) {
zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]);
}
// Partial derivatives wrt x and wrt y
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];
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);
for (int k = 0; k < zLen; k++) {
final int nK = nextIndex(k, zLen);
final int pK = previousIndex(k);
// XXX Not sure about this formula
d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
fval[pI][nJ][nK] + fval[pI][pJ][nK] -
fval[nI][nJ][pK] + fval[nI][pJ][pK] +
fval[pI][nJ][pK] - fval[pI][pJ][pK]) /
((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ;
}
}
}
// Create the interpolating splines
return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval,
dFdX, dFdY, dFdZ,
d2FdXdY, d2FdZdX, d2FdYdZ,
d3FdXdYdZ);
}
/**
* Compute the next index of an array, clipping if necessary.
* It is assumed (but not checked) that {@code i} is larger than or equal to 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;
}
/**
* Compute 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;
}
}

46
src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java Normal file
View File

@ -0,0 +1,46 @@
/*
* 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.math.analysis.interpolation;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.TrivariateRealFunction;
/**
* Interface representing a trivariate real interpolating function where the
* sample points must be specified on a regular grid.
*
* @version $Revision$ $Date$
*/
public interface TrivariateRealGridInterpolator {
/**
* Computes an interpolating function for the data set.
*
* @param xval All the x-coordinates of the interpolation points, sorted
* in increasing order.
* @param yval All the y-coordinates of the interpolation points, sorted
* in increasing order.
* @param zval All the z-coordinates of the interpolation points, sorted
* in increasing order.
* @param fval the values of the interpolation points on all the grid knots:
* {@code fval[i][j][k] = f(xval[i], yval[j], zval[k])}.
* @return a function which interpolates the data set.
* @throws MathException if arguments violate assumptions made by the
* interpolation algorithm.
*/
TrivariateRealFunction interpolate(double[] xval, double[] yval, double[] zval, double[][][] fval)
throws MathException;
}

3
src/site/xdoc/changes.xml
View File

@ -69,6 +69,9 @@ This is primarily a maintenance release, but it also includes new features and e
DummyStepInterpolator requires an additional argument for one of its constructors;
some protected fields have been removed from AbstractLeastSquaresOptimizer, AbstractScalarDifferentiableOptimizer and AbstractLinearOptimizer;
and the isOptimal(SimplexTableau) method has been removed from SimplexSolver. ">
<action dev="erans" type="add" issue="MATH-366">
Implementation of tricubic interpolation.
</action>
<action dev="erans" type="update" issue="MATH-365">
Deprecated SmoothingBicubicSplineInterpolator and SmoothingBicubicSplineInterpolatorTest.
Added BicubicSplineInterpolator and BicubicSplineInterpolatorTest.

35
src/site/xdoc/userguide/analysis.xml
View File

@ -312,13 +312,13 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
org.apache.commons.math.analysis.interpolation.BivariateRealGridInterpolator</a>
is used to find a bivariate real-valued function <code>f</code> which
for a given set of tuples
(<code>x<sub>i</sub></code>,<code>y<sub>j</sub></code>,<code>z<sub>ij</sub></code>)
yields <code>f(x<sub>i</sub>,y<sub>j</sub>)=z<sub>ij</sub></code> to the best accuracy
(<code>x<sub>i</sub></code>,<code>y<sub>j</sub></code>,<code>f<sub>ij</sub></code>)
yields <code>f(x<sub>i</sub>,y<sub>j</sub>)=f<sub>ij</sub></code> to the best accuracy
possible. The result is provided as an object implementing the
<a href="../apidocs/org/apache/commons/math/analysis/BivariateRealFunction.html">
org.apache.commons.math.analysis.BivariateRealFunction</a> interface. It can therefore
be evaluated at any point, including a point not belonging to the original set.
The array <code>x<sub>i</sub></code> and <code>y<sub>j</sub></code> must be sorted in
The arrays <code>x<sub>i</sub></code> and <code>y<sub>j</sub></code> must be sorted in
increasing order in order to define a two-dimensional regular grid.
</p>
<p>
@ -339,6 +339,35 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
smoothing of the data by computing the polynomial that best fits each of the one-dimensional curves along each
of the coordinate axes.
</p>
<p>
A <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.html">
org.apache.commons.math.analysis.interpolation.TrivariateRealGridInterpolator</a>
is used to find a bivariate real-valued function <code>f</code> which
for a given set of tuples
(<code>x<sub>i</sub></code>,<code>y<sub>j</sub></code>,<code>z<sub>k</sub></code>,
<code>f<sub>ijk</sub></code>)
yields <code>f(x<sub>i</sub>,y<sub>j</sub>,z<sub>k</sub>)=f<sub>ijk</sub></code> to the
best accuracy possible. The result is provided as an object implementing the
<a href="../apidocs/org/apache/commons/math/analysis/TrivariateRealFunction.html">
org.apache.commons.math.analysis.TrivariateRealFunction</a> interface. It can therefore
be evaluated at any point, including a point not belonging to the original set.
The arrays <code>x<sub>i</sub></code>, <code>y<sub>j</sub></code> and
<code>z<sub>k</sub></code> must be sorted in increasing order in order to define a
three-dimensional regular grid.
</p>
<p>
In <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.html">
tricubic interpolation</a>, the interpolation function is a 3rd-degree polynomial of two
variables. The coefficients are computed from the function values sampled on a regular grid,
as well as the values of the partial derivatives of the function at those grid points.
</p>
<p>
From three-dimensional data sampled on a regular grid, the
<a href="../apidocs/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.html">
org.apache.commons.math.analysis.interpolation.TricubicSplineInterpolator</a>
computes a <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.html">
tricubic interpolating function</a>.
</p>
</subsection>
<subsection name="4.4 Integration" href="integration">
<p>

187
src/test/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
View File

@ -257,4 +257,191 @@ public final class BicubicSplineInterpolatingFunctionTest {
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>
*/
@Test
public void testSplinePartialDerivatives() throws MathException {
final int N = 4;
final double[] coeff = new double[16];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
final int index = i + N * j;
coeff[i + N * j] = (i + 1) * (j + 2);
}
}
final BicubicSplineFunction f = new BicubicSplineFunction(coeff);
BivariateRealFunction derivative;
final double x = 0.435;
final double y = 0.776;
final double tol = 1e-13;
derivative = new BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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 BivariateRealFunction() {
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>
*/
@Test
public void testMatchingPartialDerivatives() throws MathException {
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
BivariateRealFunction f = new BivariateRealFunction() {
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];
BivariateRealFunction dfdX = new BivariateRealFunction() {
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];
BivariateRealFunction dfdY = new BivariateRealFunction() {
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];
BivariateRealFunction d2fdXdY = new BivariateRealFunction() {
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);
}
}
}
}

544
src/test/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunctionTest.java Normal file
View File

@ -0,0 +1,544 @@
/*
* 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.math.analysis.interpolation;
import org.apache.commons.math.MathException;
import org.apache.commons.math.DimensionMismatchException;
import org.apache.commons.math.analysis.TrivariateRealFunction;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for the bicubic function.
*
* @version $Revision: 821626 $ $Date: 2009-10-04 23:57:30 +0200 (Sun, 04 Oct 2009) $
*/
public final class TricubicSplineInterpolatingFunctionTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() throws MathException {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
double[][][] fval = new double[xval.length][yval.length][zval.length];
@SuppressWarnings("unused")
TrivariateRealFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, fval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
tcf = new TricubicSplineInterpolatingFunction(wxval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (IllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -1, -1, 2.5};
try {
tcf = new TricubicSplineInterpolatingFunction(xval, wyval, zval,
fval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (IllegalArgumentException e) {
// Expected
}
double[] wzval = new double[] {-12, -8, -9, -3, 0, 2.5};
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, wzval,
fval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (IllegalArgumentException e) {
// Expected
}
double[][][] wfval = new double[xval.length - 1][yval.length - 1][zval.length];
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
wfval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, wfval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, wfval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, wfval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
wfval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, wfval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, wfval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length][yval.length - 1][zval.length];
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
wfval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, wfval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, wfval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, wfval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
wfval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, wfval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, wfval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length][yval.length][zval.length - 1];
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
wfval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, wfval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, wfval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, wfval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
wfval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, wfval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, wfval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Test for a plane.
* <p>
* f(x, y, z) = 2 x - 3 y - 4 z + 5
* </p>
*/
@Test
public void testPlane() throws MathException {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
// Function values
TrivariateRealFunction f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return 2 * x - 3 * y - 4 * z + 5;
}
};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to x
double[][][] dFdX = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdX[i][j][k] = 2;
}
}
}
// Partial derivatives with respect to y
double[][][] dFdY = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdY[i][j][k] = -3;
}
}
}
// Partial derivatives with respect to z
double[][][] dFdZ = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdZ[i][j][k] = -4;
}
}
}
// Partial cross-derivatives
double[][][] d2FdXdY = new double[xval.length][yval.length][zval.length];
double[][][] d2FdXdZ = new double[xval.length][yval.length][zval.length];
double[][][] d2FdYdZ = new double[xval.length][yval.length][zval.length];
double[][][] d3FdXdYdZ = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdXdY[i][j][k] = 0;
d2FdXdZ[i][j][k] = 0;
d2FdYdZ[i][j][k] = 0;
d3FdXdYdZ[i][j][k] = 0;
}
}
}
TrivariateRealFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, dFdX, dFdY, dFdZ,
d2FdXdY, d2FdXdZ, d2FdYdZ,
d3FdXdYdZ);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("On sample point",
expected, result, 1e-15);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.3);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.3);
}
/**
* Sine wave.
* <p>
* f(x, y, z) = a cos [&omega; z - k<sub>y</sub> x - k<sub>y</sub> y]
* </p>
* with A = 0.2, &omega; = 0.5, k<sub>x</sub> = 2, k<sub>y</sub> = 1.
*/
@Test
public void testWave() throws MathException {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 4};
final double a = 0.2;
final double omega = 0.5;
final double kx = 2;
final double ky = 1;
// Function values
TrivariateRealFunction f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.cos(omega * z - kx * x - ky * y);
}
};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to x
double[][][] dFdX = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction dFdX_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.sin(omega * z - kx * x - ky * y) * kx;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdX[i][j][k] = dFdX_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to y
double[][][] dFdY = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction dFdY_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.sin(omega * z - kx * x - ky * y) * ky;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdY[i][j][k] = dFdY_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to z
double[][][] dFdZ = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction dFdZ_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return -a * Math.sin(omega * z - kx * x - ky * y) * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdZ[i][j][k] = dFdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial second derivatives w.r.t. (x, y)
double[][][] d2FdXdY = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction d2FdXdY_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return -a * Math.cos(omega * z - kx * x - ky * y) * kx * ky;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdXdY[i][j][k] = d2FdXdY_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial second derivatives w.r.t. (x, z)
double[][][] d2FdXdZ = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction d2FdXdZ_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.cos(omega * z - kx * x - ky * y) * kx * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdXdZ[i][j][k] = d2FdXdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial second derivatives w.r.t. (y, z)
double[][][] d2FdYdZ = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction d2FdYdZ_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.cos(omega * z - kx * x - ky * y) * ky * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdYdZ[i][j][k] = d2FdYdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial third derivatives
double[][][] d3FdXdYdZ = new double[xval.length][yval.length][zval.length];
TrivariateRealFunction d3FdXdYdZ_f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.sin(omega * z - kx * x - ky * y) * kx * ky * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d3FdXdYdZ[i][j][k] = d3FdXdYdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
TrivariateRealFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, dFdX, dFdY, dFdZ,
d2FdXdY, d2FdXdZ, d2FdYdZ,
d3FdXdYdZ);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("On sample point",
expected, result, 1e-14);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.1);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.1);
}
}

211
src/test/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatorTest.java Normal file
View File

@ -0,0 +1,211 @@
/*
* 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.math.analysis.interpolation;
import org.apache.commons.math.MathException;
import org.apache.commons.math.DimensionMismatchException;
import org.apache.commons.math.analysis.TrivariateRealFunction;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for the tricubic interpolator.
*
* @version $Revision$ $Date$
*/
public final class TricubicSplineInterpolatorTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() throws MathException {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
double[][][] fval = new double[xval.length][yval.length][zval.length];
TrivariateRealGridInterpolator interpolator = new TricubicSplineInterpolator();
@SuppressWarnings("unused")
TrivariateRealFunction p = interpolator.interpolate(xval, yval, zval, fval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
p = interpolator.interpolate(wxval, yval, zval, fval);
Assert.fail("an exception should have been thrown");
} catch (IllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -3, -1, -1};
try {
p = interpolator.interpolate(xval, wyval, zval, fval);
Assert.fail("an exception should have been thrown");
} catch (IllegalArgumentException e) {
// Expected
}
double[] wzval = new double[] {-12, -8, -5.5, -3, -4, 2.5};
try {
p = interpolator.interpolate(xval, yval, wzval, fval);
Assert.fail("an exception should have been thrown");
} catch (IllegalArgumentException e) {
// Expected
}
double[][][] wfval = new double[xval.length][yval.length + 1][zval.length];
try {
p = interpolator.interpolate(xval, yval, zval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length - 1][yval.length][zval.length];
try {
p = interpolator.interpolate(xval, yval, zval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length][yval.length][zval.length - 1];
try {
p = interpolator.interpolate(xval, yval, zval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Test of interpolator for a plane.
* <p>
* f(x, y, z) = 2 x - 3 y - z + 5
*/
@Test
public void testPlane() throws MathException {
TrivariateRealFunction f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return 2 * x - 3 * y - z + 5;
}
};
TrivariateRealGridInterpolator interpolator = new TricubicSplineInterpolator();
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
TrivariateRealFunction p = interpolator.interpolate(xval, yval, zval, fval);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("On sample point", expected, result, 1e-15);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 0.3);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 0.3);
}
/**
* Test of interpolator for a sine wave.
* <p>
* <p>
* f(x, y, z) = a cos [&omega; z - k<sub>y</sub> x - k<sub>y</sub> y]
* </p>
* with A = 0.2, &omega; = 0.5, k<sub>x</sub> = 2, k<sub>y</sub> = 1.
*/
@Test
public void testWave() throws MathException {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 4};
final double a = 0.2;
final double omega = 0.5;
final double kx = 2;
final double ky = 1;
// Function values
TrivariateRealFunction f = new TrivariateRealFunction() {
public double value(double x, double y, double z) {
return a * Math.cos(omega * z - kx * x - ky * y);
}
};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
TrivariateRealGridInterpolator interpolator = new TricubicSplineInterpolator();
TrivariateRealFunction p = interpolator.interpolate(xval, yval, zval, fval);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("On sample point",
expected, result, 1e-12);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.1);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.1);
}
}