http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java deleted file mode 100644 index 522c2a5..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolatingFunction.java +++ /dev/null @@ -1,637 +0,0 @@ -/* - * 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. - */ - public BicubicSplineFunction(double[] coeff) { - this(coeff, false); - } - - /** - * Simple constructor. - * - * @param coeff Spline coefficients. - * @param initializeDerivatives Whether to initialize the internal data - * needed for calling any of the methods that compute the partial derivatives - * this function. - */ - public BicubicSplineFunction(double[] coeff, - boolean initializeDerivatives) { - a = new double[N][N]; - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - a[i][j] = coeff[i * N + j]; - } - } - - if (initializeDerivatives) { - // Compute all partial derivatives functions. - final double[][] aX = new double[N][N]; - final double[][] aY = new double[N][N]; - final double[][] aXX = new double[N][N]; - final double[][] aYY = new double[N][N]; - final double[][] aXY = new double[N][N]; - - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - final double c = a[i][j]; - aX[i][j] = i * c; - aY[i][j] = j * c; - aXX[i][j] = (i - 1) * aX[i][j]; - aYY[i][j] = (j - 1) * aY[i][j]; - aXY[i][j] = j * aX[i][j]; - } - } - - partialDerivativeX = new BivariateFunction() { - public double value(double x, double y) { - final double x2 = x * x; - final double[] pX = {0, 1, x, x2}; - - final double y2 = y * y; - final double y3 = y2 * y; - final double[] pY = {1, y, y2, y3}; - - return apply(pX, pY, aX); - } - }; - partialDerivativeY = new BivariateFunction() { - public double value(double x, double y) { - final double x2 = x * x; - final double x3 = x2 * x; - final double[] pX = {1, x, x2, x3}; - - final double y2 = y * y; - final double[] pY = {0, 1, y, y2}; - - return apply(pX, pY, aY); - } - }; - partialDerivativeXX = new BivariateFunction() { - public double value(double x, double y) { - final double[] pX = {0, 0, 1, x}; - - final double y2 = y * y; - final double y3 = y2 * y; - final double[] pY = {1, y, y2, y3}; - - return apply(pX, pY, aXX); - } - }; - partialDerivativeYY = new BivariateFunction() { - public double value(double x, double y) { - final double x2 = x * x; - final double x3 = x2 * x; - final double[] pX = {1, x, x2, x3}; - - final double[] pY = {0, 0, 1, y}; - - return apply(pX, pY, aYY); - } - }; - partialDerivativeXY = new BivariateFunction() { - public double value(double x, double y) { - final double x2 = x * x; - final double[] pX = {0, 1, x, x2}; - - final double y2 = y * y; - final double[] pY = {0, 1, y, y2}; - - return apply(pX, pY, aXY); - } - }; - } else { - partialDerivativeX = null; - partialDerivativeY = null; - partialDerivativeXX = null; - partialDerivativeYY = null; - partialDerivativeXY = null; - } - } - - /** - * {@inheritDoc} - */ - public double value(double x, double y) { - if (x < 0 || x > 1) { - throw new OutOfRangeException(x, 0, 1); - } - if (y < 0 || y > 1) { - throw new OutOfRangeException(y, 0, 1); - } - - final double x2 = x * x; - final double x3 = x2 * x; - final double[] pX = {1, x, x2, x3}; - - final double y2 = y * y; - final double y3 = y2 * y; - final double[] pY = {1, y, y2, y3}; - - return apply(pX, pY, a); - } - - /** - * Compute the value of the bicubic polynomial. - * - * @param pX Powers of the x-coordinate. - * @param pY Powers of the y-coordinate. - * @param coeff Spline coefficients. - * @return the interpolated value. - */ - private double apply(double[] pX, double[] pY, double[][] coeff) { - double result = 0; - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - result += coeff[i][j] * pX[i] * pY[j]; - } - } - - return result; - } - - /** - * @return the partial derivative wrt {@code x}. - */ - public BivariateFunction partialDerivativeX() { - return partialDerivativeX; - } - /** - * @return the partial derivative wrt {@code y}. - */ - public BivariateFunction partialDerivativeY() { - return partialDerivativeY; - } - /** - * @return the second partial derivative wrt {@code x}. - */ - public BivariateFunction partialDerivativeXX() { - return partialDerivativeXX; - } - /** - * @return the second partial derivative wrt {@code y}. - */ - public BivariateFunction partialDerivativeYY() { - return partialDerivativeYY; - } - /** - * @return the second partial cross-derivative. - */ - public BivariateFunction partialDerivativeXY() { - return partialDerivativeXY; - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java deleted file mode 100644 index 09acd07..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicSplineInterpolator.java +++ /dev/null @@ -1,176 +0,0 @@ -/* - * 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; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/BivariateGridInterpolator.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 94d75ad..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BivariateGridInterpolator.java +++ /dev/null @@ -1,51 +0,0 @@ -/* - * 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; -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/DividedDifferenceInterpolator.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index e308160..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/DividedDifferenceInterpolator.java +++ /dev/null @@ -1,119 +0,0 @@ -/* - * 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> - * <p> - * The computational complexity is O(N^2).</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; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 9125b83..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java +++ /dev/null @@ -1,209 +0,0 @@ -/* - * 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; - - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 15ed322..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java +++ /dev/null @@ -1,239 +0,0 @@ -/* - * 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); - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/LinearInterpolator.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 7e0e69b..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/LinearInterpolator.java +++ /dev/null @@ -1,79 +0,0 @@ -/* - * 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); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/LoessInterpolator.java ---------------------------------------------------------------------- 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 deleted file mode 100644 index 7f0788d..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/LoessInterpolator.java +++ /dev/null @@ -1,473 +0,0 @@ -/* - * 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.NotPositiveException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NotFiniteNumberException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathUtils; -import org.apache.commons.math3.util.MathArrays; - -/** - * 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://www.math.tau.ac.il/~yekutiel/MA seminar/Cleveland 1979.pdf"> - * 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. - */ - 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. - */ - 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.</br> - * 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.</br> - * 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.</br> - * 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.</br> - * 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]); - } - } -}
