CurveFitter.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.fitting;

  19. import java.util.ArrayList;
  20. import java.util.List;
  21. import org.apache.commons.math3.analysis.MultivariateVectorFunction;
  22. import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
  23. import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
  24. import org.apache.commons.math3.optim.MaxEval;
  25. import org.apache.commons.math3.optim.InitialGuess;
  26. import org.apache.commons.math3.optim.PointVectorValuePair;
  27. import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
  28. import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
  29. import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
  30. import org.apache.commons.math3.optim.nonlinear.vector.Target;
  31. import org.apache.commons.math3.optim.nonlinear.vector.Weight;

  33. /**
  34.  * Fitter for parametric univariate real functions y = f(x).
  35.  * <br/>
  36.  * When a univariate real function y = f(x) does depend on some
  37.  * unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
  38.  * this class can be used to find these parameters. It does this
  39.  * by <em>fitting</em> the curve so it remains very close to a set of
  40.  * observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
  41.  * y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
  42.  * is done by finding the parameters values that minimizes the objective
  43.  * function &sum;(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
  44.  * really a least squares problem.
  45.  *
  46.  * @param <T> Function to use for the fit.
  47.  *
  48.  * @since 2.0
  49.  * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
  50.  * {@link WeightedObservedPoints} instead.
  51.  */
  52. @Deprecated
  53. public class CurveFitter<T extends ParametricUnivariateFunction> {
  54.     /** Optimizer to use for the fitting. */
  55.     private final MultivariateVectorOptimizer optimizer;
  56.     /** Observed points. */
  57.     private final List<WeightedObservedPoint> observations;

  59.     /**
  60.      * Simple constructor.
  61.      *
  62.      * @param optimizer Optimizer to use for the fitting.
  63.      * @since 3.1
  64.      */
  65.     public CurveFitter(final MultivariateVectorOptimizer optimizer) {
  66.         this.optimizer = optimizer;
  67.         observations = new ArrayList<WeightedObservedPoint>();
  68.     }

  70.     /** Add an observed (x,y) point to the sample with unit weight.
  71.      * <p>Calling this method is equivalent to call
  72.      * {@code addObservedPoint(1.0, x, y)}.</p>
  73.      * @param x abscissa of the point
  74.      * @param y observed value of the point at x, after fitting we should
  75.      * have f(x) as close as possible to this value
  76.      * @see #addObservedPoint(double, double, double)
  77.      * @see #addObservedPoint(WeightedObservedPoint)
  78.      * @see #getObservations()
  79.      */
  80.     public void addObservedPoint(double x, double y) {
  81.         addObservedPoint(1.0, x, y);
  82.     }

  84.     /** Add an observed weighted (x,y) point to the sample.
  85.      * @param weight weight of the observed point in the fit
  86.      * @param x abscissa of the point
  87.      * @param y observed value of the point at x, after fitting we should
  88.      * have f(x) as close as possible to this value
  89.      * @see #addObservedPoint(double, double)
  90.      * @see #addObservedPoint(WeightedObservedPoint)
  91.      * @see #getObservations()
  92.      */
  93.     public void addObservedPoint(double weight, double x, double y) {
  94.         observations.add(new WeightedObservedPoint(weight, x, y));
  95.     }

  97.     /** Add an observed weighted (x,y) point to the sample.
  98.      * @param observed observed point to add
  99.      * @see #addObservedPoint(double, double)
  100.      * @see #addObservedPoint(double, double, double)
  101.      * @see #getObservations()
  102.      */
  103.     public void addObservedPoint(WeightedObservedPoint observed) {
  104.         observations.add(observed);
  105.     }

  107.     /** Get the observed points.
  108.      * @return observed points
  109.      * @see #addObservedPoint(double, double)
  110.      * @see #addObservedPoint(double, double, double)
  111.      * @see #addObservedPoint(WeightedObservedPoint)
  112.      */
  113.     public WeightedObservedPoint[] getObservations() {
  114.         return observations.toArray(new WeightedObservedPoint[observations.size()]);
  115.     }

  117.     /**
  118.      * Remove all observations.
  119.      */
  120.     public void clearObservations() {
  121.         observations.clear();
  122.     }

  124.     /**
  125.      * Fit a curve.
  126.      * This method compute the coefficients of the curve that best
  127.      * fit the sample of observed points previously given through calls
  128.      * to the {@link #addObservedPoint(WeightedObservedPoint)
  129.      * addObservedPoint} method.
  130.      *
  131.      * @param f parametric function to fit.
  132.      * @param initialGuess first guess of the function parameters.
  133.      * @return the fitted parameters.
  134.      * @throws org.apache.commons.math3.exception.DimensionMismatchException
  135.      * if the start point dimension is wrong.
  136.      */
  137.     public double[] fit(T f, final double[] initialGuess) {
  138.         return fit(Integer.MAX_VALUE, f, initialGuess);
  139.     }

  141.     /**
  142.      * Fit a curve.
  143.      * This method compute the coefficients of the curve that best
  144.      * fit the sample of observed points previously given through calls
  145.      * to the {@link #addObservedPoint(WeightedObservedPoint)
  146.      * addObservedPoint} method.
  147.      *
  148.      * @param f parametric function to fit.
  149.      * @param initialGuess first guess of the function parameters.
  150.      * @param maxEval Maximum number of function evaluations.
  151.      * @return the fitted parameters.
  152.      * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
  153.      * if the number of allowed evaluations is exceeded.
  154.      * @throws org.apache.commons.math3.exception.DimensionMismatchException
  155.      * if the start point dimension is wrong.
  156.      * @since 3.0
  157.      */
  158.     public double[] fit(int maxEval, T f,
  159.                         final double[] initialGuess) {
  160.         // Prepare least squares problem.
  161.         double[] target  = new double[observations.size()];
  162.         double[] weights = new double[observations.size()];
  163.         int i = 0;
  164.         for (WeightedObservedPoint point : observations) {
  165.             target[i]  = point.getY();
  166.             weights[i] = point.getWeight();
  167.             ++i;
  168.         }

  170.         // Input to the optimizer: the model and its Jacobian.
  171.         final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);

  173.         // Perform the fit.
  174.         final PointVectorValuePair optimum
  175.             = optimizer.optimize(new MaxEval(maxEval),
  176.                                  model.getModelFunction(),
  177.                                  model.getModelFunctionJacobian(),
  178.                                  new Target(target),
  179.                                  new Weight(weights),
  180.                                  new InitialGuess(initialGuess));
  181.         // Extract the coefficients.
  182.         return optimum.getPointRef();
  183.     }

  185.     /** Vectorial function computing function theoretical values. */
  186.     private class TheoreticalValuesFunction {
  187.         /** Function to fit. */
  188.         private final ParametricUnivariateFunction f;

  190.         /**
  191.          * @param f function to fit.
  192.          */
  193.         public TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
  194.             this.f = f;
  195.         }

  197.         /**
  198.          * @return the model function values.
  199.          */
  200.         public ModelFunction getModelFunction() {
  201.             return new ModelFunction(new MultivariateVectorFunction() {
  202.                     /** {@inheritDoc} */
  203.                     public double[] value(double[] point) {
  204.                         // compute the residuals
  205.                         final double[] values = new double[observations.size()];
  206.                         int i = 0;
  207.                         for (WeightedObservedPoint observed : observations) {
  208.                             values[i++] = f.value(observed.getX(), point);
  209.                         }

  211.                         return values;
  212.                     }
  213.                 });
  214.         }

  216.         /**
  217.          * @return the model function Jacobian.
  218.          */
  219.         public ModelFunctionJacobian getModelFunctionJacobian() {
  220.             return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
  221.                     public double[][] value(double[] point) {
  222.                         final double[][] jacobian = new double[observations.size()][];
  223.                         int i = 0;
  224.                         for (WeightedObservedPoint observed : observations) {
  225.                             jacobian[i++] = f.gradient(observed.getX(), point);
  226.                         }
  227.                         return jacobian;
  228.                     }
  229.                 });
  230.         }
  231.     }
  232. }
Created with JaCoCo 0.7.2.201409121644