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

Remove deprecated interpolation and fitter classes.

This commit is contained in:
Thomas Neidhart 2015-02-25 23:02:30 +01:00
parent d389e94bee
commit 0a5cd11327
20 changed files with 0 additions and 5703 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

638
src/main/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatingFunction.java
View File

@ -1,638 +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.math4.analysis.interpolation;
import java.util.Arrays;
import org.apache.commons.math4.analysis.BivariateFunction;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NoDataException;
import org.apache.commons.math4.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.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.math4.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;
}
}

176
src/main/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolator.java
View File

@ -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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.UnivariateFunction;
import org.apache.commons.math4.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NoDataException;
import org.apache.commons.math4.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.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.math4.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;
}
}

171
src/main/java/org/apache/commons/math4/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java
View File

@ -1,171 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NoDataException;
import org.apache.commons.math4.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.exception.NotPositiveException;
import org.apache.commons.math4.exception.NullArgumentException;
import org.apache.commons.math4.fitting.PolynomialFitter;
import org.apache.commons.math4.optim.SimpleVectorValueChecker;
import org.apache.commons.math4.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer;
import org.apache.commons.math4.util.MathArrays;
import org.apache.commons.math4.util.Precision;
/**
* 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);
}
}

482
src/main/java/org/apache/commons/math4/analysis/interpolation/TricubicSplineInterpolatingFunction.java
View File

@ -1,482 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.TrivariateFunction;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NoDataException;
import org.apache.commons.math4.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.util.MathArrays;
/**
* Function that implements the
* <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
* tricubic spline interpolation</a>, as proposed in
* <quote>
* Tricubic interpolation in three dimensions<br/>
* F. Lekien and J. Marsden<br/>
* <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
* </quote>
*
* @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.
*/
public 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;
}
}

201
src/main/java/org/apache/commons/math4/analysis/interpolation/TricubicSplineInterpolator.java
View File

@ -1,201 +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.math4.analysis.interpolation;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.NoDataException;
import org.apache.commons.math4.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.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;
}
}

92
src/main/java/org/apache/commons/math4/analysis/solvers/NewtonSolver.java
View File

@ -1,92 +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.math4.analysis.solvers;
import org.apache.commons.math4.analysis.DifferentiableUnivariateFunction;
import org.apache.commons.math4.exception.TooManyEvaluationsException;
import org.apache.commons.math4.util.FastMath;
/**
* Implements <a href="http://mathworld.wolfram.com/NewtonsMethod.html">
* Newton's Method</a> for finding zeros of real univariate functions.
* <p>
* The function should be continuous but not necessarily smooth.</p>
*
* @deprecated as of 3.1, replaced by {@link NewtonRaphsonSolver}
*/
@Deprecated
public class NewtonSolver extends AbstractDifferentiableUnivariateSolver {
/** Default absolute accuracy. */
private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver.
*/
public NewtonSolver() {
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public NewtonSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Find a zero near the midpoint of {@code min} and {@code max}.
*
* @param f Function to solve.
* @param min Lower bound for the interval.
* @param max Upper bound for the interval.
* @param maxEval Maximum number of evaluations.
* @return the value where the function is zero.
* @throws org.apache.commons.math4.exception.TooManyEvaluationsException
* if the maximum evaluation count is exceeded.
* @throws org.apache.commons.math4.exception.NumberIsTooLargeException
* if {@code min >= max}.
*/
@Override
public double solve(int maxEval, final DifferentiableUnivariateFunction f,
final double min, final double max)
throws TooManyEvaluationsException {
return super.solve(maxEval, f, UnivariateSolverUtils.midpoint(min, max));
}
/**
* {@inheritDoc}
*/
@Override
protected double doSolve()
throws TooManyEvaluationsException {
final double startValue = getStartValue();
final double absoluteAccuracy = getAbsoluteAccuracy();
double x0 = startValue;
double x1;
while (true) {
x1 = x0 - (computeObjectiveValue(x0) / computeDerivativeObjectiveValue(x0));
if (FastMath.abs(x1 - x0) <= absoluteAccuracy) {
return x1;
}
x0 = x1;
}
}
}

233
src/main/java/org/apache/commons/math4/fitting/CurveFitter.java
View File

@ -1,233 +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.math4.fitting;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math4.analysis.MultivariateMatrixFunction;
import org.apache.commons.math4.analysis.MultivariateVectorFunction;
import org.apache.commons.math4.analysis.ParametricUnivariateFunction;
import org.apache.commons.math4.optim.InitialGuess;
import org.apache.commons.math4.optim.MaxEval;
import org.apache.commons.math4.optim.PointVectorValuePair;
import org.apache.commons.math4.optim.nonlinear.vector.ModelFunction;
import org.apache.commons.math4.optim.nonlinear.vector.ModelFunctionJacobian;
import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
import org.apache.commons.math4.optim.nonlinear.vector.Target;
import org.apache.commons.math4.optim.nonlinear.vector.Weight;
/**
* Fitter for parametric univariate real functions y = f(x).
* <br/>
* When a univariate real function y = f(x) does depend on some
* unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
* this class can be used to find these parameters. It does this
* by <em>fitting</em> the curve so it remains very close to a set of
* observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
* y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
* is done by finding the parameters values that minimizes the objective
* function &sum;(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
* really a least squares problem.
*
* @param <T> Function to use for the fit.
*
* @since 2.0
* @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
* {@link WeightedObservedPoints} instead.
*/
@Deprecated
public class CurveFitter<T extends ParametricUnivariateFunction> {
/** Optimizer to use for the fitting. */
private final MultivariateVectorOptimizer optimizer;
/** Observed points. */
private final List<WeightedObservedPoint> observations;
/**
* Simple constructor.
*
* @param optimizer Optimizer to use for the fitting.
* @since 3.1
*/
public CurveFitter(final MultivariateVectorOptimizer optimizer) {
this.optimizer = optimizer;
observations = new ArrayList<WeightedObservedPoint>();
}
/** Add an observed (x,y) point to the sample with unit weight.
* <p>Calling this method is equivalent to call
* {@code addObservedPoint(1.0, x, y)}.</p>
* @param x abscissa of the point
* @param y observed value of the point at x, after fitting we should
* have f(x) as close as possible to this value
* @see #addObservedPoint(double, double, double)
* @see #addObservedPoint(WeightedObservedPoint)
* @see #getObservations()
*/
public void addObservedPoint(double x, double y) {
addObservedPoint(1.0, x, y);
}
/** Add an observed weighted (x,y) point to the sample.
* @param weight weight of the observed point in the fit
* @param x abscissa of the point
* @param y observed value of the point at x, after fitting we should
* have f(x) as close as possible to this value
* @see #addObservedPoint(double, double)
* @see #addObservedPoint(WeightedObservedPoint)
* @see #getObservations()
*/
public void addObservedPoint(double weight, double x, double y) {
observations.add(new WeightedObservedPoint(weight, x, y));
}
/** Add an observed weighted (x,y) point to the sample.
* @param observed observed point to add
* @see #addObservedPoint(double, double)
* @see #addObservedPoint(double, double, double)
* @see #getObservations()
*/
public void addObservedPoint(WeightedObservedPoint observed) {
observations.add(observed);
}
/** Get the observed points.
* @return observed points
* @see #addObservedPoint(double, double)
* @see #addObservedPoint(double, double, double)
* @see #addObservedPoint(WeightedObservedPoint)
*/
public WeightedObservedPoint[] getObservations() {
return observations.toArray(new WeightedObservedPoint[observations.size()]);
}
/**
* Remove all observations.
*/
public void clearObservations() {
observations.clear();
}
/**
* Fit a curve.
* This method compute the coefficients of the curve that best
* fit the sample of observed points previously given through calls
* to the {@link #addObservedPoint(WeightedObservedPoint)
* addObservedPoint} method.
*
* @param f parametric function to fit.
* @param initialGuess first guess of the function parameters.
* @return the fitted parameters.
* @throws org.apache.commons.math4.exception.DimensionMismatchException
* if the start point dimension is wrong.
*/
public double[] fit(T f, final double[] initialGuess) {
return fit(Integer.MAX_VALUE, f, initialGuess);
}
/**
* Fit a curve.
* This method compute the coefficients of the curve that best
* fit the sample of observed points previously given through calls
* to the {@link #addObservedPoint(WeightedObservedPoint)
* addObservedPoint} method.
*
* @param f parametric function to fit.
* @param initialGuess first guess of the function parameters.
* @param maxEval Maximum number of function evaluations.
* @return the fitted parameters.
* @throws org.apache.commons.math4.exception.TooManyEvaluationsException
* if the number of allowed evaluations is exceeded.
* @throws org.apache.commons.math4.exception.DimensionMismatchException
* if the start point dimension is wrong.
* @since 3.0
*/
public double[] fit(int maxEval, T f,
final double[] initialGuess) {
// Prepare least squares problem.
double[] target = new double[observations.size()];
double[] weights = new double[observations.size()];
int i = 0;
for (WeightedObservedPoint point : observations) {
target[i] = point.getY();
weights[i] = point.getWeight();
++i;
}
// Input to the optimizer: the model and its Jacobian.
final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);
// Perform the fit.
final PointVectorValuePair optimum
= optimizer.optimize(new MaxEval(maxEval),
model.getModelFunction(),
model.getModelFunctionJacobian(),
new Target(target),
new Weight(weights),
new InitialGuess(initialGuess));
// Extract the coefficients.
return optimum.getPointRef();
}
/** Vectorial function computing function theoretical values. */
private class TheoreticalValuesFunction {
/** Function to fit. */
private final ParametricUnivariateFunction f;
/**
* @param f function to fit.
*/
public TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
this.f = f;
}
/**
* @return the model function values.
*/
public ModelFunction getModelFunction() {
return new ModelFunction(new MultivariateVectorFunction() {
/** {@inheritDoc} */
public double[] value(double[] point) {
// compute the residuals
final double[] values = new double[observations.size()];
int i = 0;
for (WeightedObservedPoint observed : observations) {
values[i++] = f.value(observed.getX(), point);
}
return values;
}
});
}
/**
* @return the model function Jacobian.
*/
public ModelFunctionJacobian getModelFunctionJacobian() {
return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
public double[][] value(double[] point) {
final double[][] jacobian = new double[observations.size()][];
int i = 0;
for (WeightedObservedPoint observed : observations) {
jacobian[i++] = f.gradient(observed.getX(), point);
}
return jacobian;
}
});
}
}
}

365
src/main/java/org/apache/commons/math4/fitting/GaussianFitter.java
View File

