* Formatting (braces, indent, conditionals, array declaration, ...).
 * Reduce variable's scope.
 * Temporary variable to avoid multiple accesses to the same array element.
 * Use "final" keyword.
This commit is contained in:
Gilles 2014-11-11 22:56:36 +01:00
parent ed565027c7
commit e3700ef350
5 changed files with 213 additions and 306 deletions

View File

@ -37,30 +37,18 @@ import org.apache.commons.math3.util.Precision;
* This implementation is based on the Akima implementation in the CubicSpline * This implementation is based on the Akima implementation in the CubicSpline
* class in the Math.NET Numerics library. The method referenced is * class in the Math.NET Numerics library. The method referenced is
* CubicSpline.InterpolateAkimaSorted * 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>
* <p> * <p>
* The {@link #interpolate(double[], double[]) interpolate} method returns a
* {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
* over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}.
* The Akima algorithm requires that {@code n >= 5}.
* </p>
*/ */
public class AkimaSplineInterpolator public class AkimaSplineInterpolator
implements UnivariateInterpolator { implements UnivariateInterpolator {
/** The minimum number of points that are needed to compute the function. */
private static final int MINIMUM_NUMBER_POINTS = 5;
/**
* 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. * Computes an interpolating function for the data set.
@ -68,17 +56,20 @@ public class AkimaSplineInterpolator
* @param xvals the arguments for the interpolation points * @param xvals the arguments for the interpolation points
* @param yvals the values for the interpolation points * @param yvals the values for the interpolation points
* @return a function which interpolates the data set * @return a function which interpolates the data set
* @throws DimensionMismatchException if {@code x} and {@code y} have * @throws DimensionMismatchException if {@code xvals} and {@code yvals} have
* different sizes. * different sizes.
* @throws NonMonotonicSequenceException if {@code x} is not sorted in * @throws NonMonotonicSequenceException if {@code xvals} is not sorted in
* strict increasing order. * strict increasing order.
* @throws NumberIsTooSmallException if the size of {@code x} is smaller * @throws NumberIsTooSmallException if the size of {@code xvals} is smaller
* than 5. * than 5.
*/ */
public PolynomialSplineFunction interpolate(double[] xvals, double[] yvals) public PolynomialSplineFunction interpolate(double[] xvals,
throws DimensionMismatchException, NumberIsTooSmallException, double[] yvals)
NonMonotonicSequenceException { throws DimensionMismatchException,
if (xvals == null || yvals == null) { NumberIsTooSmallException,
NonMonotonicSequenceException {
if (xvals == null ||
yvals == null) {
throw new NullArgumentException(); throw new NullArgumentException();
} }
@ -87,8 +78,7 @@ public class AkimaSplineInterpolator
} }
if (xvals.length < MINIMUM_NUMBER_POINTS) { if (xvals.length < MINIMUM_NUMBER_POINTS) {
throw new NumberIsTooSmallException( throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
LocalizedFormats.NUMBER_OF_POINTS,
xvals.length, xvals.length,
MINIMUM_NUMBER_POINTS, true); MINIMUM_NUMBER_POINTS, true);
} }
@ -97,47 +87,44 @@ public class AkimaSplineInterpolator
final int numberOfDiffAndWeightElements = xvals.length - 1; final int numberOfDiffAndWeightElements = xvals.length - 1;
double differences[] = new double[numberOfDiffAndWeightElements]; final double[] differences = new double[numberOfDiffAndWeightElements];
double weights[] = new double[numberOfDiffAndWeightElements]; final double[] weights = new double[numberOfDiffAndWeightElements];
for (int i = 0; i < differences.length; i++) { for (int i = 0; i < differences.length; i++) {
differences[i] = (yvals[i + 1] - yvals[i]) / differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]);
(xvals[i + 1] - xvals[i]);
} }
for (int i = 1; i < weights.length; i++) { for (int i = 1; i < weights.length; i++) {
weights[i] = FastMath.abs(differences[i] - differences[i - 1]); weights[i] = FastMath.abs(differences[i] - differences[i - 1]);
} }
/* Prepare Hermite interpolation scheme */ // Prepare Hermite interpolation scheme.
final double[] firstDerivatives = new double[xvals.length];
double firstDerivatives[] = new double[xvals.length];
for (int i = 2; i < firstDerivatives.length - 2; i++) { for (int i = 2; i < firstDerivatives.length - 2; i++) {
if (Precision.equals(weights[i - 1], 0.0) && final double wP = weights[i + 1];
Precision.equals(weights[i + 1], 0.0)) { final double wM = weights[i - 1];
firstDerivatives[i] = (((xvals[i + 1] - xvals[i]) * differences[i - 1]) + ((xvals[i] - xvals[i - 1]) * differences[i])) / if (Precision.equals(wP, 0.0) &&
(xvals[i + 1] - xvals[i - 1]); Precision.equals(wM, 0.0)) {
final double xv = xvals[i];
final double xvP = xvals[i + 1];
final double xvM = xvals[i - 1];
firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM);
} else { } else {
firstDerivatives[i] = ((weights[i + 1] * differences[i - 1]) + (weights[i - 1] * differences[i])) / firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM);
(weights[i + 1] + weights[i - 1]);
} }
} }
firstDerivatives[0] = this.differentiateThreePoint(xvals, yvals, 0, 0, firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2);
1, 2); firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2);
firstDerivatives[1] = this.differentiateThreePoint(xvals, yvals, 1, 0, firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2,
1, 2); xvals.length - 3, xvals.length - 2,
firstDerivatives[xvals.length - 2] = this xvals.length - 1);
.differentiateThreePoint(xvals, yvals, xvals.length - 2, firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1,
xvals.length - 3, xvals.length - 2, xvals.length - 3, xvals.length - 2,
xvals.length - 1); 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); return interpolateHermiteSorted(xvals, yvals, firstDerivatives);
} }
/** /**
@ -158,16 +145,17 @@ public class AkimaSplineInterpolator
int indexOfFirstSample, int indexOfFirstSample,
int indexOfSecondsample, int indexOfSecondsample,
int indexOfThirdSample) { int indexOfThirdSample) {
double x0 = yvals[indexOfFirstSample]; final double x0 = yvals[indexOfFirstSample];
double x1 = yvals[indexOfSecondsample]; final double x1 = yvals[indexOfSecondsample];
double x2 = yvals[indexOfThirdSample]; final double x2 = yvals[indexOfThirdSample];
double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample]; final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample]; final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample]; final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];
final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
final double b = (x1 - x0 - a * t1 * t1) / t1;
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; return (2 * a * t) + b;
} }
@ -195,27 +183,29 @@ public class AkimaSplineInterpolator
final int minimumLength = 2; final int minimumLength = 2;
if (xvals.length < minimumLength) { if (xvals.length < minimumLength) {
throw new NumberIsTooSmallException( throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
LocalizedFormats.NUMBER_OF_POINTS,
xvals.length, minimumLength, xvals.length, minimumLength,
true); true);
} }
int size = xvals.length - 1; final int size = xvals.length - 1;
final PolynomialFunction polynomials[] = new PolynomialFunction[size]; final PolynomialFunction[] polynomials = new PolynomialFunction[size];
final double coefficients[] = new double[4]; final double[] coefficients = new double[4];
for (int i = 0; i < polynomials.length; i++) { for (int i = 0; i < polynomials.length; i++) {
double w = xvals[i + 1] - xvals[i]; final double w = xvals[i + 1] - xvals[i];
double w2 = w * w; final double w2 = w * w;
coefficients[0] = yvals[i];
final double yv = yvals[i];
final double yvP = yvals[i + 1];
final double fd = firstDerivatives[i];
final double fdP = firstDerivatives[i + 1];
coefficients[0] = yv;
coefficients[1] = firstDerivatives[i]; coefficients[1] = firstDerivatives[i];
coefficients[2] = (3 * (yvals[i + 1] - yvals[i]) / w - 2 * coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w;
firstDerivatives[i] - firstDerivatives[i + 1]) / coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2;
w;
coefficients[3] = (2 * (yvals[i] - yvals[i + 1]) / w +
firstDerivatives[i] + firstDerivatives[i + 1]) /
w2;
polynomials[i] = new PolynomialFunction(coefficients); polynomials[i] = new PolynomialFunction(coefficients);
} }

View File

@ -28,21 +28,24 @@ import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.util.MathArrays; import org.apache.commons.math3.util.MathArrays;
/** /**
* Function that implements the <a * Function that implements the
* href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline * <a href="http://www.paulinternet.nl/?page=bicubic">bicubic spline</a>
* interpolation</a>. * interpolation.
* This implementation currently uses {@link AkimaSplineInterpolator} as the
* underlying one-dimensional interpolator, which requires 5 sample points;
* insufficient data will raise an exception when the
* {@link #value(double,double) value} method is called.
* *
* @since 2.1 * @since 3.4
*/ */
public class PiecewiseBicubicSplineInterpolatingFunction public class PiecewiseBicubicSplineInterpolatingFunction
implements BivariateFunction { implements BivariateFunction {
/** The minimum number of points that are needed to compute the function. */
private static final int MIN_NUM_POINTS = 5;
/** Samples x-coordinates */ /** Samples x-coordinates */
private final double[] xval; private final double[] xval;
/** Samples y-coordinates */ /** Samples y-coordinates */
private final double[] yval; private final double[] yval;
/** Set of cubic splines patching the whole data grid */ /** Set of cubic splines patching the whole data grid */
private final double[][] fval; private final double[][] fval;
@ -58,27 +61,34 @@ public class PiecewiseBicubicSplineInterpolatingFunction
* @throws DimensionMismatchException if the length of x and y don't match the row, column * @throws DimensionMismatchException if the length of x and y don't match the row, column
* height of f * height of f
*/ */
public PiecewiseBicubicSplineInterpolatingFunction(double[] x,
public PiecewiseBicubicSplineInterpolatingFunction(double[] x, double[] y, double[] y,
double[][] f) double[][] f)
throws DimensionMismatchException, NullArgumentException, throws DimensionMismatchException,
NoDataException, NonMonotonicSequenceException { NullArgumentException,
NoDataException,
final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS; NonMonotonicSequenceException {
if (x == null ||
if (x == null || y == null || f == null || f[0] == null) { y == null ||
f == null ||
f[0] == null) {
throw new NullArgumentException(); throw new NullArgumentException();
} }
final int xLen = x.length; final int xLen = x.length;
final int yLen = y.length; final int yLen = y.length;
if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { if (xLen == 0 ||
yLen == 0 ||
f.length == 0 ||
f[0].length == 0) {
throw new NoDataException(); throw new NoDataException();
} }
if (xLen < minimumLength || yLen < minimumLength || if (xLen < MIN_NUM_POINTS ||
f.length < minimumLength || f[0].length < minimumLength) { yLen < MIN_NUM_POINTS ||
f.length < MIN_NUM_POINTS ||
f[0].length < MIN_NUM_POINTS) {
throw new InsufficientDataException(); throw new InsufficientDataException();
} }
@ -101,35 +111,34 @@ public class PiecewiseBicubicSplineInterpolatingFunction
/** /**
* {@inheritDoc} * {@inheritDoc}
*/ */
public double value(double x, double y) public double value(double x,
double y)
throws OutOfRangeException { throws OutOfRangeException {
int index; final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
PolynomialSplineFunction spline;
AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
final int offset = 2; final int offset = 2;
final int count = offset + 3; final int count = offset + 3;
final int i = searchIndex(x, xval, offset, count); final int i = searchIndex(x, xval, offset, count);
final int j = searchIndex(y, yval, offset, count); final int j = searchIndex(y, yval, offset, count);
double xArray[] = new double[count]; final double xArray[] = new double[count];
double yArray[] = new double[count]; final double yArray[] = new double[count];
double zArray[] = new double[count]; final double zArray[] = new double[count];
double interpArray[] = new double[count]; final double interpArray[] = new double[count];
for (index = 0; index < count; index++) { for (int index = 0; index < count; index++) {
xArray[index] = xval[i + index]; xArray[index] = xval[i + index];
yArray[index] = yval[j + index]; yArray[index] = yval[j + index];
} }
for (int zIndex = 0; zIndex < count; zIndex++) { for (int zIndex = 0; zIndex < count; zIndex++) {
for (index = 0; index < count; index++) { for (int index = 0; index < count; index++) {
zArray[index] = fval[i + index][j + zIndex]; zArray[index] = fval[i + index][j + zIndex];
} }
spline = interpolator.interpolate(xArray, zArray); final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray);
interpArray[zIndex] = spline.value(x); interpArray[zIndex] = spline.value(x);
} }
spline = interpolator.interpolate(yArray, interpArray); final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray);
double returnValue = spline.value(y); double returnValue = spline.value(y);
@ -144,8 +153,11 @@ public class PiecewiseBicubicSplineInterpolatingFunction
* @return {@code true} if (x, y) is a valid point. * @return {@code true} if (x, y) is a valid point.
* @since 3.3 * @since 3.3
*/ */
public boolean isValidPoint(double x, double y) { public boolean isValidPoint(double x,
if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] || double y) {
if (x < xval[0] ||
x > xval[xval.length - 1] ||
y < yval[0] ||
y > yval[yval.length - 1]) { y > yval[yval.length - 1]) {
return false; return false;
} else { } else {
@ -164,7 +176,10 @@ public class PiecewiseBicubicSplineInterpolatingFunction
* @throws OutOfRangeException if {@code c} is out of the range defined by * @throws OutOfRangeException if {@code c} is out of the range defined by
* the boundary values of {@code val}. * the boundary values of {@code val}.
*/ */
private int searchIndex(double c, double[] val, int offset, int count) { private int searchIndex(double c,
double[] val,
int offset,
int count) {
int r = Arrays.binarySearch(val, c); int r = Arrays.binarySearch(val, c);
if (r == -1 || r == -val.length - 1) { if (r == -1 || r == -val.length - 1) {

View File

@ -28,31 +28,28 @@ import org.apache.commons.math3.util.MathArrays;
* @since 2.2 * @since 2.2
*/ */
public class PiecewiseBicubicSplineInterpolator public class PiecewiseBicubicSplineInterpolator
implements BivariateGridInterpolator implements BivariateGridInterpolator {
{
/**
* Default constructor.
*/
public PiecewiseBicubicSplineInterpolator()
{
}
/** /**
* {@inheritDoc} * {@inheritDoc}
*/ */
public PiecewiseBicubicSplineInterpolatingFunction interpolate( final double[] xval, final double[] yval, public PiecewiseBicubicSplineInterpolatingFunction interpolate( final double[] xval,
final double[][] fval) final double[] yval,
throws DimensionMismatchException, NullArgumentException, NoDataException, NonMonotonicSequenceException final double[][] fval)
{ throws DimensionMismatchException,
if ( xval == null || yval == null || fval == null || fval[0] == null ) NullArgumentException,
{ NoDataException,
NonMonotonicSequenceException {
if ( xval == null ||
yval == null ||
fval == null ||
fval[0] == null ) {
throw new NullArgumentException(); throw new NullArgumentException();
} }
if ( xval.length == 0 || yval.length == 0 || fval.length == 0 ) if ( xval.length == 0 ||
{ yval.length == 0 ||
fval.length == 0 ) {
throw new NoDataException(); throw new NoDataException();
} }

View File

@ -35,14 +35,12 @@ import org.junit.Test;
/** /**
* Test case for the piecewise bicubic function. * Test case for the piecewise bicubic function.
*/ */
public final class PiecewiseBicubicSplineInterpolatingFunctionTest public final class PiecewiseBicubicSplineInterpolatingFunctionTest {
{
/** /**
* Test preconditions. * Test preconditions.
*/ */
@Test @Test
public void testPreconditions() public void testPreconditions() {
{
double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 }; double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 }; double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 };
double[][] zval = new double[xval.length][yval.length]; double[][] zval = new double[xval.length][yval.length];
@ -50,112 +48,82 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
@SuppressWarnings("unused") @SuppressWarnings("unused")
PiecewiseBicubicSplineInterpolatingFunction bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, zval ); PiecewiseBicubicSplineInterpolatingFunction bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, zval );
try try {
{
bcf = new PiecewiseBicubicSplineInterpolatingFunction( null, yval, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( null, yval, zval );
Assert.fail( "Failed to detect x null pointer" ); Assert.fail( "Failed to detect x null pointer" );
} } catch ( NullArgumentException iae ) {
catch ( NullArgumentException iae )
{
// Expected. // Expected.
} }
try try {
{
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, null, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, null, zval );
Assert.fail( "Failed to detect y null pointer" ); Assert.fail( "Failed to detect y null pointer" );
} } catch ( NullArgumentException iae ) {
catch ( NullArgumentException iae )
{
// Expected. // Expected.
} }
try try {
{
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, null ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, null );
Assert.fail( "Failed to detect z null pointer" ); Assert.fail( "Failed to detect z null pointer" );
} } catch ( NullArgumentException iae ) {
catch ( NullArgumentException iae )
{
// Expected. // Expected.
} }
try try {
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0 }; double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect insufficient x data" ); Assert.fail( "Failed to detect insufficient x data" );
} } catch ( InsufficientDataException iae ) {
catch ( InsufficientDataException iae )
{
// Expected. // Expected.
} }
try try {
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0 }; double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect insufficient y data" ); Assert.fail( "Failed to detect insufficient y data" );
} } catch ( InsufficientDataException iae ) {
catch ( InsufficientDataException iae )
{
// Expected. // Expected.
} }
try try {
{
double zval1[][] = new double[4][4]; double zval1[][] = new double[4][4];
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, zval1 ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, zval1 );
Assert.fail( "Failed to detect insufficient z data" ); Assert.fail( "Failed to detect insufficient z data" );
} } catch ( InsufficientDataException iae ) {
catch ( InsufficientDataException iae )
{
// Expected. // Expected.
} }
try try {
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect data set array with different sizes." ); Assert.fail( "Failed to detect data set array with different sizes." );
} } catch ( DimensionMismatchException iae ) {
catch ( DimensionMismatchException iae )
{
// Expected. // Expected.
} }
try try {
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect data set array with different sizes." ); Assert.fail( "Failed to detect data set array with different sizes." );
} } catch ( DimensionMismatchException iae ) {
catch ( DimensionMismatchException iae )
{
// Expected. // Expected.
} }
// X values not sorted. // X values not sorted.
try try {
{
double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 }; double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect unsorted x arguments." ); Assert.fail( "Failed to detect unsorted x arguments." );
} } catch ( NonMonotonicSequenceException iae ) {
catch ( NonMonotonicSequenceException iae )
{
// Expected. // Expected.
} }
// Y values not sorted. // Y values not sorted.
try try {
{
double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 }; double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval ); bcf = new PiecewiseBicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect unsorted y arguments." ); Assert.fail( "Failed to detect unsorted y arguments." );
} } catch ( NonMonotonicSequenceException iae ) {
catch ( NonMonotonicSequenceException iae )
{
// Expected. // Expected.
} }
} }
@ -166,8 +134,7 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
* z = 2 x - 3 y + 5 * z = 2 x - 3 y + 5
*/ */
@Test @Test
public void testInterpolatePlane() public void testInterpolatePlane() {
{
final int numberOfElements = 10; final int numberOfElements = 10;
final double minimumX = -10; final double minimumX = -10;
final double maximumX = 10; final double maximumX = 10;
@ -178,10 +145,8 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
final double maxTolerance = 6e-14; final double maxTolerance = 6e-14;
// Function values // Function values
BivariateFunction f = new BivariateFunction() BivariateFunction f = new BivariateFunction() {
{ public double value( double x, double y ) {
public double value( double x, double y )
{
return 2 * x - 3 * y + 5; return 2 * x - 3 * y + 5;
} }
}; };
@ -196,8 +161,7 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5 * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/ */
@Test @Test
public void testInterpolationParabaloid() public void testInterpolationParabaloid() {
{
final int numberOfElements = 10; final int numberOfElements = 10;
final double minimumX = -10; final double minimumX = -10;
final double maximumX = 10; final double maximumX = 10;
@ -208,22 +172,19 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
final double maxTolerance = 6e-14; final double maxTolerance = 6e-14;
// Function values // Function values
BivariateFunction f = new BivariateFunction() BivariateFunction f = new BivariateFunction() {
{ public double value( double x, double y ) {
public double value( double x, double y )
{
return 2 * x * x - 3 * y * y + 4 * x * y - 5; return 2 * x * x - 3 * y * y + 4 * x * y - 5;
} }
}; };
testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f, testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f,
interpolationTolerance, maxTolerance ); interpolationTolerance, maxTolerance );
} }
private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY, private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY,
int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance, int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance,
double maxTolerance ) double maxTolerance ) {
{
double expected; double expected;
double actual; double actual;
double currentX; double currentX;
@ -234,11 +195,9 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
double yValues[] = new double[numberOfElements]; double yValues[] = new double[numberOfElements];
double zValues[][] = new double[numberOfElements][numberOfElements]; double zValues[][] = new double[numberOfElements][numberOfElements];
for ( int i = 0; i < numberOfElements; i++ ) for ( int i = 0; i < numberOfElements; i++ ) {
{
xValues[i] = minimumX + deltaX * (double) i; xValues[i] = minimumX + deltaX * (double) i;
for ( int j = 0; j < numberOfElements; j++ ) for ( int j = 0; j < numberOfElements; j++ ) {
{
yValues[j] = minimumY + deltaY * (double) j; yValues[j] = minimumY + deltaY * (double) j;
zValues[i][j] = f.value( xValues[i], yValues[j] ); zValues[i][j] = f.value( xValues[i], yValues[j] );
} }
@ -246,11 +205,9 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
BivariateFunction interpolation = new PiecewiseBicubicSplineInterpolatingFunction( xValues, yValues, zValues ); BivariateFunction interpolation = new PiecewiseBicubicSplineInterpolatingFunction( xValues, yValues, zValues );
for ( int i = 0; i < numberOfElements; i++ ) for ( int i = 0; i < numberOfElements; i++ ) {
{
currentX = xValues[i]; currentX = xValues[i];
for ( int j = 0; j < numberOfElements; j++ ) for ( int j = 0; j < numberOfElements; j++ ) {
{
currentY = yValues[j]; currentY = yValues[j];
expected = f.value( currentX, currentY ); expected = f.value( currentX, currentY );
actual = interpolation.value( currentX, currentY ); actual = interpolation.value( currentX, currentY );
@ -259,21 +216,18 @@ public final class PiecewiseBicubicSplineInterpolatingFunctionTest
} }
final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed. final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed.
final UniformRealDistribution distX = final UniformRealDistribution distX = new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] );
new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] ); final UniformRealDistribution distY = new UniformRealDistribution( rng, yValues[0], yValues[yValues.length - 1] );
final UniformRealDistribution distY =
new UniformRealDistribution( rng, yValues[0], yValues[yValues.length - 1] );
double sumError = 0; double sumError = 0;
for ( int i = 0; i < numberOfSamples; i++ ) for ( int i = 0; i < numberOfSamples; i++ ) {
{
currentX = distX.sample(); currentX = distX.sample();
currentY = distY.sample(); currentY = distY.sample();
expected = f.value( currentX, currentY ); expected = f.value( currentX, currentY );
actual = interpolation.value( currentX, currentY ); actual = interpolation.value( currentX, currentY );
sumError += FastMath.abs( actual - expected ); sumError += FastMath.abs( actual - expected );
assertEquals( expected, actual, maxTolerance ); assertEquals( expected, actual, maxTolerance );
} }
assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance ); assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance );
} }

