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.Arrays;
020import java.util.Comparator;
021import org.apache.commons.math3.analysis.function.Gaussian;
022import org.apache.commons.math3.exception.NullArgumentException;
023import org.apache.commons.math3.exception.NumberIsTooSmallException;
024import org.apache.commons.math3.exception.OutOfRangeException;
025import org.apache.commons.math3.exception.ZeroException;
026import org.apache.commons.math3.exception.NotStrictlyPositiveException;
027import org.apache.commons.math3.exception.util.LocalizedFormats;
028import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
029import org.apache.commons.math3.util.FastMath;
030
031/**
032 * Fits points to a {@link
033 * org.apache.commons.math3.analysis.function.Gaussian.Parametric Gaussian} function.
034 * <p>
035 * Usage example:
036 * <pre>
037 *   GaussianFitter fitter = new GaussianFitter(
038 *     new LevenbergMarquardtOptimizer());
039 *   fitter.addObservedPoint(4.0254623,  531026.0);
040 *   fitter.addObservedPoint(4.03128248, 984167.0);
041 *   fitter.addObservedPoint(4.03839603, 1887233.0);
042 *   fitter.addObservedPoint(4.04421621, 2687152.0);
043 *   fitter.addObservedPoint(4.05132976, 3461228.0);
044 *   fitter.addObservedPoint(4.05326982, 3580526.0);
045 *   fitter.addObservedPoint(4.05779662, 3439750.0);
046 *   fitter.addObservedPoint(4.0636168,  2877648.0);
047 *   fitter.addObservedPoint(4.06943698, 2175960.0);
048 *   fitter.addObservedPoint(4.07525716, 1447024.0);
049 *   fitter.addObservedPoint(4.08237071, 717104.0);
050 *   fitter.addObservedPoint(4.08366408, 620014.0);
051 *   double[] parameters = fitter.fit();
052 * </pre>
053 *
054 * @since 2.2
055 * @version $Id: GaussianFitter.java 1416643 2012-12-03 19:37:14Z tn $
056 * @deprecated As of 3.3. Please use {@link GaussianCurveFitter} and
057 * {@link WeightedObservedPoints} instead.
058 */
059@Deprecated
060public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
061    /**
062     * Constructs an instance using the specified optimizer.
063     *
064     * @param optimizer Optimizer to use for the fitting.
065     */
066    public GaussianFitter(MultivariateVectorOptimizer optimizer) {
067        super(optimizer);
068    }
069
070    /**
071     * Fits a Gaussian function to the observed points.
072     *
073     * @param initialGuess First guess values in the following order:
074     * <ul>
075     *  <li>Norm</li>
076     *  <li>Mean</li>
077     *  <li>Sigma</li>
078     * </ul>
079     * @return the parameters of the Gaussian function that best fits the
080     * observed points (in the same order as above).
081     * @since 3.0
082     */
083    public double[] fit(double[] initialGuess) {
084        final Gaussian.Parametric f = new Gaussian.Parametric() {
085                @Override
086                public double value(double x, double ... p) {
087                    double v = Double.POSITIVE_INFINITY;
088                    try {
089                        v = super.value(x, p);
090                    } catch (NotStrictlyPositiveException e) { // NOPMD
091                        // Do nothing.
092                    }
093                    return v;
094                }
095
096                @Override
097                public double[] gradient(double x, double ... p) {
098                    double[] v = { Double.POSITIVE_INFINITY,
099                                   Double.POSITIVE_INFINITY,
100                                   Double.POSITIVE_INFINITY };
101                    try {
102                        v = super.gradient(x, p);
103                    } catch (NotStrictlyPositiveException e) { // NOPMD
104                        // Do nothing.
105                    }
106                    return v;
107                }
108            };
109
110        return fit(f, initialGuess);
111    }
112
113    /**
114     * Fits a Gaussian function to the observed points.
115     *
116     * @return the parameters of the Gaussian function that best fits the
117     * observed points (in the same order as above).
118     */
119    public double[] fit() {
120        final double[] guess = (new ParameterGuesser(getObservations())).guess();
121        return fit(guess);
122    }
123
124    /**
125     * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
126     * of a {@link org.apache.commons.math3.analysis.function.Gaussian.Parametric}
127     * based on the specified observed points.
128     */
129    public static class ParameterGuesser {
130        /** Normalization factor. */
131        private final double norm;
132        /** Mean. */
133        private final double mean;
134        /** Standard deviation. */
135        private final double sigma;
136
137        /**
138         * Constructs instance with the specified observed points.
139         *
140         * @param observations Observed points from which to guess the
141         * parameters of the Gaussian.
142         * @throws NullArgumentException if {@code observations} is
143         * {@code null}.
144         * @throws NumberIsTooSmallException if there are less than 3
145         * observations.
146         */
147        public ParameterGuesser(WeightedObservedPoint[] observations) {
148            if (observations == null) {
149                throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
150            }
151            if (observations.length < 3) {
152                throw new NumberIsTooSmallException(observations.length, 3, true);
153            }
154
155            final WeightedObservedPoint[] sorted = sortObservations(observations);
156            final double[] params = basicGuess(sorted);
157
158            norm = params[0];
159            mean = params[1];
160            sigma = params[2];
161        }
162
163        /**
164         * Gets an estimation of the parameters.
165         *
166         * @return the guessed parameters, in the following order:
167         * <ul>
168         *  <li>Normalization factor</li>
169         *  <li>Mean</li>
170         *  <li>Standard deviation</li>
171         * </ul>
172         */
173        public double[] guess() {
174            return new double[] { norm, mean, sigma };
175        }
176
177        /**
178         * Sort the observations.
179         *
180         * @param unsorted Input observations.
181         * @return the input observations, sorted.
182         */
183        private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
184            final WeightedObservedPoint[] observations = unsorted.clone();
185            final Comparator<WeightedObservedPoint> cmp
186                = new Comparator<WeightedObservedPoint>() {
187                public int compare(WeightedObservedPoint p1,
188                                   WeightedObservedPoint p2) {
189                    if (p1 == null && p2 == null) {
190                        return 0;
191                    }
192                    if (p1 == null) {
193                        return -1;
194                    }
195                    if (p2 == null) {
196                        return 1;
197                    }
198                    if (p1.getX() < p2.getX()) {
199                        return -1;
200                    }
201                    if (p1.getX() > p2.getX()) {
202                        return 1;
203                    }
204                    if (p1.getY() < p2.getY()) {
205                        return -1;
206                    }
207                    if (p1.getY() > p2.getY()) {
208                        return 1;
209                    }
210                    if (p1.getWeight() < p2.getWeight()) {
211                        return -1;
212                    }
213                    if (p1.getWeight() > p2.getWeight()) {
214                        return 1;
215                    }
216                    return 0;
217                }
218            };
219
220            Arrays.sort(observations, cmp);
221            return observations;
222        }
223
224        /**
225         * Guesses the parameters based on the specified observed points.
226         *
227         * @param points Observed points, sorted.
228         * @return the guessed parameters (normalization factor, mean and
229         * sigma).
230         */
231        private double[] basicGuess(WeightedObservedPoint[] points) {
232            final int maxYIdx = findMaxY(points);
233            final double n = points[maxYIdx].getY();
234            final double m = points[maxYIdx].getX();
235
236            double fwhmApprox;
237            try {
238                final double halfY = n + ((m - n) / 2);
239                final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
240                final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
241                fwhmApprox = fwhmX2 - fwhmX1;
242            } catch (OutOfRangeException e) {
243                // TODO: Exceptions should not be used for flow control.
244                fwhmApprox = points[points.length - 1].getX() - points[0].getX();
245            }
246            final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
247
248            return new double[] { n, m, s };
249        }
250
251        /**
252         * Finds index of point in specified points with the largest Y.
253         *
254         * @param points Points to search.
255         * @return the index in specified points array.
256         */
257        private int findMaxY(WeightedObservedPoint[] points) {
258            int maxYIdx = 0;
259            for (int i = 1; i < points.length; i++) {
260                if (points[i].getY() > points[maxYIdx].getY()) {
261                    maxYIdx = i;
262                }
263            }
264            return maxYIdx;
265        }
266
267        /**
268         * Interpolates using the specified points to determine X at the
269         * specified Y.
270         *
271         * @param points Points to use for interpolation.
272         * @param startIdx Index within points from which to start the search for
273         * interpolation bounds points.
274         * @param idxStep Index step for searching interpolation bounds points.
275         * @param y Y value for which X should be determined.
276         * @return the value of X for the specified Y.
277         * @throws ZeroException if {@code idxStep} is 0.
278         * @throws OutOfRangeException if specified {@code y} is not within the
279         * range of the specified {@code points}.
280         */
281        private double interpolateXAtY(WeightedObservedPoint[] points,
282                                       int startIdx,
283                                       int idxStep,
284                                       double y)
285            throws OutOfRangeException {
286            if (idxStep == 0) {
287                throw new ZeroException();
288            }
289            final WeightedObservedPoint[] twoPoints
290                = getInterpolationPointsForY(points, startIdx, idxStep, y);
291            final WeightedObservedPoint p1 = twoPoints[0];
292            final WeightedObservedPoint p2 = twoPoints[1];
293            if (p1.getY() == y) {
294                return p1.getX();
295            }
296            if (p2.getY() == y) {
297                return p2.getX();
298            }
299            return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
300                                (p2.getY() - p1.getY()));
301        }
302
303        /**
304         * Gets the two bounding interpolation points from the specified points
305         * suitable for determining X at the specified Y.
306         *
307         * @param points Points to use for interpolation.
308         * @param startIdx Index within points from which to start search for
309         * interpolation bounds points.
310         * @param idxStep Index step for search for interpolation bounds points.
311         * @param y Y value for which X should be determined.
312         * @return the array containing two points suitable for determining X at
313         * the specified Y.
314         * @throws ZeroException if {@code idxStep} is 0.
315         * @throws OutOfRangeException if specified {@code y} is not within the
316         * range of the specified {@code points}.
317         */
318        private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
319                                                                   int startIdx,
320                                                                   int idxStep,
321                                                                   double y)
322            throws OutOfRangeException {
323            if (idxStep == 0) {
324                throw new ZeroException();
325            }
326            for (int i = startIdx;
327                 idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
328                 i += idxStep) {
329                final WeightedObservedPoint p1 = points[i];
330                final WeightedObservedPoint p2 = points[i + idxStep];
331                if (isBetween(y, p1.getY(), p2.getY())) {
332                    if (idxStep < 0) {
333                        return new WeightedObservedPoint[] { p2, p1 };
334                    } else {
335                        return new WeightedObservedPoint[] { p1, p2 };
336                    }
337                }
338            }
339
340            // Boundaries are replaced by dummy values because the raised
341            // exception is caught and the message never displayed.
342            // TODO: Exceptions should not be used for flow control.
343            throw new OutOfRangeException(y,
344                                          Double.NEGATIVE_INFINITY,
345                                          Double.POSITIVE_INFINITY);
346        }
347
348        /**
349         * Determines whether a value is between two other values.
350         *
351         * @param value Value to test whether it is between {@code boundary1}
352         * and {@code boundary2}.
353         * @param boundary1 One end of the range.
354         * @param boundary2 Other end of the range.
355         * @return {@code true} if {@code value} is between {@code boundary1} and
356         * {@code boundary2} (inclusive), {@code false} otherwise.
357         */
358        private boolean isBetween(double value,
359                                  double boundary1,
360                                  double boundary2) {
361            return (value >= boundary1 && value <= boundary2) ||
362                (value >= boundary2 && value <= boundary1);
363        }
364    }
365}