@ -1,365 +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.math4.fitting;
import java.util.Arrays;
import java.util.Comparator;
import org.apache.commons.math4.analysis.function.Gaussian;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.NullArgumentException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.exception.ZeroException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
import org.apache.commons.math4.util.FastMath;
/**
* Fits points to a {@link
* org.apache.commons.math4.analysis.function.Gaussian.Parametric Gaussian} function.
* <p>
* Usage example:
* <pre>
* GaussianFitter fitter = new GaussianFitter(
* new LevenbergMarquardtOptimizer());
* fitter.addObservedPoint(4.0254623, 531026.0);
* fitter.addObservedPoint(4.03128248, 984167.0);
* fitter.addObservedPoint(4.03839603, 1887233.0);
* fitter.addObservedPoint(4.04421621, 2687152.0);
* fitter.addObservedPoint(4.05132976, 3461228.0);
* fitter.addObservedPoint(4.05326982, 3580526.0);
* fitter.addObservedPoint(4.05779662, 3439750.0);
* fitter.addObservedPoint(4.0636168, 2877648.0);
* fitter.addObservedPoint(4.06943698, 2175960.0);
* fitter.addObservedPoint(4.07525716, 1447024.0);
* fitter.addObservedPoint(4.08237071, 717104.0);
* fitter.addObservedPoint(4.08366408, 620014.0);
* double[] parameters = fitter.fit();
* </pre>
*
* @since 2.2
* @deprecated As of 3.3. Please use {@link GaussianCurveFitter} and
* {@link WeightedObservedPoints} instead.
*/
@Deprecated
public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
/**
* Constructs an instance using the specified optimizer.
*
* @param optimizer Optimizer to use for the fitting.
*/
public GaussianFitter(MultivariateVectorOptimizer optimizer) {
super(optimizer);
}
/**
* Fits a Gaussian function to the observed points.
*
* @param initialGuess First guess values in the following order:
* <ul>
* <li>Norm</li>
* <li>Mean</li>
* <li>Sigma</li>
* </ul>
* @return the parameters of the Gaussian function that best fits the
* observed points (in the same order as above).
* @since 3.0
*/
public double[] fit(double[] initialGuess) {
final Gaussian.Parametric f = new Gaussian.Parametric() {
@Override
public double value(double x, double ... p) {
double v = Double.POSITIVE_INFINITY;
try {
v = super.value(x, p);
} catch (NotStrictlyPositiveException e) { // NOPMD
// Do nothing.
}
return v;
}
@Override
public double[] gradient(double x, double ... p) {
double[] v = { Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY };
try {
v = super.gradient(x, p);
} catch (NotStrictlyPositiveException e) { // NOPMD
// Do nothing.
}
return v;
}
};
return fit(f, initialGuess);
}
/**
* Fits a Gaussian function to the observed points.
*
* @return the parameters of the Gaussian function that best fits the
* observed points (in the same order as above).
*/
public double[] fit() {
final double[] guess = (new ParameterGuesser(getObservations())).guess();
return fit(guess);
}
/**
* Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
* of a {@link org.apache.commons.math4.analysis.function.Gaussian.Parametric}
* based on the specified observed points.
*/
public static class ParameterGuesser {
/** Normalization factor. */
private final double norm;
/** Mean. */
private final double mean;
/** Standard deviation. */
private final double sigma;
/**
* Constructs instance with the specified observed points.
*
* @param observations Observed points from which to guess the
* parameters of the Gaussian.
* @throws NullArgumentException if {@code observations} is
* {@code null}.
* @throws NumberIsTooSmallException if there are less than 3
* observations.
*/
public ParameterGuesser(WeightedObservedPoint[] observations) {
if (observations == null) {
throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
}
if (observations.length < 3) {
throw new NumberIsTooSmallException(observations.length, 3, true);
}
final WeightedObservedPoint[] sorted = sortObservations(observations);
final double[] params = basicGuess(sorted);
norm = params[0];
mean = params[1];
sigma = params[2];
}
/**
* Gets an estimation of the parameters.
*
* @return the guessed parameters, in the following order:
* <ul>
* <li>Normalization factor</li>
* <li>Mean</li>
* <li>Standard deviation</li>
* </ul>
*/
public double[] guess() {
return new double[] { norm, mean, sigma };
}
/**
* Sort the observations.
*
* @param unsorted Input observations.
* @return the input observations, sorted.
*/
private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
final WeightedObservedPoint[] observations = unsorted.clone();
final Comparator<WeightedObservedPoint> cmp
= new Comparator<WeightedObservedPoint>() {
public int compare(WeightedObservedPoint p1,
WeightedObservedPoint p2) {
if (p1 == null && p2 == null) {
return 0;
}
if (p1 == null) {
return -1;
}
if (p2 == null) {
return 1;
}
if (p1.getX() < p2.getX()) {
return -1;
}
if (p1.getX() > p2.getX()) {
return 1;
}
if (p1.getY() < p2.getY()) {
return -1;
}
if (p1.getY() > p2.getY()) {
return 1;
}
if (p1.getWeight() < p2.getWeight()) {
return -1;
}
if (p1.getWeight() > p2.getWeight()) {
return 1;
}
return 0;
}
};
Arrays.sort(observations, cmp);
return observations;
}
/**
* Guesses the parameters based on the specified observed points.
*
* @param points Observed points, sorted.
* @return the guessed parameters (normalization factor, mean and
* sigma).
*/
private double[] basicGuess(WeightedObservedPoint[] points) {
final int maxYIdx = findMaxY(points);
final double n = points[maxYIdx].getY();
final double m = points[maxYIdx].getX();
double fwhmApprox;
try {
final double halfY = n + ((m - n) / 2);
final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
fwhmApprox = fwhmX2 - fwhmX1;
} catch (OutOfRangeException e) {
// TODO: Exceptions should not be used for flow control.
fwhmApprox = points[points.length - 1].getX() - points[0].getX();
}
final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
return new double[] { n, m, s };
}
/**
* Finds index of point in specified points with the largest Y.
*
* @param points Points to search.
* @return the index in specified points array.
*/
private int findMaxY(WeightedObservedPoint[] points) {
int maxYIdx = 0;
for (int i = 1; i < points.length; i++) {
if (points[i].getY() > points[maxYIdx].getY()) {
maxYIdx = i;
}
}
return maxYIdx;
}
/**
* Interpolates using the specified points to determine X at the
* specified Y.
*
* @param points Points to use for interpolation.
* @param startIdx Index within points from which to start the search for
* interpolation bounds points.
* @param idxStep Index step for searching interpolation bounds points.
* @param y Y value for which X should be determined.
* @return the value of X for the specified Y.
* @throws ZeroException if {@code idxStep} is 0.
* @throws OutOfRangeException if specified {@code y} is not within the
* range of the specified {@code points}.
*/
private double interpolateXAtY(WeightedObservedPoint[] points,
int startIdx,
int idxStep,
double y)
throws OutOfRangeException {
if (idxStep == 0) {
throw new ZeroException();
}
final WeightedObservedPoint[] twoPoints
= getInterpolationPointsForY(points, startIdx, idxStep, y);
final WeightedObservedPoint p1 = twoPoints[0];
final WeightedObservedPoint p2 = twoPoints[1];
if (p1.getY() == y) {
return p1.getX();
}
if (p2.getY() == y) {
return p2.getX();
}
return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
(p2.getY() - p1.getY()));
}
/**
* Gets the two bounding interpolation points from the specified points
* suitable for determining X at the specified Y.
*
* @param points Points to use for interpolation.
* @param startIdx Index within points from which to start search for
* interpolation bounds points.
* @param idxStep Index step for search for interpolation bounds points.
* @param y Y value for which X should be determined.
* @return the array containing two points suitable for determining X at
* the specified Y.
* @throws ZeroException if {@code idxStep} is 0.
* @throws OutOfRangeException if specified {@code y} is not within the
* range of the specified {@code points}.
*/
private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
int startIdx,
int idxStep,
double y)
throws OutOfRangeException {
if (idxStep == 0) {
throw new ZeroException();
}
for (int i = startIdx;
idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
i += idxStep) {
final WeightedObservedPoint p1 = points[i];
final WeightedObservedPoint p2 = points[i + idxStep];
if (isBetween(y, p1.getY(), p2.getY())) {
if (idxStep < 0) {
return new WeightedObservedPoint[] { p2, p1 };
} else {
return new WeightedObservedPoint[] { p1, p2 };
}
}
}
// Boundaries are replaced by dummy values because the raised
// exception is caught and the message never displayed.
// TODO: Exceptions should not be used for flow control.
throw new OutOfRangeException(y,
Double.NEGATIVE_INFINITY,
Double.POSITIVE_INFINITY);
}
/**
* Determines whether a value is between two other values.
*
* @param value Value to test whether it is between {@code boundary1}
* and {@code boundary2}.
* @param boundary1 One end of the range.
* @param boundary2 Other end of the range.
* @return {@code true} if {@code value} is between {@code boundary1} and
* {@code boundary2} (inclusive), {@code false} otherwise.
*/
private boolean isBetween(double value,
double boundary1,
double boundary2) {
return (value >= boundary1 && value <= boundary2) ||
(value >= boundary2 && value <= boundary1);
}
}
}

384
src/main/java/org/apache/commons/math4/fitting/HarmonicFitter.java
View File

