diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/interpolation')
30 files changed, 5890 insertions, 0 deletions
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..5b89dfe --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java @@ -0,0 +1,215 @@ +/* + * 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> + * <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 + implements UnivariateInterpolator { + /** The minimum number of points that are needed to compute the function. */ + private static final int MINIMUM_NUMBER_POINTS = 5; + + /** + * 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 xvals} and {@code yvals} have + * different sizes. + * @throws NonMonotonicSequenceException if {@code xvals} is not sorted in + * strict increasing order. + * @throws NumberIsTooSmallException if the size of {@code xvals} 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; + + final double[] differences = new double[numberOfDiffAndWeightElements]; + final 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. + final double[] firstDerivatives = new double[xvals.length]; + + for (int i = 2; i < firstDerivatives.length - 2; i++) { + final double wP = weights[i + 1]; + final double wM = weights[i - 1]; + if (Precision.equals(wP, 0.0) && + 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 { + firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM); + } + } + + firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2); + firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2); + firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2, + xvals.length - 3, xvals.length - 2, + xvals.length - 1); + firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1, + xvals.length - 3, xvals.length - 2, + xvals.length - 1); + + return 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) { + final double x0 = yvals[indexOfFirstSample]; + final double x1 = yvals[indexOfSecondsample]; + final double x2 = yvals[indexOfThirdSample]; + + final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample]; + final double t1 = xvals[indexOfSecondsample] - 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; + + 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); + } + + final 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++) { + final double w = xvals[i + 1] - xvals[i]; + final double w2 = w * w; + + 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[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w; + coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2; + polynomials[i] = new PolynomialFunction(coefficients); + } + + return new PolynomialSplineFunction(xvals, polynomials); + + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java new file mode 100644 index 0000000..454d8f8 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java @@ -0,0 +1,325 @@ +/* + * 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 java.util.Arrays; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +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>. + * + * @since 3.4 + */ +public class BicubicInterpolatingFunction + 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 BicubicFunction[][] splines; + + /** + * @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. + * @throws NoDataException if any of the arrays has zero length. + */ + public BicubicInterpolatingFunction(double[] x, + double[] y, + double[][] f, + double[][] dFdX, + double[][] dFdY, + double[][] d2FdXdY) + throws DimensionMismatchException, + NoDataException, + NonMonotonicSequenceException { + 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 != 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); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + + xval = x.clone(); + yval = y.clone(); + + final int lastI = xLen - 1; + final int lastJ = yLen - 1; + splines = new BicubicFunction[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; + final double xR = xval[ip1] - xval[i]; + for (int j = 0; j < lastJ; j++) { + final int jp1 = j + 1; + final double yR = yval[jp1] - yval[j]; + final double xRyR = xR * yR; + final double[] beta = new double[] { + f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1], + dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR, + dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR, + d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR + }; + + splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta)); + } + } + } + + /** + * {@inheritDoc} + */ + public double value(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 splines[i][j].value(xN, yN); + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @return {@code true} if (x, y) is a valid point. + */ + public boolean isValidPoint(double x, double y) { + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || + y > yval[yval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @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) { + 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; + } +} + +/** + * Bicubic function. + */ +class BicubicFunction implements BivariateFunction { + /** Number of points. */ + private static final short N = 4; + /** Coefficients */ + private final double[][] a; + + /** + * Simple constructor. + * + * @param coeff Spline coefficients. + */ + BicubicFunction(double[] coeff) { + a = new double[N][N]; + for (int j = 0; j < N; j++) { + final double[] aJ = a[j]; + for (int i = 0; i < N; i++) { + aJ[i] = coeff[i * N + j]; + } + } + } + + /** + * {@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++) { + final double r = MathArrays.linearCombination(coeff[i], pY); + result += r * pX[i]; + } + + return result; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java new file mode 100644 index 0000000..4aebdff --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java @@ -0,0 +1,113 @@ +/* + * 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.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Generates a {@link BicubicInterpolatingFunction bicubic interpolating + * function}. + * <p> + * Caveat: Because the interpolation scheme requires that derivatives be + * specified at the sample points, those are approximated with finite + * differences (using the 2-points symmetric formulae). + * Since their values are undefined at the borders of the provided + * interpolation ranges, the interpolated values will be wrong at the + * edges of the patch. + * The {@code interpolate} method will return a function that overrides + * {@link BicubicInterpolatingFunction#isValidPoint(double,double)} to + * indicate points where the interpolation will be inaccurate. + * </p> + * + * @since 3.4 + */ +public class BicubicInterpolator + implements BivariateGridInterpolator { + /** + * {@inheritDoc} + */ + public BicubicInterpolatingFunction 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(); + } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + + final int xLen = xval.length; + final int yLen = yval.length; + + // Approximation to the partial derivatives using finite differences. + final double[][] dFdX = new double[xLen][yLen]; + final double[][] dFdY = new double[xLen][yLen]; + final double[][] d2FdXdY = new double[xLen][yLen]; + for (int i = 1; i < xLen - 1; i++) { + final int nI = i + 1; + final int pI = i - 1; + + final double nX = xval[nI]; + final double pX = xval[pI]; + + final double deltaX = nX - pX; + + for (int j = 1; j < yLen - 1; j++) { + final int nJ = j + 1; + final int pJ = j - 1; + + final double nY = yval[nJ]; + final double pY = yval[pJ]; + + final double deltaY = nY - pY; + + dFdX[i][j] = (fval[nI][j] - fval[pI][j]) / deltaX; + dFdY[i][j] = (fval[i][nJ] - fval[i][pJ]) / deltaY; + + final double deltaXY = deltaX * deltaY; + + d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / deltaXY; + } + } + + // Create the interpolating function. + return new BicubicInterpolatingFunction(xval, yval, fval, + dFdX, dFdY, d2FdXdY) { + /** {@inheritDoc} */ + @Override + public boolean isValidPoint(double x, double y) { + if (x < xval[1] || + x > xval[xval.length - 2] || + y < yval[1] || + y > yval[yval.length - 2]) { + return false; + } else { + return true; + } + } + }; + } +} 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 new file mode 100644 index 0000000..8f1771a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java @@ -0,0 +1,641 @@ +/* + * 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 java.util.Arrays; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +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>. Due to numerical accuracy issues this should not + * be used. + * + * @since 2.1 + * @deprecated as of 3.4 replaced by + * {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction} + */ +@Deprecated +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; + + /** + * @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. + * @throws NoDataException if any of the arrays has zero length. + */ + 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 { + 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 != 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); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + + 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; + } + } + + /** + * {@inheritDoc} + */ + public double value(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 splines[i][j].value(xN, yN); + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @return {@code true} if (x, y) is a valid point. + * @since 3.3 + */ + public boolean isValidPoint(double x, double y) { + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || + y > yval[yval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @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) { + 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. + */ + 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. + */ + 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() { + /** {@inheritDoc} */ + 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() { + /** {@inheritDoc} */ + 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() { + /** {@inheritDoc} */ + 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() { + /** {@inheritDoc} */ + 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() { + /** {@inheritDoc} */ + 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; + } +} 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 new file mode 100644 index 0000000..09acd07 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java @@ -0,0 +1,176 @@ +/* + * 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.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.NumberIsTooSmallException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Generates a bicubic interpolating function. Due to numerical accuracy issues this should not + * be used. + * + * @since 2.2 + * @deprecated as of 3.4 replaced by {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolator} + */ +@Deprecated +public class BicubicSplineInterpolator + implements BivariateGridInterpolator { + /** Whether to initialize internal data used to compute the analytical + derivatives of the splines. */ + private final boolean initializeDerivatives; + + /** + * Default constructor. + * The argument {@link #BicubicSplineInterpolator(boolean) initializeDerivatives} + * is set to {@code false}. + */ + public BicubicSplineInterpolator() { + this(false); + } + + /** + * Creates an interpolator. + * + * @param initializeDerivatives Whether to initialize the internal data + * needed for calling any of the methods that compute the partial derivatives + * of the {@link 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, + final double[][] fval) + throws NoDataException, DimensionMismatchException, + NonMonotonicSequenceException, NumberIsTooSmallException { + if (xval.length == 0 || yval.length == 0 || fval.length == 0) { + throw new NoDataException(); + } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } + + 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; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BivariateGridInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BivariateGridInterpolator.java new file mode 100644 index 0000000..94d75ad --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/BivariateGridInterpolator.java @@ -0,0 +1,51 @@ +/* + * 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.BivariateFunction; +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.NumberIsTooSmallException; + +/** + * Interface representing a bivariate real interpolating function where the + * sample points must be specified on a regular grid. + * + */ +public interface BivariateGridInterpolator { + /** + * Compute an interpolating function for the dataset. + * + * @param xval All the x-coordinates of the interpolation points, sorted + * in increasing order. + * @param yval All the y-coordinates of the interpolation points, sorted + * in increasing order. + * @param fval The values of the interpolation points on all the grid knots: + * {@code fval[i][j] = f(xval[i], yval[j])}. + * @return a function which interpolates the dataset. + * @throws NoDataException if any of the arrays has zero length. + * @throws DimensionMismatchException if the array lengths are inconsistent. + * @throws NonMonotonicSequenceException if the array is not sorted. + * @throws NumberIsTooSmallException if the number of points is too small for + * the order of the interpolation + */ + BivariateFunction interpolate(double[] xval, double[] yval, + double[][] fval) + throws NoDataException, DimensionMismatchException, + NonMonotonicSequenceException, NumberIsTooSmallException; +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/DividedDifferenceInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/DividedDifferenceInterpolator.java new file mode 100644 index 0000000..86b988c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/DividedDifferenceInterpolator.java @@ -0,0 +1,120 @@ +/* + * 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 java.io.Serializable; +import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm; +import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionNewtonForm; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; + +/** + * Implements the <a href= + * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html"> + * Divided Difference Algorithm</a> for interpolation of real univariate + * functions. For reference, see <b>Introduction to Numerical Analysis</b>, + * ISBN 038795452X, chapter 2. + * <p> + * The actual code of Neville's evaluation is in PolynomialFunctionLagrangeForm, + * this class provides an easy-to-use interface to it.</p> + * + * @since 1.2 + */ +public class DividedDifferenceInterpolator + implements UnivariateInterpolator, Serializable { + /** serializable version identifier */ + private static final long serialVersionUID = 107049519551235069L; + + /** + * Compute an interpolating function for the dataset. + * + * @param x Interpolating points array. + * @param y Interpolating values array. + * @return a function which interpolates the dataset. + * @throws DimensionMismatchException if the array lengths are different. + * @throws NumberIsTooSmallException if the number of points is less than 2. + * @throws NonMonotonicSequenceException if {@code x} is not sorted in + * strictly increasing order. + */ + public PolynomialFunctionNewtonForm interpolate(double x[], double y[]) + throws DimensionMismatchException, + NumberIsTooSmallException, + NonMonotonicSequenceException { + /** + * a[] and c[] are defined in the general formula of Newton form: + * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + + * a[n](x-c[0])(x-c[1])...(x-c[n-1]) + */ + PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true); + + /** + * When used for interpolation, the Newton form formula becomes + * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... + + * f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2]) + * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k]. + * <p> + * Note x[], y[], a[] have the same length but c[]'s size is one less.</p> + */ + final double[] c = new double[x.length-1]; + System.arraycopy(x, 0, c, 0, c.length); + + final double[] a = computeDividedDifference(x, y); + return new PolynomialFunctionNewtonForm(a, c); + } + + /** + * Return a copy of the divided difference array. + * <p> + * The divided difference array is defined recursively by <pre> + * f[x0] = f(x0) + * f[x0,x1,...,xk] = (f[x1,...,xk] - f[x0,...,x[k-1]]) / (xk - x0) + * </pre> + * <p> + * The computational complexity is \(O(n^2)\) where \(n\) is the common + * length of {@code x} and {@code y}.</p> + * + * @param x Interpolating points array. + * @param y Interpolating values array. + * @return a fresh copy of the divided difference array. + * @throws DimensionMismatchException if the array lengths are different. + * @throws NumberIsTooSmallException if the number of points is less than 2. + * @throws NonMonotonicSequenceException + * if {@code x} is not sorted in strictly increasing order. + */ + protected static double[] computeDividedDifference(final double x[], final double y[]) + throws DimensionMismatchException, + NumberIsTooSmallException, + NonMonotonicSequenceException { + PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y, true); + + final double[] divdiff = y.clone(); // initialization + + final int n = x.length; + final double[] a = new double [n]; + a[0] = divdiff[0]; + for (int i = 1; i < n; i++) { + for (int j = 0; j < n-i; j++) { + final double denominator = x[j+i] - x[j]; + divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator; + } + a[i] = divdiff[0]; + } + + return a; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java new file mode 100644 index 0000000..9125b83 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java @@ -0,0 +1,209 @@ +/* + * 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 java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.FieldElement; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathArithmeticException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.MathUtils; + +/** Polynomial interpolator using both sample values and sample derivatives. + * <p> + * The interpolation polynomials match all sample points, including both values + * and provided derivatives. There is one polynomial for each component of + * the values vector. All polynomials have the same degree. The degree of the + * polynomials depends on the number of points and number of derivatives at each + * point. For example the interpolation polynomials for n sample points without + * any derivatives all have degree n-1. The interpolation polynomials for n + * sample points with the two extreme points having value and first derivative + * and the remaining points having value only all have degree n+1. The + * interpolation polynomial for n sample points with value, first and second + * derivative for all points all have degree 3n-1. + * </p> + * + * @param <T> Type of the field elements. + * + * @since 3.2 + */ +public class FieldHermiteInterpolator<T extends FieldElement<T>> { + + /** Sample abscissae. */ + private final List<T> abscissae; + + /** Top diagonal of the divided differences array. */ + private final List<T[]> topDiagonal; + + /** Bottom diagonal of the divided differences array. */ + private final List<T[]> bottomDiagonal; + + /** Create an empty interpolator. + */ + public FieldHermiteInterpolator() { + this.abscissae = new ArrayList<T>(); + this.topDiagonal = new ArrayList<T[]>(); + this.bottomDiagonal = new ArrayList<T[]>(); + } + + /** Add a sample point. + * <p> + * This method must be called once for each sample point. It is allowed to + * mix some calls with values only with calls with values and first + * derivatives. + * </p> + * <p> + * The point abscissae for all calls <em>must</em> be different. + * </p> + * @param x abscissa of the sample point + * @param value value and derivatives of the sample point + * (if only one row is passed, it is the value, if two rows are + * passed the first one is the value and the second the derivative + * and so on) + * @exception ZeroException if the abscissa difference between added point + * and a previous point is zero (i.e. the two points are at same abscissa) + * @exception MathArithmeticException if the number of derivatives is larger + * than 20, which prevents computation of a factorial + * @throws DimensionMismatchException if derivative structures are inconsistent + * @throws NullArgumentException if x is null + */ + public void addSamplePoint(final T x, final T[] ... value) + throws ZeroException, MathArithmeticException, + DimensionMismatchException, NullArgumentException { + + MathUtils.checkNotNull(x); + T factorial = x.getField().getOne(); + for (int i = 0; i < value.length; ++i) { + + final T[] y = value[i].clone(); + if (i > 1) { + factorial = factorial.multiply(i); + final T inv = factorial.reciprocal(); + for (int j = 0; j < y.length; ++j) { + y[j] = y[j].multiply(inv); + } + } + + // update the bottom diagonal of the divided differences array + final int n = abscissae.size(); + bottomDiagonal.add(n - i, y); + T[] bottom0 = y; + for (int j = i; j < n; ++j) { + final T[] bottom1 = bottomDiagonal.get(n - (j + 1)); + if (x.equals(abscissae.get(n - (j + 1)))) { + throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x); + } + final T inv = x.subtract(abscissae.get(n - (j + 1))).reciprocal(); + for (int k = 0; k < y.length; ++k) { + bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k])); + } + bottom0 = bottom1; + } + + // update the top diagonal of the divided differences array + topDiagonal.add(bottom0.clone()); + + // update the abscissae array + abscissae.add(x); + + } + + } + + /** Interpolate value at a specified abscissa. + * @param x interpolation abscissa + * @return interpolated value + * @exception NoDataException if sample is empty + * @throws NullArgumentException if x is null + */ + public T[] value(T x) throws NoDataException, NullArgumentException { + + // safety check + MathUtils.checkNotNull(x); + if (abscissae.isEmpty()) { + throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE); + } + + final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length); + T valueCoeff = x.getField().getOne(); + for (int i = 0; i < topDiagonal.size(); ++i) { + T[] dividedDifference = topDiagonal.get(i); + for (int k = 0; k < value.length; ++k) { + value[k] = value[k].add(dividedDifference[k].multiply(valueCoeff)); + } + final T deltaX = x.subtract(abscissae.get(i)); + valueCoeff = valueCoeff.multiply(deltaX); + } + + return value; + + } + + /** Interpolate value and first derivatives at a specified abscissa. + * @param x interpolation abscissa + * @param order maximum derivation order + * @return interpolated value and derivatives (value in row 0, + * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n) + * @exception NoDataException if sample is empty + * @throws NullArgumentException if x is null + */ + public T[][] derivatives(T x, int order) throws NoDataException, NullArgumentException { + + // safety check + MathUtils.checkNotNull(x); + if (abscissae.isEmpty()) { + throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE); + } + + final T zero = x.getField().getZero(); + final T one = x.getField().getOne(); + final T[] tj = MathArrays.buildArray(x.getField(), order + 1); + tj[0] = zero; + for (int i = 0; i < order; ++i) { + tj[i + 1] = tj[i].add(one); + } + + final T[][] derivatives = + MathArrays.buildArray(x.getField(), order + 1, topDiagonal.get(0).length); + final T[] valueCoeff = MathArrays.buildArray(x.getField(), order + 1); + valueCoeff[0] = x.getField().getOne(); + for (int i = 0; i < topDiagonal.size(); ++i) { + T[] dividedDifference = topDiagonal.get(i); + final T deltaX = x.subtract(abscissae.get(i)); + for (int j = order; j >= 0; --j) { + for (int k = 0; k < derivatives[j].length; ++k) { + derivatives[j][k] = + derivatives[j][k].add(dividedDifference[k].multiply(valueCoeff[j])); + } + valueCoeff[j] = valueCoeff[j].multiply(deltaX); + if (j > 0) { + valueCoeff[j] = valueCoeff[j].add(tj[j].multiply(valueCoeff[j - 1])); + } + } + } + + return derivatives; + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java new file mode 100644 index 0000000..15ed322 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java @@ -0,0 +1,239 @@ +/* + * 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 java.util.ArrayList; +import java.util.Arrays; +import java.util.List; + +import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableVectorFunction; +import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math3.exception.MathArithmeticException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.CombinatoricsUtils; + +/** Polynomial interpolator using both sample values and sample derivatives. + * <p> + * The interpolation polynomials match all sample points, including both values + * and provided derivatives. There is one polynomial for each component of + * the values vector. All polynomials have the same degree. The degree of the + * polynomials depends on the number of points and number of derivatives at each + * point. For example the interpolation polynomials for n sample points without + * any derivatives all have degree n-1. The interpolation polynomials for n + * sample points with the two extreme points having value and first derivative + * and the remaining points having value only all have degree n+1. The + * interpolation polynomial for n sample points with value, first and second + * derivative for all points all have degree 3n-1. + * </p> + * + * @since 3.1 + */ +public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction { + + /** Sample abscissae. */ + private final List<Double> abscissae; + + /** Top diagonal of the divided differences array. */ + private final List<double[]> topDiagonal; + + /** Bottom diagonal of the divided differences array. */ + private final List<double[]> bottomDiagonal; + + /** Create an empty interpolator. + */ + public HermiteInterpolator() { + this.abscissae = new ArrayList<Double>(); + this.topDiagonal = new ArrayList<double[]>(); + this.bottomDiagonal = new ArrayList<double[]>(); + } + + /** Add a sample point. + * <p> + * This method must be called once for each sample point. It is allowed to + * mix some calls with values only with calls with values and first + * derivatives. + * </p> + * <p> + * The point abscissae for all calls <em>must</em> be different. + * </p> + * @param x abscissa of the sample point + * @param value value and derivatives of the sample point + * (if only one row is passed, it is the value, if two rows are + * passed the first one is the value and the second the derivative + * and so on) + * @exception ZeroException if the abscissa difference between added point + * and a previous point is zero (i.e. the two points are at same abscissa) + * @exception MathArithmeticException if the number of derivatives is larger + * than 20, which prevents computation of a factorial + */ + public void addSamplePoint(final double x, final double[] ... value) + throws ZeroException, MathArithmeticException { + + for (int i = 0; i < value.length; ++i) { + + final double[] y = value[i].clone(); + if (i > 1) { + double inv = 1.0 / CombinatoricsUtils.factorial(i); + for (int j = 0; j < y.length; ++j) { + y[j] *= inv; + } + } + + // update the bottom diagonal of the divided differences array + final int n = abscissae.size(); + bottomDiagonal.add(n - i, y); + double[] bottom0 = y; + for (int j = i; j < n; ++j) { + final double[] bottom1 = bottomDiagonal.get(n - (j + 1)); + final double inv = 1.0 / (x - abscissae.get(n - (j + 1))); + if (Double.isInfinite(inv)) { + throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x); + } + for (int k = 0; k < y.length; ++k) { + bottom1[k] = inv * (bottom0[k] - bottom1[k]); + } + bottom0 = bottom1; + } + + // update the top diagonal of the divided differences array + topDiagonal.add(bottom0.clone()); + + // update the abscissae array + abscissae.add(x); + + } + + } + + /** Compute the interpolation polynomials. + * @return interpolation polynomials array + * @exception NoDataException if sample is empty + */ + public PolynomialFunction[] getPolynomials() + throws NoDataException { + + // safety check + checkInterpolation(); + + // iteration initialization + final PolynomialFunction zero = polynomial(0); + PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length]; + for (int i = 0; i < polynomials.length; ++i) { + polynomials[i] = zero; + } + PolynomialFunction coeff = polynomial(1); + + // build the polynomials by iterating on the top diagonal of the divided differences array + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] tdi = topDiagonal.get(i); + for (int k = 0; k < polynomials.length; ++k) { + polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k]))); + } + coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0)); + } + + return polynomials; + + } + + /** Interpolate value at a specified abscissa. + * <p> + * Calling this method is equivalent to call the {@link PolynomialFunction#value(double) + * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials}, + * except it does not build the intermediate polynomials, so this method is faster and + * numerically more stable. + * </p> + * @param x interpolation abscissa + * @return interpolated value + * @exception NoDataException if sample is empty + */ + public double[] value(double x) + throws NoDataException { + + // safety check + checkInterpolation(); + + final double[] value = new double[topDiagonal.get(0).length]; + double valueCoeff = 1; + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] dividedDifference = topDiagonal.get(i); + for (int k = 0; k < value.length; ++k) { + value[k] += dividedDifference[k] * valueCoeff; + } + final double deltaX = x - abscissae.get(i); + valueCoeff *= deltaX; + } + + return value; + + } + + /** Interpolate value at a specified abscissa. + * <p> + * Calling this method is equivalent to call the {@link + * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials + * returned by {@link #getPolynomials() getPolynomials}, except it does not build the + * intermediate polynomials, so this method is faster and numerically more stable. + * </p> + * @param x interpolation abscissa + * @return interpolated value + * @exception NoDataException if sample is empty + */ + public DerivativeStructure[] value(final DerivativeStructure x) + throws NoDataException { + + // safety check + checkInterpolation(); + + final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length]; + Arrays.fill(value, x.getField().getZero()); + DerivativeStructure valueCoeff = x.getField().getOne(); + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] dividedDifference = topDiagonal.get(i); + for (int k = 0; k < value.length; ++k) { + value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k])); + } + final DerivativeStructure deltaX = x.subtract(abscissae.get(i)); + valueCoeff = valueCoeff.multiply(deltaX); + } + + return value; + + } + + /** Check interpolation can be performed. + * @exception NoDataException if interpolation cannot be performed + * because sample is empty + */ + private void checkInterpolation() throws NoDataException { + if (abscissae.isEmpty()) { + throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE); + } + } + + /** Create a polynomial from its coefficients. + * @param c polynomials coefficients + * @return polynomial + */ + private PolynomialFunction polynomial(double ... c) { + return new PolynomialFunction(c); + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere.java new file mode 100644 index 0000000..dc600bd --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere.java @@ -0,0 +1,385 @@ +/* + * 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 java.util.List; +import java.util.ArrayList; +import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; + +/** + * Utility class for the {@link MicrosphereProjectionInterpolator} algorithm. + * + * @since 3.6 + */ +public class InterpolatingMicrosphere { + /** Microsphere. */ + private final List<Facet> microsphere; + /** Microsphere data. */ + private final List<FacetData> microsphereData; + /** Space dimension. */ + private final int dimension; + /** Number of surface elements. */ + private final int size; + /** Maximum fraction of the facets that can be dark. */ + private final double maxDarkFraction; + /** Lowest non-zero illumination. */ + private final double darkThreshold; + /** Background value. */ + private final double background; + + /** + * Create an unitialiazed sphere. + * Sub-classes are responsible for calling the {@code add(double[]) add} + * method in order to initialize all the sphere's facets. + * + * @param dimension Dimension of the data space. + * @param size Number of surface elements of the sphere. + * @param maxDarkFraction Maximum fraction of the facets that can be dark. + * If the fraction of "non-illuminated" facets is larger, no estimation + * of the value will be performed, and the {@code background} value will + * be returned instead. + * @param darkThreshold Value of the illumination below which a facet is + * considered dark. + * @param background Value returned when the {@code maxDarkFraction} + * threshold is exceeded. + * @throws NotStrictlyPositiveException if {@code dimension <= 0} + * or {@code size <= 0}. + * @throws NotPositiveException if {@code darkThreshold < 0}. + * @throws OutOfRangeException if {@code maxDarkFraction} does not + * belong to the interval {@code [0, 1]}. + */ + protected InterpolatingMicrosphere(int dimension, + int size, + double maxDarkFraction, + double darkThreshold, + double background) { + if (dimension <= 0) { + throw new NotStrictlyPositiveException(dimension); + } + if (size <= 0) { + throw new NotStrictlyPositiveException(size); + } + if (maxDarkFraction < 0 || + maxDarkFraction > 1) { + throw new OutOfRangeException(maxDarkFraction, 0, 1); + } + if (darkThreshold < 0) { + throw new NotPositiveException(darkThreshold); + } + + this.dimension = dimension; + this.size = size; + this.maxDarkFraction = maxDarkFraction; + this.darkThreshold = darkThreshold; + this.background = background; + microsphere = new ArrayList<Facet>(size); + microsphereData = new ArrayList<FacetData>(size); + } + + /** + * Create a sphere from randomly sampled vectors. + * + * @param dimension Dimension of the data space. + * @param size Number of surface elements of the sphere. + * @param rand Unit vector generator for creating the microsphere. + * @param maxDarkFraction Maximum fraction of the facets that can be dark. + * If the fraction of "non-illuminated" facets is larger, no estimation + * of the value will be performed, and the {@code background} value will + * be returned instead. + * @param darkThreshold Value of the illumination below which a facet + * is considered dark. + * @param background Value returned when the {@code maxDarkFraction} + * threshold is exceeded. + * @throws DimensionMismatchException if the size of the generated + * vectors does not match the dimension set in the constructor. + * @throws NotStrictlyPositiveException if {@code dimension <= 0} + * or {@code size <= 0}. + * @throws NotPositiveException if {@code darkThreshold < 0}. + * @throws OutOfRangeException if {@code maxDarkFraction} does not + * belong to the interval {@code [0, 1]}. + */ + public InterpolatingMicrosphere(int dimension, + int size, + double maxDarkFraction, + double darkThreshold, + double background, + UnitSphereRandomVectorGenerator rand) { + this(dimension, size, maxDarkFraction, darkThreshold, background); + + // Generate the microsphere normals, assuming that a number of + // randomly generated normals will represent a sphere. + for (int i = 0; i < size; i++) { + add(rand.nextVector(), false); + } + } + + /** + * Copy constructor. + * + * @param other Instance to copy. + */ + protected InterpolatingMicrosphere(InterpolatingMicrosphere other) { + dimension = other.dimension; + size = other.size; + maxDarkFraction = other.maxDarkFraction; + darkThreshold = other.darkThreshold; + background = other.background; + + // Field can be shared. + microsphere = other.microsphere; + + // Field must be copied. + microsphereData = new ArrayList<FacetData>(size); + for (FacetData fd : other.microsphereData) { + microsphereData.add(new FacetData(fd.illumination(), fd.sample())); + } + } + + /** + * Perform a copy. + * + * @return a copy of this instance. + */ + public InterpolatingMicrosphere copy() { + return new InterpolatingMicrosphere(this); + } + + /** + * Get the space dimensionality. + * + * @return the number of space dimensions. + */ + public int getDimension() { + return dimension; + } + + /** + * Get the size of the sphere. + * + * @return the number of surface elements of the microspshere. + */ + public int getSize() { + return size; + } + + /** + * Estimate the value at the requested location. + * This microsphere is placed at the given {@code point}, contribution + * of the given {@code samplePoints} to each sphere facet is computed + * (illumination) and the interpolation is performed (integration of + * the illumination). + * + * @param point Interpolation point. + * @param samplePoints Sampling data points. + * @param sampleValues Sampling data values at the corresponding + * {@code samplePoints}. + * @param exponent Exponent used in the power law that computes + * the weights (distance dimming factor) of the sample data. + * @param noInterpolationTolerance When the distance between the + * {@code point} and one of the {@code samplePoints} is less than + * this value, no interpolation will be performed, and the value + * of the sample will just be returned. + * @return the estimated value at the given {@code point}. + * @throws NotPositiveException if {@code exponent < 0}. + */ + public double value(double[] point, + double[][] samplePoints, + double[] sampleValues, + double exponent, + double noInterpolationTolerance) { + if (exponent < 0) { + throw new NotPositiveException(exponent); + } + + clear(); + + // Contribution of each sample point to the illumination of the + // microsphere's facets. + final int numSamples = samplePoints.length; + for (int i = 0; i < numSamples; i++) { + // Vector between interpolation point and current sample point. + final double[] diff = MathArrays.ebeSubtract(samplePoints[i], point); + final double diffNorm = MathArrays.safeNorm(diff); + + if (FastMath.abs(diffNorm) < noInterpolationTolerance) { + // No need to interpolate, as the interpolation point is + // actually (very close to) one of the sampled points. + return sampleValues[i]; + } + + final double weight = FastMath.pow(diffNorm, -exponent); + illuminate(diff, sampleValues[i], weight); + } + + return interpolate(); + } + + /** + * Replace {@code i}-th facet of the microsphere. + * Method for initializing the microsphere facets. + * + * @param normal Facet's normal vector. + * @param copy Whether to copy the given array. + * @throws DimensionMismatchException if the length of {@code n} + * does not match the space dimension. + * @throws MaxCountExceededException if the method has been called + * more times than the size of the sphere. + */ + protected void add(double[] normal, + boolean copy) { + if (microsphere.size() >= size) { + throw new MaxCountExceededException(size); + } + if (normal.length > dimension) { + throw new DimensionMismatchException(normal.length, dimension); + } + + microsphere.add(new Facet(copy ? normal.clone() : normal)); + microsphereData.add(new FacetData(0d, 0d)); + } + + /** + * Interpolation. + * + * @return the value estimated from the current illumination of the + * microsphere. + */ + private double interpolate() { + // Number of non-illuminated facets. + int darkCount = 0; + + double value = 0; + double totalWeight = 0; + for (FacetData fd : microsphereData) { + final double iV = fd.illumination(); + if (iV != 0d) { + value += iV * fd.sample(); + totalWeight += iV; + } else { + ++darkCount; + } + } + + final double darkFraction = darkCount / (double) size; + + return darkFraction <= maxDarkFraction ? + value / totalWeight : + background; + } + + /** + * Illumination. + * + * @param sampleDirection Vector whose origin is at the interpolation + * point and tail is at the sample location. + * @param sampleValue Data value of the sample. + * @param weight Weight. + */ + private void illuminate(double[] sampleDirection, + double sampleValue, + double weight) { + for (int i = 0; i < size; i++) { + final double[] n = microsphere.get(i).getNormal(); + final double cos = MathArrays.cosAngle(n, sampleDirection); + + if (cos > 0) { + final double illumination = cos * weight; + + if (illumination > darkThreshold && + illumination > microsphereData.get(i).illumination()) { + microsphereData.set(i, new FacetData(illumination, sampleValue)); + } + } + } + } + + /** + * Reset the all the {@link Facet facets} data to zero. + */ + private void clear() { + for (int i = 0; i < size; i++) { + microsphereData.set(i, new FacetData(0d, 0d)); + } + } + + /** + * Microsphere "facet" (surface element). + */ + private static class Facet { + /** Normal vector characterizing a surface element. */ + private final double[] normal; + + /** + * @param n Normal vector characterizing a surface element + * of the microsphere. No copy is made. + */ + Facet(double[] n) { + normal = n; + } + + /** + * Return a reference to the vector normal to this facet. + * + * @return the normal vector. + */ + public double[] getNormal() { + return normal; + } + } + + /** + * Data associated with each {@link Facet}. + */ + private static class FacetData { + /** Illumination received from the sample. */ + private final double illumination; + /** Data value of the sample. */ + private final double sample; + + /** + * @param illumination Illumination. + * @param sample Data value. + */ + FacetData(double illumination, double sample) { + this.illumination = illumination; + this.sample = sample; + } + + /** + * Get the illumination. + * @return the illumination. + */ + public double illumination() { + return illumination; + } + + /** + * Get the data value. + * @return the data value. + */ + public double sample() { + return sample; + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere2D.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere2D.java new file mode 100644 index 0000000..fdc01b2 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/InterpolatingMicrosphere2D.java @@ -0,0 +1,87 @@ +/* + * 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.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +/** + * Utility class for the {@link MicrosphereProjectionInterpolator} algorithm. + * For 2D interpolation, this class constructs the microsphere as a series of + * evenly spaced facets (rather than generating random normals as in the + * base implementation). + * + * @since 3.6 + */ +public class InterpolatingMicrosphere2D extends InterpolatingMicrosphere { + /** Space dimension. */ + private static final int DIMENSION = 2; + + /** + * Create a sphere from vectors regularly sampled around a circle. + * + * @param size Number of surface elements of the sphere. + * @param maxDarkFraction Maximum fraction of the facets that can be dark. + * If the fraction of "non-illuminated" facets is larger, no estimation + * of the value will be performed, and the {@code background} value will + * be returned instead. + * @param darkThreshold Value of the illumination below which a facet is + * considered dark. + * @param background Value returned when the {@code maxDarkFraction} + * threshold is exceeded. + * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException + * if {@code size <= 0}. + * @throws org.apache.commons.math3.exception.NotPositiveException if + * {@code darkThreshold < 0}. + * @throws org.apache.commons.math3.exception.OutOfRangeException if + * {@code maxDarkFraction} does not belong to the interval {@code [0, 1]}. + */ + public InterpolatingMicrosphere2D(int size, + double maxDarkFraction, + double darkThreshold, + double background) { + super(DIMENSION, size, maxDarkFraction, darkThreshold, background); + + // Generate the microsphere normals. + for (int i = 0; i < size; i++) { + final double angle = i * MathUtils.TWO_PI / size; + + add(new double[] { FastMath.cos(angle), + FastMath.sin(angle) }, + false); + } + } + + /** + * Copy constructor. + * + * @param other Instance to copy. + */ + protected InterpolatingMicrosphere2D(InterpolatingMicrosphere2D other) { + super(other); + } + + /** + * Perform a copy. + * + * @return a copy of this instance. + */ + @Override + public InterpolatingMicrosphere2D copy() { + return new InterpolatingMicrosphere2D(this); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/LinearInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/LinearInterpolator.java new file mode 100644 index 0000000..7e0e69b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/LinearInterpolator.java @@ -0,0 +1,79 @@ +/* + * 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.util.MathArrays; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.util.LocalizedFormats; + +/** + * Implements a linear function for interpolation of real univariate functions. + * + */ +public class LinearInterpolator implements UnivariateInterpolator { + /** + * Computes a linear interpolating function for the data set. + * + * @param x the arguments for the interpolation points + * @param y 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 2. + */ + public PolynomialSplineFunction interpolate(double x[], double y[]) + throws DimensionMismatchException, + NumberIsTooSmallException, + NonMonotonicSequenceException { + if (x.length != y.length) { + throw new DimensionMismatchException(x.length, y.length); + } + + if (x.length < 2) { + throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, + x.length, 2, true); + } + + // Number of intervals. The number of data points is n + 1. + int n = x.length - 1; + + MathArrays.checkOrder(x); + + // Slope of the lines between the datapoints. + final double m[] = new double[n]; + for (int i = 0; i < n; i++) { + m[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]); + } + + final PolynomialFunction polynomials[] = new PolynomialFunction[n]; + final double coefficients[] = new double[2]; + for (int i = 0; i < n; i++) { + coefficients[0] = y[i]; + coefficients[1] = m[i]; + polynomials[i] = new PolynomialFunction(coefficients); + } + + return new PolynomialSplineFunction(x, polynomials); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/LoessInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/LoessInterpolator.java new file mode 100644 index 0000000..3ffbe93 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/LoessInterpolator.java @@ -0,0 +1,473 @@ +/* + * 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 java.io.Serializable; +import java.util.Arrays; + +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.NotFiniteNumberException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.OutOfRangeException; +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.MathUtils; + +/** + * Implements the <a href="http://en.wikipedia.org/wiki/Local_regression"> + * Local Regression Algorithm</a> (also Loess, Lowess) for interpolation of + * real univariate functions. + * <p> + * For reference, see + * <a href="http://amstat.tandfonline.com/doi/abs/10.1080/01621459.1979.10481038"> + * William S. Cleveland - Robust Locally Weighted Regression and Smoothing + * Scatterplots</a> + * </p> + * This class implements both the loess method and serves as an interpolation + * adapter to it, allowing one to build a spline on the obtained loess fit. + * + * @since 2.0 + */ +public class LoessInterpolator + implements UnivariateInterpolator, Serializable { + /** Default value of the bandwidth parameter. */ + public static final double DEFAULT_BANDWIDTH = 0.3; + /** Default value of the number of robustness iterations. */ + public static final int DEFAULT_ROBUSTNESS_ITERS = 2; + /** + * Default value for accuracy. + * @since 2.1 + */ + public static final double DEFAULT_ACCURACY = 1e-12; + /** serializable version identifier. */ + private static final long serialVersionUID = 5204927143605193821L; + /** + * The bandwidth parameter: when computing the loess fit at + * a particular point, this fraction of source points closest + * to the current point is taken into account for computing + * a least-squares regression. + * <p> + * A sensible value is usually 0.25 to 0.5.</p> + */ + private final double bandwidth; + /** + * The number of robustness iterations parameter: this many + * robustness iterations are done. + * <p> + * A sensible value is usually 0 (just the initial fit without any + * robustness iterations) to 4.</p> + */ + private final int robustnessIters; + /** + * If the median residual at a certain robustness iteration + * is less than this amount, no more iterations are done. + */ + private final double accuracy; + + /** + * Constructs a new {@link LoessInterpolator} + * with a bandwidth of {@link #DEFAULT_BANDWIDTH}, + * {@link #DEFAULT_ROBUSTNESS_ITERS} robustness iterations + * and an accuracy of {#link #DEFAULT_ACCURACY}. + * See {@link #LoessInterpolator(double, int, double)} for an explanation of + * the parameters. + */ + public LoessInterpolator() { + this.bandwidth = DEFAULT_BANDWIDTH; + this.robustnessIters = DEFAULT_ROBUSTNESS_ITERS; + this.accuracy = DEFAULT_ACCURACY; + } + + /** + * Construct a new {@link LoessInterpolator} + * with given bandwidth and number of robustness iterations. + * <p> + * Calling this constructor is equivalent to calling {link {@link + * #LoessInterpolator(double, int, double) LoessInterpolator(bandwidth, + * robustnessIters, LoessInterpolator.DEFAULT_ACCURACY)} + * </p> + * + * @param bandwidth when computing the loess fit at + * a particular point, this fraction of source points closest + * to the current point is taken into account for computing + * a least-squares regression. + * A sensible value is usually 0.25 to 0.5, the default value is + * {@link #DEFAULT_BANDWIDTH}. + * @param robustnessIters This many robustness iterations are done. + * A sensible value is usually 0 (just the initial fit without any + * robustness iterations) to 4, the default value is + * {@link #DEFAULT_ROBUSTNESS_ITERS}. + + * @see #LoessInterpolator(double, int, double) + */ + public LoessInterpolator(double bandwidth, int robustnessIters) { + this(bandwidth, robustnessIters, DEFAULT_ACCURACY); + } + + /** + * Construct a new {@link LoessInterpolator} + * with given bandwidth, number of robustness iterations and accuracy. + * + * @param bandwidth when computing the loess fit at + * a particular point, this fraction of source points closest + * to the current point is taken into account for computing + * a least-squares regression. + * A sensible value is usually 0.25 to 0.5, the default value is + * {@link #DEFAULT_BANDWIDTH}. + * @param robustnessIters This many robustness iterations are done. + * A sensible value is usually 0 (just the initial fit without any + * robustness iterations) to 4, the default value is + * {@link #DEFAULT_ROBUSTNESS_ITERS}. + * @param accuracy If the median residual at a certain robustness iteration + * is less than this amount, no more iterations are done. + * @throws OutOfRangeException if bandwidth does not lie in the interval [0,1]. + * @throws NotPositiveException if {@code robustnessIters} is negative. + * @see #LoessInterpolator(double, int) + * @since 2.1 + */ + public LoessInterpolator(double bandwidth, int robustnessIters, double accuracy) + throws OutOfRangeException, + NotPositiveException { + if (bandwidth < 0 || + bandwidth > 1) { + throw new OutOfRangeException(LocalizedFormats.BANDWIDTH, bandwidth, 0, 1); + } + this.bandwidth = bandwidth; + if (robustnessIters < 0) { + throw new NotPositiveException(LocalizedFormats.ROBUSTNESS_ITERATIONS, robustnessIters); + } + this.robustnessIters = robustnessIters; + this.accuracy = accuracy; + } + + /** + * Compute an interpolating function by performing a loess fit + * on the data at the original abscissae and then building a cubic spline + * with a + * {@link org.apache.commons.math3.analysis.interpolation.SplineInterpolator} + * on the resulting fit. + * + * @param xval the arguments for the interpolation points + * @param yval the values for the interpolation points + * @return A cubic spline built upon a loess fit to the data at the original abscissae + * @throws NonMonotonicSequenceException if {@code xval} not sorted in + * strictly increasing order. + * @throws DimensionMismatchException if {@code xval} and {@code yval} have + * different sizes. + * @throws NoDataException if {@code xval} or {@code yval} has zero size. + * @throws NotFiniteNumberException if any of the arguments and values are + * not finite real numbers. + * @throws NumberIsTooSmallException if the bandwidth is too small to + * accomodate the size of the input data (i.e. the bandwidth must be + * larger than 2/n). + */ + public final PolynomialSplineFunction interpolate(final double[] xval, + final double[] yval) + throws NonMonotonicSequenceException, + DimensionMismatchException, + NoDataException, + NotFiniteNumberException, + NumberIsTooSmallException { + return new SplineInterpolator().interpolate(xval, smooth(xval, yval)); + } + + /** + * Compute a weighted loess fit on the data at the original abscissae. + * + * @param xval Arguments for the interpolation points. + * @param yval Values for the interpolation points. + * @param weights point weights: coefficients by which the robustness weight + * of a point is multiplied. + * @return the values of the loess fit at corresponding original abscissae. + * @throws NonMonotonicSequenceException if {@code xval} not sorted in + * strictly increasing order. + * @throws DimensionMismatchException if {@code xval} and {@code yval} have + * different sizes. + * @throws NoDataException if {@code xval} or {@code yval} has zero size. + * @throws NotFiniteNumberException if any of the arguments and values are + not finite real numbers. + * @throws NumberIsTooSmallException if the bandwidth is too small to + * accomodate the size of the input data (i.e. the bandwidth must be + * larger than 2/n). + * @since 2.1 + */ + public final double[] smooth(final double[] xval, final double[] yval, + final double[] weights) + throws NonMonotonicSequenceException, + DimensionMismatchException, + NoDataException, + NotFiniteNumberException, + NumberIsTooSmallException { + if (xval.length != yval.length) { + throw new DimensionMismatchException(xval.length, yval.length); + } + + final int n = xval.length; + + if (n == 0) { + throw new NoDataException(); + } + + checkAllFiniteReal(xval); + checkAllFiniteReal(yval); + checkAllFiniteReal(weights); + + MathArrays.checkOrder(xval); + + if (n == 1) { + return new double[]{yval[0]}; + } + + if (n == 2) { + return new double[]{yval[0], yval[1]}; + } + + int bandwidthInPoints = (int) (bandwidth * n); + + if (bandwidthInPoints < 2) { + throw new NumberIsTooSmallException(LocalizedFormats.BANDWIDTH, + bandwidthInPoints, 2, true); + } + + final double[] res = new double[n]; + + final double[] residuals = new double[n]; + final double[] sortedResiduals = new double[n]; + + final double[] robustnessWeights = new double[n]; + + // Do an initial fit and 'robustnessIters' robustness iterations. + // This is equivalent to doing 'robustnessIters+1' robustness iterations + // starting with all robustness weights set to 1. + Arrays.fill(robustnessWeights, 1); + + for (int iter = 0; iter <= robustnessIters; ++iter) { + final int[] bandwidthInterval = {0, bandwidthInPoints - 1}; + // At each x, compute a local weighted linear regression + for (int i = 0; i < n; ++i) { + final double x = xval[i]; + + // Find out the interval of source points on which + // a regression is to be made. + if (i > 0) { + updateBandwidthInterval(xval, weights, i, bandwidthInterval); + } + + final int ileft = bandwidthInterval[0]; + final int iright = bandwidthInterval[1]; + + // Compute the point of the bandwidth interval that is + // farthest from x + final int edge; + if (xval[i] - xval[ileft] > xval[iright] - xval[i]) { + edge = ileft; + } else { + edge = iright; + } + + // Compute a least-squares linear fit weighted by + // the product of robustness weights and the tricube + // weight function. + // See http://en.wikipedia.org/wiki/Linear_regression + // (section "Univariate linear case") + // and http://en.wikipedia.org/wiki/Weighted_least_squares + // (section "Weighted least squares") + double sumWeights = 0; + double sumX = 0; + double sumXSquared = 0; + double sumY = 0; + double sumXY = 0; + double denom = FastMath.abs(1.0 / (xval[edge] - x)); + for (int k = ileft; k <= iright; ++k) { + final double xk = xval[k]; + final double yk = yval[k]; + final double dist = (k < i) ? x - xk : xk - x; + final double w = tricube(dist * denom) * robustnessWeights[k] * weights[k]; + final double xkw = xk * w; + sumWeights += w; + sumX += xkw; + sumXSquared += xk * xkw; + sumY += yk * w; + sumXY += yk * xkw; + } + + final double meanX = sumX / sumWeights; + final double meanY = sumY / sumWeights; + final double meanXY = sumXY / sumWeights; + final double meanXSquared = sumXSquared / sumWeights; + + final double beta; + if (FastMath.sqrt(FastMath.abs(meanXSquared - meanX * meanX)) < accuracy) { + beta = 0; + } else { + beta = (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX); + } + + final double alpha = meanY - beta * meanX; + + res[i] = beta * x + alpha; + residuals[i] = FastMath.abs(yval[i] - res[i]); + } + + // No need to recompute the robustness weights at the last + // iteration, they won't be needed anymore + if (iter == robustnessIters) { + break; + } + + // Recompute the robustness weights. + + // Find the median residual. + // An arraycopy and a sort are completely tractable here, + // because the preceding loop is a lot more expensive + System.arraycopy(residuals, 0, sortedResiduals, 0, n); + Arrays.sort(sortedResiduals); + final double medianResidual = sortedResiduals[n / 2]; + + if (FastMath.abs(medianResidual) < accuracy) { + break; + } + + for (int i = 0; i < n; ++i) { + final double arg = residuals[i] / (6 * medianResidual); + if (arg >= 1) { + robustnessWeights[i] = 0; + } else { + final double w = 1 - arg * arg; + robustnessWeights[i] = w * w; + } + } + } + + return res; + } + + /** + * Compute a loess fit on the data at the original abscissae. + * + * @param xval the arguments for the interpolation points + * @param yval the values for the interpolation points + * @return values of the loess fit at corresponding original abscissae + * @throws NonMonotonicSequenceException if {@code xval} not sorted in + * strictly increasing order. + * @throws DimensionMismatchException if {@code xval} and {@code yval} have + * different sizes. + * @throws NoDataException if {@code xval} or {@code yval} has zero size. + * @throws NotFiniteNumberException if any of the arguments and values are + * not finite real numbers. + * @throws NumberIsTooSmallException if the bandwidth is too small to + * accomodate the size of the input data (i.e. the bandwidth must be + * larger than 2/n). + */ + public final double[] smooth(final double[] xval, final double[] yval) + throws NonMonotonicSequenceException, + DimensionMismatchException, + NoDataException, + NotFiniteNumberException, + NumberIsTooSmallException { + if (xval.length != yval.length) { + throw new DimensionMismatchException(xval.length, yval.length); + } + + final double[] unitWeights = new double[xval.length]; + Arrays.fill(unitWeights, 1.0); + + return smooth(xval, yval, unitWeights); + } + + /** + * Given an index interval into xval that embraces a certain number of + * points closest to {@code xval[i-1]}, update the interval so that it + * embraces the same number of points closest to {@code xval[i]}, + * ignoring zero weights. + * + * @param xval Arguments array. + * @param weights Weights array. + * @param i Index around which the new interval should be computed. + * @param bandwidthInterval a two-element array {left, right} such that: + * {@code (left==0 or xval[i] - xval[left-1] > xval[right] - xval[i])} + * and + * {@code (right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])}. + * The array will be updated. + */ + private static void updateBandwidthInterval(final double[] xval, final double[] weights, + final int i, + final int[] bandwidthInterval) { + final int left = bandwidthInterval[0]; + final int right = bandwidthInterval[1]; + + // The right edge should be adjusted if the next point to the right + // is closer to xval[i] than the leftmost point of the current interval + int nextRight = nextNonzero(weights, right); + if (nextRight < xval.length && xval[nextRight] - xval[i] < xval[i] - xval[left]) { + int nextLeft = nextNonzero(weights, bandwidthInterval[0]); + bandwidthInterval[0] = nextLeft; + bandwidthInterval[1] = nextRight; + } + } + + /** + * Return the smallest index {@code j} such that + * {@code j > i && (j == weights.length || weights[j] != 0)}. + * + * @param weights Weights array. + * @param i Index from which to start search. + * @return the smallest compliant index. + */ + private static int nextNonzero(final double[] weights, final int i) { + int j = i + 1; + while(j < weights.length && weights[j] == 0) { + ++j; + } + return j; + } + + /** + * Compute the + * <a href="http://en.wikipedia.org/wiki/Local_regression#Weight_function">tricube</a> + * weight function + * + * @param x Argument. + * @return <code>(1 - |x|<sup>3</sup>)<sup>3</sup></code> for |x| < 1, 0 otherwise. + */ + private static double tricube(final double x) { + final double absX = FastMath.abs(x); + if (absX >= 1.0) { + return 0.0; + } + final double tmp = 1 - absX * absX * absX; + return tmp * tmp * tmp; + } + + /** + * Check that all elements of an array are finite real numbers. + * + * @param values Values array. + * @throws org.apache.commons.math3.exception.NotFiniteNumberException + * if one of the values is not a finite real number. + */ + private static void checkAllFiniteReal(final double[] values) { + for (int i = 0; i < values.length; i++) { + MathUtils.checkFinite(values[i]); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java new file mode 100644 index 0000000..58be772 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java @@ -0,0 +1,253 @@ +/* + * 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 java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + +import org.apache.commons.math3.analysis.MultivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; +import org.apache.commons.math3.util.FastMath; + +/** + * Interpolating function that implements the + * <a href="http://www.dudziak.com/microsphere.php">Microsphere Projection</a>. + * + * @deprecated Code will be removed in 4.0. Use {@link InterpolatingMicrosphere} + * and {@link MicrosphereProjectionInterpolator} instead. + */ +@Deprecated +public class MicrosphereInterpolatingFunction + implements MultivariateFunction { + /** + * Space dimension. + */ + private final int dimension; + /** + * Internal accounting data for the interpolation algorithm. + * Each element of the list corresponds to one surface element of + * the microsphere. + */ + private final List<MicrosphereSurfaceElement> microsphere; + /** + * Exponent used in the power law that computes the weights of the + * sample data. + */ + private final double brightnessExponent; + /** + * Sample data. + */ + private final Map<RealVector, Double> samples; + + /** + * Class for storing the accounting data needed to perform the + * microsphere projection. + */ + private static class MicrosphereSurfaceElement { + /** Normal vector characterizing a surface element. */ + private final RealVector normal; + /** Illumination received from the brightest sample. */ + private double brightestIllumination; + /** Brightest sample. */ + private Map.Entry<RealVector, Double> brightestSample; + + /** + * @param n Normal vector characterizing a surface element + * of the microsphere. + */ + MicrosphereSurfaceElement(double[] n) { + normal = new ArrayRealVector(n); + } + + /** + * Return the normal vector. + * @return the normal vector + */ + RealVector normal() { + return normal; + } + + /** + * Reset "illumination" and "sampleIndex". + */ + void reset() { + brightestIllumination = 0; + brightestSample = null; + } + + /** + * Store the illumination and index of the brightest sample. + * @param illuminationFromSample illumination received from sample + * @param sample current sample illuminating the element + */ + void store(final double illuminationFromSample, + final Map.Entry<RealVector, Double> sample) { + if (illuminationFromSample > this.brightestIllumination) { + this.brightestIllumination = illuminationFromSample; + this.brightestSample = sample; + } + } + + /** + * Get the illumination of the element. + * @return the illumination. + */ + double illumination() { + return brightestIllumination; + } + + /** + * Get the sample illuminating the element the most. + * @return the sample. + */ + Map.Entry<RealVector, Double> sample() { + return brightestSample; + } + } + + /** + * @param xval Arguments for the interpolation points. + * {@code xval[i][0]} is the first component of interpolation point + * {@code i}, {@code xval[i][1]} is the second component, and so on + * until {@code xval[i][d-1]}, the last component of that interpolation + * point (where {@code dimension} is thus the dimension of the sampled + * space). + * @param yval Values for the interpolation points. + * @param brightnessExponent Brightness dimming factor. + * @param microsphereElements Number of surface elements of the + * microsphere. + * @param rand Unit vector generator for creating the microsphere. + * @throws DimensionMismatchException if the lengths of {@code yval} and + * {@code xval} (equal to {@code n}, the number of interpolation points) + * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]}, + * have lengths different from {@code dimension}. + * @throws NoDataException if there an array has zero-length. + * @throws NullArgumentException if an argument is {@code null}. + */ + public MicrosphereInterpolatingFunction(double[][] xval, + double[] yval, + int brightnessExponent, + int microsphereElements, + UnitSphereRandomVectorGenerator rand) + throws DimensionMismatchException, + NoDataException, + NullArgumentException { + if (xval == null || + yval == null) { + throw new NullArgumentException(); + } + if (xval.length == 0) { + throw new NoDataException(); + } + if (xval.length != yval.length) { + throw new DimensionMismatchException(xval.length, yval.length); + } + if (xval[0] == null) { + throw new NullArgumentException(); + } + + dimension = xval[0].length; + this.brightnessExponent = brightnessExponent; + + // Copy data samples. + samples = new HashMap<RealVector, Double>(yval.length); + for (int i = 0; i < xval.length; ++i) { + final double[] xvalI = xval[i]; + if (xvalI == null) { + throw new NullArgumentException(); + } + if (xvalI.length != dimension) { + throw new DimensionMismatchException(xvalI.length, dimension); + } + + samples.put(new ArrayRealVector(xvalI), yval[i]); + } + + microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements); + // Generate the microsphere, assuming that a fairly large number of + // randomly generated normals will represent a sphere. + for (int i = 0; i < microsphereElements; i++) { + microsphere.add(new MicrosphereSurfaceElement(rand.nextVector())); + } + } + + /** + * @param point Interpolation point. + * @return the interpolated value. + * @throws DimensionMismatchException if point dimension does not math sample + */ + public double value(double[] point) throws DimensionMismatchException { + final RealVector p = new ArrayRealVector(point); + + // Reset. + for (MicrosphereSurfaceElement md : microsphere) { + md.reset(); + } + + // Compute contribution of each sample points to the microsphere elements illumination + for (Map.Entry<RealVector, Double> sd : samples.entrySet()) { + + // Vector between interpolation point and current sample point. + final RealVector diff = sd.getKey().subtract(p); + final double diffNorm = diff.getNorm(); + + if (FastMath.abs(diffNorm) < FastMath.ulp(1d)) { + // No need to interpolate, as the interpolation point is + // actually (very close to) one of the sampled points. + return sd.getValue(); + } + + for (MicrosphereSurfaceElement md : microsphere) { + final double w = FastMath.pow(diffNorm, -brightnessExponent); + md.store(cosAngle(diff, md.normal()) * w, sd); + } + + } + + // Interpolation calculation. + double value = 0; + double totalWeight = 0; + for (MicrosphereSurfaceElement md : microsphere) { + final double iV = md.illumination(); + final Map.Entry<RealVector, Double> sd = md.sample(); + if (sd != null) { + value += iV * sd.getValue(); + totalWeight += iV; + } + } + + return value / totalWeight; + } + + /** + * Compute the cosine of the angle between 2 vectors. + * + * @param v Vector. + * @param w Vector. + * @return the cosine of the angle between {@code v} and {@code w}. + */ + private double cosAngle(final RealVector v, final RealVector w) { + return v.dotProduct(w) / (v.getNorm() * w.getNorm()); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java new file mode 100644 index 0000000..d9174bc --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java @@ -0,0 +1,105 @@ +/* + * 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.MultivariateFunction; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; + +/** + * Interpolator that implements the algorithm described in + * <em>William Dudziak</em>'s + * <a href="http://www.dudziak.com/microsphere.pdf">MS thesis</a>. + * + * @since 2.1 + * @deprecated Code will be removed in 4.0. Use {@link InterpolatingMicrosphere} + * and {@link MicrosphereProjectionInterpolator} instead. + */ +@Deprecated +public class MicrosphereInterpolator + implements MultivariateInterpolator { + /** + * Default number of surface elements that composes the microsphere. + */ + public static final int DEFAULT_MICROSPHERE_ELEMENTS = 2000; + /** + * Default exponent used the weights calculation. + */ + public static final int DEFAULT_BRIGHTNESS_EXPONENT = 2; + /** + * Number of surface elements of the microsphere. + */ + private final int microsphereElements; + /** + * Exponent used in the power law that computes the weights of the + * sample data. + */ + private final int brightnessExponent; + + /** + * Create a microsphere interpolator with default settings. + * Calling this constructor is equivalent to call {@link + * #MicrosphereInterpolator(int, int) + * MicrosphereInterpolator(MicrosphereInterpolator.DEFAULT_MICROSPHERE_ELEMENTS, + * MicrosphereInterpolator.DEFAULT_BRIGHTNESS_EXPONENT)}. + */ + public MicrosphereInterpolator() { + this(DEFAULT_MICROSPHERE_ELEMENTS, DEFAULT_BRIGHTNESS_EXPONENT); + } + + /** Create a microsphere interpolator. + * @param elements Number of surface elements of the microsphere. + * @param exponent Exponent used in the power law that computes the + * weights (distance dimming factor) of the sample data. + * @throws NotPositiveException if {@code exponent < 0}. + * @throws NotStrictlyPositiveException if {@code elements <= 0}. + */ + public MicrosphereInterpolator(final int elements, + final int exponent) + throws NotPositiveException, + NotStrictlyPositiveException { + if (exponent < 0) { + throw new NotPositiveException(exponent); + } + if (elements <= 0) { + throw new NotStrictlyPositiveException(elements); + } + + microsphereElements = elements; + brightnessExponent = exponent; + } + + /** + * {@inheritDoc} + */ + public MultivariateFunction interpolate(final double[][] xval, + final double[] yval) + throws DimensionMismatchException, + NoDataException, + NullArgumentException { + final UnitSphereRandomVectorGenerator rand + = new UnitSphereRandomVectorGenerator(xval[0].length); + return new MicrosphereInterpolatingFunction(xval, yval, + brightnessExponent, + microsphereElements, + rand); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolator.java new file mode 100644 index 0000000..28f5b26 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereProjectionInterpolator.java @@ -0,0 +1,164 @@ +/* + * 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.MultivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; + +/** + * Interpolator that implements the algorithm described in + * <em>William Dudziak</em>'s + * <a href="http://www.dudziak.com/microsphere.pdf">MS thesis</a>. + * + * @since 3.6 + */ +public class MicrosphereProjectionInterpolator + implements MultivariateInterpolator { + /** Brightness exponent. */ + private final double exponent; + /** Microsphere. */ + private final InterpolatingMicrosphere microsphere; + /** Whether to share the sphere. */ + private final boolean sharedSphere; + /** Tolerance value below which no interpolation is necessary. */ + private final double noInterpolationTolerance; + + /** + * Create a microsphere interpolator. + * + * @param dimension Space dimension. + * @param elements Number of surface elements of the microsphere. + * @param exponent Exponent used in the power law that computes the + * @param maxDarkFraction Maximum fraction of the facets that can be dark. + * If the fraction of "non-illuminated" facets is larger, no estimation + * of the value will be performed, and the {@code background} value will + * be returned instead. + * @param darkThreshold Value of the illumination below which a facet is + * considered dark. + * @param background Value returned when the {@code maxDarkFraction} + * threshold is exceeded. + * @param sharedSphere Whether the sphere can be shared among the + * interpolating function instances. If {@code true}, the instances + * will share the same data, and thus will <em>not</em> be thread-safe. + * @param noInterpolationTolerance When the distance between an + * interpolated point and one of the sample points is less than this + * value, no interpolation will be performed (the value of the sample + * will be returned). + * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException + * if {@code dimension <= 0} or {@code elements <= 0}. + * @throws NotPositiveException if {@code exponent < 0}. + * @throws NotPositiveException if {@code darkThreshold < 0}. + * @throws org.apache.commons.math3.exception.OutOfRangeException if + * {@code maxDarkFraction} does not belong to the interval {@code [0, 1]}. + */ + public MicrosphereProjectionInterpolator(int dimension, + int elements, + double maxDarkFraction, + double darkThreshold, + double background, + double exponent, + boolean sharedSphere, + double noInterpolationTolerance) { + this(new InterpolatingMicrosphere(dimension, + elements, + maxDarkFraction, + darkThreshold, + background, + new UnitSphereRandomVectorGenerator(dimension)), + exponent, + sharedSphere, + noInterpolationTolerance); + } + + /** + * Create a microsphere interpolator. + * + * @param microsphere Microsphere. + * @param exponent Exponent used in the power law that computes the + * weights (distance dimming factor) of the sample data. + * @param sharedSphere Whether the sphere can be shared among the + * interpolating function instances. If {@code true}, the instances + * will share the same data, and thus will <em>not</em> be thread-safe. + * @param noInterpolationTolerance When the distance between an + * interpolated point and one of the sample points is less than this + * value, no interpolation will be performed (the value of the sample + * will be returned). + * @throws NotPositiveException if {@code exponent < 0}. + */ + public MicrosphereProjectionInterpolator(InterpolatingMicrosphere microsphere, + double exponent, + boolean sharedSphere, + double noInterpolationTolerance) + throws NotPositiveException { + if (exponent < 0) { + throw new NotPositiveException(exponent); + } + + this.microsphere = microsphere; + this.exponent = exponent; + this.sharedSphere = sharedSphere; + this.noInterpolationTolerance = noInterpolationTolerance; + } + + /** + * {@inheritDoc} + * + * @throws DimensionMismatchException if the space dimension of the + * given samples does not match the space dimension of the microsphere. + */ + public MultivariateFunction interpolate(final double[][] xval, + final double[] yval) + throws DimensionMismatchException, + NoDataException, + NullArgumentException { + if (xval == null || + yval == null) { + throw new NullArgumentException(); + } + if (xval.length == 0) { + throw new NoDataException(); + } + if (xval.length != yval.length) { + throw new DimensionMismatchException(xval.length, yval.length); + } + if (xval[0] == null) { + throw new NullArgumentException(); + } + final int dimension = microsphere.getDimension(); + if (dimension != xval[0].length) { + throw new DimensionMismatchException(xval[0].length, dimension); + } + + // Microsphere copy. + final InterpolatingMicrosphere m = sharedSphere ? microsphere : microsphere.copy(); + + return new MultivariateFunction() { + /** {inheritDoc} */ + public double value(double[] point) { + return m.value(point, + xval, + yval, + exponent, + noInterpolationTolerance); + } + }; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java new file mode 100644 index 0000000..7d76374 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java @@ -0,0 +1,51 @@ +/* + * 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.MultivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NullArgumentException; + +/** + * Interface representing a univariate real interpolating function. + * + * @since 2.1 + */ +public interface MultivariateInterpolator { + + /** + * Computes an interpolating function for the data set. + * + * @param xval the arguments for the interpolation points. + * {@code xval[i][0]} is the first component of interpolation point + * {@code i}, {@code xval[i][1]} is the second component, and so on + * until {@code xval[i][d-1]}, the last component of that interpolation + * point (where {@code d} is thus the dimension of the space). + * @param yval the values for the interpolation points + * @return a function which interpolates the data set + * @throws MathIllegalArgumentException if the arguments violate assumptions + * made by the interpolation algorithm. + * @throws DimensionMismatchException when the array dimensions are not consistent. + * @throws NoDataException if an array has zero-length. + * @throws NullArgumentException if the arguments are {@code null}. + */ + MultivariateFunction interpolate(double[][] xval, double[] yval) + throws MathIllegalArgumentException, DimensionMismatchException, + NoDataException, NullArgumentException; +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java new file mode 100644 index 0000000..6b47451 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java @@ -0,0 +1,60 @@ +/* + * 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 java.io.Serializable; + +import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; + +/** + * Implements the <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html"> + * Neville's Algorithm</a> for interpolation of real univariate functions. For + * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, + * chapter 2. + * <p> + * The actual code of Neville's algorithm is in PolynomialFunctionLagrangeForm, + * this class provides an easy-to-use interface to it.</p> + * + * @since 1.2 + */ +public class NevilleInterpolator implements UnivariateInterpolator, + Serializable { + + /** serializable version identifier */ + static final long serialVersionUID = 3003707660147873733L; + + /** + * Computes an interpolating function for the data set. + * + * @param x Interpolating points. + * @param y Interpolating values. + * @return a function which interpolates the data set + * @throws DimensionMismatchException if the array lengths are different. + * @throws NumberIsTooSmallException if the number of points is less than 2. + * @throws NonMonotonicSequenceException if two abscissae have the same + * value. + */ + public PolynomialFunctionLagrangeForm interpolate(double x[], double y[]) + throws DimensionMismatchException, + NumberIsTooSmallException, + NonMonotonicSequenceException { + return new PolynomialFunctionLagrangeForm(x, y); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java new file mode 100644 index 0000000..7dd135a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java @@ -0,0 +1,210 @@ +/* + * 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 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://www.paulinternet.nl/?page=bicubic">bicubic spline</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 3.4 + */ +public class PiecewiseBicubicSplineInterpolatingFunction + 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 */ + private final double[] xval; + /** Samples y-coordinates */ + private final double[] yval; + /** Set of cubic splines patching the whole data grid */ + 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. 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 PiecewiseBicubicSplineInterpolatingFunction(double[] x, + double[] y, + double[][] f) + throws DimensionMismatchException, + NullArgumentException, + NoDataException, + NonMonotonicSequenceException { + 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 < MIN_NUM_POINTS || + yLen < MIN_NUM_POINTS || + f.length < MIN_NUM_POINTS || + f[0].length < MIN_NUM_POINTS) { + throw new InsufficientDataException(); + } + + if (xLen != f.length) { + throw new DimensionMismatchException(xLen, f.length); + } + + if (yLen != f[0].length) { + throw new DimensionMismatchException(yLen, f[0].length); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + + xval = x.clone(); + yval = y.clone(); + fval = f.clone(); + } + + /** + * {@inheritDoc} + */ + public double value(double x, + double y) + throws OutOfRangeException { + final 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 xArray[] = new double[count]; + final double yArray[] = new double[count]; + final double zArray[] = new double[count]; + final double interpArray[] = new double[count]; + + for (int index = 0; index < count; index++) { + xArray[index] = xval[i + index]; + yArray[index] = yval[j + index]; + } + + for (int zIndex = 0; zIndex < count; zIndex++) { + for (int index = 0; index < count; index++) { + zArray[index] = fval[i + index][j + zIndex]; + } + final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray); + interpArray[zIndex] = spline.value(x); + } + + final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray); + + double returnValue = spline.value(y); + + return returnValue; + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @return {@code true} if (x, y) is a valid point. + * @since 3.3 + */ + public boolean isValidPoint(double x, + double y) { + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || + y > yval[yval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @param c Coordinate. + * @param val Coordinate samples. + * @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, which returns + // negative + // need to remove the negative sign for consistency + r = -r - offset - 1; + } else { + r -= offset; + } + + 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. + r = val.length - count; + } + + return r; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java new file mode 100644 index 0000000..826f328 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java @@ -0,0 +1,61 @@ +/* + * 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.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.util.MathArrays; + +/** + * Generates a piecewise-bicubic interpolating function. + * + * @since 2.2 + */ +public class PiecewiseBicubicSplineInterpolator + implements BivariateGridInterpolator { + + /** + * {@inheritDoc} + */ + public PiecewiseBicubicSplineInterpolatingFunction interpolate( final double[] xval, + final double[] yval, + final double[][] fval) + throws DimensionMismatchException, + NullArgumentException, + NoDataException, + NonMonotonicSequenceException { + if ( xval == null || + yval == null || + fval == null || + fval[0] == null ) { + throw new NullArgumentException(); + } + + if ( xval.length == 0 || + yval.length == 0 || + fval.length == 0 ) { + throw new NoDataException(); + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + + return new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, fval ); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java new file mode 100644 index 0000000..e1639b2 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java @@ -0,0 +1,171 @@ +/* + * 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.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.Precision; +import org.apache.commons.math3.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer; +import org.apache.commons.math3.fitting.PolynomialFitter; +import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math3.optim.SimpleVectorValueChecker; + +/** + * Generates a bicubic interpolation function. + * Prior to generating the interpolating function, the input is smoothed using + * polynomial fitting. + * + * @since 2.2 + * @deprecated To be removed in 4.0 (see MATH-1166). + */ +@Deprecated +public class SmoothingPolynomialBicubicSplineInterpolator + extends BicubicSplineInterpolator { + /** Fitter for x. */ + private final PolynomialFitter xFitter; + /** Degree of the fitting polynomial. */ + private final int xDegree; + /** Fitter for y. */ + private final PolynomialFitter yFitter; + /** Degree of the fitting polynomial. */ + private final int yDegree; + + /** + * Default constructor. The degree of the fitting polynomials is set to 3. + */ + public SmoothingPolynomialBicubicSplineInterpolator() { + this(3); + } + + /** + * @param degree Degree of the polynomial fitting functions. + * @exception NotPositiveException if degree is not positive + */ + public SmoothingPolynomialBicubicSplineInterpolator(int degree) + throws NotPositiveException { + this(degree, degree); + } + + /** + * @param xDegree Degree of the polynomial fitting functions along the + * x-dimension. + * @param yDegree Degree of the polynomial fitting functions along the + * y-dimension. + * @exception NotPositiveException if degrees are not positive + */ + public SmoothingPolynomialBicubicSplineInterpolator(int xDegree, int yDegree) + throws NotPositiveException { + if (xDegree < 0) { + throw new NotPositiveException(xDegree); + } + if (yDegree < 0) { + throw new NotPositiveException(yDegree); + } + this.xDegree = xDegree; + this.yDegree = yDegree; + + final double safeFactor = 1e2; + final SimpleVectorValueChecker checker + = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON, + safeFactor * Precision.SAFE_MIN); + xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker)); + yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker)); + } + + /** + * {@inheritDoc} + */ + @Override + public BicubicSplineInterpolatingFunction interpolate(final double[] xval, + final double[] yval, + final double[][] fval) + throws NoDataException, NullArgumentException, + DimensionMismatchException, NonMonotonicSequenceException { + if (xval.length == 0 || yval.length == 0 || fval.length == 0) { + throw new NoDataException(); + } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } + + final int xLen = xval.length; + final int yLen = yval.length; + + for (int i = 0; i < xLen; i++) { + if (fval[i].length != yLen) { + throw new DimensionMismatchException(fval[i].length, yLen); + } + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + + // For each line y[j] (0 <= j < yLen), construct a polynomial, with + // respect to variable x, fitting array fval[][j] + final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen]; + for (int j = 0; j < yLen; j++) { + xFitter.clearObservations(); + for (int i = 0; i < xLen; i++) { + xFitter.addObservedPoint(1, xval[i], fval[i][j]); + } + + // Initial guess for the fit is zero for each coefficients (of which + // there are "xDegree" + 1). + yPolyX[j] = new PolynomialFunction(xFitter.fit(new double[xDegree + 1])); + } + + // For every knot (xval[i], yval[j]) of the grid, calculate corrected + // values fval_1 + final double[][] fval_1 = new double[xLen][yLen]; + for (int j = 0; j < yLen; j++) { + final PolynomialFunction f = yPolyX[j]; + for (int i = 0; i < xLen; i++) { + fval_1[i][j] = f.value(xval[i]); + } + } + + // For each line x[i] (0 <= i < xLen), construct a polynomial, with + // respect to variable y, fitting array fval_1[i][] + final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen]; + for (int i = 0; i < xLen; i++) { + yFitter.clearObservations(); + for (int j = 0; j < yLen; j++) { + yFitter.addObservedPoint(1, yval[j], fval_1[i][j]); + } + + // Initial guess for the fit is zero for each coefficients (of which + // there are "yDegree" + 1). + xPolyY[i] = new PolynomialFunction(yFitter.fit(new double[yDegree + 1])); + } + + // For every knot (xval[i], yval[j]) of the grid, calculate corrected + // values fval_2 + final double[][] fval_2 = new double[xLen][yLen]; + for (int i = 0; i < xLen; i++) { + final PolynomialFunction f = xPolyY[i]; + for (int j = 0; j < yLen; j++) { + fval_2[i][j] = f.value(yval[j]); + } + } + + return super.interpolate(xval, yval, fval_2); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java new file mode 100644 index 0000000..f37e1b1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java @@ -0,0 +1,127 @@ +/* + * 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.NumberIsTooSmallException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.MathArrays; + +/** + * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set. + * <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, + * {@code x[0] < x[i] ... < x[n].} The x values are referred to as "knot points." + * <p> + * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest + * knot point and strictly less than the largest knot point is computed by finding the subinterval to which + * x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where + * <code>i</code> is the index of the subinterval. See {@link PolynomialSplineFunction} for more details. + * </p> + * <p> + * The interpolating polynomials satisfy: <ol> + * <li>The value of the PolynomialSplineFunction at each of the input x values equals the + * corresponding y value.</li> + * <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials + * "match up" at the knot points, as do their first and second derivatives).</li> + * </ol> + * <p> + * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires, + * <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131. + * </p> + * + */ +public class SplineInterpolator implements UnivariateInterpolator { + /** + * Computes an interpolating function for the data set. + * @param x the arguments for the interpolation points + * @param y 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 3. + */ + public PolynomialSplineFunction interpolate(double x[], double y[]) + throws DimensionMismatchException, + NumberIsTooSmallException, + NonMonotonicSequenceException { + if (x.length != y.length) { + throw new DimensionMismatchException(x.length, y.length); + } + + if (x.length < 3) { + throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, + x.length, 3, true); + } + + // Number of intervals. The number of data points is n + 1. + final int n = x.length - 1; + + MathArrays.checkOrder(x); + + // Differences between knot points + final double h[] = new double[n]; + for (int i = 0; i < n; i++) { + h[i] = x[i + 1] - x[i]; + } + + final double mu[] = new double[n]; + final double z[] = new double[n + 1]; + mu[0] = 0d; + z[0] = 0d; + double g = 0; + for (int i = 1; i < n; i++) { + g = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1]; + mu[i] = h[i] / g; + z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) / + (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g; + } + + // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) + final double b[] = new double[n]; + final double c[] = new double[n + 1]; + final double d[] = new double[n]; + + z[n] = 0d; + c[n] = 0d; + + for (int j = n -1; j >=0; j--) { + c[j] = z[j] - mu[j] * c[j + 1]; + b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d; + d[j] = (c[j + 1] - c[j]) / (3d * h[j]); + } + + final PolynomialFunction polynomials[] = new PolynomialFunction[n]; + final double coefficients[] = new double[4]; + for (int i = 0; i < n; i++) { + coefficients[0] = y[i]; + coefficients[1] = b[i]; + coefficients[2] = c[i]; + coefficients[3] = d[i]; + polynomials[i] = new PolynomialFunction(coefficients); + } + + return new PolynomialSplineFunction(x, polynomials); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java new file mode 100644 index 0000000..27e9a65 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java @@ -0,0 +1,508 @@ +/* + * 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.TrivariateFunction; +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.OutOfRangeException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Function that implements the + * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation"> + * tricubic spline interpolation</a>, as proposed in + * <blockquote> + * Tricubic interpolation in three dimensions, + * F. Lekien and J. Marsden, + * <em>Int. J. Numer. Meth. Eng</em> 2005; <b>63</b>:455-471 + * </blockquote> + * + * @since 3.4. + */ +public class TricubicInterpolatingFunction + implements TrivariateFunction { + /** + * Matrix to compute the spline coefficients from the function values + * and function derivatives values + */ + private static final double[][] AINV = { + { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 }, + {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 }, + { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 }, + { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 }, + { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 }, + { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 }, + { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 }, + { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 }, + { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 }, + { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 }, + { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 }, + { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 }, + { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 }, + { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 }, + { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 } + }; + + /** Samples x-coordinates */ + private final double[] xval; + /** Samples y-coordinates */ + private final double[] yval; + /** Samples z-coordinates */ + private final double[] zval; + /** Set of cubic splines patching the whole data grid */ + private final TricubicFunction[][][] splines; + + /** + * @param x Sample values of the x-coordinate, in increasing order. + * @param y Sample values of the y-coordinate, in increasing order. + * @param z Sample values of the y-coordinate, in increasing order. + * @param f Values of the function on every grid point. + * @param dFdX Values of the partial derivative of function with respect to x on every grid point. + * @param dFdY Values of the partial derivative of function with respect to y on every grid point. + * @param dFdZ Values of the partial derivative of function with respect to z on every grid point. + * @param d2FdXdY Values of the cross partial derivative of function on every grid point. + * @param d2FdXdZ Values of the cross partial derivative of function on every grid point. + * @param d2FdYdZ Values of the cross partial derivative of function on every grid point. + * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point. + * @throws NoDataException if any of the arrays has zero length. + * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements. + * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing. + */ + public TricubicInterpolatingFunction(double[] x, + double[] y, + double[] z, + double[][][] f, + double[][][] dFdX, + double[][][] dFdY, + double[][][] dFdZ, + double[][][] d2FdXdY, + double[][][] d2FdXdZ, + double[][][] d2FdYdZ, + double[][][] d3FdXdYdZ) + throws NoDataException, + DimensionMismatchException, + NonMonotonicSequenceException { + final int xLen = x.length; + final int yLen = y.length; + final int zLen = z.length; + + if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) { + throw new NoDataException(); + } + if (xLen != f.length) { + throw new DimensionMismatchException(xLen, f.length); + } + if (xLen != dFdX.length) { + throw new DimensionMismatchException(xLen, dFdX.length); + } + if (xLen != dFdY.length) { + throw new DimensionMismatchException(xLen, dFdY.length); + } + if (xLen != dFdZ.length) { + throw new DimensionMismatchException(xLen, dFdZ.length); + } + if (xLen != d2FdXdY.length) { + throw new DimensionMismatchException(xLen, d2FdXdY.length); + } + if (xLen != d2FdXdZ.length) { + throw new DimensionMismatchException(xLen, d2FdXdZ.length); + } + if (xLen != d2FdYdZ.length) { + throw new DimensionMismatchException(xLen, d2FdYdZ.length); + } + if (xLen != d3FdXdYdZ.length) { + throw new DimensionMismatchException(xLen, d3FdXdYdZ.length); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + MathArrays.checkOrder(z); + + xval = x.clone(); + yval = y.clone(); + zval = z.clone(); + + final int lastI = xLen - 1; + final int lastJ = yLen - 1; + final int lastK = zLen - 1; + splines = new TricubicFunction[lastI][lastJ][lastK]; + + for (int i = 0; i < lastI; i++) { + if (f[i].length != yLen) { + throw new DimensionMismatchException(f[i].length, yLen); + } + if (dFdX[i].length != yLen) { + throw new DimensionMismatchException(dFdX[i].length, yLen); + } + if (dFdY[i].length != yLen) { + throw new DimensionMismatchException(dFdY[i].length, yLen); + } + if (dFdZ[i].length != yLen) { + throw new DimensionMismatchException(dFdZ[i].length, yLen); + } + if (d2FdXdY[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdY[i].length, yLen); + } + if (d2FdXdZ[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdZ[i].length, yLen); + } + if (d2FdYdZ[i].length != yLen) { + throw new DimensionMismatchException(d2FdYdZ[i].length, yLen); + } + if (d3FdXdYdZ[i].length != yLen) { + throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen); + } + + final int ip1 = i + 1; + final double xR = xval[ip1] - xval[i]; + for (int j = 0; j < lastJ; j++) { + if (f[i][j].length != zLen) { + throw new DimensionMismatchException(f[i][j].length, zLen); + } + if (dFdX[i][j].length != zLen) { + throw new DimensionMismatchException(dFdX[i][j].length, zLen); + } + if (dFdY[i][j].length != zLen) { + throw new DimensionMismatchException(dFdY[i][j].length, zLen); + } + if (dFdZ[i][j].length != zLen) { + throw new DimensionMismatchException(dFdZ[i][j].length, zLen); + } + if (d2FdXdY[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen); + } + if (d2FdXdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen); + } + if (d2FdYdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen); + } + if (d3FdXdYdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen); + } + + final int jp1 = j + 1; + final double yR = yval[jp1] - yval[j]; + final double xRyR = xR * yR; + for (int k = 0; k < lastK; k++) { + final int kp1 = k + 1; + final double zR = zval[kp1] - zval[k]; + final double xRzR = xR * zR; + final double yRzR = yR * zR; + final double xRyRzR = xR * yRzR; + + final double[] beta = new double[] { + f[i][j][k], f[ip1][j][k], + f[i][jp1][k], f[ip1][jp1][k], + f[i][j][kp1], f[ip1][j][kp1], + f[i][jp1][kp1], f[ip1][jp1][kp1], + + dFdX[i][j][k] * xR, dFdX[ip1][j][k] * xR, + dFdX[i][jp1][k] * xR, dFdX[ip1][jp1][k] * xR, + dFdX[i][j][kp1] * xR, dFdX[ip1][j][kp1] * xR, + dFdX[i][jp1][kp1] * xR, dFdX[ip1][jp1][kp1] * xR, + + dFdY[i][j][k] * yR, dFdY[ip1][j][k] * yR, + dFdY[i][jp1][k] * yR, dFdY[ip1][jp1][k] * yR, + dFdY[i][j][kp1] * yR, dFdY[ip1][j][kp1] * yR, + dFdY[i][jp1][kp1] * yR, dFdY[ip1][jp1][kp1] * yR, + + dFdZ[i][j][k] * zR, dFdZ[ip1][j][k] * zR, + dFdZ[i][jp1][k] * zR, dFdZ[ip1][jp1][k] * zR, + dFdZ[i][j][kp1] * zR, dFdZ[ip1][j][kp1] * zR, + dFdZ[i][jp1][kp1] * zR, dFdZ[ip1][jp1][kp1] * zR, + + d2FdXdY[i][j][k] * xRyR, d2FdXdY[ip1][j][k] * xRyR, + d2FdXdY[i][jp1][k] * xRyR, d2FdXdY[ip1][jp1][k] * xRyR, + d2FdXdY[i][j][kp1] * xRyR, d2FdXdY[ip1][j][kp1] * xRyR, + d2FdXdY[i][jp1][kp1] * xRyR, d2FdXdY[ip1][jp1][kp1] * xRyR, + + d2FdXdZ[i][j][k] * xRzR, d2FdXdZ[ip1][j][k] * xRzR, + d2FdXdZ[i][jp1][k] * xRzR, d2FdXdZ[ip1][jp1][k] * xRzR, + d2FdXdZ[i][j][kp1] * xRzR, d2FdXdZ[ip1][j][kp1] * xRzR, + d2FdXdZ[i][jp1][kp1] * xRzR, d2FdXdZ[ip1][jp1][kp1] * xRzR, + + d2FdYdZ[i][j][k] * yRzR, d2FdYdZ[ip1][j][k] * yRzR, + d2FdYdZ[i][jp1][k] * yRzR, d2FdYdZ[ip1][jp1][k] * yRzR, + d2FdYdZ[i][j][kp1] * yRzR, d2FdYdZ[ip1][j][kp1] * yRzR, + d2FdYdZ[i][jp1][kp1] * yRzR, d2FdYdZ[ip1][jp1][kp1] * yRzR, + + d3FdXdYdZ[i][j][k] * xRyRzR, d3FdXdYdZ[ip1][j][k] * xRyRzR, + d3FdXdYdZ[i][jp1][k] * xRyRzR, d3FdXdYdZ[ip1][jp1][k] * xRyRzR, + d3FdXdYdZ[i][j][kp1] * xRyRzR, d3FdXdYdZ[ip1][j][kp1] * xRyRzR, + d3FdXdYdZ[i][jp1][kp1] * xRyRzR, d3FdXdYdZ[ip1][jp1][kp1] * xRyRzR, + }; + + splines[i][j][k] = new TricubicFunction(computeCoefficients(beta)); + } + } + } + } + + /** + * {@inheritDoc} + * + * @throws OutOfRangeException if any of the variables is outside its interpolation range. + */ + public double value(double x, double y, double z) + throws OutOfRangeException { + final int i = searchIndex(x, xval); + if (i == -1) { + throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); + } + final int j = searchIndex(y, yval); + if (j == -1) { + throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); + } + final int k = searchIndex(z, zval); + if (k == -1) { + throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); + } + + final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); + final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); + final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); + + return splines[i][j][k].value(xN, yN, zN); + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @param z Third coordinate. + * @return {@code true} if (x, y, z) is a valid point. + */ + public boolean isValidPoint(double x, double y, double z) { + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || + y > yval[yval.length - 1] || + z < zval[0] || + z > zval[zval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @param c Coordinate. + * @param val Coordinate samples. + * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1} + * if {@code c} is out of the range defined by the end values of {@code val}. + */ + private int searchIndex(double c, double[] val) { + if (c < val[0]) { + return -1; + } + + final int max = val.length; + for (int i = 1; i < max; i++) { + if (c <= val[i]) { + return i - 1; + } + } + + return -1; + } + + /** + * Compute the spline coefficients from the list of function values and + * function partial derivatives values at the four corners of a grid + * element. They must be specified in the following order: + * <ul> + * <li>f(0,0,0)</li> + * <li>f(1,0,0)</li> + * <li>f(0,1,0)</li> + * <li>f(1,1,0)</li> + * <li>f(0,0,1)</li> + * <li>f(1,0,1)</li> + * <li>f(0,1,1)</li> + * <li>f(1,1,1)</li> + * + * <li>f<sub>x</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>x</sub>(1,1,1)</li> + * + * <li>f<sub>y</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>y</sub>(1,1,1)</li> + * + * <li>f<sub>z</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>z</sub>(1,1,1)</li> + * + * <li>f<sub>xy</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xy</sub>(1,1,1)</li> + * + * <li>f<sub>xz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xz</sub>(1,1,1)</li> + * + * <li>f<sub>yz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>yz</sub>(1,1,1)</li> + * + * <li>f<sub>xyz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xyz</sub>(1,1,1)</li> + * </ul> + * where the subscripts indicate the partial derivative with respect to + * the corresponding variable(s). + * + * @param beta List of function values and function partial derivatives values. + * @return the spline coefficients. + */ + private double[] computeCoefficients(double[] beta) { + final int sz = 64; + final double[] a = new double[sz]; + + for (int i = 0; i < sz; i++) { + double result = 0; + final double[] row = AINV[i]; + for (int j = 0; j < sz; j++) { + result += row[j] * beta[j]; + } + a[i] = result; + } + + return a; + } +} + +/** + * 3D-spline function. + * + */ +class TricubicFunction + implements TrivariateFunction { + /** Number of points. */ + private static final short N = 4; + /** Coefficients */ + private final double[][][] a = new double[N][N][N]; + + /** + * @param aV List of spline coefficients. + */ + TricubicFunction(double[] aV) { + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + a[i][j][k] = aV[i + N * (j + N * k)]; + } + } + } + } + + /** + * @param x x-coordinate of the interpolation point. + * @param y y-coordinate of the interpolation point. + * @param z z-coordinate of the interpolation point. + * @return the interpolated value. + * @throws OutOfRangeException if {@code x}, {@code y} or + * {@code z} are not in the interval {@code [0, 1]}. + */ + public double value(double x, double y, double z) + throws OutOfRangeException { + if (x < 0 || x > 1) { + throw new OutOfRangeException(x, 0, 1); + } + if (y < 0 || y > 1) { + throw new OutOfRangeException(y, 0, 1); + } + if (z < 0 || z > 1) { + throw new OutOfRangeException(z, 0, 1); + } + + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = { 1, x, x2, x3 }; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = { 1, y, y2, y3 }; + + final double z2 = z * z; + final double z3 = z2 * z; + final double[] pZ = { 1, z, z2, z3 }; + + double result = 0; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + result += a[i][j][k] * pX[i] * pY[j] * pZ[k]; + } + } + } + + return result; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java new file mode 100644 index 0000000..ec9e3bb --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java @@ -0,0 +1,143 @@ +/* + * 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.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Generates a tricubic interpolating function. + * + * @since 3.4 + */ +public class TricubicInterpolator + implements TrivariateGridInterpolator { + /** + * {@inheritDoc} + */ + public TricubicInterpolatingFunction interpolate(final double[] xval, + final double[] yval, + final double[] zval, + final double[][][] fval) + throws NoDataException, NumberIsTooSmallException, + DimensionMismatchException, NonMonotonicSequenceException { + if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { + throw new NoDataException(); + } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + MathArrays.checkOrder(zval); + + final int xLen = xval.length; + final int yLen = yval.length; + final int zLen = zval.length; + + // Approximation to the partial derivatives using finite differences. + final double[][][] dFdX = new double[xLen][yLen][zLen]; + final double[][][] dFdY = new double[xLen][yLen][zLen]; + final double[][][] dFdZ = new double[xLen][yLen][zLen]; + final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; + final double[][][] d2FdXdZ = new double[xLen][yLen][zLen]; + final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; + final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; + + for (int i = 1; i < xLen - 1; i++) { + if (yval.length != fval[i].length) { + throw new DimensionMismatchException(yval.length, fval[i].length); + } + + final int nI = i + 1; + final int pI = i - 1; + + final double nX = xval[nI]; + final double pX = xval[pI]; + + final double deltaX = nX - pX; + + for (int j = 1; j < yLen - 1; j++) { + if (zval.length != fval[i][j].length) { + throw new DimensionMismatchException(zval.length, fval[i][j].length); + } + + final int nJ = j + 1; + final int pJ = j - 1; + + final double nY = yval[nJ]; + final double pY = yval[pJ]; + + final double deltaY = nY - pY; + final double deltaXY = deltaX * deltaY; + + for (int k = 1; k < zLen - 1; k++) { + final int nK = k + 1; + final int pK = k - 1; + + final double nZ = zval[nK]; + final double pZ = zval[pK]; + + final double deltaZ = nZ - pZ; + + dFdX[i][j][k] = (fval[nI][j][k] - fval[pI][j][k]) / deltaX; + dFdY[i][j][k] = (fval[i][nJ][k] - fval[i][pJ][k]) / deltaY; + dFdZ[i][j][k] = (fval[i][j][nK] - fval[i][j][pK]) / deltaZ; + + final double deltaXZ = deltaX * deltaZ; + final double deltaYZ = deltaY * deltaZ; + + d2FdXdY[i][j][k] = (fval[nI][nJ][k] - fval[nI][pJ][k] - fval[pI][nJ][k] + fval[pI][pJ][k]) / deltaXY; + d2FdXdZ[i][j][k] = (fval[nI][j][nK] - fval[nI][j][pK] - fval[pI][j][nK] + fval[pI][j][pK]) / deltaXZ; + d2FdYdZ[i][j][k] = (fval[i][nJ][nK] - fval[i][nJ][pK] - fval[i][pJ][nK] + fval[i][pJ][pK]) / deltaYZ; + + final double deltaXYZ = deltaXY * deltaZ; + + d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - + fval[pI][nJ][nK] + fval[pI][pJ][nK] - + fval[nI][nJ][pK] + fval[nI][pJ][pK] + + fval[pI][nJ][pK] - fval[pI][pJ][pK]) / deltaXYZ; + } + } + } + + // Create the interpolating function. + return new TricubicInterpolatingFunction(xval, yval, zval, fval, + dFdX, dFdY, dFdZ, + d2FdXdY, d2FdXdZ, d2FdYdZ, + d3FdXdYdZ) { + /** {@inheritDoc} */ + @Override + public boolean isValidPoint(double x, double y, double z) { + if (x < xval[1] || + x > xval[xval.length - 2] || + y < yval[1] || + y > yval[yval.length - 2] || + z < zval[1] || + z > zval[zval.length - 2]) { + return false; + } else { + return true; + } + } + }; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java new file mode 100644 index 0000000..96aebd3 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java @@ -0,0 +1,482 @@ +/* + * 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.TrivariateFunction; +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.OutOfRangeException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Function that implements the + * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation"> + * tricubic spline interpolation</a>, as proposed in + * <blockquote> + * Tricubic interpolation in three dimensions, + * F. Lekien and J. Marsden, + * <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471 + * </blockquote> + * + * @since 2.2 + * @deprecated To be removed in 4.0 (see MATH-1166). + */ +@Deprecated +public class TricubicSplineInterpolatingFunction + implements TrivariateFunction { + /** + * Matrix to compute the spline coefficients from the function values + * and function derivatives values + */ + private static final double[][] AINV = { + { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 }, + {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 }, + { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 }, + { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 }, + { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 }, + { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 }, + { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 }, + { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 }, + { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 }, + { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 }, + { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 }, + { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 }, + { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 }, + { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 }, + { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 } + }; + + /** Samples x-coordinates */ + private final double[] xval; + /** Samples y-coordinates */ + private final double[] yval; + /** Samples z-coordinates */ + private final double[] zval; + /** Set of cubic splines pacthing the whole data grid */ + private final TricubicSplineFunction[][][] splines; + + /** + * @param x Sample values of the x-coordinate, in increasing order. + * @param y Sample values of the y-coordinate, in increasing order. + * @param z Sample values of the y-coordinate, in increasing order. + * @param f Values of the function on every grid point. + * @param dFdX Values of the partial derivative of function with respect to x on every grid point. + * @param dFdY Values of the partial derivative of function with respect to y on every grid point. + * @param dFdZ Values of the partial derivative of function with respect to z on every grid point. + * @param d2FdXdY Values of the cross partial derivative of function on every grid point. + * @param d2FdXdZ Values of the cross partial derivative of function on every grid point. + * @param d2FdYdZ Values of the cross partial derivative of function on every grid point. + * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point. + * @throws NoDataException if any of the arrays has zero length. + * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements. + * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing. + */ + public TricubicSplineInterpolatingFunction(double[] x, + double[] y, + double[] z, + double[][][] f, + double[][][] dFdX, + double[][][] dFdY, + double[][][] dFdZ, + double[][][] d2FdXdY, + double[][][] d2FdXdZ, + double[][][] d2FdYdZ, + double[][][] d3FdXdYdZ) + throws NoDataException, + DimensionMismatchException, + NonMonotonicSequenceException { + final int xLen = x.length; + final int yLen = y.length; + final int zLen = z.length; + + if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) { + throw new NoDataException(); + } + if (xLen != f.length) { + throw new DimensionMismatchException(xLen, f.length); + } + if (xLen != dFdX.length) { + throw new DimensionMismatchException(xLen, dFdX.length); + } + if (xLen != dFdY.length) { + throw new DimensionMismatchException(xLen, dFdY.length); + } + if (xLen != dFdZ.length) { + throw new DimensionMismatchException(xLen, dFdZ.length); + } + if (xLen != d2FdXdY.length) { + throw new DimensionMismatchException(xLen, d2FdXdY.length); + } + if (xLen != d2FdXdZ.length) { + throw new DimensionMismatchException(xLen, d2FdXdZ.length); + } + if (xLen != d2FdYdZ.length) { + throw new DimensionMismatchException(xLen, d2FdYdZ.length); + } + if (xLen != d3FdXdYdZ.length) { + throw new DimensionMismatchException(xLen, d3FdXdYdZ.length); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + MathArrays.checkOrder(z); + + xval = x.clone(); + yval = y.clone(); + zval = z.clone(); + + final int lastI = xLen - 1; + final int lastJ = yLen - 1; + final int lastK = zLen - 1; + splines = new TricubicSplineFunction[lastI][lastJ][lastK]; + + for (int i = 0; i < lastI; i++) { + if (f[i].length != yLen) { + throw new DimensionMismatchException(f[i].length, yLen); + } + if (dFdX[i].length != yLen) { + throw new DimensionMismatchException(dFdX[i].length, yLen); + } + if (dFdY[i].length != yLen) { + throw new DimensionMismatchException(dFdY[i].length, yLen); + } + if (dFdZ[i].length != yLen) { + throw new DimensionMismatchException(dFdZ[i].length, yLen); + } + if (d2FdXdY[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdY[i].length, yLen); + } + if (d2FdXdZ[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdZ[i].length, yLen); + } + if (d2FdYdZ[i].length != yLen) { + throw new DimensionMismatchException(d2FdYdZ[i].length, yLen); + } + if (d3FdXdYdZ[i].length != yLen) { + throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen); + } + + final int ip1 = i + 1; + for (int j = 0; j < lastJ; j++) { + if (f[i][j].length != zLen) { + throw new DimensionMismatchException(f[i][j].length, zLen); + } + if (dFdX[i][j].length != zLen) { + throw new DimensionMismatchException(dFdX[i][j].length, zLen); + } + if (dFdY[i][j].length != zLen) { + throw new DimensionMismatchException(dFdY[i][j].length, zLen); + } + if (dFdZ[i][j].length != zLen) { + throw new DimensionMismatchException(dFdZ[i][j].length, zLen); + } + if (d2FdXdY[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen); + } + if (d2FdXdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen); + } + if (d2FdYdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen); + } + if (d3FdXdYdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen); + } + + final int jp1 = j + 1; + for (int k = 0; k < lastK; k++) { + final int kp1 = k + 1; + + final double[] beta = new double[] { + f[i][j][k], f[ip1][j][k], + f[i][jp1][k], f[ip1][jp1][k], + f[i][j][kp1], f[ip1][j][kp1], + f[i][jp1][kp1], f[ip1][jp1][kp1], + + dFdX[i][j][k], dFdX[ip1][j][k], + dFdX[i][jp1][k], dFdX[ip1][jp1][k], + dFdX[i][j][kp1], dFdX[ip1][j][kp1], + dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1], + + dFdY[i][j][k], dFdY[ip1][j][k], + dFdY[i][jp1][k], dFdY[ip1][jp1][k], + dFdY[i][j][kp1], dFdY[ip1][j][kp1], + dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1], + + dFdZ[i][j][k], dFdZ[ip1][j][k], + dFdZ[i][jp1][k], dFdZ[ip1][jp1][k], + dFdZ[i][j][kp1], dFdZ[ip1][j][kp1], + dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1], + + d2FdXdY[i][j][k], d2FdXdY[ip1][j][k], + d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k], + d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1], + d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1], + + d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k], + d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k], + d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1], + d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1], + + d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k], + d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k], + d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1], + d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1], + + d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k], + d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k], + d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1], + d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1], + }; + + splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta)); + } + } + } + } + + /** + * {@inheritDoc} + * + * @throws OutOfRangeException if any of the variables is outside its interpolation range. + */ + public double value(double x, double y, double z) + throws OutOfRangeException { + final int i = searchIndex(x, xval); + if (i == -1) { + throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); + } + final int j = searchIndex(y, yval); + if (j == -1) { + throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); + } + final int k = searchIndex(z, zval); + if (k == -1) { + throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); + } + + final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); + final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); + final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); + + return splines[i][j][k].value(xN, yN, zN); + } + + /** + * @param c Coordinate. + * @param val Coordinate samples. + * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1} + * if {@code c} is out of the range defined by the end values of {@code val}. + */ + private int searchIndex(double c, double[] val) { + if (c < val[0]) { + return -1; + } + + final int max = val.length; + for (int i = 1; i < max; i++) { + if (c <= val[i]) { + return i - 1; + } + } + + return -1; + } + + /** + * Compute the spline coefficients from the list of function values and + * function partial derivatives values at the four corners of a grid + * element. They must be specified in the following order: + * <ul> + * <li>f(0,0,0)</li> + * <li>f(1,0,0)</li> + * <li>f(0,1,0)</li> + * <li>f(1,1,0)</li> + * <li>f(0,0,1)</li> + * <li>f(1,0,1)</li> + * <li>f(0,1,1)</li> + * <li>f(1,1,1)</li> + * + * <li>f<sub>x</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>x</sub>(1,1,1)</li> + * + * <li>f<sub>y</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>y</sub>(1,1,1)</li> + * + * <li>f<sub>z</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>z</sub>(1,1,1)</li> + * + * <li>f<sub>xy</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xy</sub>(1,1,1)</li> + * + * <li>f<sub>xz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xz</sub>(1,1,1)</li> + * + * <li>f<sub>yz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>yz</sub>(1,1,1)</li> + * + * <li>f<sub>xyz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xyz</sub>(1,1,1)</li> + * </ul> + * where the subscripts indicate the partial derivative with respect to + * the corresponding variable(s). + * + * @param beta List of function values and function partial derivatives values. + * @return the spline coefficients. + */ + private double[] computeSplineCoefficients(double[] beta) { + final int sz = 64; + final double[] a = new double[sz]; + + for (int i = 0; i < sz; i++) { + double result = 0; + final double[] row = AINV[i]; + for (int j = 0; j < sz; j++) { + result += row[j] * beta[j]; + } + a[i] = result; + } + + return a; + } +} + +/** + * 3D-spline function. + * + */ +class TricubicSplineFunction + implements TrivariateFunction { + /** Number of points. */ + private static final short N = 4; + /** Coefficients */ + private final double[][][] a = new double[N][N][N]; + + /** + * @param aV List of spline coefficients. + */ + TricubicSplineFunction(double[] aV) { + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + a[i][j][k] = aV[i + N * (j + N * k)]; + } + } + } + } + + /** + * @param x x-coordinate of the interpolation point. + * @param y y-coordinate of the interpolation point. + * @param z z-coordinate of the interpolation point. + * @return the interpolated value. + * @throws OutOfRangeException if {@code x}, {@code y} or + * {@code z} are not in the interval {@code [0, 1]}. + */ + public double value(double x, double y, double z) + throws OutOfRangeException { + if (x < 0 || x > 1) { + throw new OutOfRangeException(x, 0, 1); + } + if (y < 0 || y > 1) { + throw new OutOfRangeException(y, 0, 1); + } + if (z < 0 || z > 1) { + throw new OutOfRangeException(z, 0, 1); + } + + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = { 1, x, x2, x3 }; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = { 1, y, y2, y3 }; + + final double z2 = z * z; + final double z3 = z2 * z; + final double[] pZ = { 1, z, z2, z3 }; + + double result = 0; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + result += a[i][j][k] * pX[i] * pY[j] * pZ[k]; + } + } + } + + return result; + } +} 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 new file mode 100644 index 0000000..7f43e6f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java @@ -0,0 +1,201 @@ +/* + * 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.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Generates a tricubic interpolating function. + * + * @since 2.2 + * @deprecated To be removed in 4.0 (see MATH-1166). + */ +@Deprecated +public class TricubicSplineInterpolator + implements TrivariateGridInterpolator { + /** + * {@inheritDoc} + */ + public TricubicSplineInterpolatingFunction interpolate(final double[] xval, + final double[] yval, + final double[] zval, + final double[][][] fval) + throws NoDataException, NumberIsTooSmallException, + DimensionMismatchException, NonMonotonicSequenceException { + if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { + throw new NoDataException(); + } + if (xval.length != fval.length) { + throw new DimensionMismatchException(xval.length, fval.length); + } + + MathArrays.checkOrder(xval); + MathArrays.checkOrder(yval); + MathArrays.checkOrder(zval); + + final int xLen = xval.length; + final int yLen = yval.length; + final int zLen = zval.length; + + // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets + // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k]) + // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k]) + final double[][][] fvalXY = new double[zLen][xLen][yLen]; + final double[][][] fvalZX = new double[yLen][zLen][xLen]; + for (int i = 0; i < xLen; i++) { + if (fval[i].length != yLen) { + throw new DimensionMismatchException(fval[i].length, yLen); + } + + for (int j = 0; j < yLen; j++) { + if (fval[i][j].length != zLen) { + throw new DimensionMismatchException(fval[i][j].length, zLen); + } + + for (int k = 0; k < zLen; k++) { + final double v = fval[i][j][k]; + fvalXY[k][i][j] = v; + fvalZX[j][k][i] = v; + } + } + } + + final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true); + + // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z + final BicubicSplineInterpolatingFunction[] xSplineYZ + = new BicubicSplineInterpolatingFunction[xLen]; + for (int i = 0; i < xLen; i++) { + xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]); + } + + // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x + final BicubicSplineInterpolatingFunction[] ySplineZX + = new BicubicSplineInterpolatingFunction[yLen]; + for (int j = 0; j < yLen; j++) { + ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]); + } + + // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y + final BicubicSplineInterpolatingFunction[] zSplineXY + = new BicubicSplineInterpolatingFunction[zLen]; + for (int k = 0; k < zLen; k++) { + zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]); + } + + // Partial derivatives wrt x and wrt y + final double[][][] dFdX = new double[xLen][yLen][zLen]; + final double[][][] dFdY = new double[xLen][yLen][zLen]; + final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; + for (int k = 0; k < zLen; k++) { + final BicubicSplineInterpolatingFunction f = zSplineXY[k]; + for (int i = 0; i < xLen; i++) { + final double x = xval[i]; + for (int j = 0; j < yLen; j++) { + final double y = yval[j]; + dFdX[i][j][k] = f.partialDerivativeX(x, y); + dFdY[i][j][k] = f.partialDerivativeY(x, y); + d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y); + } + } + } + + // Partial derivatives wrt y and wrt z + final double[][][] dFdZ = new double[xLen][yLen][zLen]; + final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; + for (int i = 0; i < xLen; i++) { + final BicubicSplineInterpolatingFunction f = xSplineYZ[i]; + for (int j = 0; j < yLen; j++) { + final double y = yval[j]; + for (int k = 0; k < zLen; k++) { + final double z = zval[k]; + dFdZ[i][j][k] = f.partialDerivativeY(y, z); + d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z); + } + } + } + + // Partial derivatives wrt x and wrt z + final double[][][] d2FdZdX = new double[xLen][yLen][zLen]; + for (int j = 0; j < yLen; j++) { + final BicubicSplineInterpolatingFunction f = ySplineZX[j]; + for (int k = 0; k < zLen; k++) { + final double z = zval[k]; + for (int i = 0; i < xLen; i++) { + final double x = xval[i]; + d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x); + } + } + } + + // Third partial cross-derivatives + final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; + for (int i = 0; i < xLen ; i++) { + final int nI = nextIndex(i, xLen); + final int pI = previousIndex(i); + for (int j = 0; j < yLen; j++) { + final int nJ = nextIndex(j, yLen); + final int pJ = previousIndex(j); + for (int k = 0; k < zLen; k++) { + final int nK = nextIndex(k, zLen); + final int pK = previousIndex(k); + + // XXX Not sure about this formula + d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - + fval[pI][nJ][nK] + fval[pI][pJ][nK] - + fval[nI][nJ][pK] + fval[nI][pJ][pK] + + fval[pI][nJ][pK] - fval[pI][pJ][pK]) / + ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ; + } + } + } + + // Create the interpolating splines + return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, + dFdX, dFdY, dFdZ, + d2FdXdY, d2FdZdX, d2FdYdZ, + d3FdXdYdZ); + } + + /** + * Compute the next index of an array, clipping if necessary. + * It is assumed (but not checked) that {@code i} is larger than or equal to 0. + * + * @param i Index + * @param max Upper limit of the array + * @return the next index + */ + private int nextIndex(int i, int max) { + final int index = i + 1; + return index < max ? index : index - 1; + } + /** + * Compute the previous index of an array, clipping if necessary. + * It is assumed (but not checked) that {@code i} is smaller than the size of the array. + * + * @param i Index + * @return the previous index + */ + private int previousIndex(int i) { + final int index = i - 1; + return index >= 0 ? index : 0; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java new file mode 100644 index 0000000..ec69715 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java @@ -0,0 +1,54 @@ +/* + * 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.TrivariateFunction; +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.NumberIsTooSmallException; + +/** + * Interface representing a trivariate real interpolating function where the + * sample points must be specified on a regular grid. + * + * @since 2.2 + */ +public interface TrivariateGridInterpolator { + /** + * Compute an interpolating function for the dataset. + * + * @param xval All the x-coordinates of the interpolation points, sorted + * in increasing order. + * @param yval All the y-coordinates of the interpolation points, sorted + * in increasing order. + * @param zval All the z-coordinates of the interpolation points, sorted + * in increasing order. + * @param fval the values of the interpolation points on all the grid knots: + * {@code fval[i][j][k] = f(xval[i], yval[j], zval[k])}. + * @return a function that interpolates the data set. + * @throws NoDataException if any of the arrays has zero length. + * @throws DimensionMismatchException if the array lengths are inconsistent. + * @throws NonMonotonicSequenceException if arrays are not sorted + * @throws NumberIsTooSmallException if the number of points is too small for + * the order of the interpolation + */ + TrivariateFunction interpolate(double[] xval, double[] yval, double[] zval, + double[][][] fval) + throws NoDataException, NumberIsTooSmallException, + DimensionMismatchException, NonMonotonicSequenceException; +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java new file mode 100644 index 0000000..f7a1bd1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java @@ -0,0 +1,41 @@ +/* + * 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.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; + +/** + * Interface representing a univariate real interpolating function. + * + */ +public interface UnivariateInterpolator { + /** + * Compute an interpolating function for the dataset. + * + * @param xval Arguments for the interpolation points. + * @param yval Values for the interpolation points. + * @return a function which interpolates the dataset. + * @throws MathIllegalArgumentException + * if the arguments violate assumptions made by the interpolation + * algorithm. + * @throws DimensionMismatchException if arrays lengthes do not match + */ + UnivariateFunction interpolate(double xval[], double yval[]) + throws MathIllegalArgumentException, DimensionMismatchException; +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java new file mode 100644 index 0000000..88760f5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java @@ -0,0 +1,124 @@ +/* + * 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.util.MathUtils; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; + +/** + * Adapter for classes implementing the {@link UnivariateInterpolator} + * interface. + * The data to be interpolated is assumed to be periodic. Thus values that are + * outside of the range can be passed to the interpolation function: They will + * be wrapped into the initial range before being passed to the class that + * actually computes the interpolation. + * + */ +public class UnivariatePeriodicInterpolator + implements UnivariateInterpolator { + /** Default number of extension points of the samples array. */ + public static final int DEFAULT_EXTEND = 5; + /** Interpolator. */ + private final UnivariateInterpolator interpolator; + /** Period. */ + private final double period; + /** Number of extension points. */ + private final int extend; + + /** + * Builds an interpolator. + * + * @param interpolator Interpolator. + * @param period Period. + * @param extend Number of points to be appended at the beginning and + * end of the sample arrays in order to avoid interpolation failure at + * the (periodic) boundaries of the orginal interval. The value is the + * number of sample points which the original {@code interpolator} needs + * on each side of the interpolated point. + */ + public UnivariatePeriodicInterpolator(UnivariateInterpolator interpolator, + double period, + int extend) { + this.interpolator = interpolator; + this.period = period; + this.extend = extend; + } + + /** + * Builds an interpolator. + * Uses {@link #DEFAULT_EXTEND} as the number of extension points on each side + * of the original abscissae range. + * + * @param interpolator Interpolator. + * @param period Period. + */ + public UnivariatePeriodicInterpolator(UnivariateInterpolator interpolator, + double period) { + this(interpolator, period, DEFAULT_EXTEND); + } + + /** + * {@inheritDoc} + * + * @throws NumberIsTooSmallException if the number of extension points + * is larger than the size of {@code xval}. + */ + public UnivariateFunction interpolate(double[] xval, + double[] yval) + throws NumberIsTooSmallException, NonMonotonicSequenceException { + if (xval.length < extend) { + throw new NumberIsTooSmallException(xval.length, extend, true); + } + + MathArrays.checkOrder(xval); + final double offset = xval[0]; + + final int len = xval.length + extend * 2; + final double[] x = new double[len]; + final double[] y = new double[len]; + for (int i = 0; i < xval.length; i++) { + final int index = i + extend; + x[index] = MathUtils.reduce(xval[i], period, offset); + y[index] = yval[i]; + } + + // Wrap to enable interpolation at the boundaries. + for (int i = 0; i < extend; i++) { + int index = xval.length - extend + i; + x[i] = MathUtils.reduce(xval[index], period, offset) - period; + y[i] = yval[index]; + + index = len - extend + i; + x[index] = MathUtils.reduce(xval[i], period, offset) + period; + y[index] = yval[i]; + } + + MathArrays.sortInPlace(x, y); + + final UnivariateFunction f = interpolator.interpolate(x, y); + return new UnivariateFunction() { + /** {@inheritDoc} */ + public double value(final double x) throws MathIllegalArgumentException { + return f.value(MathUtils.reduce(x, period, offset)); + } + }; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/package-info.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/package-info.java new file mode 100644 index 0000000..b4b25dd --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/package-info.java @@ -0,0 +1,22 @@ +/* + * 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. + */ +/** + * + * Univariate real functions interpolation algorithms. + * + */ +package org.apache.commons.math3.analysis.interpolation; |