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 * @since 2.0
049 * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
050 * {@link WeightedObservedPoints} instead.
051 */
052@Deprecated
053public class CurveFitter<T extends ParametricUnivariateFunction> {
054    /** Optimizer to use for the fitting. */
055    private final MultivariateVectorOptimizer optimizer;
056    /** Observed points. */
057    private final List<WeightedObservedPoint> observations;
058
059    /**
060     * Simple constructor.
061     *
062     * @param optimizer Optimizer to use for the fitting.
063     * @since 3.1
064     */
065    public CurveFitter(final MultivariateVectorOptimizer optimizer) {
066        this.optimizer = optimizer;
067        observations = new ArrayList<WeightedObservedPoint>();
068    }
069
070    /** Add an observed (x,y) point to the sample with unit weight.
071     * <p>Calling this method is equivalent to call
072     * {@code addObservedPoint(1.0, x, y)}.</p>
073     * @param x abscissa of the point
074     * @param y observed value of the point at x, after fitting we should
075     * have f(x) as close as possible to this value
076     * @see #addObservedPoint(double, double, double)
077     * @see #addObservedPoint(WeightedObservedPoint)
078     * @see #getObservations()
079     */
080    public void addObservedPoint(double x, double y) {
081        addObservedPoint(1.0, x, y);
082    }
083
084    /** Add an observed weighted (x,y) point to the sample.
085     * @param weight weight of the observed point in the fit
086     * @param x abscissa of the point
087     * @param y observed value of the point at x, after fitting we should
088     * have f(x) as close as possible to this value
089     * @see #addObservedPoint(double, double)
090     * @see #addObservedPoint(WeightedObservedPoint)
091     * @see #getObservations()
092     */
093    public void addObservedPoint(double weight, double x, double y) {
094        observations.add(new WeightedObservedPoint(weight, x, y));
095    }
096
097    /** Add an observed weighted (x,y) point to the sample.
098     * @param observed observed point to add
099     * @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    }
106
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    }
116
117    /**
118     * Remove all observations.
119     */
120    public void clearObservations() {
121        observations.clear();
122    }
123
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    }
140
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        }
169
170        // Input to the optimizer: the model and its Jacobian.
171        final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);
172
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    }
184
185    /** Vectorial function computing function theoretical values. */
186    private class TheoreticalValuesFunction {
187        /** Function to fit. */
188        private final ParametricUnivariateFunction f;
189
190        /**
191         * @param f function to fit.
192         */
193        TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
194            this.f = f;
195        }
196
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                        }
210
211                        return values;
212                    }
213                });
214        }
215
216        /**
217         * @return the model function Jacobian.
218         */
219        public ModelFunctionJacobian getModelFunctionJacobian() {
220            return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
221                    /** {@inheritDoc} */
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}