@ -1,384 +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.math4.fitting;
import org.apache.commons.math4.analysis.function.HarmonicOscillator;
import org.apache.commons.math4.exception.MathIllegalStateException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.exception.ZeroException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
import org.apache.commons.math4.util.FastMath;
/**
* Class that implements a curve fitting specialized for sinusoids.
*
* Harmonic fitting is a very simple case of curve fitting. The
* estimated coefficients are the amplitude a, the pulsation &omega; and
* the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
* searched by a least square estimator initialized with a rough guess
* based on integrals.
*
* @since 2.0
* @deprecated As of 3.3. Please use {@link HarmonicCurveFitter} and
* {@link WeightedObservedPoints} instead.
*/
@Deprecated
public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> {
/**
* Simple constructor.
* @param optimizer Optimizer to use for the fitting.
*/
public HarmonicFitter(final MultivariateVectorOptimizer optimizer) {
super(optimizer);
}
/**
* Fit an harmonic function to the observed points.
*
* @param initialGuess First guess values in the following order:
* <ul>
* <li>Amplitude</li>
* <li>Angular frequency</li>
* <li>Phase</li>
* </ul>
* @return the parameters of the harmonic function that best fits the
* observed points (in the same order as above).
*/
public double[] fit(double[] initialGuess) {
return fit(new HarmonicOscillator.Parametric(), initialGuess);
}
/**
* Fit an harmonic function to the observed points.
* An initial guess will be automatically computed.
*
* @return the parameters of the harmonic function that best fits the
* observed points (see the other {@link #fit(double[]) fit} method.
* @throws NumberIsTooSmallException if the sample is too short for the
* the first guess to be computed.
* @throws ZeroException if the first guess cannot be computed because
* the abscissa range is zero.
*/
public double[] fit() {
return fit((new ParameterGuesser(getObservations())).guess());
}
/**
* This class guesses harmonic coefficients from a sample.
* <p>The algorithm used to guess the coefficients is as follows:</p>
*
* <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
* &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
* </p>
*
* <p>From the analytical expression, we can compute two primitives :
* <pre>
* If2 (t) = &int; f<sup>2</sup> = a<sup>2</sup> &times; [t + S (t)] / 2
* If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
* where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
* </pre>
* </p>
*
* <p>We can remove S between these expressions :
* <pre>
* If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
* </pre>
* </p>
*
* <p>The preceding expression shows that If'2 (t) is a linear
* combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
* </p>
*
* <p>From the primitive, we can deduce the same form for definite
* integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
* <pre>
* If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A &times; (t<sub>i</sub> - t<sub>1</sub>) + B &times; (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
* </pre>
* </p>
*
* <p>We can find the coefficients A and B that best fit the sample
* to this linear expression by computing the definite integrals for
* each sample points.
* </p>
*
* <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A &times; x<sub>i</sub> + B &times; y<sub>i</sub>, the
* coefficients A and B that minimize a least square criterion
* &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
* <pre>
*
* &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* A = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub>
* B = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
* </pre>
* </p>
*
*
* <p>In fact, we can assume both a and &omega; are positive and
* compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
* B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
* <pre>
*
* for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
* f (t<sub>i</sub>)
* f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>)
* x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
* y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* update the sums &sum;x<sub>i</sub>x<sub>i</sub>, &sum;y<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>z<sub>i</sub> and &sum;y<sub>i</sub>z<sub>i</sub>
* end for
*
* |--------------------------
* \ | &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* a = \ | ------------------------
* \| &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
*
*
* |--------------------------
* \ | &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* &omega; = \ | ------------------------
* \| &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* </pre>
* </p>
*
* <p>Once we know &omega;, we can compute:
* <pre>
* fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
* fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
* </pre>
* </p>
*
* <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
* <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
* expressions to compute &phi;. The best estimate over the sample is
* given by averaging these expressions.
* </p>
*
* <p>Since integrals and means are involved in the preceding
* estimations, these operations run in O(n) time, where n is the
* number of measurements.</p>
*/
public static class ParameterGuesser {
/** Amplitude. */
private final double a;
/** Angular frequency. */
private final double omega;
/** Phase. */
private final double phi;
/**
* Simple constructor.
*
* @param observations Sampled observations.
* @throws NumberIsTooSmallException if the sample is too short.
* @throws ZeroException if the abscissa range is zero.
* @throws MathIllegalStateException when the guessing procedure cannot
* produce sensible results.
*/
public ParameterGuesser(WeightedObservedPoint[] observations) {
if (observations.length < 4) {
throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
observations.length, 4, true);
}
final WeightedObservedPoint[] sorted = sortObservations(observations);
final double aOmega[] = guessAOmega(sorted);
a = aOmega[0];
omega = aOmega[1];
phi = guessPhi(sorted);
}
/**
* Gets an estimation of the parameters.
*
* @return the guessed parameters, in the following order:
* <ul>
* <li>Amplitude</li>
* <li>Angular frequency</li>
* <li>Phase</li>
* </ul>
*/
public double[] guess() {
return new double[] { a, omega, phi };
}
/**
* Sort the observations with respect to the abscissa.
*
* @param unsorted Input observations.
* @return the input observations, sorted.
*/
private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
final WeightedObservedPoint[] observations = unsorted.clone();
// Since the samples are almost always already sorted, this
// method is implemented as an insertion sort that reorders the
// elements in place. Insertion sort is very efficient in this case.
WeightedObservedPoint curr = observations[0];
for (int j = 1; j < observations.length; ++j) {
WeightedObservedPoint prec = curr;
curr = observations[j];
if (curr.getX() < prec.getX()) {
// the current element should be inserted closer to the beginning
int i = j - 1;
WeightedObservedPoint mI = observations[i];
while ((i >= 0) && (curr.getX() < mI.getX())) {
observations[i + 1] = mI;
if (i-- != 0) {
mI = observations[i];
}
}
observations[i + 1] = curr;
curr = observations[j];
}
}
return observations;
}
/**
* Estimate a first guess of the amplitude and angular frequency.
* This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method
* has been called previously.
*
* @param observations Observations, sorted w.r.t. abscissa.
* @throws ZeroException if the abscissa range is zero.
* @throws MathIllegalStateException when the guessing procedure cannot
* produce sensible results.
* @return the guessed amplitude (at index 0) and circular frequency
* (at index 1).
*/
private double[] guessAOmega(WeightedObservedPoint[] observations) {
final double[] aOmega = new double[2];
// initialize the sums for the linear model between the two integrals
double sx2 = 0;
double sy2 = 0;
double sxy = 0;
double sxz = 0;
double syz = 0;
double currentX = observations[0].getX();
double currentY = observations[0].getY();
double f2Integral = 0;
double fPrime2Integral = 0;
final double startX = currentX;
for (int i = 1; i < observations.length; ++i) {
// one step forward
final double previousX = currentX;
final double previousY = currentY;
currentX = observations[i].getX();
currentY = observations[i].getY();
// update the integrals of f<sup>2</sup> and f'<sup>2</sup>
// considering a linear model for f (and therefore constant f')
final double dx = currentX - previousX;
final double dy = currentY - previousY;
final double f2StepIntegral =
dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
final double fPrime2StepIntegral = dy * dy / dx;
final double x = currentX - startX;
f2Integral += f2StepIntegral;
fPrime2Integral += fPrime2StepIntegral;
sx2 += x * x;
sy2 += f2Integral * f2Integral;
sxy += x * f2Integral;
sxz += x * fPrime2Integral;
syz += f2Integral * fPrime2Integral;
}
// compute the amplitude and pulsation coefficients
double c1 = sy2 * sxz - sxy * syz;
double c2 = sxy * sxz - sx2 * syz;
double c3 = sx2 * sy2 - sxy * sxy;
if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
final int last = observations.length - 1;
// Range of the observations, assuming that the
// observations are sorted.
final double xRange = observations[last].getX() - observations[0].getX();
if (xRange == 0) {
throw new ZeroException();
}
aOmega[1] = 2 * Math.PI / xRange;
double yMin = Double.POSITIVE_INFINITY;
double yMax = Double.NEGATIVE_INFINITY;
for (int i = 1; i < observations.length; ++i) {
final double y = observations[i].getY();
if (y < yMin) {
yMin = y;
}
if (y > yMax) {
yMax = y;
}
}
aOmega[0] = 0.5 * (yMax - yMin);
} else {
if (c2 == 0) {
// In some ill-conditioned cases (cf. MATH-844), the guesser
// procedure cannot produce sensible results.
throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
}
aOmega[0] = FastMath.sqrt(c1 / c2);
aOmega[1] = FastMath.sqrt(c2 / c3);
}
return aOmega;
}
/**
* Estimate a first guess of the phase.
*
* @param observations Observations, sorted w.r.t. abscissa.
* @return the guessed phase.
*/
private double guessPhi(WeightedObservedPoint[] observations) {
// initialize the means
double fcMean = 0;
double fsMean = 0;
double currentX = observations[0].getX();
double currentY = observations[0].getY();
for (int i = 1; i < observations.length; ++i) {
// one step forward
final double previousX = currentX;
final double previousY = currentY;
currentX = observations[i].getX();
currentY = observations[i].getY();
final double currentYPrime = (currentY - previousY) / (currentX - previousX);
double omegaX = omega * currentX;
double cosine = FastMath.cos(omegaX);
double sine = FastMath.sin(omegaX);
fcMean += omega * currentY * cosine - currentYPrime * sine;
fsMean += omega * currentY * sine + currentYPrime * cosine;
}
return FastMath.atan2(-fsMean, fcMean);
}
}
}

72
src/main/java/org/apache/commons/math4/fitting/PolynomialFitter.java
View File

@ -1,72 +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.math4.fitting;
import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
/**
* Polynomial fitting is a very simple case of {@link CurveFitter curve fitting}.
* The estimated coefficients are the polynomial coefficients (see the
* {@link #fit(double[]) fit} method).
*
* @since 2.0
* @deprecated As of 3.3. Please use {@link PolynomialCurveFitter} and
* {@link WeightedObservedPoints} instead.
*/
@Deprecated
public class PolynomialFitter extends CurveFitter<PolynomialFunction.Parametric> {
/**
* Simple constructor.
*
* @param optimizer Optimizer to use for the fitting.
*/
public PolynomialFitter(MultivariateVectorOptimizer optimizer) {
super(optimizer);
}
/**
* Get the coefficients of the polynomial fitting the weighted data points.
* The degree of the fitting polynomial is {@code guess.length - 1}.
*
* @param guess First guess for the coefficients. They must be sorted in
* increasing order of the polynomial's degree.
* @param maxEval Maximum number of evaluations of the polynomial.
* @return the coefficients of the polynomial that best fits the observed points.
* @throws org.apache.commons.math4.exception.TooManyEvaluationsException if
* the number of evaluations exceeds {@code maxEval}.
* @throws org.apache.commons.math4.exception.ConvergenceException
* if the algorithm failed to converge.
*/
public double[] fit(int maxEval, double[] guess) {
return fit(maxEval, new PolynomialFunction.Parametric(), guess);
}
/**
* Get the coefficients of the polynomial fitting the weighted data points.
* The degree of the fitting polynomial is {@code guess.length - 1}.
*
* @param guess First guess for the coefficients. They must be sorted in
* increasing order of the polynomial's degree.
* @return the coefficients of the polynomial that best fits the observed points.
* @throws org.apache.commons.math4.exception.ConvergenceException
* if the algorithm failed to converge.
*/
public double[] fit(double[] guess) {
return fit(new PolynomialFunction.Parametric(), guess);
}
}

670
src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
View File

