commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject 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...
Date Fri, 19 Jun 2009 12:36:17 GMT
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 &sum;(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 &omega; 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 &omega;. */
+    private double omega;
+
+    /** Guessed phase &phi;. */
+    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 &phi; 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 &omega; 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 &phi; 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 &omega;.
+     * @return guessed pulsation &omega;
+     */
+    public double getGuessedPulsation() {
+        return omega;
+    }
 
-  private static final long serialVersionUID = 2400399048702758814L;
+    /** Get the guessed phase &phi;.
+     * @return guessed phase &phi;
+     */
+    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 &omega; and
+ * the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
  * searched by a least square estimator initialized with a rough guess
  * based on integrals.</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 &omega; and phase &phi;. */
+    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 &omega; (index 1) and phase &phi; (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 (&omega; t + &phi;)</code>.
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class HarmonicFunction implements DifferentiableUnivariateRealFunction {
+
+    /** Amplitude a. */
+    private final double a;
+
+    /** Pulsation &omega;. */
+    private final double omega;
+
+    /** Phase &phi;. */
+    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 &omega;.
+     * @return pulsation &omega;
+     */
+    public double getPulsation() {
+        return omega;
+    }
+
+    /** Get the phase &phi;.
+     * @return phase &phi;
+     */
+    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;
-
-  }
-
 }



Mime
View raw message