MATH-1138 #comment Implemented new BiCubicSplineInterpolator, supporting Akima Spline Interpolator and updated tests
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/d8bfc8c8 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/d8bfc8c8 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/d8bfc8c8 Branch: refs/heads/master Commit: d8bfc8c8f8864f9c22e0409780d5dd3fb30497ff Parents: a3fdeb4 Author: Hank Grabowski <h...@applieddefense.com> Authored: Wed Oct 15 10:15:14 2014 -0400 Committer: Hank Grabowski <h...@applieddefense.com> Committed: Wed Oct 15 10:15:14 2014 -0400 ---------------------------------------------------------------------- .../interpolation/AkimaSplineInterpolator.java | 225 ++++++ .../BicubicSplineInterpolatingFunction.java | 620 +++------------ .../BicubicSplineInterpolator.java | 140 +--- .../TricubicSplineInterpolator.java | 35 +- .../AkimaSplineInterpolatorTest.java | 227 ++++++ .../BicubicSplineInterpolatingFunctionTest.java | 746 +++++-------------- .../BicubicSplineInterpolatorTest.java | 217 ++++-- ...PolynomialBicubicSplineInterpolatorTest.java | 10 +- .../interpolation/SplineInterpolatorTest.java | 42 +- .../TricubicSplineInterpolatorTest.java | 2 +- 10 files changed, 923 insertions(+), 1341 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/d8bfc8c8/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java new file mode 100644 index 0000000..8fb896c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java @@ -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); + + } +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/d8bfc8c8/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java index 079e9fc..c0ce3c5 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java @@ -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); + + double xArray[] = new double[count]; + double yArray[] = new double[count]; + double zArray[] = new double[count]; + double interpArray[] = new double[count]; + + 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); - final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); - final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); + double returnValue = spline.value(y); - return splines[i][j].value(xN, yN); + 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 { @@ -249,385 +154,42 @@ 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}. - */ - private int searchIndex(double c, double[] val) { - final int r = Arrays.binarySearch(val, c); - - if (r == -1 || - r == -val.length - 1) { + * @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, int offset, int count) { + int r = Arrays.binarySearch(val, c); + + 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; - } - final int last = val.length - 1; - if (r == last) { - // "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; - } - - // "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); - } - }; + // "c" in within an interpolation sub-interval, which returns + // negative + // need to remove the negative sign for consistency + r = -r - offset - 1; } else { - partialDerivativeX = null; - partialDerivativeY = null; - partialDerivativeXX = null; - partialDerivativeYY = null; - partialDerivativeXY = null; + r -= offset; } - } - /** - * {@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); + if (r < 0) { + r = 0; } - 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]; - } + 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. + r = val.length - count; } - 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; + return r; } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/d8bfc8c8/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java index 5e2c92f..36a9da2 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java @@ -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 ); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/d8bfc8c8/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java index cda6a33..0faa274 100644 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java @@ -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]; http://git-wip-us.apache.org/repos/asf/commons-math/blob/d8bfc8c8/src/test/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolatorTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolatorTest.java new file mode 100644 index 0000000..017241a --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolatorTest.java @@ -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 ); + } +}