View File

@ -31,14 +31,12 @@ import org.junit.Test;
/** /**
* Test case for the piecewise bicubic interpolator. * Test case for the piecewise bicubic interpolator.
*/ */
public final class PiecewiseBicubicSplineInterpolatorTest public final class PiecewiseBicubicSplineInterpolatorTest {
{
/** /**
* Test preconditions. * Test preconditions.
*/ */
@Test @Test
public void testPreconditions() public void testPreconditions() {
{
double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 }; double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 }; double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 };
double[][] zval = new double[xval.length][yval.length]; double[][] zval = new double[xval.length][yval.length];
@ -46,115 +44,84 @@ public final class PiecewiseBicubicSplineInterpolatorTest
@SuppressWarnings( "unused" ) @SuppressWarnings( "unused" )
BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator(); BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
try try {
{
interpolator.interpolate( null, yval, zval ); interpolator.interpolate( null, yval, zval );
Assert.fail( "Failed to detect x null pointer" ); Assert.fail( "Failed to detect x null pointer" );
} } catch ( NullArgumentException iae ) {
catch ( NullArgumentException iae )
{
// Expected. // Expected.
} }
try try {
{
interpolator.interpolate( xval, null, zval ); interpolator.interpolate( xval, null, zval );
Assert.fail( "Failed to detect y null pointer" ); Assert.fail( "Failed to detect y null pointer" );
} } catch ( NullArgumentException iae ) {
catch ( NullArgumentException iae )
{
// Expected. // Expected.
} }
try try {
{
interpolator.interpolate( xval, yval, null ); interpolator.interpolate( xval, yval, null );
Assert.fail( "Failed to detect z null pointer" ); Assert.fail( "Failed to detect z null pointer" );
} } catch ( NullArgumentException iae ) {
catch ( NullArgumentException iae )
{
// Expected. // Expected.
} }
try try {
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0 }; double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
interpolator.interpolate( xval1, yval, zval ); interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect insufficient x data" ); Assert.fail( "Failed to detect insufficient x data" );
} } catch ( InsufficientDataException iae ) {
catch ( InsufficientDataException iae )
{
// Expected. // Expected.
} }
try try {
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0 }; double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
interpolator.interpolate( xval, yval1, zval ); interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect insufficient y data" ); Assert.fail( "Failed to detect insufficient y data" );
} } catch ( InsufficientDataException iae ) {
catch ( InsufficientDataException iae )
{
// Expected. // Expected.
} }
try try {
{
double zval1[][] = new double[4][4]; double zval1[][] = new double[4][4];
interpolator.interpolate( xval, yval, zval1 ); interpolator.interpolate( xval, yval, zval1 );
Assert.fail( "Failed to detect insufficient z data" ); Assert.fail( "Failed to detect insufficient z data" );
} } catch ( InsufficientDataException iae ) {
catch ( InsufficientDataException iae )
{
// Expected. // Expected.
} }
try try {
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
interpolator.interpolate( xval1, yval, zval ); interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect data set array with different sizes." ); Assert.fail( "Failed to detect data set array with different sizes." );
} } catch ( DimensionMismatchException iae ) {
catch ( DimensionMismatchException iae )
{
// Expected. // Expected.
} }
try try {
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
interpolator.interpolate( xval, yval1, zval ); interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect data set array with different sizes." ); Assert.fail( "Failed to detect data set array with different sizes." );
} } catch ( DimensionMismatchException iae ) {
catch ( DimensionMismatchException iae )
{
// Expected. // Expected.
} }
// X values not sorted. // X values not sorted.
try try {
{
double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 }; double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
interpolator.interpolate( xval1, yval, zval ); interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect unsorted x arguments." ); Assert.fail( "Failed to detect unsorted x arguments." );
} } catch ( NonMonotonicSequenceException iae ) {
catch ( NonMonotonicSequenceException iae )
{
// Expected. // Expected.
} }
// Y values not sorted. // Y values not sorted.
try try {
{
double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 }; double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
interpolator.interpolate( xval, yval1, zval ); interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect unsorted y arguments." ); Assert.fail( "Failed to detect unsorted y arguments." );
} } catch ( NonMonotonicSequenceException iae ) {
catch ( NonMonotonicSequenceException iae )
{
// Expected. // Expected.
} }
} }
/** /**
@ -163,32 +130,26 @@ public final class PiecewiseBicubicSplineInterpolatorTest
* z = 2 x - 3 y + 5 * z = 2 x - 3 y + 5
*/ */
@Test @Test
public void testInterpolation1() public void testInterpolation1() {
{
final int sz = 21; final int sz = 21;
double[] xval = new double[sz]; double[] xval = new double[sz];
double[] yval = new double[sz]; double[] yval = new double[sz];
// Coordinate values // Coordinate values
final double delta = 1d / (sz - 1); 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; xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta; yval[i] = -20 + 30 * i * delta;
} }
// Function values // Function values
BivariateFunction f = new BivariateFunction() BivariateFunction f = new BivariateFunction() {
{ public double value( double x, double y ) {
public double value( double x, double y )
{
return 2 * x - 3 * y + 5; return 2 * x - 3 * y + 5;
} }
}; };
double[][] zval = new double[xval.length][yval.length]; double[][] zval = new double[xval.length][yval.length];
for ( int i = 0; i < xval.length; i++ ) for ( int i = 0; i < xval.length; i++ ) {
{ for ( int j = 0; j < yval.length; j++ ) {
for ( int j = 0; j < yval.length; j++ )
{
zval[i][j] = f.value(xval[i], yval[j]); zval[i][j] = f.value(xval[i], yval[j]);
} }
} }
@ -203,11 +164,9 @@ public final class PiecewiseBicubicSplineInterpolatorTest
final int numSamples = 50; final int numSamples = 50;
final double tol = 2e-14; final double tol = 2e-14;
for ( int i = 0; i < numSamples; i++ ) for ( int i = 0; i < numSamples; i++ ) {
{
x = distX.sample(); x = distX.sample();
for ( int j = 0; j < numSamples; j++ ) for ( int j = 0; j < numSamples; j++ ) {
{
y = distY.sample(); y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y)); // System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
Assert.assertEquals(f.value(x, y), p.value(x, y), tol); Assert.assertEquals(f.value(x, y), p.value(x, y), tol);
@ -222,32 +181,26 @@ public final class PiecewiseBicubicSplineInterpolatorTest
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5 * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/ */
@Test @Test
public void testInterpolation2() public void testInterpolation2() {
{
final int sz = 21; final int sz = 21;
double[] xval = new double[sz]; double[] xval = new double[sz];
double[] yval = new double[sz]; double[] yval = new double[sz];
// Coordinate values // Coordinate values
final double delta = 1d / (sz - 1); 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; xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta; yval[i] = -20 + 30 * i * delta;
} }
// Function values // Function values
BivariateFunction f = new BivariateFunction() BivariateFunction f = new BivariateFunction() {
{ public double value( double x, double y ) {
public double value( double x, double y )
{
return 2 * x * x - 3 * y * y + 4 * x * y - 5; return 2 * x * x - 3 * y * y + 4 * x * y - 5;
} }
}; };
double[][] zval = new double[xval.length][yval.length]; double[][] zval = new double[xval.length][yval.length];
for ( int i = 0; i < xval.length; i++ ) for ( int i = 0; i < xval.length; i++ ) {
{ for ( int j = 0; j < yval.length; j++ ) {
for ( int j = 0; j < yval.length; j++ )
{
zval[i][j] = f.value(xval[i], yval[j]); zval[i][j] = f.value(xval[i], yval[j]);
} }
} }
@ -262,11 +215,9 @@ public final class PiecewiseBicubicSplineInterpolatorTest
final int numSamples = 50; final int numSamples = 50;
final double tol = 5e-13; final double tol = 5e-13;
for ( int i = 0; i < numSamples; i++ ) for ( int i = 0; i < numSamples; i++ ) {
{
x = distX.sample(); x = distX.sample();
for ( int j = 0; j < numSamples; j++ ) for ( int j = 0; j < numSamples; j++ ) {
{
y = distY.sample(); y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y)); // System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
Assert.assertEquals(f.value(x, y), p.value(x, y), tol); Assert.assertEquals(f.value(x, y), p.value(x, y), tol);