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

MATH-1014 

"PolynomialCurveFitter" as replacement of "PolynomialFitter". Some tests have
been obsoleted by the refactoring (which hides the optimizer and thus avoids
some potential misuses).


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1543906 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski 2013-11-20 18:57:53 +00:00
parent 54a7796ff2
commit c1ba07bb65
4 changed files with 297 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/changes/changes.xml
View File

@ -51,6 +51,9 @@ If the output is not quite correct, check for invisible trailing spaces!
</properties> </properties>
<body> <body>
<release version="3.3" date="TBD" description="TBD"> <release version="3.3" date="TBD" description="TBD">
<action dev="erans" type="add" issue="MATH-1014">
Refactoring of curve fitters (package "o.a.c.m.fitting").
</action>
<action dev="tn" type="add" issue="MATH-970"> <action dev="tn" type="add" issue="MATH-970">
Added possibility to retrieve the best found solution of the "SimplexSolver" in case Added possibility to retrieve the best found solution of the "SimplexSolver" in case
the iteration limit has been reached. The "optimize(OptimizationData...)" method now the iteration limit has been reached. The "optimize(OptimizationData...)" method now

125
src/main/java/org/apache/commons/math3/fitting/PolynomialCurveFitter.java Normal file
View File

@ -0,0 +1,125 @@
/*
* 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.fitting;
import java.util.Collection;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.exception.MathInternalError;
import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer;
import org.apache.commons.math3.fitting.leastsquares.WithStartPoint;
import org.apache.commons.math3.fitting.leastsquares.WithMaxIterations;
import org.apache.commons.math3.linear.DiagonalMatrix;
/**
* Fits points to a {@link
* org.apache.commons.math3.analysis.polynomials.PolynomialFunction.Parametric polynomial}
* function.
* <br/>
* The size of the {@link #withStartPoint(double[]) initial guess} array defines the
* degree of the polynomial to be fitted.
* They must be sorted in increasing order of the polynomial's degree.
* The optimal values of the coefficients will be returned in the same order.
*
* @version $Id$
* @since 3.3
*/
public class PolynomialCurveFitter extends AbstractCurveFitter<LevenbergMarquardtOptimizer>
implements WithStartPoint<PolynomialCurveFitter>,
WithMaxIterations<PolynomialCurveFitter> {
/** Parametric function to be fitted. */
private static final PolynomialFunction.Parametric FUNCTION = new PolynomialFunction.Parametric();
/** Initial guess. */
private final double[] initialGuess;
/** Maximum number of iterations of the optimization algorithm. */
private final int maxIter;
/**
* Contructor used by the factory methods.
*
* @param initialGuess Initial guess.
* @param maxIter Maximum number of iterations of the optimization algorithm.
* @throws MathInternalError if {@code initialGuess} is {@code null}.
*/
private PolynomialCurveFitter(double[] initialGuess,
int maxIter) {
this.initialGuess = initialGuess;
this.maxIter = maxIter;
}
/**
* Creates a default curve fitter.
* Zero will be used as initial guess for the coefficients, and the maximum
* number of iterations of the optimization algorithm is set to
* {@link Integer#MAX_VALUE}.
*
* @param degree Degree of the polynomial to be fitted.
* @return a curve fitter.
*
* @see #withStartPoint(double[])
* @see #withMaxIterations(int)
*/
public static PolynomialCurveFitter create(int degree) {
return new PolynomialCurveFitter(new double[degree + 1], Integer.MAX_VALUE);
}
/** {@inheritDoc} */
public PolynomialCurveFitter withStartPoint(double[] start) {
return new PolynomialCurveFitter(start.clone(),
maxIter);
}
/** {@inheritDoc} */
public PolynomialCurveFitter withMaxIterations(int max) {
return new PolynomialCurveFitter(initialGuess,
max);
}
/** {@inheritDoc} */
@Override
protected LevenbergMarquardtOptimizer getOptimizer(Collection<WeightedObservedPoint> observations) {
// Prepare least-squares problem.
final int len = observations.size();
final double[] target = new double[len];
final double[] weights = new double[len];
int i = 0;
for (WeightedObservedPoint obs : observations) {
target[i] = obs.getY();
weights[i] = obs.getWeight();
++i;
}
final AbstractCurveFitter.TheoreticalValuesFunction model
= new AbstractCurveFitter.TheoreticalValuesFunction(FUNCTION,
observations);
if (initialGuess == null) {
throw new MathInternalError();
}
// Return a new optimizer set up to fit a Gaussian curve to the
// observed points.
return LevenbergMarquardtOptimizer.create()
.withMaxEvaluations(Integer.MAX_VALUE)
.withMaxIterations(maxIter)
.withStartPoint(initialGuess)
.withTarget(target)
.withWeight(new DiagonalMatrix(weights))
.withModelAndJacobian(model.getModelFunction(),
model.getModelFunctionJacobian());
}
}

3
src/main/java/org/apache/commons/math3/fitting/PolynomialFitter.java
View File

