001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.fitting;
018
019import java.util.ArrayList;
020import java.util.List;
021import org.apache.commons.math3.analysis.MultivariateVectorFunction;
022import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
023import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
024import org.apache.commons.math3.optim.MaxEval;
025import org.apache.commons.math3.optim.InitialGuess;
026import org.apache.commons.math3.optim.PointVectorValuePair;
027import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
028import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
029import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
030import org.apache.commons.math3.optim.nonlinear.vector.Target;
031import org.apache.commons.math3.optim.nonlinear.vector.Weight;
032
033/**
034 * Fitter for parametric univariate real functions y = f(x).
035 * <br/>
036 * When a univariate real function y = f(x) does depend on some
037 * unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
038 * this class can be used to find these parameters. It does this
039 * by <em>fitting</em> the curve so it remains very close to a set of
040 * observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
041 * y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
042 * is done by finding the parameters values that minimizes the objective
043 * function &sum;(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
044 * really a least squares problem.
045 *
046 * @param <T> Function to use for the fit.
047 *
048 * @version $Id: CurveFitter.java 1416643 2012-12-03 19:37:14Z tn $
049 * @since 2.0
050 * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
051 * {@link WeightedObservedPoints} instead.
052 */
053@Deprecated
054public class CurveFitter<T extends ParametricUnivariateFunction> {
055    /** Optimizer to use for the fitting. */
056    private final MultivariateVectorOptimizer optimizer;
057    /** Observed points. */
058    private final List<WeightedObservedPoint> observations;
059
060    /**
061     * Simple constructor.
062     *
063     * @param optimizer Optimizer to use for the fitting.
064     * @since 3.1
065     */
066    public CurveFitter(final MultivariateVectorOptimizer optimizer) {
067        this.optimizer = optimizer;
068        observations = new ArrayList<WeightedObservedPoint>();
069    }
070
071    /** Add an observed (x,y) point to the sample with unit weight.
072     * <p>Calling this method is equivalent to call
073     * {@code addObservedPoint(1.0, x, y)}.</p>
074     * @param x abscissa of the point
075     * @param y observed value of the point at x, after fitting we should
076     * have f(x) as close as possible to this value
077     * @see #addObservedPoint(double, double, double)
078     * @see #addObservedPoint(WeightedObservedPoint)
079     * @see #getObservations()
080     */
081    public void addObservedPoint(double x, double y) {
082        addObservedPoint(1.0, x, y);
083    }
084
085    /** Add an observed weighted (x,y) point to the sample.
086     * @param weight weight of the observed point in the fit
087     * @param x abscissa of the point
088     * @param y observed value of the point at x, after fitting we should
089     * have f(x) as close as possible to this value
090     * @see #addObservedPoint(double, double)
091     * @see #addObservedPoint(WeightedObservedPoint)
092     * @see #getObservations()
093     */
094    public void addObservedPoint(double weight, double x, double y) {
095        observations.add(new WeightedObservedPoint(weight, x, y));
096    }
097
098    /** Add an observed weighted (x,y) point to the sample.
099     * @param observed observed point to add
100     * @see #addObservedPoint(double, double)
101     * @see #addObservedPoint(double, double, double)
102     * @see #getObservations()
103     */
104    public void addObservedPoint(WeightedObservedPoint observed) {
105        observations.add(observed);
106    }
107
108    /** Get the observed points.
109     * @return observed points
110     * @see #addObservedPoint(double, double)
111     * @see #addObservedPoint(double, double, double)
112     * @see #addObservedPoint(WeightedObservedPoint)
113     */
114    public WeightedObservedPoint[] getObservations() {
115        return observations.toArray(new WeightedObservedPoint[observations.size()]);
116    }
117
118    /**
119     * Remove all observations.
120     */
121    public void clearObservations() {
122        observations.clear();
123    }
124
125    /**
126     * Fit a curve.
127     * This method compute the coefficients of the curve that best
128     * fit the sample of observed points previously given through calls
129     * to the {@link #addObservedPoint(WeightedObservedPoint)
130     * addObservedPoint} method.
131     *
132     * @param f parametric function to fit.
133     * @param initialGuess first guess of the function parameters.
134     * @return the fitted parameters.
135     * @throws org.apache.commons.math3.exception.DimensionMismatchException
136     * if the start point dimension is wrong.
137     */
138    public double[] fit(T f, final double[] initialGuess) {
139        return fit(Integer.MAX_VALUE, f, initialGuess);
140    }
141
142    /**
143     * Fit a curve.
144     * This method compute the coefficients of the curve that best
145     * fit the sample of observed points previously given through calls
146     * to the {@link #addObservedPoint(WeightedObservedPoint)
147     * addObservedPoint} method.
148     *
149     * @param f parametric function to fit.
150     * @param initialGuess first guess of the function parameters.
151     * @param maxEval Maximum number of function evaluations.
152     * @return the fitted parameters.
153     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
154     * if the number of allowed evaluations is exceeded.
155     * @throws org.apache.commons.math3.exception.DimensionMismatchException
156     * if the start point dimension is wrong.
157     * @since 3.0
158     */
159    public double[] fit(int maxEval, T f,
160                        final double[] initialGuess) {
161        // Prepare least squares problem.
162        double[] target  = new double[observations.size()];
163        double[] weights = new double[observations.size()];
164        int i = 0;
165        for (WeightedObservedPoint point : observations) {
166            target[i]  = point.getY();
167            weights[i] = point.getWeight();
168            ++i;
169        }
170
171        // Input to the optimizer: the model and its Jacobian.
172        final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);
173
174        // Perform the fit.
175        final PointVectorValuePair optimum
176            = optimizer.optimize(new MaxEval(maxEval),
177                                 model.getModelFunction(),
178                                 model.getModelFunctionJacobian(),
179                                 new Target(target),
180                                 new Weight(weights),
181                                 new InitialGuess(initialGuess));
182        // Extract the coefficients.
183        return optimum.getPointRef();
184    }
185
186    /** Vectorial function computing function theoretical values. */
187    private class TheoreticalValuesFunction {
188        /** Function to fit. */
189        private final ParametricUnivariateFunction f;
190
191        /**
192         * @param f function to fit.
193         */
194        public TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
195            this.f = f;
196        }
197
198        /**
199         * @return the model function values.
200         */
201        public ModelFunction getModelFunction() {
202            return new ModelFunction(new MultivariateVectorFunction() {
203                    /** {@inheritDoc} */
204                    public double[] value(double[] point) {
205                        // compute the residuals
206                        final double[] values = new double[observations.size()];
207                        int i = 0;
208                        for (WeightedObservedPoint observed : observations) {
209                            values[i++] = f.value(observed.getX(), point);
210                        }
211
212                        return values;
213                    }
214                });
215        }
216
217        /**
218         * @return the model function Jacobian.
219         */
220        public ModelFunctionJacobian getModelFunctionJacobian() {
221            return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
222                    public double[][] value(double[] point) {
223                        final double[][] jacobian = new double[observations.size()][];
224                        int i = 0;
225                        for (WeightedObservedPoint observed : observations) {
226                            jacobian[i++] = f.gradient(observed.getX(), point);
227                        }
228                        return jacobian;
229                    }
230                });
231        }
232    }
233}