This commit is contained in:
Gilles 2014-10-20 23:20:09 +02:00
commit 29b1936291
35 changed files with 1835 additions and 1397 deletions

12
pom.xml
View File

@ -38,9 +38,9 @@
</issueManagement>
<scm>
<connection>scm:svn:http://svn.apache.org/repos/asf/commons/proper/math/trunk</connection>
<developerConnection>scm:svn:https://svn.apache.org/repos/asf/commons/proper/math/trunk</developerConnection>
<url>http://svn.apache.org/viewvc/commons/proper/math/trunk</url>
<connection>scm:git:http://git-wip-us.apache.org/repos/asf/commons-math.git</connection>
<developerConnection>scm:git:https://git-wip-us.apache.org/repos/asf/commons-math.git</developerConnection>
<url>https://git-wip-us.apache.org/repos/asf?p=commons-math.git</url>
</scm>
<distributionManagement>
@ -203,6 +203,9 @@
<contributor>
<name>Ted Dunning</name>
</contributor>
<contributor>
<name>Ole Ersoy</name>
</contributor>
<contributor>
<name>Ajo Fod</name>
</contributor>
@ -212,6 +215,9 @@
<contributor>
<name>Ken Geis</name>
</contributor>
<contributor>
<name>Hank Grabowski</name>
</contributor>
<contributor>
<name>Bernhard Gr&#252;newaldt</name>
</contributor>

View File

@ -73,6 +73,18 @@ Users are encouraged to upgrade to this version as this release not
2. A few methods in the FastMath class are in fact slower that their
counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901).
">
<action dev="luc" type="fix" issue="MATH-1138" due-to="Hank Grabowski">
Fixed bicubic spline interpolator, using Akima splines.
</action>
<action dev="psteitz" type="add" issue="MATH-1154" >
Changed classes in the inference package that instantiate distributions to
pass null RandomGenerators to avoid initialization overhead for the default
generator.
</action>
<action dev="luc" type="add" issue="MATH-1156" >
Added all Java 8 StrictMath methods to FastMath, so FastMath remains compatible
with newer Java versions.
</action>
<action dev="tn" type="add" issue="MATH-1139" due-to="Alexey Volkov">
Added Gumbel, Laplace, Logistic and Nakagami distributions.
</action>

View File

@ -0,0 +1,225 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
import org.apache.commons.math3.util.Precision;
/**
* Computes a cubic spline interpolation for the data set using the Akima
* algorithm, as originally formulated by Hiroshi Akima in his 1970 paper
* "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures."
* J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609
* http://doi.acm.org/10.1145/321607.321609
* <p>
* This implementation is based on the Akima implementation in the CubicSpline
* class in the Math.NET Numerics library. The method referenced is
* CubicSpline.InterpolateAkimaSorted
* <p>
* The {@link #interpolate(double[], double[])} method returns a
* {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
* over the subintervals determined by the x values, x[0] < x[i] ... < x[n]. The
* Akima algorithm requires that n >= 5.
* </p>
* <p>
*/
public class AkimaSplineInterpolator
implements UnivariateInterpolator {
/**
* The minimum number of points that are needed to compute the function
*/
public static final int MINIMUM_NUMBER_POINTS = 5;
/**
* Default constructor. Builds an AkimaSplineInterpolator object
*/
public AkimaSplineInterpolator() {
}
/**
* Computes an interpolating function for the data set.
*
* @param xvals the arguments for the interpolation points
* @param yvals the values for the interpolation points
* @return a function which interpolates the data set
* @throws DimensionMismatchException if {@code x} and {@code y} have
* different sizes.
* @throws NonMonotonicSequenceException if {@code x} is not sorted in
* strict increasing order.
* @throws NumberIsTooSmallException if the size of {@code x} is smaller
* than 5.
*/
public PolynomialSplineFunction interpolate(double[] xvals, double[] yvals)
throws DimensionMismatchException, NumberIsTooSmallException,
NonMonotonicSequenceException {
if (xvals == null || yvals == null) {
throw new NullArgumentException();
}
if (xvals.length != yvals.length) {
throw new DimensionMismatchException(xvals.length, yvals.length);
}
if (xvals.length < MINIMUM_NUMBER_POINTS) {
throw new NumberIsTooSmallException(
LocalizedFormats.NUMBER_OF_POINTS,
xvals.length,
MINIMUM_NUMBER_POINTS, true);
}
MathArrays.checkOrder(xvals);
final int numberOfDiffAndWeightElements = xvals.length - 1;
double differences[] = new double[numberOfDiffAndWeightElements];
double weights[] = new double[numberOfDiffAndWeightElements];
for (int i = 0; i < differences.length; i++) {
differences[i] = (yvals[i + 1] - yvals[i]) /
(xvals[i + 1] - xvals[i]);
}
for (int i = 1; i < weights.length; i++) {
weights[i] = FastMath.abs(differences[i] - differences[i - 1]);
}
/* Prepare Hermite interpolation scheme */
double firstDerivatives[] = new double[xvals.length];
for (int i = 2; i < firstDerivatives.length - 2; i++) {
if (Precision.equals(weights[i - 1], 0.0) &&
Precision.equals(weights[i + 1], 0.0)) {
firstDerivatives[i] = (((xvals[i + 1] - xvals[i]) * differences[i - 1]) + ((xvals[i] - xvals[i - 1]) * differences[i])) /
(xvals[i + 1] - xvals[i - 1]);
} else {
firstDerivatives[i] = ((weights[i + 1] * differences[i - 1]) + (weights[i - 1] * differences[i])) /
(weights[i + 1] + weights[i - 1]);
}
}
firstDerivatives[0] = this.differentiateThreePoint(xvals, yvals, 0, 0,
1, 2);
firstDerivatives[1] = this.differentiateThreePoint(xvals, yvals, 1, 0,
1, 2);
firstDerivatives[xvals.length - 2] = this
.differentiateThreePoint(xvals, yvals, xvals.length - 2,
xvals.length - 3, xvals.length - 2,
xvals.length - 1);
firstDerivatives[xvals.length - 1] = this
.differentiateThreePoint(xvals, yvals, xvals.length - 1,
xvals.length - 3, xvals.length - 2,
xvals.length - 1);
return this.interpolateHermiteSorted(xvals, yvals, firstDerivatives);
}
/**
* Three point differentiation helper, modeled off of the same method in the
* Math.NET CubicSpline class. This is used by both the Apache Math and the
* Math.NET Akima Cubic Spline algorithms
*
* @param xvals x values to calculate the numerical derivative with
* @param yvals y values to calculate the numerical derivative with
* @param indexOfDifferentiation index of the elemnt we are calculating the derivative around
* @param indexOfFirstSample index of the first element to sample for the three point method
* @param indexOfSecondsample index of the second element to sample for the three point method
* @param indexOfThirdSample index of the third element to sample for the three point method
* @return the derivative
*/
private double differentiateThreePoint(double[] xvals, double[] yvals,
int indexOfDifferentiation,
int indexOfFirstSample,
int indexOfSecondsample,
int indexOfThirdSample) {
double x0 = yvals[indexOfFirstSample];
double x1 = yvals[indexOfSecondsample];
double x2 = yvals[indexOfThirdSample];
double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];
double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
double b = (x1 - x0 - a * t1 * t1) / t1;
return (2 * a * t) + b;
}
/**
* Creates a Hermite cubic spline interpolation from the set of (x,y) value
* pairs and their derivatives. This is modeled off of the
* InterpolateHermiteSorted method in the Math.NET CubicSpline class.
*
* @param xvals x values for interpolation
* @param yvals y values for interpolation
* @param firstDerivatives first derivative values of the function
* @return polynomial that fits the function
*/
private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals,
double[] yvals,
double[] firstDerivatives) {
if (xvals.length != yvals.length) {
throw new DimensionMismatchException(xvals.length, yvals.length);
}
if (xvals.length != firstDerivatives.length) {
throw new DimensionMismatchException(xvals.length,
firstDerivatives.length);
}
final int minimumLength = 2;
if (xvals.length < minimumLength) {
throw new NumberIsTooSmallException(
LocalizedFormats.NUMBER_OF_POINTS,
xvals.length, minimumLength,
true);
}
int size = xvals.length - 1;
final PolynomialFunction polynomials[] = new PolynomialFunction[size];
final double coefficients[] = new double[4];
for (int i = 0; i < polynomials.length; i++) {
double w = xvals[i + 1] - xvals[i];
double w2 = w * w;
coefficients[0] = yvals[i];
coefficients[1] = firstDerivatives[i];
coefficients[2] = (3 * (yvals[i + 1] - yvals[i]) / w - 2 *
firstDerivatives[i] - firstDerivatives[i + 1]) /
w;
coefficients[3] = (2 * (yvals[i] - yvals[i + 1]) / w +
firstDerivatives[i] + firstDerivatives[i + 1]) /
w2;
polynomials[i] = new PolynomialFunction(coefficients);
}
return new PolynomialSplineFunction(xvals, polynomials);
}
}

View File