@ -26,7 +26,10 @@ import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimiz
* *
* @version $Id: PolynomialFitter.java 1416643 2012-12-03 19:37:14Z tn $ * @version $Id: PolynomialFitter.java 1416643 2012-12-03 19:37:14Z tn $
* @since 2.0 * @since 2.0
* @deprecated As of 3.3. Please use {@link PolynomialCurveFitter} and
* {@link WeightedObservedPoints} instead.
*/ */
@Deprecated
public class PolynomialFitter extends CurveFitter<PolynomialFunction.Parametric> { public class PolynomialFitter extends CurveFitter<PolynomialFunction.Parametric> {
/** /**
* Simple constructor. * Simple constructor.

166
src/test/java/org/apache/commons/math3/fitting/PolynomialCurveFitterTest.java Normal file
View File

@ -0,0 +1,166 @@
/*
* 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.fitting;
import java.util.Random;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.distribution.RealDistribution;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.TestUtils;
import org.junit.Test;
import org.junit.Assert;
/**
* Test for class {@link PolynomialCurveFitter}.
*/
public class PolynomialCurveFitterTest {
@Test
public void testFit() {
final RealDistribution rng = new UniformRealDistribution(-100, 100);
rng.reseedRandomGenerator(64925784252L);
final double[] coeff = { 12.9, -3.4, 2.1 }; // 12.9 - 3.4 x + 2.1 x^2
final PolynomialFunction f = new PolynomialFunction(coeff);
// Collect data from a known polynomial.
final WeightedObservedPoints obs = new WeightedObservedPoints();
for (int i = 0; i < 100; i++) {
final double x = rng.sample();
obs.add(x, f.value(x));
}
// Start fit from initial guesses that are far from the optimal values.
final PolynomialCurveFitter fitter
= PolynomialCurveFitter.create(0).withStartPoint(new double[] { -1e-20, 3e15, -5e25 });
final double[] best = fitter.fit(obs.toList());
TestUtils.assertEquals("best != coeff", coeff, best, 1e-12);
}
@Test
public void testNoError() {
final Random randomizer = new Random(64925784252l);
for (int degree = 1; degree < 10; ++degree) {
final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
final WeightedObservedPoints obs = new WeightedObservedPoints();
for (int i = 0; i <= degree; ++i) {
obs.add(1.0, i, p.value(i));
}
final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
for (double x = -1.0; x < 1.0; x += 0.01) {
final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
(1.0 + FastMath.abs(p.value(x)));
Assert.assertEquals(0.0, error, 1.0e-6);
}
}
}
@Test
public void testSmallError() {
final Random randomizer = new Random(53882150042l);
double maxError = 0;
for (int degree = 0; degree < 10; ++degree) {
final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
final WeightedObservedPoints obs = new WeightedObservedPoints();
for (double x = -1.0; x < 1.0; x += 0.01) {
obs.add(1.0, x, p.value(x) + 0.1 * randomizer.nextGaussian());
}
final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
for (double x = -1.0; x < 1.0; x += 0.01) {
final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
(1.0 + FastMath.abs(p.value(x)));
maxError = FastMath.max(maxError, error);
Assert.assertTrue(FastMath.abs(error) < 0.1);
}
}
Assert.assertTrue(maxError > 0.01);
}
@Test
public void testRedundantSolvable() {
// Levenberg-Marquardt should handle redundant information gracefully
checkUnsolvableProblem(true);
}
@Test
public void testLargeSample() {
final Random randomizer = new Random(0x5551480dca5b369bl);
double maxError = 0;
for (int degree = 0; degree < 10; ++degree) {
final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
final WeightedObservedPoints obs = new WeightedObservedPoints();
for (int i = 0; i < 40000; ++i) {
final double x = -1.0 + i / 20000.0;
obs.add(1.0, x, p.value(x) + 0.1 * randomizer.nextGaussian());
}
final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
for (double x = -1.0; x < 1.0; x += 0.01) {
final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
(1.0 + FastMath.abs(p.value(x)));
maxError = FastMath.max(maxError, error);
Assert.assertTrue(FastMath.abs(error) < 0.01);
}
}
Assert.assertTrue(maxError > 0.001);
}
private void checkUnsolvableProblem(boolean solvable) {
final Random randomizer = new Random(1248788532l);
for (int degree = 0; degree < 10; ++degree) {
final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
final WeightedObservedPoints obs = new WeightedObservedPoints();
// reusing the same point over and over again does not bring
// information, the problem cannot be solved in this case for
// degrees greater than 1 (but one point is sufficient for
// degree 0)
for (double x = -1.0; x < 1.0; x += 0.01) {
obs.add(1.0, 0.0, p.value(0.0));
}
try {
fitter.fit(obs.toList());
Assert.assertTrue(solvable || (degree == 0));
} catch(ConvergenceException e) {
Assert.assertTrue((! solvable) && (degree > 0));
}
}
}
private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
final double[] coefficients = new double[degree + 1];
for (int i = 0; i <= degree; ++i) {
coefficients[i] = randomizer.nextGaussian();
}
return new PolynomialFunction(coefficients);
}
}