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

Fixed end of line encoding problems.

This commit is contained in:
Luc Maisonobe 2014-10-17 10:23:23 +02:00
parent bf88cac1ac
commit f44d4852ab
2 changed files with 452 additions and 452 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

450
src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java
View File

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

454
src/test/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolatorTest.java
View File

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