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