@ -1,670 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.BivariateFunction;
import org.apache.commons.math4.analysis.interpolation.BicubicSplineFunction;
import org.apache.commons.math4.analysis.interpolation.BicubicSplineInterpolatingFunction;
import org.apache.commons.math4.distribution.UniformRealDistribution;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.random.RandomGenerator;
import org.apache.commons.math4.random.Well19937c;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* Test case for the bicubic function.
*
*/
public final class BicubicSplineInterpolatingFunctionTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[][] zval = new double[xval.length][yval.length];
@SuppressWarnings("unused")
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
zval, zval, zval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
bcf = new BicubicSplineInterpolatingFunction(wxval, yval, zval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -1, -1, 2.5};
try {
bcf = new BicubicSplineInterpolatingFunction(xval, wyval, zval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[][] wzval = new double[xval.length][yval.length - 1];
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wzval = new double[xval.length - 1][yval.length];
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Test for a plane.
* <p>
* z = 2 x - 3 y + 5
*/
@Ignore@Test
public void testPlane() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x - 3 * y + 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = 2;
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = -3;
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 0;
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
double expected, result;
x = 4;
y = -3;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("On sample point",
expected, result, 1e-15);
x = 4.5;
y = -1.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.3);
x = 3.5;
y = -3.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.3);
}
/**
* Test for a paraboloid.
* <p>
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Ignore@Test
public void testParaboloid() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x * x - 3 * y * y + 4 * x * y - 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
BivariateFunction dfdX = new BivariateFunction() {
public double value(double x, double y) {
return 4 * (x + y);
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = dfdX.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
BivariateFunction dfdY = new BivariateFunction() {
public double value(double x, double y) {
return 4 * x - 6 * y;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = dfdY.value(xval[i], yval[j]);
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 4;
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
double expected, result;
x = 4;
y = -3;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("On sample point",
expected, result, 1e-15);
x = 4.5;
y = -1.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 2);
x = 3.5;
y = -3.5;
expected = f.value(x, y);
result = bcf.value(x, y);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 2);
}
/**
* Test for partial derivatives of {@link BicubicSplineFunction}.
* <p>
* f(x, y) = &Sigma;<sub>i</sub>&Sigma;<sub>j</sub> (i+1) (j+2) x<sup>i</sup> y<sup>j</sup>
*/
@Ignore@Test
public void testSplinePartialDerivatives() {
final int N = 4;
final double[] coeff = new double[16];
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
coeff[i + N * j] = (i + 1) * (j + 2);
}
}
final BicubicSplineFunction f = new BicubicSplineFunction(coeff);
BivariateFunction derivative;
final double x = 0.435;
final double y = 0.776;
final double tol = 1e-13;
derivative = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double y2 = y * y;
final double y3 = y2 * y;
final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
return yFactor * (2 + 6 * x + 12 * x2);
}
};
Assert.assertEquals("dFdX", derivative.value(x, y),
f.partialDerivativeX().value(x, y), tol);
derivative = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double x3 = x2 * x;
final double y2 = y * y;
final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
return xFactor * (3 + 8 * y + 15 * y2);
}
};
Assert.assertEquals("dFdY", derivative.value(x, y),
f.partialDerivativeY().value(x, y), tol);
derivative = new BivariateFunction() {
public double value(double x, double y) {
final double y2 = y * y;
final double y3 = y2 * y;
final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
return yFactor * (6 + 24 * x);
}
};
Assert.assertEquals("d2FdX2", derivative.value(x, y),
f.partialDerivativeXX().value(x, y), tol);
derivative = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double x3 = x2 * x;
final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
return xFactor * (8 + 30 * y);
}
};
Assert.assertEquals("d2FdY2", derivative.value(x, y),
f.partialDerivativeYY().value(x, y), tol);
derivative = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double y2 = y * y;
final double yFactor = 3 + 8 * y + 15 * y2;
return yFactor * (2 + 6 * x + 12 * x2);
}
};
Assert.assertEquals("d2FdXdY", derivative.value(x, y),
f.partialDerivativeXY().value(x, y), tol);
}
/**
* Test that the partial derivatives computed from a
* {@link BicubicSplineInterpolatingFunction} match the input data.
* <p>
* f(x, y) = 5
* - 3 x + 2 y
* - x y + 2 x<sup>2</sup> - 3 y<sup>2</sup>
* + 4 x<sup>2</sup> y - x y<sup>2</sup> - 3 x<sup>3</sup> + y<sup>3</sup>
*/
@Ignore@Test
public void testMatchingPartialDerivatives() {
final int sz = 21;
double[] val = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
val[i] = i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double x3 = x2 * x;
final double y2 = y * y;
final double y3 = y2 * y;
return 5
- 3 * x + 2 * y
- x * y + 2 * x2 - 3 * y2
+ 4 * x2 * y - x * y2 - 3 * x3 + y3;
}
};
double[][] fval = new double[sz][sz];
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
fval[i][j] = f.value(val[i], val[j]);
}
}
// Partial derivatives with respect to x
double[][] dFdX = new double[sz][sz];
BivariateFunction dfdX = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double y2 = y * y;
return - 3 - y + 4 * x + 8 * x * y - y2 - 9 * x2;
}
};
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
dFdX[i][j] = dfdX.value(val[i], val[j]);
}
}
// Partial derivatives with respect to y
double[][] dFdY = new double[sz][sz];
BivariateFunction dfdY = new BivariateFunction() {
public double value(double x, double y) {
final double x2 = x * x;
final double y2 = y * y;
return 2 - x - 6 * y + 4 * x2 - 2 * x * y + 3 * y2;
}
};
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
dFdY[i][j] = dfdY.value(val[i], val[j]);
}
}
// Partial cross-derivatives
double[][] d2FdXdY = new double[sz][sz];
BivariateFunction d2fdXdY = new BivariateFunction() {
public double value(double x, double y) {
return -1 + 8 * x - 2 * y;
}
};
for (int i = 0; i < sz; i++) {
for (int j = 0; j < sz; j++) {
d2FdXdY[i][j] = d2fdXdY.value(val[i], val[j]);
}
}
BicubicSplineInterpolatingFunction bcf
= new BicubicSplineInterpolatingFunction(val, val, fval, dFdX, dFdY, d2FdXdY);
double x, y;
double expected, result;
final double tol = 1e-12;
for (int i = 0; i < sz; i++) {
x = val[i];
for (int j = 0; j < sz; j++) {
y = val[j];
expected = dfdX.value(x, y);
result = bcf.partialDerivativeX(x, y);
Assert.assertEquals(x + " " + y + " dFdX", expected, result, tol);
expected = dfdY.value(x, y);
result = bcf.partialDerivativeY(x, y);
Assert.assertEquals(x + " " + y + " dFdY", expected, result, tol);
expected = d2fdXdY.value(x, y);
result = bcf.partialDerivativeXY(x, y);
Assert.assertEquals(x + " " + y + " d2FdXdY", expected, result, tol);
}
}
}
/**
* Interpolating a plane.
* <p>
* z = 2 x - 3 y + 5
*/
@Test
public void testInterpolation1() {
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x - 3 * y + 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = 2;
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = -3;
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 0;
}
}
final BivariateFunction bcf
= new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final int numSamples = 50;
final double tol = 6;
for (int i = 0; i < numSamples; i++) {
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol);
}
// System.out.println();
}
}
/**
* Interpolating a paraboloid.
* <p>
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Test
public void testInterpolation2() {
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x * x - 3 * y * y + 4 * x * y - 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to x
double[][] dZdX = new double[xval.length][yval.length];
BivariateFunction dfdX = new BivariateFunction() {
public double value(double x, double y) {
return 4 * (x + y);
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdX[i][j] = dfdX.value(xval[i], yval[j]);
}
}
// Partial derivatives with respect to y
double[][] dZdY = new double[xval.length][yval.length];
BivariateFunction dfdY = new BivariateFunction() {
public double value(double x, double y) {
return 4 * x - 6 * y;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdY[i][j] = dfdY.value(xval[i], yval[j]);
}
}
// Partial cross-derivatives
double[][] dZdXdY = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
dZdXdY[i][j] = 4;
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final double tol = 224;
for (int i = 0; i < sz; i++) {
x = distX.sample();
for (int j = 0; j < sz; j++) {
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol);
}
// System.out.println();
}
}
@Test
public void testIsValidPoint() {
final double xMin = -12;
final double xMax = 34;
final double yMin = 5;
final double yMax = 67;
final double[] xval = new double[] { xMin, xMax };
final double[] yval = new double[] { yMin, yMax };
final double[][] f = new double[][] { { 1, 2 },
{ 3, 4 } };
final double[][] dFdX = f;
final double[][] dFdY = f;
final double[][] dFdXdY = f;
final BicubicSplineInterpolatingFunction bcf
= new BicubicSplineInterpolatingFunction(xval, yval, f,
dFdX, dFdY, dFdXdY);
double x, y;
x = xMin;
y = yMin;
Assert.assertTrue(bcf.isValidPoint(x, y));
// Ensure that no exception is thrown.
bcf.value(x, y);
x = xMax;
y = yMax;
Assert.assertTrue(bcf.isValidPoint(x, y));
// Ensure that no exception is thrown.
bcf.value(x, y);
final double xRange = xMax - xMin;
final double yRange = yMax - yMin;
x = xMin + xRange / 3.4;
y = yMin + yRange / 1.2;
Assert.assertTrue(bcf.isValidPoint(x, y));
// Ensure that no exception is thrown.
bcf.value(x, y);
final double small = 1e-8;
x = xMin - small;
y = yMax;
Assert.assertFalse(bcf.isValidPoint(x, y));
// Ensure that an exception would have been thrown.
try {
bcf.value(x, y);
Assert.fail("OutOfRangeException expected");
} catch (OutOfRangeException expected) {}
x = xMin;
y = yMax + small;
Assert.assertFalse(bcf.isValidPoint(x, y));
// Ensure that an exception would have been thrown.
try {
bcf.value(x, y);
Assert.fail("OutOfRangeException expected");
} catch (OutOfRangeException expected) {}
}
}

186
src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatorTest.java
View File

@ -1,186 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.BivariateFunction;
import org.apache.commons.math4.analysis.interpolation.BicubicSplineInterpolator;
import org.apache.commons.math4.analysis.interpolation.BivariateGridInterpolator;
import org.apache.commons.math4.distribution.UniformRealDistribution;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
import org.apache.commons.math4.random.RandomGenerator;
import org.apache.commons.math4.random.Well19937c;
import org.junit.Assert;
import org.junit.Test;
/**
* Test case for the bicubic interpolator.
*
*/
public final class BicubicSplineInterpolatorTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[][] zval = new double[xval.length][yval.length];
BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
@SuppressWarnings("unused")
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
p = interpolator.interpolate(wxval, yval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -3, -1, -1};
try {
p = interpolator.interpolate(xval, wyval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[][] wzval = new double[xval.length][yval.length + 1];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wzval = new double[xval.length - 1][yval.length];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Interpolating a plane.
* <p>
* z = 2 x - 3 y + 5
*/
@Test
public void testInterpolation1() {
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x - 3 * y + 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final int numSamples = 50;
final double tol = 6;
for (int i = 0; i < numSamples; i++) {
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
Assert.assertEquals(f.value(x, y), p.value(x, y), tol);
}
// System.out.println();
}
}
/**
* Interpolating a paraboloid.
* <p>
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Test
public void testInterpolation2() {
final int sz = 21;
double[] xval = new double[sz];
double[] yval = new double[sz];
// Coordinate values
final double delta = 1d / (sz - 1);
for (int i = 0; i < sz; i++) {
xval[i] = -1 + 15 * i * delta;
yval[i] = -20 + 30 * i * delta;
}
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x * x - 3 * y * y + 4 * x * y - 5;
}
};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
double x, y;
final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
final UniformRealDistribution distX
= new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
final UniformRealDistribution distY
= new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
final int numSamples = 50;
final double tol = 251;
for (int i = 0; i < numSamples; i++) {
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
y = distY.sample();
// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
Assert.assertEquals(f.value(x, y), p.value(x, y), tol);
}
// System.out.println();
}
}
}

181
src/test/java/org/apache/commons/math4/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolatorTest.java
View File

@ -1,181 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.BivariateFunction;
import org.apache.commons.math4.analysis.interpolation.BivariateGridInterpolator;
import org.apache.commons.math4.analysis.interpolation.SmoothingPolynomialBicubicSplineInterpolator;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Test case for the smoothing bicubic interpolator.
*
*/
public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
double[] yval = new double[] {-4, -3, -1, 2.5, 3};
double[][] zval = new double[xval.length][yval.length];
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(0);
@SuppressWarnings("unused")
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
p = interpolator.interpolate(wxval, yval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -3, -1, -1};
try {
p = interpolator.interpolate(xval, wyval, zval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[][] wzval = new double[xval.length][yval.length + 1];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wzval = new double[xval.length - 1][yval.length];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wzval = new double[xval.length][yval.length - 1];
try {
p = interpolator.interpolate(xval, yval, wzval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Test of interpolator for a plane.
* <p>
* z = 2 x - 3 y + 5
*/
@Test
public void testPlane() {
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x - 3 * y + 5
+ ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
}
};
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);
double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5, 3.5};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
double x, y;
double expected, result;
x = 4;
y = -3;
expected = f.value(x, y);
result = p.value(x, y);
Assert.assertEquals("On sample point", expected, result, 2);
x = 4.5;
y = -1.5;
expected = f.value(x, y);
result = p.value(x, y);
Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);
x = 3.5;
y = -3.5;
expected = f.value(x, y);
result = p.value(x, y);
Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
/**
* Test of interpolator for a paraboloid.
* <p>
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Test
public void testParaboloid() {
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
return 2 * x * x - 3 * y * y + 4 * x * y - 5
+ ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
}
};
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);
double[] xval = new double[] {3, 4, 5, 6.5, 7.5, 8};
double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
double[][] zval = new double[xval.length][yval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
zval[i][j] = f.value(xval[i], yval[j]);
}
}
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
double x, y;
double expected, result;
x = 5;
y = 0.5;
expected = f.value(x, y);
result = p.value(x, y);
Assert.assertEquals("On sample point", expected, result, 2);
x = 4.5;
y = -1.5;
expected = f.value(x, y);
result = p.value(x, y);
Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);
x = 3.5;
y = -3.5;
expected = f.value(x, y);
result = p.value(x, y);
Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
}

545
src/test/java/org/apache/commons/math4/analysis/interpolation/TricubicSplineInterpolatingFunctionTest.java
View File

@ -1,545 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.TrivariateFunction;
import org.apache.commons.math4.analysis.interpolation.TricubicSplineInterpolatingFunction;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Test case for the bicubic function.
*
*/
public final class TricubicSplineInterpolatingFunctionTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
double[][][] fval = new double[xval.length][yval.length][zval.length];
@SuppressWarnings("unused")
TrivariateFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, fval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
tcf = new TricubicSplineInterpolatingFunction(wxval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -1, -1, 2.5};
try {
tcf = new TricubicSplineInterpolatingFunction(xval, wyval, zval,
fval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wzval = new double[] {-12, -8, -9, -3, 0, 2.5};
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, wzval,
fval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[][][] wfval = new double[xval.length - 1][yval.length - 1][zval.length];
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
wfval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, wfval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, wfval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, wfval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
wfval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, wfval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, wfval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length][yval.length - 1][zval.length];
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
wfval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, wfval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, wfval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, wfval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
wfval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, wfval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, wfval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length][yval.length][zval.length - 1];
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
wfval, fval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, wfval, fval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, wfval, fval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, wfval,
fval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
wfval, fval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, wfval, fval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, wfval, fval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
try {
tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, fval, fval, fval,
fval, fval, fval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Test for a plane.
* <p>
* f(x, y, z) = 2 x - 3 y - 4 z + 5
* </p>
*/
@Test
public void testPlane() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
// Function values
TrivariateFunction f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return 2 * x - 3 * y - 4 * z + 5;
}
};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to x
double[][][] dFdX = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdX[i][j][k] = 2;
}
}
}
// Partial derivatives with respect to y
double[][][] dFdY = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdY[i][j][k] = -3;
}
}
}
// Partial derivatives with respect to z
double[][][] dFdZ = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdZ[i][j][k] = -4;
}
}
}
// Partial cross-derivatives
double[][][] d2FdXdY = new double[xval.length][yval.length][zval.length];
double[][][] d2FdXdZ = new double[xval.length][yval.length][zval.length];
double[][][] d2FdYdZ = new double[xval.length][yval.length][zval.length];
double[][][] d3FdXdYdZ = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdXdY[i][j][k] = 0;
d2FdXdZ[i][j][k] = 0;
d2FdYdZ[i][j][k] = 0;
d3FdXdYdZ[i][j][k] = 0;
}
}
}
TrivariateFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, dFdX, dFdY, dFdZ,
d2FdXdY, d2FdXdZ, d2FdYdZ,
d3FdXdYdZ);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("On sample point",
expected, result, 1e-15);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.3);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.3);
}
/**
* Sine wave.
* <p>
* f(x, y, z) = a cos [&omega; z - k<sub>y</sub> x - k<sub>y</sub> y]
* </p>
* with A = 0.2, &omega; = 0.5, k<sub>x</sub> = 2, k<sub>y</sub> = 1.
*/
@Test
public void testWave() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 4};
final double a = 0.2;
final double omega = 0.5;
final double kx = 2;
final double ky = 1;
// Function values
TrivariateFunction f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.cos(omega * z - kx * x - ky * y);
}
};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to x
double[][][] dFdX = new double[xval.length][yval.length][zval.length];
TrivariateFunction dFdX_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.sin(omega * z - kx * x - ky * y) * kx;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdX[i][j][k] = dFdX_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to y
double[][][] dFdY = new double[xval.length][yval.length][zval.length];
TrivariateFunction dFdY_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.sin(omega * z - kx * x - ky * y) * ky;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdY[i][j][k] = dFdY_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial derivatives with respect to z
double[][][] dFdZ = new double[xval.length][yval.length][zval.length];
TrivariateFunction dFdZ_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return -a * FastMath.sin(omega * z - kx * x - ky * y) * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
dFdZ[i][j][k] = dFdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial second derivatives w.r.t. (x, y)
double[][][] d2FdXdY = new double[xval.length][yval.length][zval.length];
TrivariateFunction d2FdXdY_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return -a * FastMath.cos(omega * z - kx * x - ky * y) * kx * ky;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdXdY[i][j][k] = d2FdXdY_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial second derivatives w.r.t. (x, z)
double[][][] d2FdXdZ = new double[xval.length][yval.length][zval.length];
TrivariateFunction d2FdXdZ_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.cos(omega * z - kx * x - ky * y) * kx * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdXdZ[i][j][k] = d2FdXdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial second derivatives w.r.t. (y, z)
double[][][] d2FdYdZ = new double[xval.length][yval.length][zval.length];
TrivariateFunction d2FdYdZ_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.cos(omega * z - kx * x - ky * y) * ky * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d2FdYdZ[i][j][k] = d2FdYdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
// Partial third derivatives
double[][][] d3FdXdYdZ = new double[xval.length][yval.length][zval.length];
TrivariateFunction d3FdXdYdZ_f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.sin(omega * z - kx * x - ky * y) * kx * ky * omega;
}
};
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
d3FdXdYdZ[i][j][k] = d3FdXdYdZ_f.value(xval[i], yval[j], zval[k]);
}
}
}
TrivariateFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval,
fval, dFdX, dFdY, dFdZ,
d2FdXdY, d2FdXdZ, d2FdYdZ,
d3FdXdYdZ);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("On sample point",
expected, result, 1e-14);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.1);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = tcf.value(x, y, z);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.1);
}
}

