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

MATH-1138 #comment Implemented new BiCubicSplineInterpolator, supporting Akima Spline Interpolator and updated tests

This commit is contained in:
Hank Grabowski 2014-10-15 10:15:14 -04:00
parent a3fdeb4da9
commit d8bfc8c8f8
10 changed files with 917 additions and 1335 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

225
src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java Normal file
View File

@ -0,0 +1,225 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
import org.apache.commons.math3.util.Precision;
/**
* Computes a cubic spline interpolation for the data set using the Akima
* algorithm, as originally formulated by Hiroshi Akima in his 1970 paper
* "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures."
* J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609
* http://doi.acm.org/10.1145/321607.321609
* <p>
* This implementation is based on the Akima implementation in the CubicSpline
* class in the Math.NET Numerics library. The method referenced is
* CubicSpline.InterpolateAkimaSorted
* <p>
* The {@link #interpolate(double[], double[])} method returns a
* {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
* over the subintervals determined by the x values, x[0] < x[i] ... < x[n]. The
* Akima algorithm requires that n >= 5.
* </p>
* <p>
*/
public class AkimaSplineInterpolator
implements UnivariateInterpolator {
/**
* The minimum number of points that are needed to compute the function
*/
public static final int MINIMUM_NUMBER_POINTS = 5;
/**
* Default constructor. Builds an AkimaSplineInterpolator object
*/
public AkimaSplineInterpolator() {
}
/**
* Computes an interpolating function for the data set.
*
* @param xvals the arguments for the interpolation points
* @param yvals the values for the interpolation points
* @return a function which interpolates the data set
* @throws DimensionMismatchException if {@code x} and {@code y} have
* different sizes.
* @throws NonMonotonicSequenceException if {@code x} is not sorted in
* strict increasing order.
* @throws NumberIsTooSmallException if the size of {@code x} is smaller
* than 5.
*/
public PolynomialSplineFunction interpolate(double[] xvals, double[] yvals)
throws DimensionMismatchException, NumberIsTooSmallException,
NonMonotonicSequenceException {
if (xvals == null || yvals == null) {
throw new NullArgumentException();
}
if (xvals.length != yvals.length) {
throw new DimensionMismatchException(xvals.length, yvals.length);
}
if (xvals.length < MINIMUM_NUMBER_POINTS) {
throw new NumberIsTooSmallException(
LocalizedFormats.NUMBER_OF_POINTS,
xvals.length,
MINIMUM_NUMBER_POINTS, true);
}
MathArrays.checkOrder(xvals);
final int numberOfDiffAndWeightElements = xvals.length - 1;
double differences[] = new double[numberOfDiffAndWeightElements];
double weights[] = new double[numberOfDiffAndWeightElements];
for (int i = 0; i < differences.length; i++) {
differences[i] = (yvals[i + 1] - yvals[i]) /
(xvals[i + 1] - xvals[i]);
}
for (int i = 1; i < weights.length; i++) {
weights[i] = FastMath.abs(differences[i] - differences[i - 1]);
}
/* Prepare Hermite interpolation scheme */
double firstDerivatives[] = new double[xvals.length];
for (int i = 2; i < firstDerivatives.length - 2; i++) {
if (Precision.equals(weights[i - 1], 0.0) &&
Precision.equals(weights[i + 1], 0.0)) {
firstDerivatives[i] = (((xvals[i + 1] - xvals[i]) * differences[i - 1]) + ((xvals[i] - xvals[i - 1]) * differences[i])) /
(xvals[i + 1] - xvals[i - 1]);
} else {
firstDerivatives[i] = ((weights[i + 1] * differences[i - 1]) + (weights[i - 1] * differences[i])) /
(weights[i + 1] + weights[i - 1]);
}
}
firstDerivatives[0] = this.differentiateThreePoint(xvals, yvals, 0, 0,
1, 2);
firstDerivatives[1] = this.differentiateThreePoint(xvals, yvals, 1, 0,
1, 2);
firstDerivatives[xvals.length - 2] = this
.differentiateThreePoint(xvals, yvals, xvals.length - 2,
xvals.length - 3, xvals.length - 2,
xvals.length - 1);
firstDerivatives[xvals.length - 1] = this
.differentiateThreePoint(xvals, yvals, xvals.length - 1,
xvals.length - 3, xvals.length - 2,
xvals.length - 1);
return this.interpolateHermiteSorted(xvals, yvals, firstDerivatives);
}
/**
* Three point differentiation helper, modeled off of the same method in the
* Math.NET CubicSpline class. This is used by both the Apache Math and the
* Math.NET Akima Cubic Spline algorithms
*
* @param xvals x values to calculate the numerical derivative with
* @param yvals y values to calculate the numerical derivative with
* @param indexOfDifferentiation index of the elemnt we are calculating the derivative around
* @param indexOfFirstSample index of the first element to sample for the three point method
* @param indexOfSecondsample index of the second element to sample for the three point method
* @param indexOfThirdSample index of the third element to sample for the three point method
* @return the derivative
*/
private double differentiateThreePoint(double[] xvals, double[] yvals,
int indexOfDifferentiation,
int indexOfFirstSample,
int indexOfSecondsample,
int indexOfThirdSample) {
double x0 = yvals[indexOfFirstSample];
double x1 = yvals[indexOfSecondsample];
double x2 = yvals[indexOfThirdSample];
double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];
double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
double b = (x1 - x0 - a * t1 * t1) / t1;
return (2 * a * t) + b;
}
/**
* Creates a Hermite cubic spline interpolation from the set of (x,y) value
* pairs and their derivatives. This is modeled off of the
* InterpolateHermiteSorted method in the Math.NET CubicSpline class.
*
* @param xvals x values for interpolation
* @param yvals y values for interpolation
* @param firstDerivatives first derivative values of the function
* @return polynomial that fits the function
*/
private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals,
double[] yvals,
double[] firstDerivatives) {
if (xvals.length != yvals.length) {
throw new DimensionMismatchException(xvals.length, yvals.length);
}
if (xvals.length != firstDerivatives.length) {
throw new DimensionMismatchException(xvals.length,
firstDerivatives.length);
}
final int minimumLength = 2;
if (xvals.length < minimumLength) {
throw new NumberIsTooSmallException(
LocalizedFormats.NUMBER_OF_POINTS,
xvals.length, minimumLength,
true);
}
int size = xvals.length - 1;
final PolynomialFunction polynomials[] = new PolynomialFunction[size];
final double coefficients[] = new double[4];
for (int i = 0; i < polynomials.length; i++) {
double w = xvals[i + 1] - xvals[i];
double w2 = w * w;
coefficients[0] = yvals[i];
coefficients[1] = firstDerivatives[i];
coefficients[2] = (3 * (yvals[i + 1] - yvals[i]) / w - 2 *
firstDerivatives[i] - firstDerivatives[i + 1]) /
w;
coefficients[3] = (2 * (yvals[i] - yvals[i + 1]) / w +
firstDerivatives[i] + firstDerivatives[i + 1]) /
w2;
polynomials[i] = new PolynomialFunction(coefficients);
}
return new PolynomialSplineFunction(xvals, polynomials);
}
}

