You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/06/19 14:36:17 UTC
svn commit: r786479 - in /commons/proper/math/trunk/src:
java/org/apache/commons/math/optimization/fitting/
mantissa/src/org/spaceroots/mantissa/fitting/
mantissa/tests-src/org/spaceroots/mantissa/fitting/ site/xdoc/
site/xdoc/userguide/ test/org/apach...
Author: luc
Date: Fri Jun 19 12:36:16 2009
New Revision: 786479
URL: http://svn.apache.org/viewvc?rev=786479&view=rev
Log:
merged curve fitting from mantissa into commons-math
Added:
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/CurveFitter.java
- copied, changed from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/AbstractCurveFitter.java
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java
- copied, changed from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicCoefficientsGuesser.java
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java
- copied, changed from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicFitter.java
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java (with props)
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java (with props)
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/PolynomialFitter.java
- copied, changed from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialFitter.java
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java (with props)
commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/package.html
- copied unchanged from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/package.html
commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/
commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/HarmonicFitterTest.java
- copied, changed from r784713, commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/HarmonicFitterTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/PolynomialFitterTest.java
- copied, changed from r784713, commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/PolynomialFitterTest.java
Removed:
commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/
commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/
Modified:
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/site/xdoc/userguide/index.xml
commons/proper/math/trunk/src/site/xdoc/userguide/optimization.xml
Copied: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/CurveFitter.java (from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/AbstractCurveFitter.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/CurveFitter.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/CurveFitter.java&p1=commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/AbstractCurveFitter.java&r1=784713&r2=786479&rev=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/AbstractCurveFitter.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/CurveFitter.java Fri Jun 19 12:36:16 2009
@@ -1,222 +1,191 @@
-// 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.
+/*
+ * 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.spaceroots.mantissa.fitting;
+package org.apache.commons.math.optimization.fitting;
-import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
-import org.spaceroots.mantissa.estimation.*;
-
-/** This class is the base class for all curve fitting classes in the package.
+import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
+import org.apache.commons.math.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
+import org.apache.commons.math.optimization.OptimizationException;
+import org.apache.commons.math.optimization.VectorialPointValuePair;
+
+/** Fitter for parametric univariate real functions y = f(x).
+ * <p>When a univariate real function y = f(x) does depend on some
+ * unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
+ * this class can be used to find these parameters. It does this
+ * by <em>fitting</em> the curve so it remains very close to a set of
+ * observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
+ * y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
+ * is done by finding the parameters values that minimizes the objective
+ * function ∑(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
+ * really a least squares problem.</p>
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class CurveFitter {
- * <p>This class handles all common features of curve fitting like the
- * sample points handling. It declares two methods ({@link
- * #valueAt} and {@link #partial}) which should be implemented by
- * sub-classes to define the precise shape of the curve they
- * represent.</p>
+ /** Optimizer to use for the fitting. */
+ private final DifferentiableMultivariateVectorialOptimizer optimizer;
- * @version $Id$
- * @author L. Maisonobe
+ /** Observed points. */
+ private final List<WeightedObservedPoint> observations;
- */
+ /** Simple constructor.
+ * @param optimizer optimizer to use for the fitting
+ */
+ public CurveFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) {
+ this.optimizer = optimizer;
+ observations = new ArrayList<WeightedObservedPoint>();
+ }
-public abstract class AbstractCurveFitter
- implements EstimationProblem, Serializable {
+ /** Add an observed (x,y) point to the sample with unit weight.
+ * <p>Calling this method is equivalent to call
+ * <code>addObservedPoint(1.0, x, y)</code>.</p>
+ * @param x abscissa of the point
+ * @param y observed value of the point at x, after fitting we should
+ * have f(x) as close as possible to this value
+ * @see #addObservedPoint(double, double, double)
+ * @see #addObservedPoint(WeightedObservedPoint)
+ * @see #getObservations()
+ */
+ public void addObservedPoint(double x, double y) {
+ addObservedPoint(1.0, x, y);
+ }
- /** Simple constructor.
- * @param n number of coefficients in the underlying function
- * @param estimator estimator to use for the fitting
- */
- protected AbstractCurveFitter(int n, Estimator estimator) {
-
- coefficients = new EstimatedParameter[n];
- measurements = new ArrayList();
- this.estimator = estimator;
- }
-
- /** Simple constructor.
- * @param coefficients first estimate of the coefficients. A
- * reference to this array is hold by the newly created object. Its
- * elements will be adjusted during the fitting process and they will
- * be set to the adjusted coefficients at the end.
- * @param estimator estimator to use for the fitting
- */
- protected AbstractCurveFitter(EstimatedParameter[] coefficients,
- Estimator estimator) {
-
- this.coefficients = coefficients;
- measurements = new ArrayList();
- this.estimator = estimator;
- }
-
- /** Add a weighted (x,y) pair to the sample.
- * @param weight weight of this pair in the fit
- * @param x abscissa
- * @param y ordinate, we have <code>y = f (x)</code>
- */
- public void addWeightedPair(double weight, double x, double y) {
- measurements.add(new FitMeasurement(weight, x, y));
- }
-
- /** Perform the fitting.
-
- * <p>This method compute the coefficients of the curve that best
- * fit the sample of weighted pairs previously given through calls
- * to the {@link #addWeightedPair addWeightedPair} method.</p>
-
- * @return coefficients of the curve
- * @exception EstimationException if the fitting is not possible
- * (for example if the sample has to few independant points)
-
- */
- public double[] fit()
- throws EstimationException {
- // perform the fit
- estimator.estimate(this);
-
- // extract the coefficients
- double[] fittedCoefficients = new double[coefficients.length];
- for (int i = 0; i < coefficients.length; ++i) {
- fittedCoefficients[i] = coefficients[i].getEstimate();
- }
-
- return fittedCoefficients;
-
- }
-
- public WeightedMeasurement[] getMeasurements() {
- return (WeightedMeasurement[]) measurements.toArray(new FitMeasurement[measurements.size()]);
- }
-
- /** Get the unbound parameters of the problem.
- * For a curve fitting, none of the function coefficient is bound.
- * @return unbound parameters
- */
- public EstimatedParameter[] getUnboundParameters() {
- return (EstimatedParameter[]) coefficients.clone();
- }
-
- /** Get all the parameters of the problem.
- * @return parameters
- */
- public EstimatedParameter[] getAllParameters() {
- return (EstimatedParameter[]) coefficients.clone();
- }
-
- /** Utility method to sort the measurements with respect to the abscissa.
-
- * <p>This method is provided as a utility for derived classes. As
- * an example, the {@link HarmonicFitter} class needs it in order to
- * compute a first guess of the coefficients to initialize the
- * estimation algorithm.</p>
-
- */
- protected void sortMeasurements() {
-
- // Since the samples are almost always already sorted, this
- // method is implemented as an insertion sort that reorders the
- // elements in place. Insertion sort is very efficient in this case.
- FitMeasurement curr = (FitMeasurement) measurements.get(0);
- for (int j = 1; j < measurements.size (); ++j) {
- FitMeasurement prec = curr;
- curr = (FitMeasurement) measurements.get(j);
- if (curr.x < prec.x) {
- // the current element should be inserted closer to the beginning
- int i = j - 1;
- FitMeasurement mI = (FitMeasurement) measurements.get(i);
- while ((i >= 0) && (curr.x < mI.x)) {
- measurements.set(i + 1, mI);
- if (i-- != 0) {
- mI = (FitMeasurement) measurements.get(i);
- } else {
- mI = null;
- }
- }
- measurements.set(i + 1, curr);
- curr = (FitMeasurement) measurements.get(j);
- }
- }
-
- }
-
- /** Get the value of the function at x according to the current parameters value.
- * @param x abscissa at which the theoretical value is requested
- * @return theoretical value at x
- */
- public abstract double valueAt(double x);
-
- /** Get the derivative of the function at x with respect to parameter p.
- * @param x abscissa at which the partial derivative is requested
- * @param p parameter with respect to which the derivative is requested
- * @return partial derivative
- */
- public abstract double partial(double x, EstimatedParameter p);
-
- /** This class represents the fit measurements.
- * One measurement is a weighted pair (x, y), where <code>y = f
- * (x)</code> is the value of the function at x abscissa. This class
- * is an inner class because the methods related to the computation
- * of f values and derivative are proveded by the fitter
- * implementations.
- */
- public class FitMeasurement
- extends WeightedMeasurement {
+ /** Add an observed weighted (x,y) point to the sample.
+ * @param weight weight of the observed point in the fit
+ * @param x abscissa of the point
+ * @param y observed value of the point at x, after fitting we should
+ * have f(x) as close as possible to this value
+ * @see #addObservedPoint(double, double)
+ * @see #addObservedPoint(WeightedObservedPoint)
+ * @see #getObservations()
+ */
+ public void addObservedPoint(double weight, double x, double y) {
+ observations.add(new WeightedObservedPoint(weight, x, y));
+ }
- /** Simple constructor.
- * @param weight weight of the measurement in the fitting process
- * @param x abscissa of the measurement
- * @param y ordinate of the measurement
+ /** Add an observed weighted (x,y) point to the sample.
+ * @param observed observed point to add
+ * @see #addObservedPoint(double, double)
+ * @see #addObservedPoint(double, double, double)
+ * @see #getObservations()
*/
- public FitMeasurement(double weight, double x, double y) {
- super(weight, y);
- this.x = x;
+ public void addObservedPoint(WeightedObservedPoint observed) {
+ observations.add(observed);
}
- /** Get the value of the fitted function at x.
- * @return theoretical value at the measurement abscissa
+ /** Get the observed points.
+ * @return observed points
+ * @see #addObservedPoint(double, double)
+ * @see #addObservedPoint(double, double, double)
+ * @see #addObservedPoint(WeightedObservedPoint)
*/
- public double getTheoreticalValue() {
- return valueAt(x);
+ public WeightedObservedPoint[] getObservations() {
+ return observations.toArray(new WeightedObservedPoint[observations.size()]);
}
- /** Partial derivative with respect to one function coefficient.
- * @param p parameter with respect to which the derivative is requested
- * @return partial derivative
+ /** Fit a curve.
+ * <p>This method compute the coefficients of the curve that best
+ * fit the sample of weighted pairs previously given through calls
+ * to the {@link #addWeightedPair addWeightedPair} method.</p>
+ * @param f parametric function to fit
+ * @param initialGuess first guess of the function parameters
+ * @return fitted parameters
+ * @exception FunctionEvaluationException if the objective function throws one during
+ * the search
+ * @exception OptimizationException if the algorithm failed to converge
+ * @exception IllegalArgumentException if the start point dimension is wrong
*/
- public double getPartial(EstimatedParameter p) {
- return partial(x, p);
+ public double[] fit(final ParametricRealFunction f,
+ final double[] initialGuess)
+ throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
+
+ // prepare least squares problem
+ double[] target = new double[observations.size()];
+ double[] weights = new double[observations.size()];
+ int i = 0;
+ for (WeightedObservedPoint point : observations) {
+ target[i] = point.getY();
+ weights[i] = point.getWeight();
+ ++i;
+ }
+
+ // perform the fit
+ VectorialPointValuePair optimum =
+ optimizer.optimize(new TheoreticalValuesFunction(f), target, weights, initialGuess);
+
+ // extract the coefficients
+ return optimum.getPointRef();
+
}
- /** Abscissa of the measurement. */
- public final double x;
+ /** Vectorial function computing function theoretical values. */
+ private class TheoreticalValuesFunction
+ implements DifferentiableMultivariateVectorialFunction {
+
+ /** Function to fit. */
+ private final ParametricRealFunction f;
+
+ /** Simple constructor.
+ * @param f function to fit.
+ */
+ public TheoreticalValuesFunction(final ParametricRealFunction f) {
+ this.f = f;
+ }
- private static final long serialVersionUID = -2682582852369995960L;
+ /** {@inheritDoc} */
+ public MultivariateMatrixFunction jacobian() {
+ return new MultivariateMatrixFunction() {
+ public double[][] value(double[] point)
+ throws FunctionEvaluationException, IllegalArgumentException {
+
+ final double[][] jacobian = new double[observations.size()][];
+
+ int i = 0;
+ for (WeightedObservedPoint observed : observations) {
+ jacobian[i++] = f.gradient(observed.getX(), point);
+ }
- }
+ return jacobian;
- /** Coefficients of the function */
- protected EstimatedParameter[] coefficients;
+ }
+ };
+ }
+
+ /** {@inheritDoc} */
+ public double[] value(double[] point)
+ throws FunctionEvaluationException, IllegalArgumentException {
+
+ // compute the residuals
+ final double[] values = new double[observations.size()];
+ int i = 0;
+ for (WeightedObservedPoint observed : observations) {
+ values[i++] = f.value(observed.getX(), point);
+ }
- /** Measurements vector */
- protected List measurements;
+ return values;
- /** Estimator for the fitting problem. */
- private Estimator estimator;
+ }
+
+ }
}
Copied: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java (from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicCoefficientsGuesser.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java&p1=commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicCoefficientsGuesser.java&r1=784713&r2=786479&rev=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicCoefficientsGuesser.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java Fri Jun 19 12:36:16 2009
@@ -1,32 +1,23 @@
-// 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.spaceroots.mantissa.fitting;
-
-import java.io.Serializable;
-
-import org.spaceroots.mantissa.functions.FunctionException;
-import org.spaceroots.mantissa.functions.ExhaustedSampleException;
-import org.spaceroots.mantissa.functions.vectorial.SampledFunctionIterator;
-import org.spaceroots.mantissa.functions.vectorial.VectorialValuedPair;
+/*
+ * 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.
+ */
-import org.spaceroots.mantissa.quadrature.vectorial.EnhancedSimpsonIntegratorSampler;
+package org.apache.commons.math.optimization.fitting;
-import org.spaceroots.mantissa.estimation.EstimationException;
+import org.apache.commons.math.optimization.OptimizationException;
/** This class guesses harmonic coefficients from a sample.
@@ -127,143 +118,183 @@
* estimations, these operations run in O(n) time, where n is the
* number of measurements.</p>
- * @version $Id$
- * @author L. Maisonobe
+ * @version $Revision$ $Date$
+ * @since 2.0
*/
+public class HarmonicCoefficientsGuesser {
-public class HarmonicCoefficientsGuesser
- implements Serializable{
+ /** Sampled observations. */
+ private final WeightedObservedPoint[] observations;
- public HarmonicCoefficientsGuesser(AbstractCurveFitter.FitMeasurement[] measurements) {
- this.measurements =
- (AbstractCurveFitter.FitMeasurement[]) measurements.clone();
- a = Double.NaN;
- omega = Double.NaN;
- }
-
- /** Estimate a first guess of the coefficients.
-
- * @exception ExhaustedSampleException if the sample is exhausted.
-
- * @exception FunctionException if the integrator throws one.
-
- * @exception EstimationException if the sample is too short or if
- * the first guess cannot be computed (when the elements under the
- * square roots are negative).
- * */
- public void guess()
- throws ExhaustedSampleException, FunctionException, EstimationException {
- guessAOmega();
- guessPhi();
- }
-
- /** Estimate a first guess of the a and ω coefficients.
-
- * @exception ExhaustedSampleException if the sample is exhausted.
-
- * @exception FunctionException if the integrator throws one.
-
- * @exception EstimationException if the sample is too short or if
- * the first guess cannot be computed (when the elements under the
- * square roots are negative).
-
- */
- private void guessAOmega()
- throws ExhaustedSampleException, FunctionException, EstimationException {
-
- // initialize the sums for the linear model between the two integrals
- double sx2 = 0.0;
- double sy2 = 0.0;
- double sxy = 0.0;
- double sxz = 0.0;
- double syz = 0.0;
-
- // build the integrals sampler
- F2FP2Iterator iter = new F2FP2Iterator(measurements);
- SampledFunctionIterator sampler =
- new EnhancedSimpsonIntegratorSampler(iter);
- VectorialValuedPair p0 = sampler.nextSamplePoint();
- double p0X = p0.x;
- double[] p0Y = p0.y;
-
- // get the points for the linear model
- while (sampler.hasNext()) {
-
- VectorialValuedPair point = sampler.nextSamplePoint();
- double pX = point.x;
- double[] pY = point.y;
-
- double x = pX - p0X;
- double y = pY[0] - p0Y[0];
- double z = pY[1] - p0Y[1];
-
- sx2 += x * x;
- sy2 += y * y;
- sxy += x * y;
- sxz += x * z;
- syz += y * z;
+ /** Guessed amplitude. */
+ private double a;
+ /** Guessed pulsation ω. */
+ private double omega;
+
+ /** Guessed phase φ. */
+ private double phi;
+
+ /** Simple constructor.
+ * @param observations sampled observations
+ */
+ public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) {
+ this.observations = observations.clone();
+ a = Double.NaN;
+ omega = Double.NaN;
}
- // compute the amplitude and pulsation coefficients
- double c1 = sy2 * sxz - sxy * syz;
- double c2 = sxy * sxz - sx2 * syz;
- double c3 = sx2 * sy2 - sxy * sxy;
- if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) {
- throw new EstimationException("unable to guess a first estimate");
+ /** Estimate a first guess of the coefficients.
+ * @exception OptimizationException if the sample is too short or if
+ * the first guess cannot be computed (when the elements under the
+ * square roots are negative).
+ * */
+ public void guess() throws OptimizationException {
+ sortObservations();
+ guessAOmega();
+ guessPhi();
}
- a = Math.sqrt(c1 / c2);
- omega = Math.sqrt(c2 / c3);
-
- }
-
- /** Estimate a first guess of the φ coefficient.
-
- * @exception ExhaustedSampleException if the sample is exhausted.
-
- * @exception FunctionException if the sampler throws one.
- */
- private void guessPhi()
- throws ExhaustedSampleException, FunctionException {
+ /** Sort the observations with respect to the abscissa.
+ */
+ private void sortObservations() {
+
+ // Since the samples are almost always already sorted, this
+ // method is implemented as an insertion sort that reorders the
+ // elements in place. Insertion sort is very efficient in this case.
+ WeightedObservedPoint curr = observations[0];
+ for (int j = 1; j < observations.length; ++j) {
+ WeightedObservedPoint prec = curr;
+ curr = observations[j];
+ if (curr.getX() < prec.getX()) {
+ // the current element should be inserted closer to the beginning
+ int i = j - 1;
+ WeightedObservedPoint mI = observations[i];
+ while ((i >= 0) && (curr.getX() < mI.getX())) {
+ observations[i + 1] = mI;
+ if (i-- != 0) {
+ mI = observations[i];
+ } else {
+ mI = null;
+ }
+ }
+ observations[i + 1] = curr;
+ curr = observations[j];
+ }
+ }
- SampledFunctionIterator iter = new FFPIterator(measurements);
+ }
- // initialize the means
- double fcMean = 0.0;
- double fsMean = 0.0;
+ /** Estimate a first guess of the a and ω coefficients.
+ * @exception OptimizationException if the sample is too short or if
+ * the first guess cannot be computed (when the elements under the
+ * square roots are negative).
+ */
+ private void guessAOmega() throws OptimizationException {
+
+ // initialize the sums for the linear model between the two integrals
+ double sx2 = 0.0;
+ double sy2 = 0.0;
+ double sxy = 0.0;
+ double sxz = 0.0;
+ double syz = 0.0;
+
+ double currentX = observations[0].getX();
+ double currentY = observations[0].getY();
+ double f2Integral = 0;
+ double fPrime2Integral = 0;
+ final double startX = currentX;
+ for (int i = 1; i < observations.length; ++i) {
+
+ // one step forward
+ final double previousX = currentX;
+ final double previousY = currentY;
+ currentX = observations[i].getX();
+ currentY = observations[i].getY();
+
+ // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
+ // considering a linear model for f (and therefore constant f')
+ final double dx = currentX - previousX;
+ final double dy = currentY - previousY;
+ final double f2StepIntegral =
+ dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
+ final double fPrime2StepIntegral = dy * dy / dx;
+
+ final double x = currentX - startX;
+ f2Integral += f2StepIntegral;
+ fPrime2Integral += fPrime2StepIntegral;
+
+ sx2 += x * x;
+ sy2 += f2Integral * f2Integral;
+ sxy += x * f2Integral;
+ sxz += x * fPrime2Integral;
+ syz += f2Integral * fPrime2Integral;
+
+ }
+
+ // compute the amplitude and pulsation coefficients
+ double c1 = sy2 * sxz - sxy * syz;
+ double c2 = sxy * sxz - sx2 * syz;
+ double c3 = sx2 * sy2 - sxy * sxy;
+ if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) {
+ throw new OptimizationException("unable to first guess the harmonic coefficients");
+ }
+ a = Math.sqrt(c1 / c2);
+ omega = Math.sqrt(c2 / c3);
- while (iter.hasNext()) {
- VectorialValuedPair point = iter.nextSamplePoint();
- double omegaX = omega * point.x;
- double cosine = Math.cos(omegaX);
- double sine = Math.sin(omegaX);
- fcMean += omega * point.y[0] * cosine - point.y[1] * sine;
- fsMean += omega * point.y[0] * sine + point.y[1] * cosine;
}
- phi = Math.atan2(-fsMean, fcMean);
+ /** Estimate a first guess of the φ coefficient.
+ */
+ private void guessPhi() {
+
+ // initialize the means
+ double fcMean = 0.0;
+ double fsMean = 0.0;
+
+ double currentX = observations[0].getX();
+ double currentY = observations[0].getY();
+ for (int i = 1; i < observations.length; ++i) {
+
+ // one step forward
+ final double previousX = currentX;
+ final double previousY = currentY;
+ currentX = observations[i].getX();
+ currentY = observations[i].getY();
+ final double currentYPrime = (currentY - previousY) / (currentX - previousX);
+
+ double omegaX = omega * currentX;
+ double cosine = Math.cos(omegaX);
+ double sine = Math.sin(omegaX);
+ fcMean += omega * currentY * cosine - currentYPrime * sine;
+ fsMean += omega * currentY * sine + currentYPrime * cosine;
- }
+ }
- public double getOmega() {
- return omega;
- }
+ phi = Math.atan2(-fsMean, fcMean);
- public double getA() {
- return a;
- }
+ }
- public double getPhi() {
- return phi;
- }
+ /** Get the guessed amplitude a.
+ * @return guessed amplitude a;
+ */
+ public double getGuessedAmplitude() {
+ return a;
+ }
- private AbstractCurveFitter.FitMeasurement[] measurements;
- private double a;
- private double omega;
- private double phi;
+ /** Get the guessed pulsation ω.
+ * @return guessed pulsation ω
+ */
+ public double getGuessedPulsation() {
+ return omega;
+ }
- private static final long serialVersionUID = 2400399048702758814L;
+ /** Get the guessed phase φ.
+ * @return guessed phase φ
+ */
+ public double getGuessedPhase() {
+ return phi;
+ }
}
Copied: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java (from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicFitter.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java&p1=commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicFitter.java&r1=784713&r2=786479&rev=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicFitter.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java Fri Jun 19 12:36:16 2009
@@ -1,168 +1,132 @@
-// 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.spaceroots.mantissa.fitting;
-
-import org.spaceroots.mantissa.estimation.EstimatedParameter;
-import org.spaceroots.mantissa.estimation.EstimationException;
-import org.spaceroots.mantissa.estimation.Estimator;
-import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
-import org.spaceroots.mantissa.functions.ExhaustedSampleException;
-import org.spaceroots.mantissa.functions.FunctionException;
+/*
+ * 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.
+ */
-/** This class implements a curve fitting specialized for sinusoids.
+package org.apache.commons.math.optimization.fitting;
+
+import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
+import org.apache.commons.math.optimization.OptimizationException;
+/** This class implements a curve fitting specialized for sinusoids.
* <p>Harmonic fitting is a very simple case of curve fitting. The
- * estimated coefficients are the amplitude a, the pulsation omega and
- * the phase phi: <code>f (t) = a cos (omega t + phi)</code>. They are
+ * estimated coefficients are the amplitude a, the pulsation ω and
+ * the phase φ: <code>f (t) = a cos (ω t + φ)</code>. They are
* searched by a least square estimator initialized with a rough guess
* based on integrals.</p>
-
- * <p>This class <emph>is by no means optimized</emph>, neither versus
- * space nor versus time performance.</p>
-
- * @version $Id$
- * @author L. Maisonobe
-
+ * @version $Revision$ $Date$
+ * @since 2.0
*/
+public class HarmonicFitter {
-public class HarmonicFitter
- extends AbstractCurveFitter {
-
- /** Simple constructor.
- * @param estimator estimator to use for the fitting
- */
- public HarmonicFitter(Estimator estimator) {
- super(3, estimator);
- coefficients[0] = new EstimatedParameter("a", 2.0 * Math.PI);
- coefficients[1] = new EstimatedParameter("omega", 0.0);
- coefficients[2] = new EstimatedParameter("phi", 0.0);
- firstGuessNeeded = true;
- }
-
- /**
- * Simple constructor.
-
- * <p>This constructor can be used when a first estimate of the
- * coefficients is already known.</p>
-
- * @param coefficients first estimate of the coefficients.
- * A reference to this array is hold by the newly created
- * object. Its elements will be adjusted during the fitting process
- * and they will be set to the adjusted coefficients at the end.
- * @param estimator estimator to use for the fitting
-
- */
- public HarmonicFitter(EstimatedParameter[] coefficients,
- Estimator estimator) {
- super(coefficients, estimator);
- firstGuessNeeded = false;
- }
-
- public double[] fit()
- throws EstimationException {
- if (firstGuessNeeded) {
- if (measurements.size() < 4) {
- throw new EstimationException("sample must contain at least {0} points",
- new String[] {
- Integer.toString(4)
- });
- }
-
- sortMeasurements();
-
- try {
- HarmonicCoefficientsGuesser guesser =
- new HarmonicCoefficientsGuesser((FitMeasurement[]) getMeasurements());
- guesser.guess();
-
- coefficients[0].setEstimate(guesser.getA());
- coefficients[1].setEstimate(guesser.getOmega());
- coefficients[2].setEstimate(guesser.getPhi());
- } catch(ExhaustedSampleException e) {
- throw new EstimationException(e);
- } catch(FunctionException e) {
- throw new EstimationException(e);
- }
+ /** Fitter for the coefficients. */
+ private final CurveFitter fitter;
- firstGuessNeeded = false;
+ /** Values for amplitude, pulsation ω and phase φ. */
+ private double[] parameters;
+ /** Simple constructor.
+ * @param optimizer optimizer to use for the fitting
+ */
+ public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) {
+ this.fitter = new CurveFitter(optimizer);
+ parameters = null;
}
- return super.fit();
+ /** Simple constructor.
+ * <p>This constructor can be used when a first guess of the
+ * coefficients is already known.</p>
+ * @param optimizer optimizer to use for the fitting
+ * @param initialGuess guessed values for amplitude (index 0),
+ * pulsation ω (index 1) and phase φ (index 2)
+ */
+ public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer,
+ final double[] initialGuess) {
+ this.fitter = new CurveFitter(optimizer);
+ this.parameters = initialGuess.clone();
+ }
- }
+ /** Add an observed weighted (x,y) point to the sample.
+ * @param weight weight of the observed point in the fit
+ * @param x abscissa of the point
+ * @param y observed value of the point at x, after fitting we should
+ * have P(x) as close as possible to this value
+ */
+ public void addObservedPoint(double weight, double x, double y) {
+ fitter.addObservedPoint(weight, x, y);
+ }
- /** Get the current amplitude coefficient estimate.
- * Get a, where <code>f (t) = a cos (omega t + phi)</code>
- * @return current amplitude coefficient estimate
- */
- public double getAmplitude() {
- return coefficients[0].getEstimate();
- }
-
- /** Get the current pulsation coefficient estimate.
- * Get omega, where <code>f (t) = a cos (omega t + phi)</code>
- * @return current pulsation coefficient estimate
- */
- public double getPulsation() {
- return coefficients[1].getEstimate();
- }
-
- /** Get the current phase coefficient estimate.
- * Get phi, where <code>f (t) = a cos (omega t + phi)</code>
- * @return current phase coefficient estimate
- */
- public double getPhase() {
- return coefficients[2].getEstimate();
- }
-
- /** Get the value of the function at x according to the current parameters value.
- * @param x abscissa at which the theoretical value is requested
- * @return theoretical value at x
- */
- public double valueAt(double x) {
- double a = coefficients[0].getEstimate();
- double omega = coefficients[1].getEstimate();
- double phi = coefficients[2].getEstimate();
- return a * Math.cos(omega * x + phi);
- }
-
- /** Get the derivative of the function at x with respect to parameter p.
- * @param x abscissa at which the partial derivative is requested
- * @param p parameter with respect to which the derivative is requested
- * @return partial derivative
- */
- public double partial(double x, EstimatedParameter p) {
- double a = coefficients[0].getEstimate();
- double omega = coefficients[1].getEstimate();
- double phi = coefficients[2].getEstimate();
- if (p == coefficients[0]) {
- return Math.cos(omega * x + phi);
- } else if (p == coefficients[1]) {
- return -a * x * Math.sin(omega * x + phi);
- } else {
- return -a * Math.sin(omega * x + phi);
+ /** Fit an harmonic function to the observed points.
+ * @return harmonic function best fitting the observed points
+ * @throws OptimizationException if the sample is too short or if
+ * the first guess cannot be computed
+ */
+ public HarmonicFunction fit() throws OptimizationException {
+ try {
+
+ // shall we compute the first guess of the parameters ourselves ?
+ if (parameters == null) {
+ final WeightedObservedPoint[] observations = fitter.getObservations();
+ if (observations.length < 4) {
+ throw new OptimizationException("sample contains {0} observed points, at least {1} are required",
+ observations.length, 4);
+ }
+
+ HarmonicCoefficientsGuesser guesser = new HarmonicCoefficientsGuesser(observations);
+ guesser.guess();
+ parameters = new double[] {
+ guesser.getGuessedAmplitude(),
+ guesser.getGuessedPulsation(),
+ guesser.getGuessedPhase()
+ };
+
+ }
+
+ double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters);
+ return new HarmonicFunction(fitted[0], fitted[1], fitted[2]);
+
+ } catch (FunctionEvaluationException fee) {
+ // this should never happen
+ throw MathRuntimeException.createInternalError(fee);
+ }
}
- }
- /** Indicator of the need to compute a first guess of the coefficients. */
- private boolean firstGuessNeeded;
+ /** Parametric harmonic function. */
+ private static class ParametricHarmonicFunction implements ParametricRealFunction {
+
+ /** {@inheritDoc} */
+ public double value(double x, double[] parameters) {
+ final double a = parameters[0];
+ final double omega = parameters[1];
+ final double phi = parameters[2];
+ return a * Math.cos(omega * x + phi);
+ }
+
+ /** {@inheritDoc} */
+ public double[] gradient(double x, double[] parameters) {
+ final double a = parameters[0];
+ final double omega = parameters[1];
+ final double phi = parameters[2];
+ final double alpha = omega * x + phi;
+ final double cosAlpha = Math.cos(alpha);
+ final double sinAlpha = Math.sin(alpha);
+ return new double[] { cosAlpha, -a * x * sinAlpha, -a * sinAlpha };
+ }
- private static final long serialVersionUID = -8722683066277473450L;
+ }
}
Added: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java?rev=786479&view=auto
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java (added)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java Fri Jun 19 12:36:16 2009
@@ -0,0 +1,79 @@
+/*
+ * 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.math.optimization.fitting;
+
+import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
+
+/** Harmonic function of the form <code>f (t) = a cos (ω t + φ)</code>.
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class HarmonicFunction implements DifferentiableUnivariateRealFunction {
+
+ /** Amplitude a. */
+ private final double a;
+
+ /** Pulsation ω. */
+ private final double omega;
+
+ /** Phase φ. */
+ private final double phi;
+
+ /** Simple constructor.
+ * @param a amplitude
+ * @param omega pulsation
+ * @param phi phase
+ */
+ public HarmonicFunction(double a, double omega, double phi) {
+ this.a = a;
+ this.omega = omega;
+ this.phi = phi;
+ }
+
+ /** {@inheritDoc} */
+ public double value(double x) {
+ return a * Math.cos(omega * x + phi);
+ }
+
+ /** {@inheritDoc} */
+ public HarmonicFunction derivative() {
+ return new HarmonicFunction(a * omega, omega, phi + Math.PI / 2);
+ }
+
+ /** Get the amplitude a.
+ * @return amplitude a;
+ */
+ public double getAmplitude() {
+ return a;
+ }
+
+ /** Get the pulsation ω.
+ * @return pulsation ω
+ */
+ public double getPulsation() {
+ return omega;
+ }
+
+ /** Get the phase φ.
+ * @return phase φ
+ */
+ public double getPhase() {
+ return phi;
+ }
+
+}
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Added: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java?rev=786479&view=auto
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java (added)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java Fri Jun 19 12:36:16 2009
@@ -0,0 +1,50 @@
+/*
+ * 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.math.optimization.fitting;
+
+import org.apache.commons.math.FunctionEvaluationException;
+
+/**
+ * An interface representing a real function that depends on one independent
+ * variable plus some extra parameters.
+ *
+ * @version $Revision$ $Date$
+ */
+public interface ParametricRealFunction {
+
+ /**
+ * Compute the value of the function.
+ * @param x the point for which the function value should be computed
+ * @param parameters function parameters
+ * @return the value
+ * @throws FunctionEvaluationException if the function evaluation fails
+ */
+ public double value(double x, double[] parameters)
+ throws FunctionEvaluationException;
+
+ /**
+ * Compute the gradient of the function with respect to its parameters.
+ * @param x the point for which the function value should be computed
+ * @param parameters function parameters
+ * @return the value
+ * @throws FunctionEvaluationException if the function evaluation fails
+ */
+ public double[] gradient(double x, double[] parameters)
+ throws FunctionEvaluationException;
+
+}
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Copied: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/PolynomialFitter.java (from r784713, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialFitter.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/PolynomialFitter.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/PolynomialFitter.java&p1=commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialFitter.java&r1=784713&r2=786479&rev=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialFitter.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/PolynomialFitter.java Fri Jun 19 12:36:16 2009
@@ -1,107 +1,103 @@
-// 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.spaceroots.mantissa.fitting;
-
-import org.spaceroots.mantissa.estimation.EstimatedParameter;
-import org.spaceroots.mantissa.estimation.Estimator;
-import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
+/*
+ * 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.
+ */
-/** This class implements a curve fitting specialized for polynomials.
+package org.apache.commons.math.optimization.fitting;
+
+import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
+import org.apache.commons.math.optimization.OptimizationException;
+/** This class implements a curve fitting specialized for polynomials.
* <p>Polynomial fitting is a very simple case of curve fitting. The
- * estimated coefficients are the polynom coefficients. They are
+ * estimated coefficients are the polynomial coefficients. They are
* searched by a least square estimator.</p>
-
- * <p>This class <emph>is by no means optimized</emph>, neither in
- * space nor in time performance.</p>
-
- * @see PolynomialCoefficient
-
- * @version $Id$
- * @author L. Maisonobe
-
+ * @version $Revision$ $Date$
+ * @since 2.0
*/
-public class PolynomialFitter
- extends AbstractCurveFitter {
+public class PolynomialFitter {
- /** Simple constructor.
+ /** Fitter for the coefficients. */
+ private final CurveFitter fitter;
- * <p>The polynomial fitter built this way are complete polynoms,
- * ie. a n-degree polynom has n+1 coefficients. In order to build
- * fitter for sparse polynoms (for example <code>a x^20 - b
- * x^30</code>, on should first build the coefficients array and
- * provide it to {@link
- * #PolynomialFitter(PolynomialCoefficient[], int, double, double,
- * double)}.</p>
- * @param degree maximal degree of the polynom
- * @param estimator estimator to use for the fitting
- */
- public PolynomialFitter(int degree, Estimator estimator) {
- super(degree + 1, estimator);
- for (int i = 0; i < coefficients.length; ++i) {
- coefficients[i] = new PolynomialCoefficient(i);
- }
- }
+ /** Polynomial degree. */
+ private final int degree;
- /** Simple constructor.
+ /** Simple constructor.
+ * <p>The polynomial fitter built this way are complete polynomials,
+ * ie. a n-degree polynomial has n+1 coefficients.</p>
+ * @param degree maximal degree of the polynomial
+ * @param optimizer optimizer to use for the fitting
+ */
+ public PolynomialFitter(int degree, final DifferentiableMultivariateVectorialOptimizer optimizer) {
+ this.fitter = new CurveFitter(optimizer);
+ this.degree = degree;
+ }
- * <p>This constructor can be used either when a first estimate of
- * the coefficients is already known (which is of little interest
- * because the fit cost is the same whether a first guess is known or
- * not) or when one needs to handle sparse polynoms like <code>a
- * x^20 - b x^30</code>.</p>
-
- * @param coefficients first estimate of the coefficients.
- * A reference to this array is hold by the newly created
- * object. Its elements will be adjusted during the fitting process
- * and they will be set to the adjusted coefficients at the end.
- * @param estimator estimator to use for the fitting
- */
- public PolynomialFitter(PolynomialCoefficient[] coefficients,
- Estimator estimator) {
- super(coefficients, estimator);
- }
-
- /** Get the value of the function at x according to the current parameters value.
- * @param x abscissa at which the theoretical value is requested
- * @return theoretical value at x
- */
- public double valueAt(double x) {
- double y = coefficients[coefficients.length - 1].getEstimate();
- for (int i = coefficients.length - 2; i >= 0; --i) {
- y = y * x + coefficients[i].getEstimate();
+ /** Add an observed weighted (x,y) point to the sample.
+ * @param weight weight of the observed point in the fit
+ * @param x abscissa of the point
+ * @param y observed value of the point at x, after fitting we should
+ * have P(x) as close as possible to this value
+ */
+ public void addObservedPoint(double weight, double x, double y) {
+ fitter.addObservedPoint(weight, x, y);
}
- return y;
- }
- /** Get the derivative of the function at x with respect to parameter p.
- * @param x abscissa at which the partial derivative is requested
- * @param p parameter with respect to which the derivative is requested
- * @return partial derivative
- */
- public double partial(double x, EstimatedParameter p) {
- if (p instanceof PolynomialCoefficient) {
- return Math.pow(x, ((PolynomialCoefficient) p).degree);
+ /** Get the polynomial fitting the weighted (x, y) points.
+ * @return polynomial function best fitting the observed points
+ * @exception OptimizationException if the algorithm failed to converge
+ */
+ public PolynomialFunction fit()
+ throws OptimizationException {
+ try {
+ return new PolynomialFunction(fitter.fit(new ParametricPolynomial(), new double[degree + 1]));
+ } catch (FunctionEvaluationException fee) {
+ // this should never happen
+ throw MathRuntimeException.createInternalError(fee);
+ }
}
- throw new RuntimeException("internal error");
- }
- private static final long serialVersionUID = -744904084649890769L;
+ /** Dedicated parametric polynomial class. */
+ private static class ParametricPolynomial implements ParametricRealFunction {
+
+ /** {@inheritDoc} */
+ public double[] gradient(double x, double[] parameters)
+ throws FunctionEvaluationException {
+ final double[] gradient = new double[parameters.length];
+ double xn = 1.0;
+ for (int i = 0; i < parameters.length; ++i) {
+ gradient[i] = xn;
+ xn *= x;
+ }
+ return gradient;
+ }
+
+ /** {@inheritDoc} */
+ public double value(final double x, final double[] parameters) {
+ double y = 0;
+ for (int i = parameters.length - 1; i >= 0; --i) {
+ y = y * x + parameters[i];
+ }
+ return y;
+ }
+
+ }
}
Added: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java?rev=786479&view=auto
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java (added)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java Fri Jun 19 12:36:16 2009
@@ -0,0 +1,75 @@
+/*
+ * 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.math.optimization.fitting;
+
+import java.io.Serializable;
+
+/** This class is a simple container for weighted observed point in
+ * {@link CurveFitter curve fitting}.
+ * <p>Instances of this class are guaranteed to be immutable.</p>
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class WeightedObservedPoint implements Serializable {
+
+ /** Serializable version id. */
+ private static final long serialVersionUID = 5306874947404636157L;
+
+ /** Weight of the measurement in the fitting process. */
+ private final double weight;
+
+ /** Abscissa of the point. */
+ private final double x;
+
+ /** Observed value of the function at x. */
+ private final double y;
+
+ /** Simple constructor.
+ * @param weight weight of the measurement in the fitting process
+ * @param x abscissa of the measurement
+ * @param y ordinate of the measurement
+ */
+ public WeightedObservedPoint(final double weight, final double x, final double y) {
+ this.weight = weight;
+ this.x = x;
+ this.y = y;
+ }
+
+ /** Get the weight of the measurement in the fitting process.
+ * @return weight of the measurement in the fitting process
+ */
+ public double getWeight() {
+ return weight;
+ }
+
+ /** Get the abscissa of the point.
+ * @return abscissa of the point
+ */
+ public double getX() {
+ return x;
+ }
+
+ /** Get the observed value of the function at x.
+ * @return observed value of the function at x
+ */
+ public double getY() {
+ return y;
+ }
+
+}
+
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=786479&r1=786478&r2=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Jun 19 12:36:16 2009
@@ -39,6 +39,9 @@
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
+ <action dev="luc" type="add" >
+ Added curve fitting with a general case and two specific cases (polynomial and harmonic).
+ </action>
<action dev="psteitz" type="update" issue="MATH-276" due-to="Mark Anderson">
Optimized Complex isNaN(), isInfinite() by moving computation to constructor.
</action>
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/index.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/index.xml?rev=786479&r1=786478&r2=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/index.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/index.xml Fri Jun 19 12:36:16 2009
@@ -120,6 +120,7 @@
<li><a href="analysis.html#a12.3_Linear_Programming">12.3 Linear Programming</a></li>
<li><a href="optimization.html#a12.4_Direct_Methods">12.4 Direct Methods</a></li>
<li><a href="analysis.html#a12.5_General_Case">12.5 General Case</a></li>
+ <li><a href="analysis.html#a12.6_Curve_Fitting">12.6 Curve Fitting</a></li>
</ul></li>
<li><a href="ode.html">13. Ordinary Differential Equations Integration</a>
<ul>
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/optimization.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/optimization.xml?rev=786479&r1=786478&r2=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/optimization.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/optimization.xml Fri Jun 19 12:36:16 2009
@@ -40,6 +40,7 @@
using direct search methods (i.e. not using derivatives),</li>
<li>the general package handles multivariate scalar or vector functions
using derivatives.</li>
+ <li>the fitting package handles curve fitting by univariate real functions</li>
</ul>
</p>
<p>
@@ -212,6 +213,55 @@
solver).
</p>
</subsection>
+ <subsection name="12.6 Curve Fitting" href="fitting">
+ <p>
+ The fitting package deals with curve fitting for univariate real functions.
+ When a univariate real function y = f(x) does depend on some unknown parameters
+ p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>, curve fitting can be used to
+ find these parameters. It does this by <em>fitting</em> the curve so it remains
+ very close to a set of observed points (x<sub>0</sub>, y<sub>0</sub>),
+ (x<sub>1</sub>, y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This
+ fitting is done by finding the parameters values that minimizes the objective
+ function sum(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is really a least
+ squares problem.
+ </p>
+ <p>
+ For all provided curve fitters, the operating principle is the same. Users must first
+ create an instance of the fitter, then add the observed points and once the complete
+ sample of observed points has been added they must call the <code>fit</code> method
+ which will compute the parameters that best fit the sample. A weight is associated
+ with each observed point, this allows to take into account uncertainty on some points
+ when they come from loosy measurements for example. If no such information exist and
+ all points should be treated the same, it is safe to put 1.0 as the weight for all points.
+ </p>
+ <p>
+ The <a
+ href="../apidocs/org/apache/commons/math/optimization/fitting/CurveFitter.html">
+ CurveFitter</a> class provides curve fitting for general curves. Users must
+ provide their own implementation of the curve template as a class implementing
+ the <a
+ href="../apidocs/org/apache/commons/math/optimization/fitting/ParametricRealFunction.html">
+ ParametricRealFunction</a> interface and they must provide the initial guess of the
+ parameters. The more specialized <a
+ href="../apidocs/org/apache/commons/math/optimization/fitting/PolynomialFitter.html">
+ PolynomialFitter</a> and <a
+ href="../apidocs/org/apache/commons/math/optimization/fitting/HarmonicFitter.html">
+ HarmonicFitter</a> classes require neither an implementation of the parametric real function
+ not an initial guess as they are able to compute them by themselves.
+ </p>
+ <p>
+ An example of fitting a polynomial is given here:
+ </p>
+ <source>PolynomialFitter fitter = new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
+fitter.addObservedPoint(-1.00, 2.021170021833143);
+fitter.addObservedPoint(-0.99 2.221135431136975);
+fitter.addObservedPoint(-0.98 2.09985277659314);
+fitter.addObservedPoint(-0.97 2.0211192647627025);
+// lots of lines ommitted
+fitter.addObservedPoint( 0.99, -2.4345814727089854);
+PolynomialFunction fitted = fitter.fit();
+ </source>
+ </subsection>
</section>
</body>
</document>
Copied: commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/HarmonicFitterTest.java (from r784713, commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/HarmonicFitterTest.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/HarmonicFitterTest.java?p2=commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/HarmonicFitterTest.java&p1=commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/HarmonicFitterTest.java&r1=784713&r2=786479&rev=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/HarmonicFitterTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/HarmonicFitterTest.java Fri Jun 19 12:36:16 2009
@@ -15,161 +15,118 @@
// specific language governing permissions and limitations
// under the License.
-package org.spaceroots.mantissa.fitting;
+package org.apache.commons.math.optimization.fitting;
-import java.util.Random;
-import junit.framework.*;
-
-import org.spaceroots.mantissa.estimation.EstimationException;
-import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
-import org.spaceroots.mantissa.estimation.WeightedMeasurement;
-
-public class HarmonicFitterTest
- extends TestCase {
-
- public HarmonicFitterTest(String name) {
- super(name);
- }
-
- public void testNoError()
- throws EstimationException {
- HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
-
- HarmonicFitter fitter =
- new HarmonicFitter(new LevenbergMarquardtEstimator());
- for (double x = 0.0; x < 1.3; x += 0.01) {
- fitter.addWeightedPair(1.0, x, f.valueAt(x));
- }
-
- double[] coeffs = fitter.fit();
-
- HarmonicFunction fitted = new HarmonicFunction(coeffs[0],
- coeffs[1],
- coeffs[2]);
- assertTrue(Math.abs(coeffs[0] - f.getA()) < 1.0e-13);
- assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 1.0e-13);
- assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.0e-13);
-
- for (double x = -1.0; x < 1.0; x += 0.01) {
- assertTrue(Math.abs(f.valueAt(x) - fitted.valueAt(x)) < 1.0e-13);
- }
-
- }
-
- public void test1PercentError()
- throws EstimationException {
- Random randomizer = new Random(64925784252l);
- HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
-
- HarmonicFitter fitter =
- new HarmonicFitter(new LevenbergMarquardtEstimator());
- for (double x = 0.0; x < 10.0; x += 0.1) {
- fitter.addWeightedPair(1.0, x,
- f.valueAt(x) + 0.01 * randomizer.nextGaussian());
- }
-
- double[] coeffs = fitter.fit();
-
- new HarmonicFunction(coeffs[0], coeffs[1], coeffs[2]);
- assertTrue(Math.abs(coeffs[0] - f.getA()) < 7.6e-4);
- assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 2.7e-3);
- assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.3e-2);
-
- WeightedMeasurement[] measurements = fitter.getMeasurements();
- for (int i = 0; i < measurements.length; ++i) {
- WeightedMeasurement m = measurements[i];
- assertTrue(Math.abs(measurements[i].getMeasuredValue()
- - m.getTheoreticalValue()) < 0.04);
- }
-
- }
-
- public void testUnsorted()
- throws EstimationException {
- Random randomizer = new Random(64925784252l);
- HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
-
- HarmonicFitter fitter =
- new HarmonicFitter(new LevenbergMarquardtEstimator());
-
- // build a regularly spaced array of measurements
- int size = 100;
- double[] xTab = new double[size];
- double[] yTab = new double[size];
- for (int i = 0; i < size; ++i) {
- xTab[i] = 0.1 * i;
- yTab[i] = f.valueAt (xTab[i]) + 0.01 * randomizer.nextGaussian();
- }
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
- // shake it
- for (int i = 0; i < size; ++i) {
- int i1 = randomizer.nextInt(size);
- int i2 = randomizer.nextInt(size);
- double xTmp = xTab[i1];
- double yTmp = yTab[i1];
- xTab[i1] = xTab[i2];
- yTab[i1] = yTab[i2];
- xTab[i2] = xTmp;
- yTab[i2] = yTmp;
- }
-
- // pass it to the fitter
- for (int i = 0; i < size; ++i) {
- fitter.addWeightedPair(1.0, xTab[i], yTab[i]);
- }
-
- double[] coeffs = fitter.fit();
-
- new HarmonicFunction(coeffs[0], coeffs[1], coeffs[2]);
- assertTrue(Math.abs(coeffs[0] - f.getA()) < 7.6e-4);
- assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 3.5e-3);
- assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.5e-2);
-
- WeightedMeasurement[] measurements = fitter.getMeasurements();
- for (int i = 0; i < measurements.length; ++i) {
- WeightedMeasurement m = measurements[i];
- assertTrue(Math.abs(m.getMeasuredValue() - m.getTheoreticalValue())
- < 0.04);
- }
- }
-
- public static Test suite() {
- return new TestSuite(HarmonicFitterTest.class);
- }
-
- /** Center an angle with respect to another one. */
- private static double center(double a, double ref) {
- double twoPi = Math.PI + Math.PI;
- return a - twoPi * Math.floor((a + Math.PI - ref) / twoPi);
- }
-
- private static class HarmonicFunction {
- public HarmonicFunction(double a, double omega, double phi) {
- this.a = a;
- this.omega = omega;
- this.phi = phi;
- }
-
- public double valueAt(double x) {
- return a * Math.cos(omega * x + phi);
- }
-
- public double getA() {
- return a;
- }
+import java.util.Random;
- public double getOmega() {
- return omega;
- }
+import org.apache.commons.math.optimization.OptimizationException;
+import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
+import org.apache.commons.math.util.MathUtils;
+import org.junit.Test;
+
+public class HarmonicFitterTest {
+
+ @Test
+ public void testNoError() throws OptimizationException {
+ HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
+
+ HarmonicFitter fitter =
+ new HarmonicFitter(new LevenbergMarquardtOptimizer());
+ for (double x = 0.0; x < 1.3; x += 0.01) {
+ fitter.addObservedPoint(1.0, x, f.value(x));
+ }
+
+ HarmonicFunction fitted = fitter.fit();
+ assertEquals(f.getAmplitude(), fitted.getAmplitude(), 1.0e-13);
+ assertEquals(f.getPulsation(), fitted.getPulsation(), 1.0e-13);
+ assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.0e-13);
+
+ for (double x = -1.0; x < 1.0; x += 0.01) {
+ assertTrue(Math.abs(f.value(x) - fitted.value(x)) < 1.0e-13);
+ }
+
+ }
+
+ @Test
+ public void test1PercentError() throws OptimizationException {
+ Random randomizer = new Random(64925784252l);
+ HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
+
+ HarmonicFitter fitter =
+ new HarmonicFitter(new LevenbergMarquardtOptimizer());
+ for (double x = 0.0; x < 10.0; x += 0.1) {
+ fitter.addObservedPoint(1.0, x,
+ f.value(x) + 0.01 * randomizer.nextGaussian());
+ }
+
+ HarmonicFunction fitted = fitter.fit();
+ assertEquals(f.getAmplitude(), fitted.getAmplitude(), 7.6e-4);
+ assertEquals(f.getPulsation(), fitted.getPulsation(), 2.7e-3);
+ assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.3e-2);
+
+ }
+
+ @Test
+ public void testInitialGuess() throws OptimizationException {
+ Random randomizer = new Random(45314242l);
+ HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
+
+ HarmonicFitter fitter =
+ new HarmonicFitter(new LevenbergMarquardtOptimizer(), new double[] { 0.15, 3.6, 4.5 });
+ for (double x = 0.0; x < 10.0; x += 0.1) {
+ fitter.addObservedPoint(1.0, x,
+ f.value(x) + 0.01 * randomizer.nextGaussian());
+ }
+
+ HarmonicFunction fitted = fitter.fit();
+ assertEquals(f.getAmplitude(), fitted.getAmplitude(), 1.2e-3);
+ assertEquals(f.getPulsation(), fitted.getPulsation(), 3.3e-3);
+ assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.7e-2);
+
+ }
+
+ @Test
+ public void testUnsorted() throws OptimizationException {
+ Random randomizer = new Random(64925784252l);
+ HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
+
+ HarmonicFitter fitter =
+ new HarmonicFitter(new LevenbergMarquardtOptimizer());
+
+ // build a regularly spaced array of measurements
+ int size = 100;
+ double[] xTab = new double[size];
+ double[] yTab = new double[size];
+ for (int i = 0; i < size; ++i) {
+ xTab[i] = 0.1 * i;
+ yTab[i] = f.value(xTab[i]) + 0.01 * randomizer.nextGaussian();
+ }
+
+ // shake it
+ for (int i = 0; i < size; ++i) {
+ int i1 = randomizer.nextInt(size);
+ int i2 = randomizer.nextInt(size);
+ double xTmp = xTab[i1];
+ double yTmp = yTab[i1];
+ xTab[i1] = xTab[i2];
+ yTab[i1] = yTab[i2];
+ xTab[i2] = xTmp;
+ yTab[i2] = yTmp;
+ }
+
+ // pass it to the fitter
+ for (int i = 0; i < size; ++i) {
+ fitter.addObservedPoint(1.0, xTab[i], yTab[i]);
+ }
+
+ HarmonicFunction fitted = fitter.fit();
+ assertEquals(f.getAmplitude(), fitted.getAmplitude(), 7.6e-4);
+ assertEquals(f.getPulsation(), fitted.getPulsation(), 3.5e-3);
+ assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.5e-2);
- public double getPhi() {
- return phi;
}
- private double a;
- private double omega;
- private double phi;
-
- }
-
}
Copied: commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/PolynomialFitterTest.java (from r784713, commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/PolynomialFitterTest.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/PolynomialFitterTest.java?p2=commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/PolynomialFitterTest.java&p1=commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/PolynomialFitterTest.java&r1=784713&r2=786479&rev=786479&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/fitting/PolynomialFitterTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/fitting/PolynomialFitterTest.java Fri Jun 19 12:36:16 2009
@@ -15,150 +15,120 @@
// specific language governing permissions and limitations
// under the License.
-package org.spaceroots.mantissa.fitting;
+package org.apache.commons.math.optimization.fitting;
-import java.util.Random;
-import junit.framework.*;
-
-import org.spaceroots.mantissa.estimation.EstimationException;
-import org.spaceroots.mantissa.estimation.Estimator;
-import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
-import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
-
-public class PolynomialFitterTest
- extends TestCase {
-
- public PolynomialFitterTest(String name) {
- super(name);
- }
-
- public void testNoError()
- throws EstimationException {
- Random randomizer = new Random(64925784252l);
- for (int degree = 0; degree < 10; ++degree) {
- Polynom p = new Polynom(degree);
- for (int i = 0; i <= degree; ++i) {
- p.initCoeff (i, randomizer.nextGaussian());
- }
-
- PolynomialFitter fitter =
- new PolynomialFitter(degree, new LevenbergMarquardtEstimator());
- for (int i = 0; i <= degree; ++i) {
- fitter.addWeightedPair(1.0, i, p.valueAt(i));
- }
-
- Polynom fitted = new Polynom(fitter.fit());
-
- for (double x = -1.0; x < 1.0; x += 0.01) {
- double error = Math.abs(p.valueAt(x) - fitted.valueAt(x))
- / (1.0 + Math.abs(p.valueAt(x)));
- assertTrue(Math.abs(error) < 1.0e-5);
- }
-
- }
-
- }
-
- public void testSmallError()
- throws EstimationException {
- Random randomizer = new Random(53882150042l);
- for (int degree = 0; degree < 10; ++degree) {
- Polynom p = new Polynom(degree);
- for (int i = 0; i <= degree; ++i) {
- p.initCoeff(i, randomizer.nextGaussian());
- }
-
- PolynomialFitter fitter =
- new PolynomialFitter(degree, new LevenbergMarquardtEstimator());
- for (double x = -1.0; x < 1.0; x += 0.01) {
- fitter.addWeightedPair(1.0, x,
- p.valueAt(x) + 0.1 * randomizer.nextGaussian());
- }
-
- Polynom fitted = new Polynom(fitter.fit());
-
- for (double x = -1.0; x < 1.0; x += 0.01) {
- double error = Math.abs(p.valueAt(x) - fitted.valueAt(x))
- / (1.0 + Math.abs(p.valueAt(x)));
- assertTrue(Math.abs(error) < 0.1);
- }
- }
-
- }
-
- public void testRedundantSolvable() {
- // Levenberg-Marquardt should handle redundant information gracefully
- checkUnsolvableProblem(new LevenbergMarquardtEstimator(), true);
- }
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
- public void testRedundantUnsolvable() {
- // Gauss-Newton should not be able to solve redundant information
- checkUnsolvableProblem(new GaussNewtonEstimator(10, 1.0e-7, 1.0e-7,
- 1.0e-10),
- false);
- }
-
- private void checkUnsolvableProblem(Estimator estimator,
- boolean solvable) {
- Random randomizer = new Random(1248788532l);
- for (int degree = 0; degree < 10; ++degree) {
- Polynom p = new Polynom(degree);
- for (int i = 0; i <= degree; ++i) {
- p.initCoeff(i, randomizer.nextGaussian());
- }
-
- PolynomialFitter fitter = new PolynomialFitter(degree, estimator);
-
- // reusing the same point over and over again does not bring
- // information, the problem cannot be solved in this case for
- // degrees greater than 1 (but one point is sufficient for
- // degree 0)
- for (double x = -1.0; x < 1.0; x += 0.01) {
- fitter.addWeightedPair(1.0, 0.0, p.valueAt(0.0));
- }
-
- try {
- fitter.fit();
- assertTrue(solvable || (degree == 0));
- } catch(EstimationException e) {
- assertTrue((! solvable) && (degree > 0));
- }
-
- }
-
- }
-
- public static Test suite() {
- return new TestSuite(PolynomialFitterTest.class);
- }
-
- private static class Polynom {
-
- public Polynom(int degree) {
- coeffs = new double[degree + 1];
- for (int i = 0; i < coeffs.length; ++i) {
- coeffs[i] = 0.0;
- }
- }
-
- public Polynom(double[]coeffs) {
- this.coeffs = coeffs;
- }
-
- public void initCoeff(int i, double c) {
- coeffs[i] = c;
- }
+import java.util.Random;
- public double valueAt(double x) {
- double y = coeffs[coeffs.length - 1];
- for (int i = coeffs.length - 2; i >= 0; --i) {
- y = y * x + coeffs[i];
- }
- return y;
+import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
+import org.apache.commons.math.optimization.OptimizationException;
+import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
+import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
+import org.junit.Test;
+
+public class PolynomialFitterTest {
+
+ @Test
+ public void testNoError() throws OptimizationException {
+ Random randomizer = new Random(64925784252l);
+ for (int degree = 1; degree < 10; ++degree) {
+ PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
+
+ PolynomialFitter fitter =
+ new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
+ for (int i = 0; i <= degree; ++i) {
+ fitter.addObservedPoint(1.0, i, p.value(i));
+ }
+
+ PolynomialFunction fitted = fitter.fit();
+
+ for (double x = -1.0; x < 1.0; x += 0.01) {
+ double error = Math.abs(p.value(x) - fitted.value(x)) /
+ (1.0 + Math.abs(p.value(x)));
+ assertEquals(0.0, error, 1.0e-6);
+ }
+
+ }
+
+ }
+
+ @Test
+ public void testSmallError() throws OptimizationException {
+ Random randomizer = new Random(53882150042l);
+ double maxError = 0;
+ for (int degree = 0; degree < 10; ++degree) {
+ PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
+
+ PolynomialFitter fitter =
+ new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
+ for (double x = -1.0; x < 1.0; x += 0.01) {
+ fitter.addObservedPoint(1.0, x,
+ p.value(x) + 0.1 * randomizer.nextGaussian());
+ }
+
+ PolynomialFunction fitted = fitter.fit();
+
+ for (double x = -1.0; x < 1.0; x += 0.01) {
+ double error = Math.abs(p.value(x) - fitted.value(x)) /
+ (1.0 + Math.abs(p.value(x)));
+ maxError = Math.max(maxError, error);
+ assertTrue(Math.abs(error) < 0.1);
+ }
+ }
+ assertTrue(maxError > 0.01);
+
+ }
+
+ @Test
+ public void testRedundantSolvable() {
+ // Levenberg-Marquardt should handle redundant information gracefully
+ checkUnsolvableProblem(new LevenbergMarquardtOptimizer(), true);
+ }
+
+ @Test
+ public void testRedundantUnsolvable() {
+ // Gauss-Newton should not be able to solve redundant information
+ DifferentiableMultivariateVectorialOptimizer optimizer =
+ new GaussNewtonOptimizer(true);
+ checkUnsolvableProblem(optimizer, false);
+ }
+
+ private void checkUnsolvableProblem(DifferentiableMultivariateVectorialOptimizer optimizer,
+ boolean solvable) {
+ Random randomizer = new Random(1248788532l);
+ for (int degree = 0; degree < 10; ++degree) {
+ PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
+
+ PolynomialFitter fitter = new PolynomialFitter(degree, optimizer);
+
+ // reusing the same point over and over again does not bring
+ // information, the problem cannot be solved in this case for
+ // degrees greater than 1 (but one point is sufficient for
+ // degree 0)
+ for (double x = -1.0; x < 1.0; x += 0.01) {
+ fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
+ }
+
+ try {
+ fitter.fit();
+ assertTrue(solvable || (degree == 0));
+ } catch(OptimizationException e) {
+ 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);
}
- private double[] coeffs;
-
- }
-
}