214
src/test/java/org/apache/commons/math4/analysis/interpolation/TricubicSplineInterpolatorTest.java
View File

@ -1,214 +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.math4.analysis.interpolation;
import org.apache.commons.math4.analysis.TrivariateFunction;
import org.apache.commons.math4.analysis.interpolation.TricubicSplineInterpolator;
import org.apache.commons.math4.analysis.interpolation.TrivariateGridInterpolator;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* Test case for the tricubic interpolator.
*
*/
public final class TricubicSplineInterpolatorTest {
/**
* Test preconditions.
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
double[][][] fval = new double[xval.length][yval.length][zval.length];
TrivariateGridInterpolator interpolator = new TricubicSplineInterpolator();
@SuppressWarnings("unused")
TrivariateFunction p = interpolator.interpolate(xval, yval, zval, fval);
double[] wxval = new double[] {3, 2, 5, 6.5};
try {
p = interpolator.interpolate(wxval, yval, zval, fval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wyval = new double[] {-4, -3, -1, -1};
try {
p = interpolator.interpolate(xval, wyval, zval, fval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[] wzval = new double[] {-12, -8, -5.5, -3, -4, 2.5};
try {
p = interpolator.interpolate(xval, yval, wzval, fval);
Assert.fail("an exception should have been thrown");
} catch (MathIllegalArgumentException e) {
// Expected
}
double[][][] wfval = new double[xval.length][yval.length + 1][zval.length];
try {
p = interpolator.interpolate(xval, yval, zval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length - 1][yval.length][zval.length];
try {
p = interpolator.interpolate(xval, yval, zval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
wfval = new double[xval.length][yval.length][zval.length - 1];
try {
p = interpolator.interpolate(xval, yval, zval, wfval);
Assert.fail("an exception should have been thrown");
} catch (DimensionMismatchException e) {
// Expected
}
}
/**
* Test of interpolator for a plane.
* <p>
* f(x, y, z) = 2 x - 3 y - z + 5
*/
@Ignore@Test
public void testPlane() {
TrivariateFunction f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return 2 * x - 3 * y - z + 5;
}
};
TrivariateGridInterpolator interpolator = new TricubicSplineInterpolator();
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
TrivariateFunction p = interpolator.interpolate(xval, yval, zval, fval);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("On sample point", expected, result, 1e-15);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 0.3);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 0.3);
}
/**
* Test of interpolator for a sine wave.
* <p>
* <p>
* f(x, y, z) = a cos [&omega; z - k<sub>y</sub> x - k<sub>y</sub> y]
* </p>
* with A = 0.2, &omega; = 0.5, k<sub>x</sub> = 2, k<sub>y</sub> = 1.
*/
@Ignore@Test
public void testWave() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
double[] zval = new double[] {-12, -8, -5.5, -3, 0, 4};
final double a = 0.2;
final double omega = 0.5;
final double kx = 2;
final double ky = 1;
// Function values
TrivariateFunction f = new TrivariateFunction() {
public double value(double x, double y, double z) {
return a * FastMath.cos(omega * z - kx * x - ky * y);
}
};
double[][][] fval = new double[xval.length][yval.length][zval.length];
for (int i = 0; i < xval.length; i++) {
for (int j = 0; j < yval.length; j++) {
for (int k = 0; k < zval.length; k++) {
fval[i][j][k] = f.value(xval[i], yval[j], zval[k]);
}
}
}
TrivariateGridInterpolator interpolator = new TricubicSplineInterpolator();
TrivariateFunction p = interpolator.interpolate(xval, yval, zval, fval);
double x, y, z;
double expected, result;
x = 4;
y = -3;
z = 0;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("On sample point",
expected, result, 1e-12);
x = 4.5;
y = -1.5;
z = -4.25;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("Half-way between sample points (middle of the patch)",
expected, result, 0.1);
x = 3.5;
y = -3.5;
z = -10;
expected = f.value(x, y, z);
result = p.value(x, y, z);
Assert.assertEquals("Half-way between sample points (border of the patch)",
expected, result, 0.1);
}
}