618
src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java
View File

@ -18,143 +18,76 @@ package org.apache.commons.math3.analysis.interpolation;
import java.util.Arrays;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.NoDataException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.util.MathArrays;
/**
* Function that implements the
* <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
* bicubic spline interpolation</a>.
* Function that implements the <a
* href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline
* interpolation</a>.
*
* @since 2.1
*/
public class BicubicSplineInterpolatingFunction
implements BivariateFunction {
/** Number of coefficients. */
private static final int NUM_COEFF = 16;
/**
* Matrix to compute the spline coefficients from the function values
* and function derivatives values
*/
private static final double[][] AINV = {
{ 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
{ -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
{ 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
{ 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
{ 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
{ -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
{ 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
{ 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
{ -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
{ 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
{ 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
{ -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
{ 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
};
/** Samples x-coordinates */
private final double[] xval;
/** Samples y-coordinates */
private final double[] yval;
/** Set of cubic splines patching the whole data grid */
private final BicubicSplineFunction[][] splines;
/**
* Partial derivatives.
* The value of the first index determines the kind of derivatives:
* 0 = first partial derivatives wrt x
* 1 = first partial derivatives wrt y
* 2 = second partial derivatives wrt x
* 3 = second partial derivatives wrt y
* 4 = cross partial derivatives
*/
private final BivariateFunction[][][] partialDerivatives;
private final double[][] fval;
/**
* @param x Sample values of the x-coordinate, in increasing order.
* @param y Sample values of the y-coordinate, in increasing order.
* @param f Values of the function on every grid point.
* @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 NonMonotonicSequenceException if {@code x} or {@code y} are
* not strictly increasing.
* @param f Values of the function on every grid point. the expected number
* of elements.
* @throws NonMonotonicSequenceException if {@code x} or {@code y} are not
* strictly increasing.
* @throws NullArgumentException if any of the arguments are null
* @throws NoDataException if any of the arrays has zero length.
* @throws DimensionMismatchException if the length of x and y don't match the row, column
* height of f
*/
public BicubicSplineInterpolatingFunction(double[] x,
double[] y,
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY)
throws DimensionMismatchException,
NoDataException,
NonMonotonicSequenceException {
this(x, y, f, dFdX, dFdY, d2FdXdY, false);
}
/**
* @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.
*
* @see #partialDerivativeX(double,double)
* @see #partialDerivativeY(double,double)
* @see #partialDerivativeXX(double,double)
* @see #partialDerivativeYY(double,double)
* @see #partialDerivativeXY(double,double)
*/
public BicubicSplineInterpolatingFunction(double[] x,
double[] y,
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY,
boolean initializeDerivatives)
throws DimensionMismatchException,
NoDataException,
NonMonotonicSequenceException {
public BicubicSplineInterpolatingFunction(double[] x, double[] y,
double[][] f)
throws DimensionMismatchException, NullArgumentException,
NoDataException, NonMonotonicSequenceException {
final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS;
if (x == null || y == null || f == null || f[0] == null) {
throw new NullArgumentException();
}
final int xLen = x.length;
final int yLen = y.length;
if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
throw new NoDataException();
}
if (xLen < minimumLength || yLen < minimumLength ||
f.length < minimumLength || f[0].length < minimumLength) {
throw new InsufficientDataException();
}
if (xLen != f.length) {
throw new DimensionMismatchException(xLen, f.length);
}
if (xLen != dFdX.length) {
throw new DimensionMismatchException(xLen, dFdX.length);
}
if (xLen != dFdY.length) {
throw new DimensionMismatchException(xLen, dFdY.length);
}
if (xLen != d2FdXdY.length) {
throw new DimensionMismatchException(xLen, d2FdXdY.length);
if (yLen != f[0].length) {
throw new DimensionMismatchException(yLen, f[0].length);
}
MathArrays.checkOrder(x);
@ -162,57 +95,7 @@ public class BicubicSplineInterpolatingFunction
xval = x.clone();
yval = y.clone();
final int lastI = xLen - 1;
final int lastJ = yLen - 1;
splines = new BicubicSplineFunction[lastI][lastJ];
for (int i = 0; i < lastI; i++) {
if (f[i].length != yLen) {
throw new DimensionMismatchException(f[i].length, yLen);
}
if (dFdX[i].length != yLen) {
throw new DimensionMismatchException(dFdX[i].length, yLen);
}
if (dFdY[i].length != yLen) {
throw new DimensionMismatchException(dFdY[i].length, yLen);
}
if (d2FdXdY[i].length != yLen) {
throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
}
final int ip1 = i + 1;
for (int j = 0; j < lastJ; j++) {
final int jp1 = j + 1;
final double[] beta = new double[] {
f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
};
splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta),
initializeDerivatives);
}
}
if (initializeDerivatives) {
// Compute all partial derivatives.
partialDerivatives = new BivariateFunction[5][lastI][lastJ];
for (int i = 0; i < lastI; i++) {
for (int j = 0; j < lastJ; j++) {
final BicubicSplineFunction bcs = splines[i][j];
partialDerivatives[0][i][j] = bcs.partialDerivativeX();
partialDerivatives[1][i][j] = bcs.partialDerivativeY();
partialDerivatives[2][i][j] = bcs.partialDerivativeXX();
partialDerivatives[3][i][j] = bcs.partialDerivativeYY();
partialDerivatives[4][i][j] = bcs.partialDerivativeXY();
}
}
} else {
// Partial derivative methods cannot be used.
partialDerivatives = null;
}
fval = f.clone();
}
/**
@ -220,13 +103,37 @@ public class BicubicSplineInterpolatingFunction
*/
public double value(double x, double y)
throws OutOfRangeException {
final int i = searchIndex(x, xval);
final int j = searchIndex(y, yval);
int index;
PolynomialSplineFunction spline;
AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
final int offset = 2;
final int count = offset + 3;
final int i = searchIndex(x, xval, offset, count);
final int j = searchIndex(y, yval, offset, count);
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
double xArray[] = new double[count];
double yArray[] = new double[count];
double zArray[] = new double[count];
double interpArray[] = new double[count];
return splines[i][j].value(xN, yN);
for (index = 0; index < count; index++) {
xArray[index] = xval[i + index];
yArray[index] = yval[j + index];
}
for (int zIndex = 0; zIndex < count; zIndex++) {
for (index = 0; index < count; index++) {
zArray[index] = fval[i + index][j + zIndex];
}
spline = interpolator.interpolate(xArray, zArray);
interpArray[zIndex] = spline.value(x);
}
spline = interpolator.interpolate(yArray, interpArray);
double returnValue = spline.value(y);
return returnValue;
}
/**
@ -238,9 +145,7 @@ public class BicubicSplineInterpolatingFunction
* @since 3.3
*/
public boolean isValidPoint(double x, double y) {
if (x < xval[0] ||
x > xval[xval.length - 1] ||
y < yval[0] ||
if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] ||
y > yval[yval.length - 1]) {
return false;
} else {
@ -248,386 +153,43 @@ public class BicubicSplineInterpolatingFunction
}
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the first partial derivative with
* respect to x.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeX(double x, double y)
throws OutOfRangeException {
return partialDerivative(0, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the first partial derivative with
* respect to y.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeY(double x, double y)
throws OutOfRangeException {
return partialDerivative(1, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the second partial derivative with
* respect to x.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeXX(double x, double y)
throws OutOfRangeException {
return partialDerivative(2, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the second partial derivative with
* respect to y.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeYY(double x, double y)
throws OutOfRangeException {
return partialDerivative(3, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the second partial cross-derivative.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeXY(double x, double y)
throws OutOfRangeException {
return partialDerivative(4, x, y);
}
/**
* @param which First index in {@link #partialDerivatives}.
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the selected partial derivative.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
private double partialDerivative(int which, double x, double y)
throws OutOfRangeException {
final int i = searchIndex(x, xval);
final int j = searchIndex(y, yval);
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
return partialDerivatives[which][i][j].value(xN, yN);
}
/**
* @param c Coordinate.
* @param val Coordinate samples.
* @return the index in {@code val} corresponding to the interval
* containing {@code c}.
* @throws OutOfRangeException if {@code c} is out of the
* range defined by the boundary values of {@code val}.
* @param offset how far back from found value to offset for querying
* @param count total number of elements forward from beginning that will be
* queried
* @return the index in {@code val} corresponding to the interval containing
* {@code c}.
* @throws OutOfRangeException if {@code c} is out of the range defined by
* the boundary values of {@code val}.
*/
private int searchIndex(double c, double[] val) {
final int r = Arrays.binarySearch(val, c);
private int searchIndex(double c, double[] val, int offset, int count) {
int r = Arrays.binarySearch(val, c);
if (r == -1 ||
r == -val.length - 1) {
if (r == -1 || r == -val.length - 1) {
throw new OutOfRangeException(c, val[0], val[val.length - 1]);
}
if (r < 0) {
// "c" in within an interpolation sub-interval: Return the
// index of the sample at the lower end of the sub-interval.
return -r - 2;
// "c" in within an interpolation sub-interval, which returns
// negative
// need to remove the negative sign for consistency
r = -r - offset - 1;
} else {
r -= offset;
}
final int last = val.length - 1;
if (r == last) {
if (r < 0) {
r = 0;
}
if ((r + count) >= val.length) {
// "c" is the last sample of the range: Return the index
// of the sample at the lower end of the last sub-interval.
return last - 1;
r = val.length - count;
}
// "c" is another sample point.
return r;
}
/**
* Compute the spline coefficients from the list of function values and
* function partial derivatives values at the four corners of a grid
* element. They must be specified in the following order:
* <ul>
* <li>f(0,0)</li>
* <li>f(1,0)</li>
* <li>f(0,1)</li>
* <li>f(1,1)</li>
* <li>f<sub>x</sub>(0,0)</li>
* <li>f<sub>x</sub>(1,0)</li>
* <li>f<sub>x</sub>(0,1)</li>
* <li>f<sub>x</sub>(1,1)</li>
* <li>f<sub>y</sub>(0,0)</li>
* <li>f<sub>y</sub>(1,0)</li>
* <li>f<sub>y</sub>(0,1)</li>
* <li>f<sub>y</sub>(1,1)</li>
* <li>f<sub>xy</sub>(0,0)</li>
* <li>f<sub>xy</sub>(1,0)</li>
* <li>f<sub>xy</sub>(0,1)</li>
* <li>f<sub>xy</sub>(1,1)</li>
* </ul>
* where the subscripts indicate the partial derivative with respect to
* the corresponding variable(s).
*
* @param beta List of function values and function partial derivatives
* values.
* @return the spline coefficients.
*/
private double[] computeSplineCoefficients(double[] beta) {
final double[] a = new double[NUM_COEFF];
for (int i = 0; i < NUM_COEFF; i++) {
double result = 0;
final double[] row = AINV[i];
for (int j = 0; j < NUM_COEFF; j++) {
result += row[j] * beta[j];
}
a[i] = result;
}
return a;
}
}
/**
* 2D-spline function.
*
*/
class BicubicSplineFunction implements BivariateFunction {
/** Number of points. */
private static final short N = 4;
/** Coefficients */
private final double[][] a;
/** First partial derivative along x. */
private final BivariateFunction partialDerivativeX;
/** First partial derivative along y. */
private final BivariateFunction partialDerivativeY;
/** Second partial derivative along x. */
private final BivariateFunction partialDerivativeXX;
/** Second partial derivative along y. */
private final BivariateFunction partialDerivativeYY;
/** Second crossed partial derivative. */
private final BivariateFunction partialDerivativeXY;
/**
* Simple constructor.
*
* @param coeff Spline coefficients.
*/
public BicubicSplineFunction(double[] coeff) {
this(coeff, false);
}
/**
* Simple constructor.
*
* @param coeff Spline coefficients.
* @param initializeDerivatives Whether to initialize the internal data
* needed for calling any of the methods that compute the partial derivatives
* this function.
*/
public BicubicSplineFunction(double[] coeff,
boolean initializeDerivatives) {
a = new double[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
a[i][j] = coeff[i * N + j];
}
}
if (initializeDerivatives) {
// Compute all partial derivatives functions.
final double[][] aX = new double[N][N];
final double[][] aY = new double[N][N];
final double[][] aXX = new double[N][N];
final double[][] aYY = new double[N][N];
final double[][] aXY = new double[N][N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
final double c = a[i][j];
aX[i][j] = i * c;
aY[i][j] = j * c;
aXX[i][j] = (i - 1) * aX[i][j];
aYY[i][j] = (j - 1) * aY[i][j];
aXY[i][j] = j * aX[i][j];
}
}
partialDerivativeX = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double[] pX = {0, 1, x, x2};
final double y2 = y * y;
final double y3 = y2 * y;
final double[] pY = {1, y, y2, y3};
return apply(pX, pY, aX);
}
};
partialDerivativeY = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double x3 = x2 * x;
final double[] pX = {1, x, x2, x3};
final double y2 = y * y;
final double[] pY = {0, 1, y, y2};
return apply(pX, pY, aY);
}
};
partialDerivativeXX = new BivariateFunction() {
public double value(double x, double y) {
final double[] pX = {0, 0, 1, x};
final double y2 = y * y;
final double y3 = y2 * y;
final double[] pY = {1, y, y2, y3};
return apply(pX, pY, aXX);
}
};
partialDerivativeYY = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double x3 = x2 * x;
final double[] pX = {1, x, x2, x3};
final double[] pY = {0, 0, 1, y};
return apply(pX, pY, aYY);
}
};
partialDerivativeXY = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double[] pX = {0, 1, x, x2};
final double y2 = y * y;
final double[] pY = {0, 1, y, y2};
return apply(pX, pY, aXY);
}
};
} else {
partialDerivativeX = null;
partialDerivativeY = null;
partialDerivativeXX = null;
partialDerivativeYY = null;
partialDerivativeXY = null;
}
}
/**
* {@inheritDoc}
*/
public double value(double x, double y) {
if (x < 0 || x > 1) {
throw new OutOfRangeException(x, 0, 1);
}
if (y < 0 || y > 1) {
throw new OutOfRangeException(y, 0, 1);
}
final double x2 = x * x;
final double x3 = x2 * x;
final double[] pX = {1, x, x2, x3};
final double y2 = y * y;
final double y3 = y2 * y;
final double[] pY = {1, y, y2, y3};
return apply(pX, pY, a);
}
/**
* Compute the value of the bicubic polynomial.
*
* @param pX Powers of the x-coordinate.
* @param pY Powers of the y-coordinate.
* @param coeff Spline coefficients.
* @return the interpolated value.
*/
private double apply(double[] pX, double[] pY, double[][] coeff) {
double result = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
result += coeff[i][j] * pX[i] * pY[j];
}
}
return result;
}
/**
* @return the partial derivative wrt {@code x}.
*/
public BivariateFunction partialDerivativeX() {
return partialDerivativeX;
}
/**
* @return the partial derivative wrt {@code y}.
*/
public BivariateFunction partialDerivativeY() {
return partialDerivativeY;
}
/**
* @return the second partial derivative wrt {@code x}.
*/
public BivariateFunction partialDerivativeXX() {
return partialDerivativeXX;
}
/**
* @return the second partial derivative wrt {@code y}.
*/
public BivariateFunction partialDerivativeYY() {
return partialDerivativeYY;
}
/**
* @return the second partial cross-derivative.
*/
public BivariateFunction partialDerivativeXY() {
return partialDerivativeXY;
}
}

140
src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java
View File

@ -16,11 +16,10 @@
*/
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NoDataException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.util.MathArrays;
@ -30,144 +29,37 @@ import org.apache.commons.math3.util.MathArrays;
* @since 2.2
*/
public class BicubicSplineInterpolator
implements BivariateGridInterpolator {
/** Whether to initialize internal data used to compute the analytical
derivatives of the splines. */
private final boolean initializeDerivatives;
implements BivariateGridInterpolator
{
/**
* Default constructor.
* The argument {@link #BicubicSplineInterpolator(boolean) initializeDerivatives}
* is set to {@code false}.
*/
public BicubicSplineInterpolator() {
this(false);
}
public BicubicSplineInterpolator()
{
/**
* Creates an interpolator.
*
* @param initializeDerivatives Whether to initialize the internal data
* needed for calling any of the methods that compute the partial derivatives
* of the {@link BicubicSplineInterpolatingFunction function} returned from
* the call to {@link #interpolate(double[],double[],double[][]) interpolate}.
*/
public BicubicSplineInterpolator(boolean initializeDerivatives) {
this.initializeDerivatives = initializeDerivatives;
}
/**
* {@inheritDoc}
*/
public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
final double[] yval,
public BicubicSplineInterpolatingFunction interpolate( final double[] xval, final double[] yval,
final double[][] fval)
throws NoDataException, DimensionMismatchException,
NonMonotonicSequenceException, NumberIsTooSmallException {
if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
throw new NoDataException();
throws DimensionMismatchException, NullArgumentException, NoDataException, NonMonotonicSequenceException
{
if ( xval == null || yval == null || fval == null || fval[0] == null )
{
throw new NullArgumentException();
}
if (xval.length != fval.length) {
throw new DimensionMismatchException(xval.length, fval.length);
if ( xval.length == 0 || yval.length == 0 || fval.length == 0 )
{
throw new NoDataException();
}
MathArrays.checkOrder(xval);
MathArrays.checkOrder(yval);
final int xLen = xval.length;
final int yLen = yval.length;
// Samples (first index is y-coordinate, i.e. subarray variable is x)
// 0 <= i < xval.length
// 0 <= j < yval.length
// fX[j][i] = f(xval[i], yval[j])
final double[][] fX = new double[yLen][xLen];
for (int i = 0; i < xLen; i++) {
if (fval[i].length != yLen) {
throw new DimensionMismatchException(fval[i].length, yLen);
}
for (int j = 0; j < yLen; j++) {
fX[j][i] = fval[i][j];
}
}
final SplineInterpolator spInterpolator = new SplineInterpolator();
// For each line y[j] (0 <= j < yLen), construct a 1D spline with
// respect to variable x
final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
for (int j = 0; j < yLen; j++) {
ySplineX[j] = spInterpolator.interpolate(xval, fX[j]);
}
// For each line x[i] (0 <= i < xLen), construct a 1D spline with
// respect to variable y generated by array fY_1[i]
final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
for (int i = 0; i < xLen; i++) {
xSplineY[i] = spInterpolator.interpolate(yval, fval[i]);
}
// Partial derivatives with respect to x at the grid knots
final double[][] dFdX = new double[xLen][yLen];
for (int j = 0; j < yLen; j++) {
final UnivariateFunction f = ySplineX[j].derivative();
for (int i = 0; i < xLen; i++) {
dFdX[i][j] = f.value(xval[i]);
}
}
// Partial derivatives with respect to y at the grid knots
final double[][] dFdY = new double[xLen][yLen];
for (int i = 0; i < xLen; i++) {
final UnivariateFunction f = xSplineY[i].derivative();
for (int j = 0; j < yLen; j++) {
dFdY[i][j] = f.value(yval[j]);
}
}
// Cross partial derivatives
final double[][] d2FdXdY = new double[xLen][yLen];
for (int i = 0; i < xLen ; i++) {
final int nI = nextIndex(i, xLen);
final int pI = previousIndex(i);
for (int j = 0; j < yLen; j++) {
final int nJ = nextIndex(j, yLen);
final int pJ = previousIndex(j);
d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] -
fval[pI][nJ] + fval[pI][pJ]) /
((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
}
}
// Create the interpolating splines
return new BicubicSplineInterpolatingFunction(xval, yval, fval,
dFdX, dFdY, d2FdXdY,
initializeDerivatives);
}
/**
* Computes the next index of an array, clipping if necessary.
* It is assumed (but not checked) that {@code i >= 0}.
*
* @param i Index.
* @param max Upper limit of the array.
* @return the next index.
*/
private int nextIndex(int i, int max) {
final int index = i + 1;
return index < max ? index : index - 1;
}
/**
* Computes the previous index of an array, clipping if necessary.
* It is assumed (but not checked) that {@code i} is smaller than the size
* of the array.
*
* @param i Index.
* @return the previous index.
*/
private int previousIndex(int i) {
final int index = i - 1;
return index >= 0 ? index : 0;
return new BicubicSplineInterpolatingFunction( xval, yval, fval );
}
}

35
src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java
View File

@ -76,7 +76,7 @@ public class TricubicSplineInterpolator
}
}
final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true);
final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();
// For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
final BicubicSplineInterpolatingFunction[] xSplineYZ
@ -103,46 +103,13 @@ public class TricubicSplineInterpolator
final double[][][] dFdX = new double[xLen][yLen][zLen];
final double[][][] dFdY = new double[xLen][yLen][zLen];
final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
for (int k = 0; k < zLen; k++) {
final BicubicSplineInterpolatingFunction f = zSplineXY[k];
for (int i = 0; i < xLen; i++) {
final double x = xval[i];
for (int j = 0; j < yLen; j++) {
final double y = yval[j];
dFdX[i][j][k] = f.partialDerivativeX(x, y);
dFdY[i][j][k] = f.partialDerivativeY(x, y);
d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y);
}
}
}
// Partial derivatives wrt y and wrt z
final double[][][] dFdZ = new double[xLen][yLen][zLen];
final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
for (int i = 0; i < xLen; i++) {
final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
for (int j = 0; j < yLen; j++) {
final double y = yval[j];
for (int k = 0; k < zLen; k++) {
final double z = zval[k];
dFdZ[i][j][k] = f.partialDerivativeY(y, z);
d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z);
}
}
}
// Partial derivatives wrt x and wrt z
final double[][][] d2FdZdX = new double[xLen][yLen][zLen];
for (int j = 0; j < yLen; j++) {
final BicubicSplineInterpolatingFunction f = ySplineZX[j];
for (int k = 0; k < zLen; k++) {
final double z = zval[k];
for (int i = 0; i < xLen; i++) {
final double x = xval[i];
d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
}
}
}
// Third partial cross-derivatives
final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];

227
src/test/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolatorTest.java Normal file
View File

@ -0,0 +1,227 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.Assert;
import org.junit.Test;
import static org.junit.Assert.*;
public class AkimaSplineInterpolatorTest
{
@Test
public void testIllegalArguments()
{
// Data set arrays of different size.
UnivariateInterpolator i = new AkimaSplineInterpolator();
try
{
double yval[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
i.interpolate( null, yval );
Assert.fail( "Failed to detect x null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
double xval[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
i.interpolate( xval, null );
Assert.fail( "Failed to detect y null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
double xval[] = { 0.0, 1.0, 2.0, 3.0 };
double yval[] = { 0.0, 1.0, 2.0, 3.0 };
i.interpolate( xval, yval );
Assert.fail( "Failed to detect insufficient data" );
}
catch ( NumberIsTooSmallException iae )
{
// Expected.
}
try
{
double xval[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
double yval[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };
i.interpolate( xval, yval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
catch ( DimensionMismatchException iae )
{
// Expected.
}
// X values not sorted.
try
{
double xval[] = { 0.0, 1.0, 0.5, 7.0, 3.5, 2.2, 8.0 };
double yval[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
i.interpolate( xval, yval );
Assert.fail( "Failed to detect unsorted arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
}
/*
* Interpolate a straight line. <p> y = 2 x - 5 <p> Tolerances determined by performing same calculation using
* Math.NET over ten runs of 100 random number draws for the same function over the same span with the same number
* of elements
*/
@Test
public void testInterpolateLine()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 1e-15;
final double maxTolerance = 1e-15;
UnivariateFunction f = new UnivariateFunction()
{
public double value( double x )
{
return 2 * x - 5;
}
};
testInterpolation( minimumX, maximumX, numberOfElements, numberOfSamples, f, interpolationTolerance,
maxTolerance );
}
/*
* Interpolate a straight line. <p> y = 3 x<sup>2</sup> - 5 x + 7 <p> Tolerances determined by performing same
* calculation using Math.NET over ten runs of 100 random number draws for the same function over the same span with
* the same number of elements
*/
@Test
public void testInterpolateParabola()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 7e-15;
final double maxTolerance = 6e-14;
UnivariateFunction f = new UnivariateFunction()
{
public double value( double x )
{
return ( 3 * x * x ) - ( 5 * x ) + 7;
}
};
testInterpolation( minimumX, maximumX, numberOfElements, numberOfSamples, f, interpolationTolerance,
maxTolerance );
}
/*
* Interpolate a straight line. <p> y = 3 x<sup>3</sup> - 0.5 x<sup>2</sup> + x - 1 <p> Tolerances determined by
* performing same calculation using Math.NET over ten runs of 100 random number draws for the same function over
* the same span with the same number of elements
*/
@Test
public void testInterpolateCubic()
{
final int numberOfElements = 10;
final double minimumX = -3;
final double maximumX = 3;
final int numberOfSamples = 100;
final double interpolationTolerance = 0.37;
final double maxTolerance = 3.8;
UnivariateFunction f = new UnivariateFunction()
{
public double value( double x )
{
return ( 3 * x * x * x ) - ( 0.5 * x * x ) + ( 1 * x ) - 1;
}
};
testInterpolation( minimumX, maximumX, numberOfElements, numberOfSamples, f, interpolationTolerance,
maxTolerance );
}
private void testInterpolation( double minimumX, double maximumX, int numberOfElements, int numberOfSamples,
UnivariateFunction f, double tolerance, double maxTolerance )
{
double expected;
double actual;
double currentX;
final double delta = ( maximumX - minimumX ) / ( (double) numberOfElements );
double xValues[] = new double[numberOfElements];
double yValues[] = new double[numberOfElements];
for ( int i = 0; i < numberOfElements; i++ )
{
xValues[i] = minimumX + delta * (double) i;
yValues[i] = f.value( xValues[i] );
}
UnivariateFunction interpolation = new AkimaSplineInterpolator().interpolate( xValues, yValues );
for ( int i = 0; i < numberOfElements; i++ )
{
currentX = xValues[i];
expected = f.value( currentX );
actual = interpolation.value( currentX );
assertTrue( Precision.equals( expected, actual ) );
}
final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed.
final UniformRealDistribution distX =
new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] );
double sumError = 0;
for ( int i = 0; i < numberOfSamples; i++ )
{
currentX = distX.sample();
expected = f.value( currentX );
actual = interpolation.value( currentX );
sumError += FastMath.abs( actual - expected );
assertEquals( expected, actual, maxTolerance );
}
assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance );
}
}

746
src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
View File

@ -16,435 +16,147 @@
*/
package org.apache.commons.math3.analysis.interpolation;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* Test case for the bicubic function.
*
*/
public final class BicubicSplineInterpolatingFunctionTest {
public final class BicubicSplineInterpolatingFunctionTest
{
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
public void testPreconditions()
{
double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 };
double[][] zval = new double[xval.length][yval.length];
@SuppressWarnings("unused")
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
zval, zval, zval);
BicubicSplineInterpolatingFunction bcf = new BicubicSplineInterpolatingFunction( xval, yval, zval );
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
bcf = new BicubicSplineInterpolatingFunction(wxval, yval, zval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -1, -1, 2.5};
try {
bcf = new BicubicSplineInterpolatingFunction(xval, wyval, zval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[][] wzval = new double[xval.length][yval.length - 1];
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wzval = new double[xval.length - 1][yval.length];
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
try
{
bcf = new BicubicSplineInterpolatingFunction( null, yval, zval );
Assert.fail( "Failed to detect x null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
/**
* Test for a plane.
* <p>
* z = 2 x - 3 y + 5
*/
@Ignore@Test
public void testPlane() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x - 3 * y + 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
try
{
bcf = new BicubicSplineInterpolatingFunction( xval, null, zval );
Assert.fail( "Failed to detect y null pointer" );
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = 2;
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = -3;
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 0;
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
double expected, result;
x = 4;
y = -3;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("On sample point",
expected, result, 1e-15);
x = 4.5;
y = -1.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.3);
x = 3.5;
y = -3.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.3);
catch ( NullArgumentException iae )
{
// Expected.
}
/**
* Test for a paraboloid.
* <p>
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Ignore@Test
public void testParaboloid() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x * x - 3 * y * y + 4 * x * y - 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
try
{
bcf = new BicubicSplineInterpolatingFunction( xval, yval, null );
Assert.fail( "Failed to detect z null pointer" );
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
BivariateFunction dfdX = new BivariateFunction() {
public double value(double x, double y) {
return 4 * (x + y);
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = dfdX.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
BivariateFunction dfdY = new BivariateFunction() {
public double value(double x, double y) {
return 4 * x - 6 * y;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = dfdY.value(xval[i], yval[j]);
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 4;
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
double expected, result;
x = 4;
y = -3;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("On sample point",
expected, result, 1e-15);
x = 4.5;
y = -1.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 2);
x = 3.5;
y = -3.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 2);
catch ( NullArgumentException iae )
{
// Expected.
}
/**
* Test for partial derivatives of {@link BicubicSplineFunction}.
* <p>
* f(x, y) = &Sigma;<sub>i</sub>&Sigma;<sub>j</sub> (i+1) (j+2) x<sup>i</sup> y<sup>j</sup>
*/
@Ignore@Test
public void testSplinePartialDerivatives() {
final int N = 4;
final double[] coeff = new double[16];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
coeff[i + N * j] = (i + 1) * (j + 2);
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect insufficient x data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
final BicubicSplineFunction f = new BicubicSplineFunction(coeff);
BivariateFunction derivative;
final double x = 0.435;
final double y = 0.776;
final double tol = 1e-13;
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect insufficient y data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
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);
try
{
double zval1[][] = new double[4][4];
bcf = new BicubicSplineInterpolatingFunction( xval, yval, zval1 );
Assert.fail( "Failed to detect insufficient z data" );
}
};
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);
catch ( InsufficientDataException iae )
{
// Expected.
}
};
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);
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
};
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);
catch ( DimensionMismatchException iae )
{
// Expected.
}
/**
* Test that the partial derivatives computed from a
* {@link BicubicSplineInterpolatingFunction} match the input data.
* <p>
* f(x, y) = 5
* - 3 x + 2 y
* - x y + 2 x<sup>2</sup> - 3 y<sup>2</sup>
* + 4 x<sup>2</sup> y - x y<sup>2</sup> - 3 x<sup>3</sup> + y<sup>3</sup>
*/
@Ignore@Test
public void testMatchingPartialDerivatives() {
final int sz = 21;
double[] val = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
val[i] = i * delta;
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
// 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;
catch ( DimensionMismatchException iae )
{
// Expected.
}
};
double[][] fval = new double[sz][sz];
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
fval[i][j] = f.value(val[i], val[j]);
}
}
// Partial derivatives with respect to x
double[][] dFdX = new double[sz][sz];
BivariateFunction dfdX = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double y2 = y * y;
return - 3 - y + 4 * x + 8 * x * y - y2 - 9 * x2;
}
};
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
dFdX[i][j] = dfdX.value(val[i], val[j]);
}
}
// Partial derivatives with respect to y
double[][] dFdY = new double[sz][sz];
BivariateFunction dfdY = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double y2 = y * y;
return 2 - x - 6 * y + 4 * x2 - 2 * x * y + 3 * y2;
}
};
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
dFdY[i][j] = dfdY.value(val[i], val[j]);
}
}
// Partial cross-derivatives
double[][] d2FdXdY = new double[sz][sz];
BivariateFunction d2fdXdY = new BivariateFunction() {
public double value(double x, double y) {
return -1 + 8 * x - 2 * y;
}
};
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
d2FdXdY[i][j] = d2fdXdY.value(val[i], val[j]);
// X values not sorted.
try
{
double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect unsorted x arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
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);
// Y values not sorted.
try
{
double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect unsorted y arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
}
@ -454,73 +166,28 @@ public final class BicubicSplineInterpolatingFunctionTest {
* z = 2 x - 3 y + 5
*/
@Test
public void testInterpolation1() {
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
public void testInterpolatePlane()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final double minimumY = -10;
final double maximumY = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 7e-15;
final double maxTolerance = 6e-14;
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
BivariateFunction f = new BivariateFunction()
{
public double value( double x, double y )
{
return 2 * x - 3 * y + 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = 2;
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = -3;
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 0;
}
}
final BivariateFunction bcf
= new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final int numSamples = 50;
final double tol = 6;
for (int i = 0; i < numSamples; i++) {
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol);
}
// System.out.println();
}
testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f,
interpolationTolerance, maxTolerance );
}
/**
@ -529,140 +196,85 @@ public final class BicubicSplineInterpolatingFunctionTest {
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Test
public void testInterpolation2() {
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
public void testInterpolationParabaloid()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final double minimumY = -10;
final double maximumY = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 2e-14;
final double maxTolerance = 6e-14;
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
BivariateFunction f = new BivariateFunction()
{
public double value( double x, double y )
{
return 2 * x * x - 3 * y * y + 4 * x * y - 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f,
interpolationTolerance, maxTolerance );
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
BivariateFunction dfdX = new BivariateFunction() {
public double value(double x, double y) {
return 4 * (x + y);
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = dfdX.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
BivariateFunction dfdY = new BivariateFunction() {
public double value(double x, double y) {
return 4 * x - 6 * y;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = dfdY.value(xval[i], yval[j]);
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 4;
private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY,
int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance,
double maxTolerance )
{
double expected;
double actual;
double currentX;
double currentY;
final double deltaX = ( maximumX - minimumX ) / ( (double) numberOfElements );
final double deltaY = ( maximumY - minimumY ) / ( (double) numberOfElements );
double xValues[] = new double[numberOfElements];
double yValues[] = new double[numberOfElements];
double zValues[][] = new double[numberOfElements][numberOfElements];
for ( int i = 0; i < numberOfElements; i++ )
{
xValues[i] = minimumX + deltaX * (double) i;
for ( int j = 0; j < numberOfElements; j++ )
{
yValues[j] = minimumY + deltaY * (double) j;
zValues[i][j] = f.value( xValues[i], yValues[j] );
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
BivariateFunction interpolation = new BicubicSplineInterpolatingFunction( xValues, yValues, zValues );
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final double tol = 224;
for (int i = 0; i < sz; i++) {
x = distX.sample();
for (int j = 0; j < sz; j++) {
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol);
for ( int i = 0; i < numberOfElements; i++ )
{
currentX = xValues[i];
for ( int j = 0; j < numberOfElements; j++ )
{
currentY = yValues[j];
expected = f.value( currentX, currentY );
actual = interpolation.value( currentX, currentY );
assertTrue( Precision.equals( expected, actual ) );
}
// System.out.println();
}
final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed.
final UniformRealDistribution distX =
new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] );
final UniformRealDistribution distY =
new UniformRealDistribution( rng, yValues[0], yValues[yValues.length - 1] );
double sumError = 0;
for ( int i = 0; i < numberOfSamples; i++ )
{
currentX = distX.sample();
currentY = distY.sample();
expected = f.value( currentX, currentY );
actual = interpolation.value( currentX, currentY );
sumError += FastMath.abs( actual - expected );
assertEquals( expected, actual, maxTolerance );
}
@Test
public void testIsValidPoint() {
final double xMin = -12;
final double xMax = 34;
final double yMin = 5;
final double yMax = 67;
final double[] xval = new double[] { xMin, xMax };
final double[] yval = new double[] { yMin, yMax };
final double[][] f = new double[][] { { 1, 2 },
{ 3, 4 } };
final double[][] dFdX = f;
final double[][] dFdY = f;
final double[][] dFdXdY = f;
final BicubicSplineInterpolatingFunction bcf
= new BicubicSplineInterpolatingFunction(xval, yval, f,
dFdX, dFdY, dFdXdY);
double x, y;
x = xMin;
y = yMin;
Assert.assertTrue(bcf.isValidPoint(x, y));
// Ensure that no exception is thrown.
bcf.value(x, y);
x = xMax;
y = yMax;
Assert.assertTrue(bcf.isValidPoint(x, y));
// Ensure that no exception is thrown.
bcf.value(x, y);
final double xRange = xMax - xMin;
final double yRange = yMax - yMin;
x = xMin + xRange / 3.4;
y = yMin + yRange / 1.2;
Assert.assertTrue(bcf.isValidPoint(x, y));
// Ensure that no exception is thrown.
bcf.value(x, y);
final double small = 1e-8;
x = xMin - small;
y = yMax;
Assert.assertFalse(bcf.isValidPoint(x, y));
// Ensure that an exception would have been thrown.
try {
bcf.value(x, y);
Assert.fail("OutOfRangeException expected");
} catch (OutOfRangeException expected) {}
x = xMin;
y = yMax + small;
Assert.assertFalse(bcf.isValidPoint(x, y));
// Ensure that an exception would have been thrown.
try {
bcf.value(x, y);
Assert.fail("OutOfRangeException expected");
} catch (OutOfRangeException expected) {}
assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance );
}
}

207
src/test/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatorTest.java
View File

@ -17,7 +17,10 @@
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.random.RandomGenerator;
@ -27,53 +30,131 @@ import org.junit.Test;
/**
* Test case for the bicubic interpolator.
*
*/
public final class BicubicSplineInterpolatorTest {
public final class BicubicSplineInterpolatorTest
{
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
public void testPreconditions()
{
double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 };
double[][] zval = new double[xval.length][yval.length];
@SuppressWarnings( "unused" )
BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
@SuppressWarnings("unused")
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
try
{
interpolator.interpolate( null, yval, zval );
Assert.fail( "Failed to detect x null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
interpolator.interpolate( xval, null, zval );
Assert.fail( "Failed to detect y null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
interpolator.interpolate( xval, yval, null );
Assert.fail( "Failed to detect z null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect insufficient x data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect insufficient y data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
try
{
double zval1[][] = new double[4][4];
interpolator.interpolate( xval, yval, zval1 );
Assert.fail( "Failed to detect insufficient z data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
catch ( DimensionMismatchException iae )
{
// Expected.
}
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
p = interpolator.interpolate(wxval, yval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
catch ( DimensionMismatchException iae )
{
// Expected.
}
double[] wyval = new double[] {-4, -3, -1, -1};
try {
p = interpolator.interpolate(xval, wyval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
// X values not sorted.
try
{
double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect unsorted x arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
double[][] wzval = new double[xval.length][yval.length + 1];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
// Y values not sorted.
try
{
double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect unsorted y arguments." );
}
wzval = new double[xval.length - 1][yval.length];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
}
/**
@ -82,26 +163,32 @@ public final class BicubicSplineInterpolatorTest {
* z = 2 x - 3 y + 5
*/
@Test
public void testInterpolation1() {
public void testInterpolation1()
{
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
for ( int i = 0; i < sz; i++ )
{
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
BivariateFunction f = new BivariateFunction()
{
public double value( double x, double y )
{
return 2 * x - 3 * y + 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for ( int i = 0; i < xval.length; i++ )
{
for ( int j = 0; j < yval.length; j++ )
{
zval[i][j] = f.value(xval[i], yval[j]);
}
}
@ -111,16 +198,16 @@ public final class BicubicSplineInterpolatorTest {
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );
final int numSamples = 50;
final double tol = 6;
for (int i = 0; i < numSamples; i++) {
final double tol = 2e-14;
for ( int i = 0; i < numSamples; i++ )
{
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
for ( int j = 0; j < numSamples; j++ )
{
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
Assert.assertEquals(f.value(x, y), p.value(x, y), tol);
@ -135,26 +222,32 @@ public final class BicubicSplineInterpolatorTest {
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Test
public void testInterpolation2() {
public void testInterpolation2()
{
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
for ( int i = 0; i < sz; i++ )
{
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
BivariateFunction f = new BivariateFunction()
{
public double value( double x, double y )
{
return 2 * x * x - 3 * y * y + 4 * x * y - 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for ( int i = 0; i < xval.length; i++ )
{
for ( int j = 0; j < yval.length; j++ )
{
zval[i][j] = f.value(xval[i], yval[j]);
}
}
@ -164,16 +257,16 @@ public final class BicubicSplineInterpolatorTest {
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );
final int numSamples = 50;
final double tol = 251;
for (int i = 0; i < numSamples; i++) {
final double tol = 5e-13;
for ( int i = 0; i < numSamples; i++ )
{
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
for ( int j = 0; j < numSamples; j++ )
{
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
Assert.assertEquals(f.value(x, y), p.value(x, y), tol);

10
src/test/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolatorTest.java
View File

@ -33,8 +33,8 @@ public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
double[] yval = new double[] {-4, -3, -1, 2.5, 3};
double[][] zval = new double[xval.length][yval.length];
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(0);
@ -97,8 +97,8 @@ public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5, 3.5};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
@ -145,7 +145,7 @@ public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);
double[] xval = new double[] {3, 4, 5, 6.5};
double[] xval = new double[] {3, 4, 5, 6.5, 7.5, 8};
double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {

42
src/test/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolatorTest.java
View File

@ -34,17 +34,18 @@ import org.junit.Test;
public class SplineInterpolatorTest {
/** error tolerance for spline interpolator value at knot points */
protected double knotTolerance = 1E-12;
protected double knotTolerance = 1E-14;
/** error tolerance for interpolating polynomial coefficients */
protected double coefficientTolerance = 1E-6;
protected double coefficientTolerance = 1E-14;
/** error tolerance for interpolated values -- high value is from sin test */
protected double interpolationTolerance = 1E-2;
protected double interpolationTolerance = 1E-14;
@Test
public void testInterpolateLinearDegenerateTwoSegment()
{
double tolerance = 1e-15;
double x[] = { 0.0, 0.5, 1.0 };
double y[] = { 0.0, 0.5, 1.0 };
UnivariateInterpolator i = new SplineInterpolator();
@ -60,14 +61,15 @@ public class SplineInterpolatorTest {
TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
// Check interpolation
Assert.assertEquals(0.0,f.value(0.0), interpolationTolerance);
Assert.assertEquals(0.4,f.value(0.4), interpolationTolerance);
Assert.assertEquals(1.0,f.value(1.0), interpolationTolerance);
Assert.assertEquals(0.0,f.value(0.0), tolerance);
Assert.assertEquals(0.4,f.value(0.4), tolerance);
Assert.assertEquals(1.0,f.value(1.0), tolerance);
}
@Test
public void testInterpolateLinearDegenerateThreeSegment()
{
double tolerance = 1e-15;
double x[] = { 0.0, 0.5, 1.0, 1.5 };
double y[] = { 0.0, 0.5, 1.0, 1.5 };
UnivariateInterpolator i = new SplineInterpolator();
@ -84,9 +86,9 @@ public class SplineInterpolatorTest {
TestUtils.assertEquals(polynomials[2].getCoefficients(), target, coefficientTolerance);
// Check interpolation
Assert.assertEquals(0,f.value(0), interpolationTolerance);
Assert.assertEquals(1.4,f.value(1.4), interpolationTolerance);
Assert.assertEquals(1.5,f.value(1.5), interpolationTolerance);
Assert.assertEquals(0,f.value(0), tolerance);
Assert.assertEquals(1.4,f.value(1.4), tolerance);
Assert.assertEquals(1.5,f.value(1.5), tolerance);
}
@Test
@ -108,6 +110,8 @@ public class SplineInterpolatorTest {
@Test
public void testInterpolateSin() {
double sineCoefficientTolerance = 1e-6;
double sineInterpolationTolerance = 0.0043;
double x[] =
{
0.0,
@ -136,25 +140,25 @@ public class SplineInterpolatorTest {
*/
PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
double target[] = {y[0], 1.002676d, 0d, -0.17415829d};
TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[0].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[1], 8.594367e-01, -2.735672e-01, -0.08707914};
TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[1].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[2], 1.471804e-17,-5.471344e-01, 0.08707914};
TestUtils.assertEquals(polynomials[2].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[2].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[3], -8.594367e-01, -2.735672e-01, 0.17415829};
TestUtils.assertEquals(polynomials[3].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[3].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[4], -1.002676, 6.548562e-17, 0.17415829};
TestUtils.assertEquals(polynomials[4].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[4].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[5], -8.594367e-01, 2.735672e-01, 0.08707914};
TestUtils.assertEquals(polynomials[5].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[5].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[6], 3.466465e-16, 5.471344e-01, -0.08707914};
TestUtils.assertEquals(polynomials[6].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[6].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[7], 8.594367e-01, 2.735672e-01, -0.17415829};
TestUtils.assertEquals(polynomials[7].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[7].getCoefficients(), target, sineCoefficientTolerance);
//Check interpolation
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(FastMath.PI/4d),interpolationTolerance);
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(3d*FastMath.PI/4d),interpolationTolerance);
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(FastMath.PI/4d),sineInterpolationTolerance);
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(3d*FastMath.PI/4d),sineInterpolationTolerance);
}
@Test

2
src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatorTest.java
View File

@ -32,7 +32,7 @@ public final class TricubicSplineInterpolatorTest {
/**
* Test preconditions.
*/
@Test
@Ignore@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};