@ -18,143 +18,76 @@ package org.apache.commons.math3.analysis.interpolation;
import java.util.Arrays;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.NoDataException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.util.MathArrays;
/**
* Function that implements the
* <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
* bicubic spline interpolation</a>.
* Function that implements the <a
* href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline
* interpolation</a>.
*
* @since 2.1
*/
public class BicubicSplineInterpolatingFunction
implements BivariateFunction {
/** Number of coefficients. */
private static final int NUM_COEFF = 16;
/**
* Matrix to compute the spline coefficients from the function values
* and function derivatives values
*/
private static final double[][] AINV = {
{ 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
{ -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
{ 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
{ 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
{ 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
{ 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
{ -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
{ 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
{ 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
{ -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
{ 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
{ 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
{ -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
{ 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
};
/** Samples x-coordinates */
private final double[] xval;
/** Samples y-coordinates */
private final double[] yval;
/** Set of cubic splines patching the whole data grid */
private final BicubicSplineFunction[][] splines;
/**
* Partial derivatives.
* The value of the first index determines the kind of derivatives:
* 0 = first partial derivatives wrt x
* 1 = first partial derivatives wrt y
* 2 = second partial derivatives wrt x
* 3 = second partial derivatives wrt y
* 4 = cross partial derivatives
*/
private final BivariateFunction[][][] partialDerivatives;
private final double[][] fval;
/**
* @param x Sample values of the x-coordinate, in increasing order.
* @param y Sample values of the y-coordinate, in increasing order.
* @param f Values of the function on every grid point.
* @param dFdX Values of the partial derivative of function with respect
* to x on every grid point.
* @param dFdY Values of the partial derivative of function with respect
* to y on every grid point.
* @param d2FdXdY Values of the cross partial derivative of function on
* every grid point.
* @throws DimensionMismatchException if the various arrays do not contain
* the expected number of elements.
* @throws NonMonotonicSequenceException if {@code x} or {@code y} are
* not strictly increasing.
* @param f Values of the function on every grid point. the expected number
* of elements.
* @throws NonMonotonicSequenceException if {@code x} or {@code y} are not
* strictly increasing.
* @throws NullArgumentException if any of the arguments are null
* @throws NoDataException if any of the arrays has zero length.
* @throws DimensionMismatchException if the length of x and y don't match the row, column
* height of f
*/
public BicubicSplineInterpolatingFunction(double[] x,
double[] y,
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY)
throws DimensionMismatchException,
NoDataException,
NonMonotonicSequenceException {
this(x, y, f, dFdX, dFdY, d2FdXdY, false);
}
/**
* @param x Sample values of the x-coordinate, in increasing order.
* @param y Sample values of the y-coordinate, in increasing order.
* @param f Values of the function on every grid point.
* @param dFdX Values of the partial derivative of function with respect
* to x on every grid point.
* @param dFdY Values of the partial derivative of function with respect
* to y on every grid point.
* @param d2FdXdY Values of the cross partial derivative of function on
* every grid point.
* @param initializeDerivatives Whether to initialize the internal data
* needed for calling any of the methods that compute the partial derivatives
* this function.
* @throws DimensionMismatchException if the various arrays do not contain
* the expected number of elements.
* @throws NonMonotonicSequenceException if {@code x} or {@code y} are
* not strictly increasing.
* @throws NoDataException if any of the arrays has zero length.
*
* @see #partialDerivativeX(double,double)
* @see #partialDerivativeY(double,double)
* @see #partialDerivativeXX(double,double)
* @see #partialDerivativeYY(double,double)
* @see #partialDerivativeXY(double,double)
*/
public BicubicSplineInterpolatingFunction(double[] x,
double[] y,
double[][] f,
double[][] dFdX,
double[][] dFdY,
double[][] d2FdXdY,
boolean initializeDerivatives)
throws DimensionMismatchException,
NoDataException,
NonMonotonicSequenceException {
public BicubicSplineInterpolatingFunction(double[] x, double[] y,
double[][] f)
throws DimensionMismatchException, NullArgumentException,
NoDataException, NonMonotonicSequenceException {
final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS;
if (x == null || y == null || f == null || f[0] == null) {
throw new NullArgumentException();
}
final int xLen = x.length;
final int yLen = y.length;
if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
throw new NoDataException();
}
if (xLen < minimumLength || yLen < minimumLength ||
f.length < minimumLength || f[0].length < minimumLength) {
throw new InsufficientDataException();
}
if (xLen != f.length) {
throw new DimensionMismatchException(xLen, f.length);
}
if (xLen != dFdX.length) {
throw new DimensionMismatchException(xLen, dFdX.length);
}
if (xLen != dFdY.length) {
throw new DimensionMismatchException(xLen, dFdY.length);
}
if (xLen != d2FdXdY.length) {
throw new DimensionMismatchException(xLen, d2FdXdY.length);
if (yLen != f[0].length) {
throw new DimensionMismatchException(yLen, f[0].length);
}
MathArrays.checkOrder(x);
@ -162,57 +95,7 @@ public class BicubicSplineInterpolatingFunction
xval = x.clone();
yval = y.clone();
final int lastI = xLen - 1;
final int lastJ = yLen - 1;
splines = new BicubicSplineFunction[lastI][lastJ];
for (int i = 0; i < lastI; i++) {
if (f[i].length != yLen) {
throw new DimensionMismatchException(f[i].length, yLen);
}
if (dFdX[i].length != yLen) {
throw new DimensionMismatchException(dFdX[i].length, yLen);
}
if (dFdY[i].length != yLen) {
throw new DimensionMismatchException(dFdY[i].length, yLen);
}
if (d2FdXdY[i].length != yLen) {
throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
}
final int ip1 = i + 1;
for (int j = 0; j < lastJ; j++) {
final int jp1 = j + 1;
final double[] beta = new double[] {
f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
};
splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta),
initializeDerivatives);
}
}
if (initializeDerivatives) {
// Compute all partial derivatives.
partialDerivatives = new BivariateFunction[5][lastI][lastJ];
for (int i = 0; i < lastI; i++) {
for (int j = 0; j < lastJ; j++) {
final BicubicSplineFunction bcs = splines[i][j];
partialDerivatives[0][i][j] = bcs.partialDerivativeX();
partialDerivatives[1][i][j] = bcs.partialDerivativeY();
partialDerivatives[2][i][j] = bcs.partialDerivativeXX();
partialDerivatives[3][i][j] = bcs.partialDerivativeYY();
partialDerivatives[4][i][j] = bcs.partialDerivativeXY();
}
}
} else {
// Partial derivative methods cannot be used.
partialDerivatives = null;
}
fval = f.clone();
}
/**
@ -220,13 +103,37 @@ public class BicubicSplineInterpolatingFunction
*/
public double value(double x, double y)
throws OutOfRangeException {
final int i = searchIndex(x, xval);
final int j = searchIndex(y, yval);
int index;
PolynomialSplineFunction spline;
AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
final int offset = 2;
final int count = offset + 3;
final int i = searchIndex(x, xval, offset, count);
final int j = searchIndex(y, yval, offset, count);
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
double xArray[] = new double[count];
double yArray[] = new double[count];
double zArray[] = new double[count];
double interpArray[] = new double[count];
return splines[i][j].value(xN, yN);
for (index = 0; index < count; index++) {
xArray[index] = xval[i + index];
yArray[index] = yval[j + index];
}
for (int zIndex = 0; zIndex < count; zIndex++) {
for (index = 0; index < count; index++) {
zArray[index] = fval[i + index][j + zIndex];
}
spline = interpolator.interpolate(xArray, zArray);
interpArray[zIndex] = spline.value(x);
}
spline = interpolator.interpolate(yArray, interpArray);
double returnValue = spline.value(y);
return returnValue;
}
/**
@ -238,9 +145,7 @@ public class BicubicSplineInterpolatingFunction
* @since 3.3
*/
public boolean isValidPoint(double x, double y) {
if (x < xval[0] ||
x > xval[xval.length - 1] ||
y < yval[0] ||
if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] ||
y > yval[yval.length - 1]) {
return false;
} else {
@ -248,386 +153,43 @@ public class BicubicSplineInterpolatingFunction
}
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the first partial derivative with
* respect to x.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeX(double x, double y)
throws OutOfRangeException {
return partialDerivative(0, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the first partial derivative with
* respect to y.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeY(double x, double y)
throws OutOfRangeException {
return partialDerivative(1, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the second partial derivative with
* respect to x.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeXX(double x, double y)
throws OutOfRangeException {
return partialDerivative(2, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the second partial derivative with
* respect to y.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeYY(double x, double y)
throws OutOfRangeException {
return partialDerivative(3, x, y);
}
/**
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the second partial cross-derivative.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
public double partialDerivativeXY(double x, double y)
throws OutOfRangeException {
return partialDerivative(4, x, y);
}
/**
* @param which First index in {@link #partialDerivatives}.
* @param x x-coordinate.
* @param y y-coordinate.
* @return the value at point (x, y) of the selected partial derivative.
* @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
* the range defined by the boundary values of {@code xval} (resp.
* {@code yval}).
* @throws NullPointerException if the internal data were not initialized
* (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
* double[][],double[][],double[][],boolean) constructor}).
*/
private double partialDerivative(int which, double x, double y)
throws OutOfRangeException {
final int i = searchIndex(x, xval);
final int j = searchIndex(y, yval);
final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
return partialDerivatives[which][i][j].value(xN, yN);
}
/**
* @param c Coordinate.
* @param val Coordinate samples.
* @return the index in {@code val} corresponding to the interval
* containing {@code c}.
* @throws OutOfRangeException if {@code c} is out of the
* range defined by the boundary values of {@code val}.
* @param offset how far back from found value to offset for querying
* @param count total number of elements forward from beginning that will be
* queried
* @return the index in {@code val} corresponding to the interval containing
* {@code c}.
* @throws OutOfRangeException if {@code c} is out of the range defined by
* the boundary values of {@code val}.
*/
private int searchIndex(double c, double[] val) {
final int r = Arrays.binarySearch(val, c);
private int searchIndex(double c, double[] val, int offset, int count) {
int r = Arrays.binarySearch(val, c);
if (r == -1 ||
r == -val.length - 1) {
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;
// "c" in within an interpolation sub-interval, which returns
// negative
// need to remove the negative sign for consistency
r = -r - offset - 1;
} else {
r -= offset;
}
final int last = val.length - 1;
if (r == last) {
if (r < 0) {
r = 0;
}
if ((r + count) >= val.length) {
// "c" is the last sample of the range: Return the index
// of the sample at the lower end of the last sub-interval.
return last - 1;
r = val.length - count;
}
// "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;
}
}

View File

@ -16,12 +16,10 @@
*/
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NoDataException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.util.MathArrays;
/**
@ -30,144 +28,37 @@ import org.apache.commons.math3.util.MathArrays;
* @since 2.2
*/
public class BicubicSplineInterpolator
implements BivariateGridInterpolator {
/** Whether to initialize internal data used to compute the analytical
derivatives of the splines. */
private final boolean initializeDerivatives;
implements BivariateGridInterpolator
{
/**
* Default constructor.
* The argument {@link #BicubicSplineInterpolator(boolean) initializeDerivatives}
* is set to {@code false}.
*/
public BicubicSplineInterpolator() {
this(false);
}
public BicubicSplineInterpolator()
{
/**
* Creates an interpolator.
*
* @param initializeDerivatives Whether to initialize the internal data
* needed for calling any of the methods that compute the partial derivatives
* of the {@link BicubicSplineInterpolatingFunction function} returned from
* the call to {@link #interpolate(double[],double[],double[][]) interpolate}.
*/
public BicubicSplineInterpolator(boolean initializeDerivatives) {
this.initializeDerivatives = initializeDerivatives;
}
/**
* {@inheritDoc}
*/
public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
final double[] yval,
public BicubicSplineInterpolatingFunction interpolate( final double[] xval, final double[] yval,
final double[][] fval)
throws NoDataException, DimensionMismatchException,
NonMonotonicSequenceException, NumberIsTooSmallException {
if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
throw new NoDataException();
throws DimensionMismatchException, NullArgumentException, NoDataException, NonMonotonicSequenceException
{
if ( xval == null || yval == null || fval == null || fval[0] == null )
{
throw new NullArgumentException();
}
if (xval.length != fval.length) {
throw new DimensionMismatchException(xval.length, fval.length);
if ( xval.length == 0 || yval.length == 0 || fval.length == 0 )
{
throw new NoDataException();
}
MathArrays.checkOrder(xval);
MathArrays.checkOrder(yval);
final int xLen = xval.length;
final int yLen = yval.length;
// Samples (first index is y-coordinate, i.e. subarray variable is x)
// 0 <= i < xval.length
// 0 <= j < yval.length
// fX[j][i] = f(xval[i], yval[j])
final double[][] fX = new double[yLen][xLen];
for (int i = 0; i < xLen; i++) {
if (fval[i].length != yLen) {
throw new DimensionMismatchException(fval[i].length, yLen);
}
for (int j = 0; j < yLen; j++) {
fX[j][i] = fval[i][j];
}
}
final SplineInterpolator spInterpolator = new SplineInterpolator();
// For each line y[j] (0 <= j < yLen), construct a 1D spline with
// respect to variable x
final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
for (int j = 0; j < yLen; j++) {
ySplineX[j] = spInterpolator.interpolate(xval, fX[j]);
}
// For each line x[i] (0 <= i < xLen), construct a 1D spline with
// respect to variable y generated by array fY_1[i]
final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
for (int i = 0; i < xLen; i++) {
xSplineY[i] = spInterpolator.interpolate(yval, fval[i]);
}
// Partial derivatives with respect to x at the grid knots
final double[][] dFdX = new double[xLen][yLen];
for (int j = 0; j < yLen; j++) {
final UnivariateFunction f = ySplineX[j].derivative();
for (int i = 0; i < xLen; i++) {
dFdX[i][j] = f.value(xval[i]);
}
}
// Partial derivatives with respect to y at the grid knots
final double[][] dFdY = new double[xLen][yLen];
for (int i = 0; i < xLen; i++) {
final UnivariateFunction f = xSplineY[i].derivative();
for (int j = 0; j < yLen; j++) {
dFdY[i][j] = f.value(yval[j]);
}
}
// Cross partial derivatives
final double[][] d2FdXdY = new double[xLen][yLen];
for (int i = 0; i < xLen ; i++) {
final int nI = nextIndex(i, xLen);
final int pI = previousIndex(i);
for (int j = 0; j < yLen; j++) {
final int nJ = nextIndex(j, yLen);
final int pJ = previousIndex(j);
d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] -
fval[pI][nJ] + fval[pI][pJ]) /
((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
}
}
// Create the interpolating splines
return new BicubicSplineInterpolatingFunction(xval, yval, fval,
dFdX, dFdY, d2FdXdY,
initializeDerivatives);
}
/**
* Computes the next index of an array, clipping if necessary.
* It is assumed (but not checked) that {@code i >= 0}.
*
* @param i Index.
* @param max Upper limit of the array.
* @return the next index.
*/
private int nextIndex(int i, int max) {
final int index = i + 1;
return index < max ? index : index - 1;
}
/**
* Computes the previous index of an array, clipping if necessary.
* It is assumed (but not checked) that {@code i} is smaller than the size
* of the array.
*
* @param i Index.
* @return the previous index.
*/
private int previousIndex(int i) {
final int index = i - 1;
return index >= 0 ? index : 0;
return new BicubicSplineInterpolatingFunction( xval, yval, fval );
}
}

View File

@ -76,7 +76,7 @@ public class TricubicSplineInterpolator
}
}
final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true);
final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();
// For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
final BicubicSplineInterpolatingFunction[] xSplineYZ
@ -103,46 +103,13 @@ public class TricubicSplineInterpolator
final double[][][] dFdX = new double[xLen][yLen][zLen];
final double[][][] dFdY = new double[xLen][yLen][zLen];
final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
for (int k = 0; k < zLen; k++) {
final BicubicSplineInterpolatingFunction f = zSplineXY[k];
for (int i = 0; i < xLen; i++) {
final double x = xval[i];
for (int j = 0; j < yLen; j++) {
final double y = yval[j];
dFdX[i][j][k] = f.partialDerivativeX(x, y);
dFdY[i][j][k] = f.partialDerivativeY(x, y);
d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y);
}
}
}
// Partial derivatives wrt y and wrt z
final double[][][] dFdZ = new double[xLen][yLen][zLen];
final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
for (int i = 0; i < xLen; i++) {
final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
for (int j = 0; j < yLen; j++) {
final double y = yval[j];
for (int k = 0; k < zLen; k++) {
final double z = zval[k];
dFdZ[i][j][k] = f.partialDerivativeY(y, z);
d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z);
}
}
}
// Partial derivatives wrt x and wrt z
final double[][][] d2FdZdX = new double[xLen][yLen][zLen];
for (int j = 0; j < yLen; j++) {
final BicubicSplineInterpolatingFunction f = ySplineZX[j];
for (int k = 0; k < zLen; k++) {
final double z = zval[k];
for (int i = 0; i < xLen; i++) {
final double x = xval[i];
d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
}
}
}
// Third partial cross-derivatives
final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];

View File

@ -22,8 +22,9 @@ import org.apache.commons.math3.exception.util.ExceptionContext;
import org.apache.commons.math3.exception.util.ExceptionContextProvider;
/**
* Base class for all exceptions that signal a mismatch between the
* current state and the user's expectations.
* Base class for all exceptions that signal that the process
* throwing the exception is in a state that does not comply with
* the set of states that it is designed to be in.
*
* @since 2.2
*/

View File

@ -293,6 +293,7 @@ public enum LocalizedFormats implements Localizable {
OVERFLOW_IN_FRACTION("overflow in fraction {0}/{1}, cannot negate"),
OVERFLOW_IN_ADDITION("overflow in addition: {0} + {1}"),
OVERFLOW_IN_SUBTRACTION("overflow in subtraction: {0} - {1}"),
OVERFLOW_IN_MULTIPLICATION("overflow in multiplication: {0} * {1}"),
PERCENTILE_IMPLEMENTATION_CANNOT_ACCESS_METHOD("cannot access {0} method in percentile implementation {1}"),
PERCENTILE_IMPLEMENTATION_UNSUPPORTED_METHOD("percentile implementation {0} does not support {1}"),
PERMUTATION_EXCEEDS_N("permutation size ({0}) exceeds permuation domain ({1})"), /* keep */

View File

@ -119,7 +119,19 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
// set up integration control objects
stepStart = equations.getTime();
stepSize = forward ? step : -step;
if (forward) {
if (stepStart + step >= t) {
stepSize = t - stepStart;
} else {
stepSize = step;
}
} else {
if (stepStart - step <= t) {
stepSize = t - stepStart;
} else {
stepSize = -step;
}
}
initIntegration(equations.getTime(), y0, t);
// main integration loop

View File

@ -128,7 +128,7 @@ public class Percentile extends AbstractUnivariateStatistic implements Serializa
* can be reset with {@link #withEstimationType(EstimationType)}</li>
* <li>default NaN strategy: {@link NaNStrategy#REMOVED},
* can be reset with {@link #withNaNStrategy(NaNStrategy)}</li>
* <li>a KthSelector that makes use of {@link PivotingStrategy#MEDIAN_OF_3},
* <li>a KthSelector that makes use of {@link MedianOf3PivotingStrategy},
* can be reset with {@link #withKthSelector(KthSelector)}</li>
* </ul>
*/

View File

@ -119,7 +119,8 @@ public class BinomialTest {
throw new NullArgumentException();
}
final BinomialDistribution distribution = new BinomialDistribution(numberOfTrials, probability);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final BinomialDistribution distribution = new BinomialDistribution(null, numberOfTrials, probability);
switch (alternativeHypothesis) {
case GREATER_THAN:
return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);

View File

@ -155,8 +155,9 @@ public class ChiSquareTest {
throws NotPositiveException, NotStrictlyPositiveException,
DimensionMismatchException, MaxCountExceededException {
ChiSquaredDistribution distribution =
new ChiSquaredDistribution(expected.length - 1.0);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final ChiSquaredDistribution distribution =
new ChiSquaredDistribution(null, expected.length - 1.0);
return 1.0 - distribution.cumulativeProbability(chiSquare(expected, observed));
}
@ -311,8 +312,8 @@ public class ChiSquareTest {
checkArray(counts);
double df = ((double) counts.length -1) * ((double) counts[0].length - 1);
ChiSquaredDistribution distribution;
distribution = new ChiSquaredDistribution(df);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final ChiSquaredDistribution distribution = new ChiSquaredDistribution(df);
return 1 - distribution.cumulativeProbability(chiSquare(counts));
}
@ -507,8 +508,9 @@ public class ChiSquareTest {
throws DimensionMismatchException, NotPositiveException, ZeroException,
MaxCountExceededException {
ChiSquaredDistribution distribution;
distribution = new ChiSquaredDistribution((double) observed1.length - 1);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final ChiSquaredDistribution distribution =
new ChiSquaredDistribution(null, (double) observed1.length - 1);
return 1 - distribution.cumulativeProbability(
chiSquareDataSetsComparison(observed1, observed2));

View File

@ -152,10 +152,10 @@ public class GTest {
throws NotPositiveException, NotStrictlyPositiveException,
DimensionMismatchException, MaxCountExceededException {
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final ChiSquaredDistribution distribution =
new ChiSquaredDistribution(expected.length - 1.0);
return 1.0 - distribution.cumulativeProbability(
g(expected, observed));
new ChiSquaredDistribution(null, expected.length - 1.0);
return 1.0 - distribution.cumulativeProbability(g(expected, observed));
}
/**
@ -183,10 +183,10 @@ public class GTest {
throws NotPositiveException, NotStrictlyPositiveException,
DimensionMismatchException, MaxCountExceededException {
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final ChiSquaredDistribution distribution =
new ChiSquaredDistribution(expected.length - 2.0);
return 1.0 - distribution.cumulativeProbability(
g(expected, observed));
new ChiSquaredDistribution(null, expected.length - 2.0);
return 1.0 - distribution.cumulativeProbability(g(expected, observed));
}
/**
@ -472,8 +472,10 @@ public class GTest {
final long[] observed2)
throws DimensionMismatchException, NotPositiveException, ZeroException,
MaxCountExceededException {
final ChiSquaredDistribution distribution = new ChiSquaredDistribution(
(double) observed1.length - 1);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final ChiSquaredDistribution distribution =
new ChiSquaredDistribution(null, (double) observed1.length - 1);
return 1 - distribution.cumulativeProbability(
gDataSetsComparison(observed1, observed2));
}

View File

@ -181,7 +181,8 @@ public class MannWhitneyUTest {
final double z = (Umin - EU) / FastMath.sqrt(VarU);
// No try-catch or advertised exception because args are valid
final NormalDistribution standardNormal = new NormalDistribution(0, 1);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final NormalDistribution standardNormal = new NormalDistribution(null, 0, 1);
return 2 * standardNormal.cumulativeProbability(z);
}

View File

@ -124,9 +124,10 @@ public class OneWayAnova {
throws NullArgumentException, DimensionMismatchException,
ConvergenceException, MaxCountExceededException {
AnovaStats a = anovaStats(categoryData);
final AnovaStats a = anovaStats(categoryData);
// No try-catch or advertised exception because args are valid
FDistribution fdist = new FDistribution(a.dfbg, a.dfwg);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final FDistribution fdist = new FDistribution(null, a.dfbg, a.dfwg);
return 1.0 - fdist.cumulativeProbability(a.F);
}
@ -167,7 +168,8 @@ public class OneWayAnova {
ConvergenceException, MaxCountExceededException {
final AnovaStats a = anovaStats(categoryData, allowOneElementData);
final FDistribution fdist = new FDistribution(a.dfbg, a.dfwg);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final FDistribution fdist = new FDistribution(null, a.dfbg, a.dfwg);
return 1.0 - fdist.cumulativeProbability(a.F);
}

View File

@ -1056,8 +1056,9 @@ public class TTest {
final double v, final double n)
throws MaxCountExceededException, MathIllegalArgumentException {
double t = FastMath.abs(t(m, mu, v, n));
TDistribution distribution = new TDistribution(n - 1);
final double t = FastMath.abs(t(m, mu, v, n));
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final TDistribution distribution = new TDistribution(null, n - 1);
return 2.0 * distribution.cumulativeProbability(-t);
}
@ -1086,7 +1087,8 @@ public class TTest {
final double t = FastMath.abs(t(m1, m2, v1, v2, n1, n2));
final double degreesOfFreedom = df(v1, v2, n1, n2);
TDistribution distribution = new TDistribution(degreesOfFreedom);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final TDistribution distribution = new TDistribution(null, degreesOfFreedom);
return 2.0 * distribution.cumulativeProbability(-t);
}
@ -1115,7 +1117,8 @@ public class TTest {
final double t = FastMath.abs(homoscedasticT(m1, m2, v1, v2, n1, n2));
final double degreesOfFreedom = n1 + n2 - 2;
TDistribution distribution = new TDistribution(degreesOfFreedom);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final TDistribution distribution = new TDistribution(null, degreesOfFreedom);
return 2.0 * distribution.cumulativeProbability(-t);
}

View File

@ -253,7 +253,8 @@ public class WilcoxonSignedRankTest {
final double z = (Wmin - ES - 0.5) / FastMath.sqrt(VarS);
// No try-catch or advertised exception because args are valid
final NormalDistribution standardNormal = new NormalDistribution(0, 1);
// pass a null rng to avoid unneeded overhead as we will not sample from this distribution
final NormalDistribution standardNormal = new NormalDistribution(null, 0, 1);
return 2*standardNormal.cumulativeProbability(z);
}

View File

@ -18,6 +18,9 @@ package org.apache.commons.math3.util;
import java.io.PrintStream;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
/**
* Faster, more accurate, portable alternative to {@link Math} and
* {@link StrictMath} for large scale computation.
@ -804,6 +807,24 @@ public class FastMath {
return nextAfter(a, Float.POSITIVE_INFINITY);
}
/** Compute next number towards negative infinity.
* @param a number to which neighbor should be computed
* @return neighbor of a towards negative infinity
* @since 3.4
*/
public static double nextDown(final double a) {
return nextAfter(a, Double.NEGATIVE_INFINITY);
}
/** Compute next number towards negative infinity.
* @param a number to which neighbor should be computed
* @return neighbor of a towards negative infinity
* @since 3.4
*/
public static float nextDown(final float a) {
return nextAfter(a, Float.NEGATIVE_INFINITY);
}
/** Returns a pseudo-random number between 0.0 and 1.0.
* <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
* @return a random number between 0.0 and 1.0
@ -3634,6 +3655,319 @@ public class FastMath {
return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
}
/** Convert a long to interger, detecting overflows
* @param n number to convert to int
* @return integer with same valie as n if no overflows occur
* @exception MathArithmeticException if n cannot fit into an int
* @since 3.4
*/
public static int toIntExact(final long n) throws MathArithmeticException {
if (n < Integer.MIN_VALUE || n > Integer.MAX_VALUE) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW);
}
return (int) n;
}
/** Increment a number, detecting overflows.
* @param n number to increment
* @return n+1 if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static int incrementExact(final int n) throws MathArithmeticException {
if (n == Integer.MAX_VALUE) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, n, 1);
}
return n + 1;
}
/** Increment a number, detecting overflows.
* @param n number to increment
* @return n+1 if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static long incrementExact(final long n) throws MathArithmeticException {
if (n == Long.MAX_VALUE) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, n, 1);
}
return n + 1;
}
/** Decrement a number, detecting overflows.
* @param n number to decrement
* @return n-1 if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static int decrementExact(final int n) throws MathArithmeticException {
if (n == Integer.MIN_VALUE) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
}
return n - 1;
}
/** Decrement a number, detecting overflows.
* @param n number to decrement
* @return n-1 if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static long decrementExact(final long n) throws MathArithmeticException {
if (n == Long.MIN_VALUE) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
}
return n - 1;
}
/** Add two numbers, detecting overflows.
* @param a first number to add
* @param b second number to add
* @return a+b if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static int addExact(final int a, final int b) throws MathArithmeticException {
// compute sum
final int sum = a + b;
// check for overflow
if ((a ^ b) >= 0 && (sum ^ b) < 0) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, b);
}
return sum;
}
/** Add two numbers, detecting overflows.
* @param a first number to add
* @param b second number to add
* @return a+b if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static long addExact(final long a, final long b) throws MathArithmeticException {
// compute sum
final long sum = a + b;
// check for overflow
if ((a ^ b) >= 0 && (sum ^ b) < 0) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, b);
}
return sum;
}
/** Subtract two numbers, detecting overflows.
* @param a first number
* @param b second number to subtract from a
* @return a-b if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static int subtractExact(final int a, final int b) {
// compute subtraction
final int sub = a - b;
// check for overflow
if ((a ^ b) < 0 && (sub ^ b) >= 0) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, a, b);
}
return sub;
}
/** Subtract two numbers, detecting overflows.
* @param a first number
* @param b second number to subtract from a
* @return a-b if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static long subtractExact(final long a, final long b) {
// compute subtraction
final long sub = a - b;
// check for overflow
if ((a ^ b) < 0 && (sub ^ b) >= 0) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, a, b);
}
return sub;
}
/** Multiply two numbers, detecting overflows.
* @param a first number to multiply
* @param b second number to multiply
* @return a*b if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static int multiplyExact(final int a, final int b) {
if (((b > 0) && (a > Integer.MAX_VALUE / b || a < Integer.MIN_VALUE / b)) ||
((b < -1) && (a > Integer.MIN_VALUE / b || a < Integer.MAX_VALUE / b)) ||
((b == -1) && (a == Integer.MIN_VALUE))) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
}
return a * b;
}
/** Multiply two numbers, detecting overflows.
* @param a first number to multiply
* @param b second number to multiply
* @return a*b if no overflows occur
* @exception MathArithmeticException if an overflow occurs
* @since 3.4
*/
public static long multiplyExact(final long a, final long b) {
if (((b > 0l) && (a > Long.MAX_VALUE / b || a < Long.MIN_VALUE / b)) ||
((b < -1l) && (a > Long.MIN_VALUE / b || a < Long.MAX_VALUE / b)) ||
((b == -1l) && (a == Long.MIN_VALUE))) {
throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
}
return a * b;
}
/** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
* <p>
* This methods returns the same value as integer division when
* a and b are same signs, but returns a different value when
* they are opposite (i.e. q is negative).
* </p>
* @param a dividend
* @param b divisor
* @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
* @exception MathArithmeticException if b == 0
* @see #floorMod(int, int)
* @since 3.4
*/
public static int floorDiv(final int a, final int b) throws MathArithmeticException {
if (b == 0) {
throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
}
final int m = a % b;
if ((a ^ b) >= 0 || m == 0) {
// a an b have same sign, or division is exact
return a / b;
} else {
// a and b have opposite signs and division is not exact
return (a / b) - 1;
}
}
/** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
* <p>
* This methods returns the same value as integer division when
* a and b are same signs, but returns a different value when
* they are opposite (i.e. q is negative).
* </p>
* @param a dividend
* @param b divisor
* @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
* @exception MathArithmeticException if b == 0
* @see #floorMod(long, long)
* @since 3.4
*/
public static long floorDiv(final long a, final long b) throws MathArithmeticException {
if (b == 0l) {
throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
}
final long m = a % b;
if ((a ^ b) >= 0l || m == 0l) {
// a an b have same sign, or division is exact
return a / b;
} else {
// a and b have opposite signs and division is not exact
return (a / b) - 1l;
}
}
/** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
* <p>
* This methods returns the same value as integer modulo when
* a and b are same signs, but returns a different value when
* they are opposite (i.e. q is negative).
* </p>
* @param a dividend
* @param b divisor
* @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
* @exception MathArithmeticException if b == 0
* @see #floorDiv(int, int)
* @since 3.4
*/
public static int floorMod(final int a, final int b) throws MathArithmeticException {
if (b == 0) {
throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
}
final int m = a % b;
if ((a ^ b) >= 0 || m == 0) {
// a an b have same sign, or division is exact
return m;
} else {
// a and b have opposite signs and division is not exact
return b + m;
}
}
/** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
* <p>
* This methods returns the same value as integer modulo when
* a and b are same signs, but returns a different value when
* they are opposite (i.e. q is negative).
* </p>
* @param a dividend
* @param b divisor
* @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
* @exception MathArithmeticException if b == 0
* @see #floorDiv(long, long)
* @since 3.4
*/
public static long floorMod(final long a, final long b) {
if (b == 0l) {
throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
}
final long m = a % b;
if ((a ^ b) >= 0l || m == 0l) {
// a an b have same sign, or division is exact
return m;
} else {
// a and b have opposite signs and division is not exact
return b + m;
}
}
/**
* Returns the first argument with the sign of the second argument.
* A NaN {@code sign} argument is treated as positive.

View File

@ -36,7 +36,7 @@ public class KthSelector implements Serializable {
/** Minimum selection size for insertion sort rather than selection. */
private static final int MIN_SELECT_SIZE = 15;
/** A {@link PivotingStrategy} used for pivoting */
/** A {@link PivotingStrategyInterface} used for pivoting */
private final PivotingStrategyInterface pivotingStrategy;
/**

View File

@ -265,6 +265,7 @@ OVERFLOW = d\u00e9passement de capacit\u00e9
OVERFLOW_IN_FRACTION = d\u00e9passement de capacit\u00e9 pour la fraction {0}/{1}, son signe ne peut \u00eatre chang\u00e9
OVERFLOW_IN_ADDITION = d\u00e9passement de capacit\u00e9 pour l''addition : {0} + {1}
OVERFLOW_IN_SUBTRACTION = d\u00e9passement de capacit\u00e9 pour la soustraction : {0} - {1}
OVERFLOW_IN_MULTIPLICATION = d\u00e9passement de capacit\u00e9 pour la multiplication : {0} * {1}
PERCENTILE_IMPLEMENTATION_CANNOT_ACCESS_METHOD = acc\u00e8s impossible \u00e0 la m\u00e9thode {0} dans la mise en \u0153uvre du pourcentage {1}
PERCENTILE_IMPLEMENTATION_UNSUPPORTED_METHOD = l''implantation de pourcentage {0} ne dispose pas de la m\u00e9thode {1}
PERMUTATION_EXCEEDS_N = la taille de la permutation ({0}) d\u00e9passe le domaine de la permutation ({1})

View File

@ -41,7 +41,7 @@
href="http://commons.apache.org/math/javadocs/api-2.2/index.html"/>
<item name="Issue Tracking" href="/issue-tracking.html"/>
<item name="Source Repository (current)"
href="http://svn.apache.org/viewvc/commons/proper/math/trunk"/>
href="http://git-wip-us.apache.org/repos/asf/commons-math.git"/>
<item name="Wiki"
href="http://wiki.apache.org/commons/Math"/>
<item name="Developers Guide" href="/developers.html"/>

View File

@ -50,7 +50,7 @@
under the heading "Repository Checkout" on the
<a href="https://git-wip-us.apache.org/">Git at the ASF page</a>.
The git url for the current development sources of Commons Math
is <source>http://apacheid@git-wip-us.apache.org/repos/asf/commons-math.git</source>
is <source>http://git-wip-us.apache.org/repos/asf/commons-math.git</source>
for anonymous read-only access and
<source>https://apacheid@git-wip-us.apache.org/repos/asf/commons-math.git</source>
(where apacheid should be replaced by each committer Apache ID) for committers
@ -92,9 +92,9 @@
</p>
<li>
Generating patches: The requested format for generating patches is
the Unified Diff format, which can be easily generated using the svn
the Unified Diff format, which can be easily generated using the git
client or various IDEs.
<source>svn diff > patch </source>
<source>git diff -p > patch </source>
Run this command from the top-level project directory (where pom.xml
resides).
</li>
@ -132,7 +132,7 @@
</li>
<li>Submit code as attachments to the JIRA ticket. Please use one
ticket for each feature, adding multiple patches to the ticket
as necessary. Use the svn diff command to generate your patches as
as necessary. Use the git diff command to generate your patches as
diffs. Please do not submit modified copies of existing java files. Be
patient (but not <strong>too</strong> patient) with committers reviewing
patches. Post a *nudge* message to commons-dev with a reference to the
@ -154,12 +154,18 @@
<p>
Committers should configure the <source>user.name</source>,
<source>user.email</source> and <source>core.autocrlf</source>
git repository or global options with <source>git config</source>.
The first two options set the identity and mail of the committer.
The third option deals with line endings to achieve consistency
in line endings. Windows users should configure this option to
<source>true</source> while OS X and Linux users should configure
it to <source>input</source>.
git repository or global settings with <source>git config</source>.
The first two settings define the identity and mail of the committer.
The third setting deals with line endings to achieve consistency
in line endings. Windows users should configure this setting to
<source>true</source> (thus forcing git to convert CR/LF line endings
in the workspace while maintaining LF only line endings in the repository)
while OS X and Linux users should configure it to <source>input</source>
(thus forcing git to only strip accidental CR/LF when committing into
the repository, but never when cheking out files from the repository). See <a
href="http://www.git-scm.com/book/en/Customizing-Git-Git-Configuration">Customizing
Git - Git Configuration</a> in the git book for explanation about how to
configure these settings and more.
</p>
</subsection>
<subsection name='Documentation'>

View File

@ -71,12 +71,9 @@
</p>
</subsection>
<!--
<subsection name="Nightly Builds">
<subsection name="No Nightly Builds">
<p>
<a href="http://people.apache.org/builds/commons/nightly/commons-math/">
Nightly builds</a> are built once a day from the current SVN HEAD.
This is (nearly) the latest code and so should be treated with
caution!
There are no nightly builds for Apache Commons Math anymore.
</p>
</subsection>
-->

View File

@ -0,0 +1,227 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.Assert;
import org.junit.Test;
import static org.junit.Assert.*;
public class AkimaSplineInterpolatorTest
{
@Test
public void testIllegalArguments()
{
// Data set arrays of different size.
UnivariateInterpolator i = new AkimaSplineInterpolator();
try
{
double yval[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
i.interpolate( null, yval );
Assert.fail( "Failed to detect x null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
double xval[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
i.interpolate( xval, null );
Assert.fail( "Failed to detect y null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
double xval[] = { 0.0, 1.0, 2.0, 3.0 };
double yval[] = { 0.0, 1.0, 2.0, 3.0 };
i.interpolate( xval, yval );
Assert.fail( "Failed to detect insufficient data" );
}
catch ( NumberIsTooSmallException iae )
{
// Expected.
}
try
{
double xval[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
double yval[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };
i.interpolate( xval, yval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
catch ( DimensionMismatchException iae )
{
// Expected.
}
// X values not sorted.
try
{
double xval[] = { 0.0, 1.0, 0.5, 7.0, 3.5, 2.2, 8.0 };
double yval[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
i.interpolate( xval, yval );
Assert.fail( "Failed to detect unsorted arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
}
/*
* Interpolate a straight line. <p> y = 2 x - 5 <p> Tolerances determined by performing same calculation using
* Math.NET over ten runs of 100 random number draws for the same function over the same span with the same number
* of elements
*/
@Test
public void testInterpolateLine()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 1e-15;
final double maxTolerance = 1e-15;
UnivariateFunction f = new UnivariateFunction()
{
public double value( double x )
{
return 2 * x - 5;
}
};
testInterpolation( minimumX, maximumX, numberOfElements, numberOfSamples, f, interpolationTolerance,
maxTolerance );
}
/*
* Interpolate a straight line. <p> y = 3 x<sup>2</sup> - 5 x + 7 <p> Tolerances determined by performing same
* calculation using Math.NET over ten runs of 100 random number draws for the same function over the same span with
* the same number of elements
*/
@Test
public void testInterpolateParabola()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 7e-15;
final double maxTolerance = 6e-14;
UnivariateFunction f = new UnivariateFunction()
{
public double value( double x )
{
return ( 3 * x * x ) - ( 5 * x ) + 7;
}
};
testInterpolation( minimumX, maximumX, numberOfElements, numberOfSamples, f, interpolationTolerance,
maxTolerance );
}
/*
* Interpolate a straight line. <p> y = 3 x<sup>3</sup> - 0.5 x<sup>2</sup> + x - 1 <p> Tolerances determined by
* performing same calculation using Math.NET over ten runs of 100 random number draws for the same function over
* the same span with the same number of elements
*/
@Test
public void testInterpolateCubic()
{
final int numberOfElements = 10;
final double minimumX = -3;
final double maximumX = 3;
final int numberOfSamples = 100;
final double interpolationTolerance = 0.37;
final double maxTolerance = 3.8;
UnivariateFunction f = new UnivariateFunction()
{
public double value( double x )
{
return ( 3 * x * x * x ) - ( 0.5 * x * x ) + ( 1 * x ) - 1;
}
};
testInterpolation( minimumX, maximumX, numberOfElements, numberOfSamples, f, interpolationTolerance,
maxTolerance );
}
private void testInterpolation( double minimumX, double maximumX, int numberOfElements, int numberOfSamples,
UnivariateFunction f, double tolerance, double maxTolerance )
{
double expected;
double actual;
double currentX;
final double delta = ( maximumX - minimumX ) / ( (double) numberOfElements );
double xValues[] = new double[numberOfElements];
double yValues[] = new double[numberOfElements];
for ( int i = 0; i < numberOfElements; i++ )
{
xValues[i] = minimumX + delta * (double) i;
yValues[i] = f.value( xValues[i] );
}
UnivariateFunction interpolation = new AkimaSplineInterpolator().interpolate( xValues, yValues );
for ( int i = 0; i < numberOfElements; i++ )
{
currentX = xValues[i];
expected = f.value( currentX );
actual = interpolation.value( currentX );
assertTrue( Precision.equals( expected, actual ) );
}
final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed.
final UniformRealDistribution distX =
new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] );
double sumError = 0;
for ( int i = 0; i < numberOfSamples; i++ )
{
currentX = distX.sample();
expected = f.value( currentX );
actual = interpolation.value( currentX );
sumError += FastMath.abs( actual - expected );
assertEquals( expected, actual, maxTolerance );
}
assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance );
}
}

View File

@ -16,435 +16,147 @@
*/
package org.apache.commons.math3.analysis.interpolation;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* Test case for the bicubic function.
*
*/
public final class BicubicSplineInterpolatingFunctionTest {
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};
public void testPreconditions()
{
double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 };
double[][] zval = new double[xval.length][yval.length];
@SuppressWarnings("unused")
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
zval, zval, zval);
BicubicSplineInterpolatingFunction bcf = new BicubicSplineInterpolatingFunction( xval, yval, 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
try
{
bcf = new BicubicSplineInterpolatingFunction( null, yval, zval );
Assert.fail( "Failed to detect x null pointer" );
}
catch ( NullArgumentException iae )
{
// 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]);
}
try
{
bcf = new BicubicSplineInterpolatingFunction( xval, null, zval );
Assert.fail( "Failed to detect y null pointer" );
}
// 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);
catch ( NullArgumentException iae )
{
// Expected.
}
/**
* 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]);
}
try
{
bcf = new BicubicSplineInterpolatingFunction( xval, yval, null );
Assert.fail( "Failed to detect z null pointer" );
}
// 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);
catch ( NullArgumentException iae )
{
// Expected.
}
/**
* 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);
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect insufficient x data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
final BicubicSplineFunction f = new BicubicSplineFunction(coeff);
BivariateFunction derivative;
final double x = 0.435;
final double y = 0.776;
final double tol = 1e-13;
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect insufficient y data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
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);
try
{
double zval1[][] = new double[4][4];
bcf = new BicubicSplineInterpolatingFunction( xval, yval, zval1 );
Assert.fail( "Failed to detect insufficient z data" );
}
};
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);
catch ( InsufficientDataException iae )
{
// Expected.
}
};
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);
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
};
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);
catch ( DimensionMismatchException iae )
{
// Expected.
}
/**
* 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;
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
// 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;
catch ( DimensionMismatchException iae )
{
// Expected.
}
};
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]);
// X values not sorted.
try
{
double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
bcf = new BicubicSplineInterpolatingFunction( xval1, yval, zval );
Assert.fail( "Failed to detect unsorted x arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
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);
// Y values not sorted.
try
{
double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
bcf = new BicubicSplineInterpolatingFunction( xval, yval1, zval );
Assert.fail( "Failed to detect unsorted y arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
}
@ -454,73 +166,28 @@ public final class BicubicSplineInterpolatingFunctionTest {
* 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;
}
public void testInterpolatePlane()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final double minimumY = -10;
final double maximumY = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 7e-15;
final double maxTolerance = 6e-14;
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
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();
}
testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f,
interpolationTolerance, maxTolerance );
}
/**
@ -529,140 +196,85 @@ public final class BicubicSplineInterpolatingFunctionTest {
* 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;
}
public void testInterpolationParabaloid()
{
final int numberOfElements = 10;
final double minimumX = -10;
final double maximumX = 10;
final double minimumY = -10;
final double maximumY = 10;
final int numberOfSamples = 100;
final double interpolationTolerance = 2e-14;
final double maxTolerance = 6e-14;
// Function values
BivariateFunction f = new BivariateFunction() {
public double value(double x, double y) {
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]);
}
testInterpolation( minimumX, maximumX, minimumY, maximumY, numberOfElements, numberOfSamples, f,
interpolationTolerance, maxTolerance );
}
// 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;
private void testInterpolation( double minimumX, double maximumX, double minimumY, double maximumY,
int numberOfElements, int numberOfSamples, BivariateFunction f, double tolerance,
double maxTolerance )
{
double expected;
double actual;
double currentX;
double currentY;
final double deltaX = ( maximumX - minimumX ) / ( (double) numberOfElements );
final double deltaY = ( maximumY - minimumY ) / ( (double) numberOfElements );
double xValues[] = new double[numberOfElements];
double yValues[] = new double[numberOfElements];
double zValues[][] = new double[numberOfElements][numberOfElements];
for ( int i = 0; i < numberOfElements; i++ )
{
xValues[i] = minimumX + deltaX * (double) i;
for ( int j = 0; j < numberOfElements; j++ )
{
yValues[j] = minimumY + deltaY * (double) j;
zValues[i][j] = f.value( xValues[i], yValues[j] );
}
}
BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
dZdX, dZdY, dZdXdY);
double x, y;
BivariateFunction interpolation = new BicubicSplineInterpolatingFunction( xValues, yValues, zValues );
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);
for ( int i = 0; i < numberOfElements; i++ )
{
currentX = xValues[i];
for ( int j = 0; j < numberOfElements; j++ )
{
currentY = yValues[j];
expected = f.value( currentX, currentY );
actual = interpolation.value( currentX, currentY );
assertTrue( Precision.equals( expected, actual ) );
}
// System.out.println();
}
final RandomGenerator rng = new Well19937c( 1234567L ); // "tol" depends on the seed.
final UniformRealDistribution distX =
new UniformRealDistribution( rng, xValues[0], xValues[xValues.length - 1] );
final UniformRealDistribution distY =
new UniformRealDistribution( rng, yValues[0], yValues[yValues.length - 1] );
double sumError = 0;
for ( int i = 0; i < numberOfSamples; i++ )
{
currentX = distX.sample();
currentY = distY.sample();
expected = f.value( currentX, currentY );
actual = interpolation.value( currentX, currentY );
sumError += FastMath.abs( actual - expected );
assertEquals( expected, actual, maxTolerance );
}
@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) {}
assertEquals( 0.0, ( sumError / (double) numberOfSamples ), tolerance );
}
}

View File

@ -17,7 +17,10 @@
package org.apache.commons.math3.analysis.interpolation;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.random.RandomGenerator;
@ -27,53 +30,131 @@ import org.junit.Test;
/**
* Test case for the bicubic interpolator.
*
*/
public final class BicubicSplineInterpolatorTest {
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};
public void testPreconditions()
{
double[] xval = new double[] { 3, 4, 5, 6.5, 7.5 };
double[] yval = new double[] { -4, -3, -1, 2.5, 3.5 };
double[][] zval = new double[xval.length][yval.length];
@SuppressWarnings( "unused" )
BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
@SuppressWarnings("unused")
BivariateFunction p = interpolator.interpolate(xval, yval, zval);
try
{
interpolator.interpolate( null, yval, zval );
Assert.fail( "Failed to detect x null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
interpolator.interpolate( xval, null, zval );
Assert.fail( "Failed to detect y null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
interpolator.interpolate( xval, yval, null );
Assert.fail( "Failed to detect z null pointer" );
}
catch ( NullArgumentException iae )
{
// Expected.
}
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0 };
interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect insufficient x data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0 };
interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect insufficient y data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
try
{
double zval1[][] = new double[4][4];
interpolator.interpolate( xval, yval, zval1 );
Assert.fail( "Failed to detect insufficient z data" );
}
catch ( InsufficientDataException iae )
{
// Expected.
}
try
{
double xval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
catch ( DimensionMismatchException iae )
{
// Expected.
}
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
try
{
double yval1[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect data set array with different sizes." );
}
catch ( DimensionMismatchException iae )
{
// 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
// X values not sorted.
try
{
double xval1[] = { 0.0, 1.0, 0.5, 7.0, 3.5 };
interpolator.interpolate( xval1, yval, zval );
Assert.fail( "Failed to detect unsorted x arguments." );
}
catch ( NonMonotonicSequenceException iae )
{
// 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
// Y values not sorted.
try
{
double yval1[] = { 0.0, 1.0, 1.5, 0.0, 3.0 };
interpolator.interpolate( xval, yval1, zval );
Assert.fail( "Failed to detect unsorted y arguments." );
}
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
catch ( NonMonotonicSequenceException iae )
{
// Expected.
}
}
/**
@ -82,26 +163,32 @@ public final class BicubicSplineInterpolatorTest {
* z = 2 x - 3 y + 5
*/
@Test
public void testInterpolation1() {
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++) {
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) {
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++) {
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]);
}
}
@ -111,16 +198,16 @@ public final class BicubicSplineInterpolatorTest {
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 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++) {
final double tol = 2e-14;
for ( int i = 0; i < numSamples; i++ )
{
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
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);
@ -135,26 +222,32 @@ public final class BicubicSplineInterpolatorTest {
* z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
*/
@Test
public void testInterpolation2() {
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++) {
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) {
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++) {
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]);
}
}
@ -164,16 +257,16 @@ public final class BicubicSplineInterpolatorTest {
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 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++) {
final double tol = 5e-13;
for ( int i = 0; i < numSamples; i++ )
{
x = distX.sample();
for (int j = 0; j < numSamples; j++) {
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);

View File

@ -33,8 +33,8 @@ public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
*/
@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};
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);
@ -97,8 +97,8 @@ public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2, 2.5};
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++) {
@ -145,7 +145,7 @@ public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);
double[] xval = new double[] {3, 4, 5, 6.5};
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++) {

View File

@ -34,17 +34,18 @@ import org.junit.Test;
public class SplineInterpolatorTest {
/** error tolerance for spline interpolator value at knot points */
protected double knotTolerance = 1E-12;
protected double knotTolerance = 1E-14;
/** error tolerance for interpolating polynomial coefficients */
protected double coefficientTolerance = 1E-6;
protected double coefficientTolerance = 1E-14;
/** error tolerance for interpolated values -- high value is from sin test */
protected double interpolationTolerance = 1E-2;
protected double interpolationTolerance = 1E-14;
@Test
public void testInterpolateLinearDegenerateTwoSegment()
{
double tolerance = 1e-15;
double x[] = { 0.0, 0.5, 1.0 };
double y[] = { 0.0, 0.5, 1.0 };
UnivariateInterpolator i = new SplineInterpolator();
@ -60,14 +61,15 @@ public class SplineInterpolatorTest {
TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
// Check interpolation
Assert.assertEquals(0.0,f.value(0.0), interpolationTolerance);
Assert.assertEquals(0.4,f.value(0.4), interpolationTolerance);
Assert.assertEquals(1.0,f.value(1.0), interpolationTolerance);
Assert.assertEquals(0.0,f.value(0.0), tolerance);
Assert.assertEquals(0.4,f.value(0.4), tolerance);
Assert.assertEquals(1.0,f.value(1.0), tolerance);
}
@Test
public void testInterpolateLinearDegenerateThreeSegment()
{
double tolerance = 1e-15;
double x[] = { 0.0, 0.5, 1.0, 1.5 };
double y[] = { 0.0, 0.5, 1.0, 1.5 };
UnivariateInterpolator i = new SplineInterpolator();
@ -84,9 +86,9 @@ public class SplineInterpolatorTest {
TestUtils.assertEquals(polynomials[2].getCoefficients(), target, coefficientTolerance);
// Check interpolation
Assert.assertEquals(0,f.value(0), interpolationTolerance);
Assert.assertEquals(1.4,f.value(1.4), interpolationTolerance);
Assert.assertEquals(1.5,f.value(1.5), interpolationTolerance);
Assert.assertEquals(0,f.value(0), tolerance);
Assert.assertEquals(1.4,f.value(1.4), tolerance);
Assert.assertEquals(1.5,f.value(1.5), tolerance);
}
@Test
@ -108,6 +110,8 @@ public class SplineInterpolatorTest {
@Test
public void testInterpolateSin() {
double sineCoefficientTolerance = 1e-6;
double sineInterpolationTolerance = 0.0043;
double x[] =
{
0.0,
@ -136,25 +140,25 @@ public class SplineInterpolatorTest {
*/
PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
double target[] = {y[0], 1.002676d, 0d, -0.17415829d};
TestUtils.assertEquals(polynomials[0].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[0].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[1], 8.594367e-01, -2.735672e-01, -0.08707914};
TestUtils.assertEquals(polynomials[1].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[1].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[2], 1.471804e-17,-5.471344e-01, 0.08707914};
TestUtils.assertEquals(polynomials[2].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[2].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[3], -8.594367e-01, -2.735672e-01, 0.17415829};
TestUtils.assertEquals(polynomials[3].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[3].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[4], -1.002676, 6.548562e-17, 0.17415829};
TestUtils.assertEquals(polynomials[4].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[4].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[5], -8.594367e-01, 2.735672e-01, 0.08707914};
TestUtils.assertEquals(polynomials[5].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[5].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[6], 3.466465e-16, 5.471344e-01, -0.08707914};
TestUtils.assertEquals(polynomials[6].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[6].getCoefficients(), target, sineCoefficientTolerance);
target = new double[]{y[7], 8.594367e-01, 2.735672e-01, -0.17415829};
TestUtils.assertEquals(polynomials[7].getCoefficients(), target, coefficientTolerance);
TestUtils.assertEquals(polynomials[7].getCoefficients(), target, sineCoefficientTolerance);
//Check interpolation
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(FastMath.PI/4d),interpolationTolerance);
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(3d*FastMath.PI/4d),interpolationTolerance);
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(FastMath.PI/4d),sineInterpolationTolerance);
Assert.assertEquals(FastMath.sqrt(2d) / 2d,f.value(3d*FastMath.PI/4d),sineInterpolationTolerance);
}
@Test

View File

@ -32,7 +32,7 @@ public final class TricubicSplineInterpolatorTest {
/**
* Test preconditions.
*/
@Test
@Ignore@Test
public void testPreconditions() {
double[] xval = new double[] {3, 4, 5, 6.5};
double[] yval = new double[] {-4, -3, -1, 2.5};

View File

@ -27,6 +27,25 @@ public class AbstractIntegerDistributionTest {
protected final DiceDistribution diceDistribution = new DiceDistribution();
protected final double p = diceDistribution.probability(1);
@Test
public void testInverseCumulativeProbabilityMethod()
{
double precision = 0.000000000000001;
Assert.assertEquals(1, diceDistribution.inverseCumulativeProbability(0));
Assert.assertEquals(1, diceDistribution.inverseCumulativeProbability((1d-Double.MIN_VALUE)/6d));
Assert.assertEquals(2, diceDistribution.inverseCumulativeProbability((1d+precision)/6d));
Assert.assertEquals(2, diceDistribution.inverseCumulativeProbability((2d-Double.MIN_VALUE)/6d));
Assert.assertEquals(3, diceDistribution.inverseCumulativeProbability((2d+precision)/6d));
Assert.assertEquals(3, diceDistribution.inverseCumulativeProbability((3d-Double.MIN_VALUE)/6d));
Assert.assertEquals(4, diceDistribution.inverseCumulativeProbability((3d+precision)/6d));
Assert.assertEquals(4, diceDistribution.inverseCumulativeProbability((4d-Double.MIN_VALUE)/6d));
Assert.assertEquals(5, diceDistribution.inverseCumulativeProbability((4d+precision)/6d));
Assert.assertEquals(5, diceDistribution.inverseCumulativeProbability((5d-precision)/6d));//Can't use Double.MIN
Assert.assertEquals(6, diceDistribution.inverseCumulativeProbability((5d+precision)/6d));
Assert.assertEquals(6, diceDistribution.inverseCumulativeProbability((6d-precision)/6d));//Can't use Double.MIN
Assert.assertEquals(6, diceDistribution.inverseCumulativeProbability((6d)/6d));
}
@Test
public void testCumulativeProbabilitiesSingleArguments() {
for (int i = 1; i < 7; i++) {
@ -90,7 +109,7 @@ public class AbstractIntegerDistributionTest {
}
public double getNumericalVariance() {
return 12.5 - 3.5 * 3.5; // E(X^2) - E(X)^2
return 70/24; // E(X^2) - E(X)^2
}
public int getSupportLowerBound() {

View File

@ -23,14 +23,13 @@ import java.util.Locale;
import java.util.ResourceBundle;
import org.junit.Assert;
import org.junit.Test;
public class LocalizedFormatsTest {
@Test
public void testMessageNumber() {
Assert.assertEquals(319, LocalizedFormats.values().length);
Assert.assertEquals(320, LocalizedFormats.values().length);
}
@Test

View File

@ -26,7 +26,7 @@ import org.junit.Test;
/** Unit tests for {@link EvaluationRmsChecker}. */
public class EvaluationRmsCheckerTest {
/** check {@link EvaluationRmsChecker#converged(int, Evaluation, Evaluation)}. */
/** check {@link ConvergenceChecker#converged(int, Object, Object)}. */
@Test
public void testConverged() {
//setup

View File

@ -310,4 +310,28 @@ public class ClassicalRungeKuttaIntegratorTest {
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
@Test
public void testTooLargeFirstStep() {
RungeKuttaIntegrator integ = new ClassicalRungeKuttaIntegrator(0.5);
final double start = 0.0;
final double end = 0.001;
FirstOrderDifferentialEquations equations = new FirstOrderDifferentialEquations() {
public int getDimension() {
return 1;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
Assert.assertTrue(t >= FastMath.nextAfter(start, Double.NEGATIVE_INFINITY));
Assert.assertTrue(t <= FastMath.nextAfter(end, Double.POSITIVE_INFINITY));
yDot[0] = -100.0 * y[0];
}
};
integ.integrate(equations, start, new double[] { 1.0 }, end, new double[1]);
}
}

View File

@ -24,6 +24,7 @@ import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.junit.Assert;
import org.junit.Test;
import org.junit.runner.RunWith;
@ -159,12 +160,23 @@ public class FastMathStrictComparisonTest {
}
private static void callMethods(Method mathMethod, Method fastMethod,
Object[] params, int[] entries) throws IllegalAccessException,
InvocationTargetException {
Object[] params, int[] entries) throws IllegalAccessException {
try {
Object expected = mathMethod.invoke(mathMethod, params);
Object actual = fastMethod.invoke(mathMethod, params);
if (!expected.equals(actual)) {
Object expected;
try {
expected = mathMethod.invoke(mathMethod, params);
} catch (InvocationTargetException ite) {
expected = ite.getCause();
}
Object actual;
try {
actual = fastMethod.invoke(mathMethod, params);
} catch (InvocationTargetException ite) {
actual = ite.getCause();
}
if (expected instanceof ArithmeticException) {
Assert.assertEquals(MathArithmeticException.class, actual.getClass());
} else if (!expected.equals(actual)) {
reportFailedResults(mathMethod, params, expected, actual, entries);
}
} catch (IllegalArgumentException e) {

View File

@ -19,13 +19,16 @@ package org.apache.commons.math3.util;
import java.lang.reflect.Method;
import java.lang.reflect.Modifier;
import java.lang.reflect.Type;
import java.math.BigInteger;
import org.apache.commons.math3.TestUtils;
import org.apache.commons.math3.dfp.Dfp;
import org.apache.commons.math3.dfp.DfpField;
import org.apache.commons.math3.dfp.DfpMath;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.apache.commons.math3.random.MersenneTwister;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well1024a;
import org.junit.Assert;
import org.junit.Before;
import org.junit.Ignore;
@ -958,13 +961,13 @@ public class FastMathTest {
@Test
public void testNextAfter() {
// 0x402fffffffffffff 0x404123456789abcd -> 4030000000000000
Assert.assertEquals(16.0, FastMath.nextAfter(15.999999999999998, 34.27555555555555), 0.0);
Assert.assertEquals(16.0, FastMath.nextUp(15.999999999999998), 0.0);
// 0xc02fffffffffffff 0x404123456789abcd -> c02ffffffffffffe
Assert.assertEquals(-15.999999999999996, FastMath.nextAfter(-15.999999999999998, 34.27555555555555), 0.0);
// 0x402fffffffffffff 0x400123456789abcd -> 402ffffffffffffe
Assert.assertEquals(15.999999999999996, FastMath.nextAfter(15.999999999999998, 2.142222222222222), 0.0);
Assert.assertEquals(15.999999999999996, FastMath.nextDown(15.999999999999998), 0.0);
// 0xc02fffffffffffff 0x400123456789abcd -> c02ffffffffffffe
Assert.assertEquals(-15.999999999999996, FastMath.nextAfter(-15.999999999999998, 2.142222222222222), 0.0);
@ -1037,8 +1040,8 @@ public class FastMathTest {
Assert.assertEquals(-Float.MAX_VALUE,FastMath.nextAfter(Float.NEGATIVE_INFINITY, 0F), 0F);
Assert.assertEquals(Float.MAX_VALUE,FastMath.nextAfter(Float.POSITIVE_INFINITY, 0F), 0F);
Assert.assertEquals(Float.NaN,FastMath.nextAfter(Float.NaN, 0F), 0F);
Assert.assertEquals(Float.POSITIVE_INFINITY,FastMath.nextAfter(Float.MAX_VALUE, Float.POSITIVE_INFINITY), 0F);
Assert.assertEquals(Float.NEGATIVE_INFINITY,FastMath.nextAfter(-Float.MAX_VALUE, Float.NEGATIVE_INFINITY), 0F);
Assert.assertEquals(Float.POSITIVE_INFINITY,FastMath.nextUp(Float.MAX_VALUE), 0F);
Assert.assertEquals(Float.NEGATIVE_INFINITY,FastMath.nextDown(-Float.MAX_VALUE), 0F);
Assert.assertEquals(Float.MIN_VALUE, FastMath.nextAfter(0F, 1F), 0F);
Assert.assertEquals(-Float.MIN_VALUE, FastMath.nextAfter(0F, -1F), 0F);
Assert.assertEquals(0F, FastMath.nextAfter(Float.MIN_VALUE, -1F), 0F);
@ -1182,4 +1185,422 @@ public class FastMathTest {
}
}
@Test
public void testIncrementExactInt() {
int[] specialValues = new int[] {
Integer.MIN_VALUE, Integer.MIN_VALUE + 1, Integer.MIN_VALUE + 2,
Integer.MAX_VALUE, Integer.MAX_VALUE - 1, Integer.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Integer.MIN_VALUE / 2), 0 - (Integer.MIN_VALUE / 2), 1 - (Integer.MIN_VALUE / 2),
-1 + (Integer.MAX_VALUE / 2), 0 + (Integer.MAX_VALUE / 2), 1 + (Integer.MAX_VALUE / 2),
};
for (int a : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdSum = bdA.add(BigInteger.ONE);
if (bdSum.compareTo(BigInteger.valueOf(Integer.MIN_VALUE)) < 0 ||
bdSum.compareTo(BigInteger.valueOf(Integer.MAX_VALUE)) > 0) {
try {
FastMath.incrementExact(a);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdSum, BigInteger.valueOf(FastMath.incrementExact(a)));
}
}
}
@Test
public void testDecrementExactInt() {
int[] specialValues = new int[] {
Integer.MIN_VALUE, Integer.MIN_VALUE + 1, Integer.MIN_VALUE + 2,
Integer.MAX_VALUE, Integer.MAX_VALUE - 1, Integer.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Integer.MIN_VALUE / 2), 0 - (Integer.MIN_VALUE / 2), 1 - (Integer.MIN_VALUE / 2),
-1 + (Integer.MAX_VALUE / 2), 0 + (Integer.MAX_VALUE / 2), 1 + (Integer.MAX_VALUE / 2),
};
for (int a : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdSub = bdA.subtract(BigInteger.ONE);
if (bdSub.compareTo(BigInteger.valueOf(Integer.MIN_VALUE)) < 0 ||
bdSub.compareTo(BigInteger.valueOf(Integer.MAX_VALUE)) > 0) {
try {
FastMath.decrementExact(a);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdSub, BigInteger.valueOf(FastMath.decrementExact(a)));
}
}
}
@Test
public void testAddExactInt() {
int[] specialValues = new int[] {
Integer.MIN_VALUE, Integer.MIN_VALUE + 1, Integer.MIN_VALUE + 2,
Integer.MAX_VALUE, Integer.MAX_VALUE - 1, Integer.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Integer.MIN_VALUE / 2), 0 - (Integer.MIN_VALUE / 2), 1 - (Integer.MIN_VALUE / 2),
-1 + (Integer.MAX_VALUE / 2), 0 + (Integer.MAX_VALUE / 2), 1 + (Integer.MAX_VALUE / 2),
};
for (int a : specialValues) {
for (int b : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdB = BigInteger.valueOf(b);
BigInteger bdSum = bdA.add(bdB);
if (bdSum.compareTo(BigInteger.valueOf(Integer.MIN_VALUE)) < 0 ||
bdSum.compareTo(BigInteger.valueOf(Integer.MAX_VALUE)) > 0) {
try {
FastMath.addExact(a, b);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdSum, BigInteger.valueOf(FastMath.addExact(a, b)));
}
}
}
}
@Test
public void testAddExactLong() {
long[] specialValues = new long[] {
Long.MIN_VALUE, Long.MIN_VALUE + 1, Long.MIN_VALUE + 2,
Long.MAX_VALUE, Long.MAX_VALUE - 1, Long.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Long.MIN_VALUE / 2), 0 - (Long.MIN_VALUE / 2), 1 - (Long.MIN_VALUE / 2),
-1 + (Long.MAX_VALUE / 2), 0 + (Long.MAX_VALUE / 2), 1 + (Long.MAX_VALUE / 2),
};
for (long a : specialValues) {
for (long b : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdB = BigInteger.valueOf(b);
BigInteger bdSum = bdA.add(bdB);
if (bdSum.compareTo(BigInteger.valueOf(Long.MIN_VALUE)) < 0 ||
bdSum.compareTo(BigInteger.valueOf(Long.MAX_VALUE)) > 0) {
try {
FastMath.addExact(a, b);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdSum, BigInteger.valueOf(FastMath.addExact(a, b)));
}
}
}
}
@Test
public void testSubtractExactInt() {
int[] specialValues = new int[] {
Integer.MIN_VALUE, Integer.MIN_VALUE + 1, Integer.MIN_VALUE + 2,
Integer.MAX_VALUE, Integer.MAX_VALUE - 1, Integer.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Integer.MIN_VALUE / 2), 0 - (Integer.MIN_VALUE / 2), 1 - (Integer.MIN_VALUE / 2),
-1 + (Integer.MAX_VALUE / 2), 0 + (Integer.MAX_VALUE / 2), 1 + (Integer.MAX_VALUE / 2),
};
for (int a : specialValues) {
for (int b : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdB = BigInteger.valueOf(b);
BigInteger bdSub = bdA.subtract(bdB);
if (bdSub.compareTo(BigInteger.valueOf(Integer.MIN_VALUE)) < 0 ||
bdSub.compareTo(BigInteger.valueOf(Integer.MAX_VALUE)) > 0) {
try {
FastMath.subtractExact(a, b);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdSub, BigInteger.valueOf(FastMath.subtractExact(a, b)));
}
}
}
}
@Test
public void testSubtractExactLong() {
long[] specialValues = new long[] {
Long.MIN_VALUE, Long.MIN_VALUE + 1, Long.MIN_VALUE + 2,
Long.MAX_VALUE, Long.MAX_VALUE - 1, Long.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Long.MIN_VALUE / 2), 0 - (Long.MIN_VALUE / 2), 1 - (Long.MIN_VALUE / 2),
-1 + (Long.MAX_VALUE / 2), 0 + (Long.MAX_VALUE / 2), 1 + (Long.MAX_VALUE / 2),
};
for (long a : specialValues) {
for (long b : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdB = BigInteger.valueOf(b);
BigInteger bdSub = bdA.subtract(bdB);
if (bdSub.compareTo(BigInteger.valueOf(Long.MIN_VALUE)) < 0 ||
bdSub.compareTo(BigInteger.valueOf(Long.MAX_VALUE)) > 0) {
try {
FastMath.subtractExact(a, b);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdSub, BigInteger.valueOf(FastMath.subtractExact(a, b)));
}
}
}
}
@Test
public void testMultiplyExactInt() {
int[] specialValues = new int[] {
Integer.MIN_VALUE, Integer.MIN_VALUE + 1, Integer.MIN_VALUE + 2,
Integer.MAX_VALUE, Integer.MAX_VALUE - 1, Integer.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Integer.MIN_VALUE / 2), 0 - (Integer.MIN_VALUE / 2), 1 - (Integer.MIN_VALUE / 2),
-1 + (Integer.MAX_VALUE / 2), 0 + (Integer.MAX_VALUE / 2), 1 + (Integer.MAX_VALUE / 2),
};
for (int a : specialValues) {
for (int b : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdB = BigInteger.valueOf(b);
BigInteger bdMul = bdA.multiply(bdB);
if (bdMul.compareTo(BigInteger.valueOf(Integer.MIN_VALUE)) < 0 ||
bdMul.compareTo(BigInteger.valueOf(Integer.MAX_VALUE)) > 0) {
try {
FastMath.multiplyExact(a, b);
Assert.fail("an exception should have been thrown " + a + b);
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdMul, BigInteger.valueOf(FastMath.multiplyExact(a, b)));
}
}
}
}
@Test
public void testMultiplyExactLong() {
long[] specialValues = new long[] {
Long.MIN_VALUE, Long.MIN_VALUE + 1, Long.MIN_VALUE + 2,
Long.MAX_VALUE, Long.MAX_VALUE - 1, Long.MAX_VALUE - 2,
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
-1 - (Long.MIN_VALUE / 2), 0 - (Long.MIN_VALUE / 2), 1 - (Long.MIN_VALUE / 2),
-1 + (Long.MAX_VALUE / 2), 0 + (Long.MAX_VALUE / 2), 1 + (Long.MAX_VALUE / 2),
};
for (long a : specialValues) {
for (long b : specialValues) {
BigInteger bdA = BigInteger.valueOf(a);
BigInteger bdB = BigInteger.valueOf(b);
BigInteger bdMul = bdA.multiply(bdB);
if (bdMul.compareTo(BigInteger.valueOf(Long.MIN_VALUE)) < 0 ||
bdMul.compareTo(BigInteger.valueOf(Long.MAX_VALUE)) > 0) {
try {
FastMath.multiplyExact(a, b);
Assert.fail("an exception should have been thrown " + a + b);
} catch (MathArithmeticException mae) {
// expected
}
} else {
Assert.assertEquals(bdMul, BigInteger.valueOf(FastMath.multiplyExact(a, b)));
}
}
}
}
@Test(expected=MathArithmeticException.class)
public void testToIntExactTooLow() {
FastMath.toIntExact(-1l + Integer.MIN_VALUE);
}
@Test(expected=MathArithmeticException.class)
public void testToIntExactTooHigh() {
FastMath.toIntExact(+1l + Integer.MAX_VALUE);
}
@Test
public void testToIntExact() {
for (int n = -1000; n < 1000; ++n) {
Assert.assertEquals(n, FastMath.toIntExact(0l + n));
}
Assert.assertEquals(Integer.MIN_VALUE, FastMath.toIntExact(0l + Integer.MIN_VALUE));
Assert.assertEquals(Integer.MAX_VALUE, FastMath.toIntExact(0l + Integer.MAX_VALUE));
}
@Test
public void testFloorDivInt() {
Assert.assertEquals(+1, FastMath.floorDiv(+4, +3));
Assert.assertEquals(-2, FastMath.floorDiv(-4, +3));
Assert.assertEquals(-2, FastMath.floorDiv(+4, -3));
Assert.assertEquals(+1, FastMath.floorDiv(-4, -3));
try {
FastMath.floorDiv(1, 0);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
for (int a = -100; a <= 100; ++a) {
for (int b = -100; b <= 100; ++b) {
if (b != 0) {
Assert.assertEquals(poorManFloorDiv(a, b), FastMath.floorDiv(a, b));
}
}
}
}
@Test
public void testFloorModInt() {
Assert.assertEquals(+1, FastMath.floorMod(+4, +3));
Assert.assertEquals(+2, FastMath.floorMod(-4, +3));
Assert.assertEquals(-2, FastMath.floorMod(+4, -3));
Assert.assertEquals(-1, FastMath.floorMod(-4, -3));
try {
FastMath.floorMod(1, 0);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
for (int a = -100; a <= 100; ++a) {
for (int b = -100; b <= 100; ++b) {
if (b != 0) {
Assert.assertEquals(poorManFloorMod(a, b), FastMath.floorMod(a, b));
}
}
}
}
@Test
public void testFloorDivModInt() {
RandomGenerator generator = new Well1024a(0x7ccab45edeaab90al);
for (int i = 0; i < 10000; ++i) {
int a = generator.nextInt();
int b = generator.nextInt();
if (b == 0) {
try {
FastMath.floorDiv(a, b);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
int d = FastMath.floorDiv(a, b);
int m = FastMath.floorMod(a, b);
Assert.assertEquals(FastMath.toIntExact(poorManFloorDiv(a, b)), d);
Assert.assertEquals(FastMath.toIntExact(poorManFloorMod(a, b)), m);
Assert.assertEquals(a, d * b + m);
if (b < 0) {
Assert.assertTrue(m <= 0);
Assert.assertTrue(-m < -b);
} else {
Assert.assertTrue(m >= 0);
Assert.assertTrue(m < b);
}
}
}
}
@Test
public void testFloorDivLong() {
Assert.assertEquals(+1l, FastMath.floorDiv(+4l, +3l));
Assert.assertEquals(-2l, FastMath.floorDiv(-4l, +3l));
Assert.assertEquals(-2l, FastMath.floorDiv(+4l, -3l));
Assert.assertEquals(+1l, FastMath.floorDiv(-4l, -3l));
try {
FastMath.floorDiv(1l, 0l);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
for (long a = -100l; a <= 100l; ++a) {
for (long b = -100l; b <= 100l; ++b) {
if (b != 0) {
Assert.assertEquals(poorManFloorDiv(a, b), FastMath.floorDiv(a, b));
}
}
}
}
@Test
public void testFloorModLong() {
Assert.assertEquals(+1l, FastMath.floorMod(+4l, +3l));
Assert.assertEquals(+2l, FastMath.floorMod(-4l, +3l));
Assert.assertEquals(-2l, FastMath.floorMod(+4l, -3l));
Assert.assertEquals(-1l, FastMath.floorMod(-4l, -3l));
try {
FastMath.floorMod(1l, 0l);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
for (long a = -100l; a <= 100l; ++a) {
for (long b = -100l; b <= 100l; ++b) {
if (b != 0) {
Assert.assertEquals(poorManFloorMod(a, b), FastMath.floorMod(a, b));
}
}
}
}
@Test
public void testFloorDivModLong() {
RandomGenerator generator = new Well1024a(0xb87b9bc14c96ccd5l);
for (int i = 0; i < 10000; ++i) {
long a = generator.nextLong();
long b = generator.nextLong();
if (b == 0) {
try {
FastMath.floorDiv(a, b);
Assert.fail("an exception should have been thrown");
} catch (MathArithmeticException mae) {
// expected
}
} else {
long d = FastMath.floorDiv(a, b);
long m = FastMath.floorMod(a, b);
Assert.assertEquals(poorManFloorDiv(a, b), d);
Assert.assertEquals(poorManFloorMod(a, b), m);
Assert.assertEquals(a, d * b + m);
if (b < 0) {
Assert.assertTrue(m <= 0);
Assert.assertTrue(-m < -b);
} else {
Assert.assertTrue(m >= 0);
Assert.assertTrue(m < b);
}
}
}
}
private long poorManFloorDiv(long a, long b) {
// find q0, r0 such that a = q0 b + r0
BigInteger q0 = BigInteger.valueOf(a / b);
BigInteger r0 = BigInteger.valueOf(a % b);
BigInteger fd = BigInteger.valueOf(Integer.MIN_VALUE);
BigInteger bigB = BigInteger.valueOf(b);
for (int k = -2; k < 2; ++k) {
// find another pair q, r such that a = q b + r
BigInteger bigK = BigInteger.valueOf(k);
BigInteger q = q0.subtract(bigK);
BigInteger r = r0.add(bigK.multiply(bigB));
if (r.abs().compareTo(bigB.abs()) < 0 &&
(r.longValue() == 0l || ((r.longValue() ^ b) & 0x8000000000000000l) == 0)) {
if (fd.compareTo(q) < 0) {
fd = q;
}
}
}
return fd.longValue();
}
private long poorManFloorMod(long a, long b) {
return a - b * poorManFloorDiv(a, b);
}
}