111
src/test/java/org/apache/commons/math4/analysis/solvers/NewtonSolverTest.java
View File

@ -1,111 +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.math4.analysis.solvers;
import org.apache.commons.math4.analysis.DifferentiableUnivariateFunction;
import org.apache.commons.math4.analysis.QuinticFunction;
import org.apache.commons.math4.analysis.UnivariateFunction;
import org.apache.commons.math4.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math4.analysis.differentiation.UnivariateDifferentiableFunction;
import org.apache.commons.math4.analysis.function.Sin;
import org.apache.commons.math4.analysis.solvers.NewtonSolver;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* @deprecated
*/
@Deprecated
public final class NewtonSolverTest {
/**
*
*/
@Test
public void testSinZero() {
DifferentiableUnivariateFunction f = new Sin();
double result;
NewtonSolver solver = new NewtonSolver();
result = solver.solve(100, f, 3, 4);
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 1, 4);
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() > 0);
}
/**
*
*/
@Test
public void testQuinticZero() {
final UnivariateDifferentiableFunction q = new QuinticFunction();
DifferentiableUnivariateFunction f = new DifferentiableUnivariateFunction() {
public double value(double x) {
return q.value(x);
}
public UnivariateFunction derivative() {
return new UnivariateFunction() {
public double value(double x) {
return q.value(new DerivativeStructure(1, 1, 0, x)).getPartialDerivative(1);
}
};
}
};
double result;
NewtonSolver solver = new NewtonSolver();
result = solver.solve(100, f, -0.2, 0.2);
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, -0.1, 0.3);
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, -0.3, 0.45);
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.3, 0.7);
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.2, 0.6);
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.05, 0.95);
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.85, 1.25);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.8, 1.2);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.85, 1.75);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.55, 1.45);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(100, f, 0.85, 5);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
}
}

143
src/test/java/org/apache/commons/math4/fitting/CurveFitterTest.java
View File

@ -1,143 +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.math4.fitting;
import org.apache.commons.math4.analysis.ParametricUnivariateFunction;
import org.apache.commons.math4.fitting.CurveFitter;
import org.apache.commons.math4.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@Deprecated
public class CurveFitterTest {
@Test
public void testMath303() {
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
CurveFitter<ParametricUnivariateFunction> fitter = new CurveFitter<ParametricUnivariateFunction>(optimizer);
fitter.addObservedPoint(2.805d, 0.6934785852953367d);
fitter.addObservedPoint(2.74333333333333d, 0.6306772025518496d);
fitter.addObservedPoint(1.655d, 0.9474675497289684);
fitter.addObservedPoint(1.725d, 0.9013594835804194d);
ParametricUnivariateFunction sif = new SimpleInverseFunction();
double[] initialguess1 = new double[1];
initialguess1[0] = 1.0d;
Assert.assertEquals(1, fitter.fit(sif, initialguess1).length);
double[] initialguess2 = new double[2];
initialguess2[0] = 1.0d;
initialguess2[1] = .5d;
Assert.assertEquals(2, fitter.fit(sif, initialguess2).length);
}
@Test
public void testMath304() {
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
CurveFitter<ParametricUnivariateFunction> fitter = new CurveFitter<ParametricUnivariateFunction>(optimizer);
fitter.addObservedPoint(2.805d, 0.6934785852953367d);
fitter.addObservedPoint(2.74333333333333d, 0.6306772025518496d);
fitter.addObservedPoint(1.655d, 0.9474675497289684);
fitter.addObservedPoint(1.725d, 0.9013594835804194d);
ParametricUnivariateFunction sif = new SimpleInverseFunction();
double[] initialguess1 = new double[1];
initialguess1[0] = 1.0d;
Assert.assertEquals(1.6357215104109237, fitter.fit(sif, initialguess1)[0], 1.0e-14);
double[] initialguess2 = new double[1];
initialguess2[0] = 10.0d;
Assert.assertEquals(1.6357215104109237, fitter.fit(sif, initialguess1)[0], 1.0e-14);
}
@Test
public void testMath372() {
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
CurveFitter<ParametricUnivariateFunction> curveFitter = new CurveFitter<ParametricUnivariateFunction>(optimizer);
curveFitter.addObservedPoint( 15, 4443);
curveFitter.addObservedPoint( 31, 8493);
curveFitter.addObservedPoint( 62, 17586);
curveFitter.addObservedPoint(125, 30582);
curveFitter.addObservedPoint(250, 45087);
curveFitter.addObservedPoint(500, 50683);
ParametricUnivariateFunction f = new ParametricUnivariateFunction() {
public double value(double x, double ... parameters) {
double a = parameters[0];
double b = parameters[1];
double c = parameters[2];
double d = parameters[3];
return d + ((a - d) / (1 + FastMath.pow(x / c, b)));
}
public double[] gradient(double x, double ... parameters) {
double a = parameters[0];
double b = parameters[1];
double c = parameters[2];
double d = parameters[3];
double[] gradients = new double[4];
double den = 1 + FastMath.pow(x / c, b);
// derivative with respect to a
gradients[0] = 1 / den;
// derivative with respect to b
// in the reported (invalid) issue, there was a sign error here
gradients[1] = -((a - d) * FastMath.pow(x / c, b) * FastMath.log(x / c)) / (den * den);
// derivative with respect to c
gradients[2] = (b * FastMath.pow(x / c, b - 1) * (x / (c * c)) * (a - d)) / (den * den);
// derivative with respect to d
gradients[3] = 1 - (1 / den);
return gradients;
}
};
double[] initialGuess = new double[] { 1500, 0.95, 65, 35000 };
double[] estimatedParameters = curveFitter.fit(f, initialGuess);
Assert.assertEquals( 2411.00, estimatedParameters[0], 500.00);
Assert.assertEquals( 1.62, estimatedParameters[1], 0.04);
Assert.assertEquals( 111.22, estimatedParameters[2], 0.30);
Assert.assertEquals(55347.47, estimatedParameters[3], 300.00);
Assert.assertTrue(optimizer.getRMS() < 600.0);
}
private static class SimpleInverseFunction implements ParametricUnivariateFunction {
public double value(double x, double ... parameters) {
return parameters[0] / x + (parameters.length < 2 ? 0 : parameters[1]);
}
public double[] gradient(double x, double ... doubles) {
double[] gradientVector = new double[doubles.length];
gradientVector[0] = 1 / x;
if (doubles.length >= 2) {
gradientVector[1] = 1;
}
return gradientVector;
}
}
}

364
src/test/java/org/apache/commons/math4/fitting/GaussianFitterTest.java
View File

