View Javadoc
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  
18  package org.apache.commons.math3.optimization.fitting;
19  
20  import java.util.Arrays;
21  import java.util.Comparator;
22  
23  import org.apache.commons.math3.analysis.function.Gaussian;
24  import org.apache.commons.math3.exception.NullArgumentException;
25  import org.apache.commons.math3.exception.NumberIsTooSmallException;
26  import org.apache.commons.math3.exception.OutOfRangeException;
27  import org.apache.commons.math3.exception.ZeroException;
28  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
29  import org.apache.commons.math3.exception.util.LocalizedFormats;
30  import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
31  import org.apache.commons.math3.util.FastMath;
32  
33  /**
34   * Fits points to a {@link
35   * org.apache.commons.math3.analysis.function.Gaussian.Parametric Gaussian} function.
36   * <p>
37   * Usage example:
38   * <pre>
39   *   GaussianFitter fitter = new GaussianFitter(
40   *     new LevenbergMarquardtOptimizer());
41   *   fitter.addObservedPoint(4.0254623,  531026.0);
42   *   fitter.addObservedPoint(4.03128248, 984167.0);
43   *   fitter.addObservedPoint(4.03839603, 1887233.0);
44   *   fitter.addObservedPoint(4.04421621, 2687152.0);
45   *   fitter.addObservedPoint(4.05132976, 3461228.0);
46   *   fitter.addObservedPoint(4.05326982, 3580526.0);
47   *   fitter.addObservedPoint(4.05779662, 3439750.0);
48   *   fitter.addObservedPoint(4.0636168,  2877648.0);
49   *   fitter.addObservedPoint(4.06943698, 2175960.0);
50   *   fitter.addObservedPoint(4.07525716, 1447024.0);
51   *   fitter.addObservedPoint(4.08237071, 717104.0);
52   *   fitter.addObservedPoint(4.08366408, 620014.0);
53   *   double[] parameters = fitter.fit();
54   * </pre>
55   *
56   * @since 2.2
57   * @deprecated As of 3.1 (to be removed in 4.0).
58   */
59  @Deprecated
60  public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
61      /**
62       * Constructs an instance using the specified optimizer.
63       *
64       * @param optimizer Optimizer to use for the fitting.
65       */
66      public GaussianFitter(DifferentiableMultivariateVectorOptimizer optimizer) {
67          super(optimizer);
68      }
69  
70      /**
71       * Fits a Gaussian function to the observed points.
72       *
73       * @param initialGuess First guess values in the following order:
74       * <ul>
75       *  <li>Norm</li>
76       *  <li>Mean</li>
77       *  <li>Sigma</li>
78       * </ul>
79       * @return the parameters of the Gaussian function that best fits the
80       * observed points (in the same order as above).
81       * @since 3.0
82       */
83      public double[] fit(double[] initialGuess) {
84          final Gaussian.Parametric f = new Gaussian.Parametric() {
85                  @Override
86                  public double value(double x, double ... p) {
87                      double v = Double.POSITIVE_INFINITY;
88                      try {
89                          v = super.value(x, p);
90                      } catch (NotStrictlyPositiveException e) { // NOPMD
91                          // Do nothing.
92                      }
93                      return v;
94                  }
95  
96                  @Override
97                  public double[] gradient(double x, double ... p) {
98                      double[] v = { Double.POSITIVE_INFINITY,
99                                     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 }