@ -1,364 +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.math4.fitting;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
import org.apache.commons.math4.fitting.GaussianFitter;
import org.apache.commons.math4.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;
import org.junit.Assert;
import org.junit.Test;
/**
* Tests {@link GaussianFitter}.
*
* @since 2.2
*/
@Deprecated
public class GaussianFitterTest {
/** Good data. */
protected static final double[][] DATASET1 = new double[][] {
{4.0254623, 531026.0},
{4.02804905, 664002.0},
{4.02934242, 787079.0},
{4.03128248, 984167.0},
{4.03386923, 1294546.0},
{4.03580929, 1560230.0},
{4.03839603, 1887233.0},
{4.0396894, 2113240.0},
{4.04162946, 2375211.0},
{4.04421621, 2687152.0},
{4.04550958, 2862644.0},
{4.04744964, 3078898.0},
{4.05003639, 3327238.0},
{4.05132976, 3461228.0},
{4.05326982, 3580526.0},
{4.05585657, 3576946.0},
{4.05779662, 3439750.0},
{4.06038337, 3220296.0},
{4.06167674, 3070073.0},
{4.0636168, 2877648.0},
{4.06620355, 2595848.0},
{4.06749692, 2390157.0},
{4.06943698, 2175960.0},
{4.07202373, 1895104.0},
{4.0733171, 1687576.0},
{4.07525716, 1447024.0},
{4.0778439, 1130879.0},
{4.07978396, 904900.0},
{4.08237071, 717104.0},
{4.08366408, 620014.0}
};
/** Poor data: right of peak not symmetric with left of peak. */
protected static final double[][] DATASET2 = new double[][] {
{-20.15, 1523.0},
{-19.65, 1566.0},
{-19.15, 1592.0},
{-18.65, 1927.0},
{-18.15, 3089.0},
{-17.65, 6068.0},
{-17.15, 14239.0},
{-16.65, 34124.0},
{-16.15, 64097.0},
{-15.65, 110352.0},
{-15.15, 164742.0},
{-14.65, 209499.0},
{-14.15, 267274.0},
{-13.65, 283290.0},
{-13.15, 275363.0},
{-12.65, 258014.0},
{-12.15, 225000.0},
{-11.65, 200000.0},
{-11.15, 190000.0},
{-10.65, 185000.0},
{-10.15, 180000.0},
{ -9.65, 179000.0},
{ -9.15, 178000.0},
{ -8.65, 177000.0},
{ -8.15, 176000.0},
{ -7.65, 175000.0},
{ -7.15, 174000.0},
{ -6.65, 173000.0},
{ -6.15, 172000.0},
{ -5.65, 171000.0},
{ -5.15, 170000.0}
};
/** Poor data: long tails. */
protected static final double[][] DATASET3 = new double[][] {
{-90.15, 1513.0},
{-80.15, 1514.0},
{-70.15, 1513.0},
{-60.15, 1514.0},
{-50.15, 1513.0},
{-40.15, 1514.0},
{-30.15, 1513.0},
{-20.15, 1523.0},
{-19.65, 1566.0},
{-19.15, 1592.0},
{-18.65, 1927.0},
{-18.15, 3089.0},
{-17.65, 6068.0},
{-17.15, 14239.0},
{-16.65, 34124.0},
{-16.15, 64097.0},
{-15.65, 110352.0},
{-15.15, 164742.0},
{-14.65, 209499.0},
{-14.15, 267274.0},
{-13.65, 283290.0},
{-13.15, 275363.0},
{-12.65, 258014.0},
{-12.15, 214073.0},
{-11.65, 182244.0},
{-11.15, 136419.0},
{-10.65, 97823.0},
{-10.15, 58930.0},
{ -9.65, 35404.0},
{ -9.15, 16120.0},
{ -8.65, 9823.0},
{ -8.15, 5064.0},
{ -7.65, 2575.0},
{ -7.15, 1642.0},
{ -6.65, 1101.0},
{ -6.15, 812.0},
{ -5.65, 690.0},
{ -5.15, 565.0},
{ 5.15, 564.0},
{ 15.15, 565.0},
{ 25.15, 564.0},
{ 35.15, 565.0},
{ 45.15, 564.0},
{ 55.15, 565.0},
{ 65.15, 564.0},
{ 75.15, 565.0}
};
/** Poor data: right of peak is missing. */
protected static final double[][] DATASET4 = new double[][] {
{-20.15, 1523.0},
{-19.65, 1566.0},
{-19.15, 1592.0},
{-18.65, 1927.0},
{-18.15, 3089.0},
{-17.65, 6068.0},
{-17.15, 14239.0},
{-16.65, 34124.0},
{-16.15, 64097.0},
{-15.65, 110352.0},
{-15.15, 164742.0},
{-14.65, 209499.0},
{-14.15, 267274.0},
{-13.65, 283290.0}
};
/** Good data, but few points. */
protected static final double[][] DATASET5 = new double[][] {
{4.0254623, 531026.0},
{4.03128248, 984167.0},
{4.03839603, 1887233.0},
{4.04421621, 2687152.0},
{4.05132976, 3461228.0},
{4.05326982, 3580526.0},
{4.05779662, 3439750.0},
{4.0636168, 2877648.0},
{4.06943698, 2175960.0},
{4.07525716, 1447024.0},
{4.08237071, 717104.0},
{4.08366408, 620014.0}
};
/**
* Basic.
*/
@Test
public void testFit01() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
addDatasetToGaussianFitter(DATASET1, fitter);
double[] parameters = fitter.fit();
Assert.assertEquals(3496978.1837704973, parameters[0], 1e-4);
Assert.assertEquals(4.054933085999146, parameters[1], 1e-4);
Assert.assertEquals(0.015039355620304326, parameters[2], 1e-4);
}
/**
* Zero points is not enough observed points.
*/
@Test(expected=MathIllegalArgumentException.class)
public void testFit02() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
fitter.fit();
}
/**
* Two points is not enough observed points.
*/
@Test(expected=MathIllegalArgumentException.class)
public void testFit03() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
addDatasetToGaussianFitter(new double[][] {
{4.0254623, 531026.0},
{4.02804905, 664002.0}},
fitter);
fitter.fit();
}
/**
* Poor data: right of peak not symmetric with left of peak.
*/
@Test
public void testFit04() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
addDatasetToGaussianFitter(DATASET2, fitter);
double[] parameters = fitter.fit();
Assert.assertEquals(233003.2967252038, parameters[0], 1e-4);
Assert.assertEquals(-10.654887521095983, parameters[1], 1e-4);
Assert.assertEquals(4.335937353196641, parameters[2], 1e-4);
}
/**
* Poor data: long tails.
*/
@Test
public void testFit05() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
addDatasetToGaussianFitter(DATASET3, fitter);
double[] parameters = fitter.fit();
Assert.assertEquals(283863.81929180305, parameters[0], 1e-4);
Assert.assertEquals(-13.29641995105174, parameters[1], 1e-4);
Assert.assertEquals(1.7297330293549908, parameters[2], 1e-4);
}
/**
* Poor data: right of peak is missing.
*/
@Test
public void testFit06() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
addDatasetToGaussianFitter(DATASET4, fitter);
double[] parameters = fitter.fit();
Assert.assertEquals(285250.66754309234, parameters[0], 1e-4);
Assert.assertEquals(-13.528375695228455, parameters[1], 1e-4);
Assert.assertEquals(1.5204344894331614, parameters[2], 1e-4);
}
/**
* Basic with smaller dataset.
*/
@Test
public void testFit07() {
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
addDatasetToGaussianFitter(DATASET5, fitter);
double[] parameters = fitter.fit();
Assert.assertEquals(3514384.729342235, parameters[0], 1e-4);
Assert.assertEquals(4.054970307455625, parameters[1], 1e-4);
Assert.assertEquals(0.015029412832160017, parameters[2], 1e-4);
}
@Test
public void testMath519() {
// The optimizer will try negative sigma values but "GaussianFitter"
// will catch the raised exceptions and return NaN values instead.
final double[] data = {
1.1143831578403364E-29,
4.95281403484594E-28,
1.1171347211930288E-26,
1.7044813962636277E-25,
1.9784716574832164E-24,
1.8630236407866774E-23,
1.4820532905097742E-22,
1.0241963854632831E-21,
6.275077366673128E-21,
3.461808994532493E-20,
1.7407124684715706E-19,
8.056687953553974E-19,
3.460193945992071E-18,
1.3883326374011525E-17,
5.233894983671116E-17,
1.8630791465263745E-16,
6.288759227922111E-16,
2.0204433920597856E-15,
6.198768938576155E-15,
1.821419346860626E-14,
5.139176445538471E-14,
1.3956427429045787E-13,
3.655705706448139E-13,
9.253753324779779E-13,
2.267636001476696E-12,
5.3880460095836855E-12,
1.2431632654852931E-11
};
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
for (int i = 0; i < data.length; i++) {
fitter.addObservedPoint(i, data[i]);
}
final double[] p = fitter.fit();
Assert.assertEquals(53.1572792, p[1], 1e-7);
Assert.assertEquals(5.75214622, p[2], 1e-8);
}
@Test
public void testMath798() {
final GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
// When the data points are not commented out below, the fit stalls.
// This is expected however, since the whole dataset hardly looks like
// a Gaussian.
// When commented out, the fit proceeds fine.
fitter.addObservedPoint(0.23, 395.0);
//fitter.addObservedPoint(0.68, 0.0);
fitter.addObservedPoint(1.14, 376.0);
//fitter.addObservedPoint(1.59, 0.0);
fitter.addObservedPoint(2.05, 163.0);
//fitter.addObservedPoint(2.50, 0.0);
fitter.addObservedPoint(2.95, 49.0);
//fitter.addObservedPoint(3.41, 0.0);
fitter.addObservedPoint(3.86, 16.0);
//fitter.addObservedPoint(4.32, 0.0);
fitter.addObservedPoint(4.77, 1.0);
final double[] p = fitter.fit();
// Values are copied from a previous run of this test.
Assert.assertEquals(420.8397296167364, p[0], 1e-12);
Assert.assertEquals(0.603770729862231, p[1], 1e-15);
Assert.assertEquals(1.0786447936766612, p[2], 1e-14);
}
/**
* Adds the specified points to specified <code>GaussianFitter</code>
* instance.
*
* @param points data points where first dimension is a point index and
* second dimension is an array of length two representing the point
* with the first value corresponding to X and the second value
* corresponding to Y
* @param fitter fitter to which the points in <code>points</code> should be
* added as observed points
*/
protected static void addDatasetToGaussianFitter(double[][] points,
GaussianFitter fitter) {
for (int i = 0; i < points.length; i++) {
fitter.addObservedPoint(points[i][0], points[i][1]);
}
}
}

187
src/test/java/org/apache/commons/math4/fitting/HarmonicFitterTest.java
View File

@ -1,187 +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.math4.fitting;
import java.util.Random;
import org.apache.commons.math4.analysis.function.HarmonicOscillator;
import org.apache.commons.math4.exception.MathIllegalStateException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.fitting.HarmonicFitter;
import org.apache.commons.math4.fitting.WeightedObservedPoint;
import org.apache.commons.math4.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathUtils;
import org.junit.Test;
import org.junit.Assert;
@Deprecated
public class HarmonicFitterTest {
@Test(expected=NumberIsTooSmallException.class)
public void testPreconditions1() {
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
fitter.fit();
}
@Test
public void testNoError() {
final double a = 0.2;
final double w = 3.4;
final double p = 4.1;
HarmonicOscillator f = new HarmonicOscillator(a, w, p);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
for (double x = 0.0; x < 1.3; x += 0.01) {
fitter.addObservedPoint(1, x, f.value(x));
}
final double[] fitted = fitter.fit();
Assert.assertEquals(a, fitted[0], 1.0e-13);
Assert.assertEquals(w, fitted[1], 1.0e-13);
Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1e-13);
HarmonicOscillator ff = new HarmonicOscillator(fitted[0], fitted[1], fitted[2]);
for (double x = -1.0; x < 1.0; x += 0.01) {
Assert.assertTrue(FastMath.abs(f.value(x) - ff.value(x)) < 1e-13);
}
}
@Test
public void test1PercentError() {
Random randomizer = new Random(64925784252l);
final double a = 0.2;
final double w = 3.4;
final double p = 4.1;
HarmonicOscillator f = new HarmonicOscillator(a, w, p);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addObservedPoint(1, x,
f.value(x) + 0.01 * randomizer.nextGaussian());
}
final double[] fitted = fitter.fit();
Assert.assertEquals(a, fitted[0], 7.6e-4);
Assert.assertEquals(w, fitted[1], 2.7e-3);
Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.3e-2);
}
@Test
public void testTinyVariationsData() {
Random randomizer = new Random(64925784252l);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addObservedPoint(1, x, 1e-7 * randomizer.nextGaussian());
}
fitter.fit();
// This test serves to cover the part of the code of "guessAOmega"
// when the algorithm using integrals fails.
}
@Test
public void testInitialGuess() {
Random randomizer = new Random(45314242l);
final double a = 0.2;
final double w = 3.4;
final double p = 4.1;
HarmonicOscillator f = new HarmonicOscillator(a, w, p);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addObservedPoint(1, x,
f.value(x) + 0.01 * randomizer.nextGaussian());
}
final double[] fitted = fitter.fit(new double[] { 0.15, 3.6, 4.5 });
Assert.assertEquals(a, fitted[0], 1.2e-3);
Assert.assertEquals(w, fitted[1], 3.3e-3);
Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.7e-2);
}
@Test
public void testUnsorted() {
Random randomizer = new Random(64925784252l);
final double a = 0.2;
final double w = 3.4;
final double p = 4.1;
HarmonicOscillator f = new HarmonicOscillator(a, w, p);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
// build a regularly spaced array of measurements
int size = 100;
double[] xTab = new double[size];
double[] yTab = new double[size];
for (int i = 0; i < size; ++i) {
xTab[i] = 0.1 * i;
yTab[i] = f.value(xTab[i]) + 0.01 * randomizer.nextGaussian();
}
// shake it
for (int i = 0; i < size; ++i) {
int i1 = randomizer.nextInt(size);
int i2 = randomizer.nextInt(size);
double xTmp = xTab[i1];
double yTmp = yTab[i1];
xTab[i1] = xTab[i2];
yTab[i1] = yTab[i2];
xTab[i2] = xTmp;
yTab[i2] = yTmp;
}
// pass it to the fitter
for (int i = 0; i < size; ++i) {
fitter.addObservedPoint(1, xTab[i], yTab[i]);
}
final double[] fitted = fitter.fit();
Assert.assertEquals(a, fitted[0], 7.6e-4);
Assert.assertEquals(w, fitted[1], 3.5e-3);
Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.5e-2);
}
@Test(expected=MathIllegalStateException.class)
public void testMath844() {
final double[] y = { 0, 1, 2, 3, 2, 1,
0, -1, -2, -3, -2, -1,
0, 1, 2, 3, 2, 1,
0, -1, -2, -3, -2, -1,
0, 1, 2, 3, 2, 1, 0 };
final int len = y.length;
final WeightedObservedPoint[] points = new WeightedObservedPoint[len];
for (int i = 0; i < len; i++) {
points[i] = new WeightedObservedPoint(1, i, y[i]);
}
// The guesser fails because the function is far from an harmonic
// function: It is a triangular periodic function with amplitude 3
// and period 12, and all sample points are taken at integer abscissae
// so function values all belong to the integer subset {-3, -2, -1, 0,
// 1, 2, 3}.
new HarmonicFitter.ParameterGuesser(points);
}
}

288
src/test/java/org/apache/commons/math4/fitting/PolynomialFitterTest.java
View File

@ -1,288 +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.math4.fitting;
import java.util.Random;
import org.apache.commons.math4.TestUtils;
import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math4.analysis.polynomials.PolynomialFunction.Parametric;
import org.apache.commons.math4.distribution.RealDistribution;
import org.apache.commons.math4.distribution.UniformRealDistribution;
import org.apache.commons.math4.exception.ConvergenceException;
import org.apache.commons.math4.exception.TooManyEvaluationsException;
import org.apache.commons.math4.fitting.CurveFitter;
import org.apache.commons.math4.fitting.PolynomialFitter;
import org.apache.commons.math4.optim.SimpleVectorValueChecker;
import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
import org.apache.commons.math4.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer;
import org.apache.commons.math4.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;
import org.apache.commons.math4.util.FastMath;
import org.junit.Test;
import org.junit.Assert;
/**
* Test for class {@link CurveFitter} where the function to fit is a
* polynomial.
*/
@Deprecated
public class PolynomialFitterTest {
@Test
public void testFit() {
final RealDistribution rng = new UniformRealDistribution(-100, 100);
rng.reseedRandomGenerator(64925784252L);
final LevenbergMarquardtOptimizer optim = new LevenbergMarquardtOptimizer();
final PolynomialFitter fitter = new PolynomialFitter(optim);
final double[] coeff = { 12.9, -3.4, 2.1 }; // 12.9 - 3.4 x + 2.1 x^2
final PolynomialFunction f = new PolynomialFunction(coeff);
// Collect data from a known polynomial.
for (int i = 0; i < 100; i++) {
final double x = rng.sample();
fitter.addObservedPoint(x, f.value(x));
}
// Start fit from initial guesses that are far from the optimal values.
final double[] best = fitter.fit(new double[] { -1e-20, 3e15, -5e25 });
TestUtils.assertEquals("best != coeff", coeff, best, 1e-12);
}
@Test
public void testNoError() {
Random randomizer = new Random(64925784252l);
for (int degree = 1; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
for (int i = 0; i <= degree; ++i) {
fitter.addObservedPoint(1.0, i, p.value(i));
}
final double[] init = new double[degree + 1];
PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = FastMath.abs(p.value(x) - fitted.value(x)) /
(1.0 + FastMath.abs(p.value(x)));
Assert.assertEquals(0.0, error, 1.0e-6);
}
}
}
@Test
public void testSmallError() {
Random randomizer = new Random(53882150042l);
double maxError = 0;
for (int degree = 0; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addObservedPoint(1.0, x,
p.value(x) + 0.1 * randomizer.nextGaussian());
}
final double[] init = new double[degree + 1];
PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = FastMath.abs(p.value(x) - fitted.value(x)) /
(1.0 + FastMath.abs(p.value(x)));
maxError = FastMath.max(maxError, error);
Assert.assertTrue(FastMath.abs(error) < 0.1);
}
}
Assert.assertTrue(maxError > 0.01);
}
@Test
public void testMath798() {
final double tol = 1e-14;
final SimpleVectorValueChecker checker = new SimpleVectorValueChecker(tol, tol);
final double[] init = new double[] { 0, 0 };
final int maxEval = 3;
final double[] lm = doMath798(new LevenbergMarquardtOptimizer(checker), maxEval, init);
final double[] gn = doMath798(new GaussNewtonOptimizer(checker), maxEval, init);
for (int i = 0; i <= 1; i++) {
Assert.assertEquals(lm[i], gn[i], tol);
}
}
/**
* This test shows that the user can set the maximum number of iterations
* to avoid running for too long.
* But in the test case, the real problem is that the tolerance is way too
* stringent.
*/
@Test(expected=TooManyEvaluationsException.class)
public void testMath798WithToleranceTooLow() {
final double tol = 1e-100;
final SimpleVectorValueChecker checker = new SimpleVectorValueChecker(tol, tol);
final double[] init = new double[] { 0, 0 };
final int maxEval = 10000; // Trying hard to fit.
@SuppressWarnings("unused")
final double[] gn = doMath798(new GaussNewtonOptimizer(checker), maxEval, init);
}
/**
* This test shows that the user can set the maximum number of iterations
* to avoid running for too long.
* Even if the real problem is that the tolerance is way too stringent, it
* is possible to get the best solution so far, i.e. a checker will return
* the point when the maximum iteration count has been reached.
*/
@Test
public void testMath798WithToleranceTooLowButNoException() {
final double tol = 1e-100;
final double[] init = new double[] { 0, 0 };
final int maxEval = 10000; // Trying hard to fit.
final SimpleVectorValueChecker checker = new SimpleVectorValueChecker(tol, tol, maxEval);
final double[] lm = doMath798(new LevenbergMarquardtOptimizer(checker), maxEval, init);
final double[] gn = doMath798(new GaussNewtonOptimizer(checker), maxEval, init);
for (int i = 0; i <= 1; i++) {
Assert.assertEquals(lm[i], gn[i], 1e-15);
}
}
/**
* @param optimizer Optimizer.
* @param maxEval Maximum number of function evaluations.
* @param init First guess.
* @return the solution found by the given optimizer.
*/
private double[] doMath798(MultivariateVectorOptimizer optimizer,
int maxEval,
double[] init) {
final CurveFitter<Parametric> fitter = new CurveFitter<Parametric>(optimizer);
fitter.addObservedPoint(-0.2, -7.12442E-13);
fitter.addObservedPoint(-0.199, -4.33397E-13);
fitter.addObservedPoint(-0.198, -2.823E-13);
fitter.addObservedPoint(-0.197, -1.40405E-13);
fitter.addObservedPoint(-0.196, -7.80821E-15);
fitter.addObservedPoint(-0.195, 6.20484E-14);
fitter.addObservedPoint(-0.194, 7.24673E-14);
fitter.addObservedPoint(-0.193, 1.47152E-13);
fitter.addObservedPoint(-0.192, 1.9629E-13);
fitter.addObservedPoint(-0.191, 2.12038E-13);
fitter.addObservedPoint(-0.19, 2.46906E-13);
fitter.addObservedPoint(-0.189, 2.77495E-13);
fitter.addObservedPoint(-0.188, 2.51281E-13);
fitter.addObservedPoint(-0.187, 2.64001E-13);
fitter.addObservedPoint(-0.186, 2.8882E-13);
fitter.addObservedPoint(-0.185, 3.13604E-13);
fitter.addObservedPoint(-0.184, 3.14248E-13);
fitter.addObservedPoint(-0.183, 3.1172E-13);
fitter.addObservedPoint(-0.182, 3.12912E-13);
fitter.addObservedPoint(-0.181, 3.06761E-13);
fitter.addObservedPoint(-0.18, 2.8559E-13);
fitter.addObservedPoint(-0.179, 2.86806E-13);
fitter.addObservedPoint(-0.178, 2.985E-13);
fitter.addObservedPoint(-0.177, 2.67148E-13);
fitter.addObservedPoint(-0.176, 2.94173E-13);
fitter.addObservedPoint(-0.175, 3.27528E-13);
fitter.addObservedPoint(-0.174, 3.33858E-13);
fitter.addObservedPoint(-0.173, 2.97511E-13);
fitter.addObservedPoint(-0.172, 2.8615E-13);
fitter.addObservedPoint(-0.171, 2.84624E-13);
final double[] coeff = fitter.fit(maxEval,
new PolynomialFunction.Parametric(),
init);
return coeff;
}
@Test
public void testRedundantSolvable() {
// Levenberg-Marquardt should handle redundant information gracefully
checkUnsolvableProblem(new LevenbergMarquardtOptimizer(), true);
}
@Test
public void testRedundantUnsolvable() {
// Gauss-Newton should not be able to solve redundant information
checkUnsolvableProblem(new GaussNewtonOptimizer(true, new SimpleVectorValueChecker(1e-15, 1e-15)), false);
}
@Test
public void testLargeSample() {
Random randomizer = new Random(0x5551480dca5b369bl);
double maxError = 0;
for (int degree = 0; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
for (int i = 0; i < 40000; ++i) {
double x = -1.0 + i / 20000.0;
fitter.addObservedPoint(1.0, x,
p.value(x) + 0.1 * randomizer.nextGaussian());
}
final double[] init = new double[degree + 1];
PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = FastMath.abs(p.value(x) - fitted.value(x)) /
(1.0 + FastMath.abs(p.value(x)));
maxError = FastMath.max(maxError, error);
Assert.assertTrue(FastMath.abs(error) < 0.01);
}
}
Assert.assertTrue(maxError > 0.001);
}
private void checkUnsolvableProblem(MultivariateVectorOptimizer optimizer,
boolean solvable) {
Random randomizer = new Random(1248788532l);
for (int degree = 0; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter = new PolynomialFitter(optimizer);
// reusing the same point over and over again does not bring
// information, the problem cannot be solved in this case for
// degrees greater than 1 (but one point is sufficient for
// degree 0)
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
}
try {
final double[] init = new double[degree + 1];
fitter.fit(init);
Assert.assertTrue(solvable || (degree == 0));
} catch(ConvergenceException e) {
Assert.assertTrue((! solvable) && (degree > 0));
}
}
}
private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
final double[] coefficients = new double[degree + 1];
for (int i = 0; i <= degree; ++i) {
coefficients[i] = randomizer.nextGaussian();
}
return new PolynomialFunction(coefficients);
}
}