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  package org.apache.commons.math3.fitting;
18  
19  import java.util.Arrays;
20  import java.util.Comparator;
21  import org.apache.commons.math3.analysis.function.Gaussian;
22  import org.apache.commons.math3.exception.NullArgumentException;
23  import org.apache.commons.math3.exception.NumberIsTooSmallException;
24  import org.apache.commons.math3.exception.OutOfRangeException;
25  import org.apache.commons.math3.exception.ZeroException;
26  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
27  import org.apache.commons.math3.exception.util.LocalizedFormats;
28  import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
29  import org.apache.commons.math3.util.FastMath;
30  
31  /**
32   * Fits points to a {@link
33   * org.apache.commons.math3.analysis.function.Gaussian.Parametric Gaussian} function.
34   * <p>
35   * Usage example:
36   * <pre>
37   *   GaussianFitter fitter = new GaussianFitter(
38   *     new LevenbergMarquardtOptimizer());
39   *   fitter.addObservedPoint(4.0254623,  531026.0);
40   *   fitter.addObservedPoint(4.03128248, 984167.0);
41   *   fitter.addObservedPoint(4.03839603, 1887233.0);
42   *   fitter.addObservedPoint(4.04421621, 2687152.0);
43   *   fitter.addObservedPoint(4.05132976, 3461228.0);
44   *   fitter.addObservedPoint(4.05326982, 3580526.0);
45   *   fitter.addObservedPoint(4.05779662, 3439750.0);
46   *   fitter.addObservedPoint(4.0636168,  2877648.0);
47   *   fitter.addObservedPoint(4.06943698, 2175960.0);
48   *   fitter.addObservedPoint(4.07525716, 1447024.0);
49   *   fitter.addObservedPoint(4.08237071, 717104.0);
50   *   fitter.addObservedPoint(4.08366408, 620014.0);
51   *   double[] parameters = fitter.fit();
52   * </pre>
53   *
54   * @since 2.2
55   * @deprecated As of 3.3. Please use {@link GaussianCurveFitter} and
56   * {@link WeightedObservedPoints} instead.
57   */
58  @Deprecated
59  public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
60      /**
61       * Constructs an instance using the specified optimizer.
62       *
63       * @param optimizer Optimizer to use for the fitting.
64       */
65      public GaussianFitter(MultivariateVectorOptimizer optimizer) {
66          super(optimizer);
67      }
68  
69      /**
70       * Fits a Gaussian function to the observed points.
71       *
72       * @param initialGuess First guess values in the following order:
73       * <ul>
74       *  <li>Norm</li>
75       *  <li>Mean</li>
76       *  <li>Sigma</li>
77       * </ul>
78       * @return the parameters of the Gaussian function that best fits the
79       * observed points (in the same order as above).
80       * @since 3.0
81       */
82      public double[] fit(double[] initialGuess) {
83          final Gaussian.Parametric f = new Gaussian.Parametric() {
84                  @Override
85                  public double value(double x, double ... p) {
86                      double v = Double.POSITIVE_INFINITY;
87                      try {
88                          v = super.value(x, p);
89                      } catch (NotStrictlyPositiveException e) { // NOPMD
90                          // Do nothing.
91                      }
92                      return v;
93                  }
94  
95                  @Override
96                  public double[] gradient(double x, double ... p) {
97                      double[] v = { Double.POSITIVE_INFINITY,
98                                     Double.POSITIVE_INFINITY,
99                                     Double.POSITIVE_INFINITY };
100                     try {
101                         v = super.gradient(x, p);
102                     } catch (NotStrictlyPositiveException e) { // NOPMD
103                         // Do nothing.
104                     }
105                     return v;
106                 }
107             };
108 
109         return fit(f, initialGuess);
110     }
111 
112     /**
113      * Fits a Gaussian function to the observed points.
114      *
115      * @return the parameters of the Gaussian function that best fits the
116      * observed points (in the same order as above).
117      */
118     public double[] fit() {
119         final double[] guess = (new ParameterGuesser(getObservations())).guess();
120         return fit(guess);
121     }
122 
123     /**
124      * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
125      * of a {@link org.apache.commons.math3.analysis.function.Gaussian.Parametric}
126      * based on the specified observed points.
127      */
128     public static class ParameterGuesser {
129         /** Normalization factor. */
130         private final double norm;
131         /** Mean. */
132         private final double mean;
133         /** Standard deviation. */
134         private final double sigma;
135 
136         /**
137          * Constructs instance with the specified observed points.
138          *
139          * @param observations Observed points from which to guess the
140          * parameters of the Gaussian.
141          * @throws NullArgumentException if {@code observations} is
142          * {@code null}.
143          * @throws NumberIsTooSmallException if there are less than 3
144          * observations.
145          */
146         public ParameterGuesser(WeightedObservedPoint[] observations) {
147             if (observations == null) {
148                 throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
149             }
150             if (observations.length < 3) {
151                 throw new NumberIsTooSmallException(observations.length, 3, true);
152             }
153 
154             final WeightedObservedPoint[] sorted = sortObservations(observations);
155             final double[] params = basicGuess(sorted);
156 
157             norm = params[0];
158             mean = params[1];
159             sigma = params[2];
160         }
161 
162         /**
163          * Gets an estimation of the parameters.
164          *
165          * @return the guessed parameters, in the following order:
166          * <ul>
167          *  <li>Normalization factor</li>
168          *  <li>Mean</li>
169          *  <li>Standard deviation</li>
170          * </ul>
171          */
172         public double[] guess() {
173             return new double[] { norm, mean, sigma };
174         }
175 
176         /**
177          * Sort the observations.
178          *
179          * @param unsorted Input observations.
180          * @return the input observations, sorted.
181          */
182         private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
183             final WeightedObservedPoint[] observations = unsorted.clone();
184             final Comparator<WeightedObservedPoint> cmp
185                 = new Comparator<WeightedObservedPoint>() {
186                 public int compare(WeightedObservedPoint p1,
187                                    WeightedObservedPoint p2) {
188                     if (p1 == null && p2 == null) {
189                         return 0;
190                     }
191                     if (p1 == null) {
192                         return -1;
193                     }
194                     if (p2 == null) {
195                         return 1;
196                     }
197                     if (p1.getX() < p2.getX()) {
198                         return -1;
199                     }
200                     if (p1.getX() > p2.getX()) {
201                         return 1;
202                     }
203                     if (p1.getY() < p2.getY()) {
204                         return -1;
205                     }
206                     if (p1.getY() > p2.getY()) {
207                         return 1;
208                     }
209                     if (p1.getWeight() < p2.getWeight()) {
210                         return -1;
211                     }
212                     if (p1.getWeight() > p2.getWeight()) {
213                         return 1;
214                     }
215                     return 0;
216                 }
217             };
218 
219             Arrays.sort(observations, cmp);
220             return observations;
221         }
222 
223         /**
224          * Guesses the parameters based on the specified observed points.
225          *
226          * @param points Observed points, sorted.
227          * @return the guessed parameters (normalization factor, mean and
228          * sigma).
229          */
230         private double[] basicGuess(WeightedObservedPoint[] points) {
231             final int maxYIdx = findMaxY(points);
232             final double n = points[maxYIdx].getY();
233             final double m = points[maxYIdx].getX();
234 
235             double fwhmApprox;
236             try {
237                 final double halfY = n + ((m - n) / 2);
238                 final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
239                 final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
240                 fwhmApprox = fwhmX2 - fwhmX1;
241             } catch (OutOfRangeException e) {
242                 // TODO: Exceptions should not be used for flow control.
243                 fwhmApprox = points[points.length - 1].getX() - points[0].getX();
244             }
245             final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
246 
247             return new double[] { n, m, s };
248         }
249 
250         /**
251          * Finds index of point in specified points with the largest Y.
252          *
253          * @param points Points to search.
254          * @return the index in specified points array.
255          */
256         private int findMaxY(WeightedObservedPoint[] points) {
257             int maxYIdx = 0;
258             for (int i = 1; i < points.length; i++) {
259                 if (points[i].getY() > points[maxYIdx].getY()) {
260                     maxYIdx = i;
261                 }
262             }
263             return maxYIdx;
264         }
265 
266         /**
267          * Interpolates using the specified points to determine X at the
268          * specified Y.
269          *
270          * @param points Points to use for interpolation.
271          * @param startIdx Index within points from which to start the search for
272          * interpolation bounds points.
273          * @param idxStep Index step for searching interpolation bounds points.
274          * @param y Y value for which X should be determined.
275          * @return the value of X for the specified Y.
276          * @throws ZeroException if {@code idxStep} is 0.
277          * @throws OutOfRangeException if specified {@code y} is not within the
278          * range of the specified {@code points}.
279          */
280         private double interpolateXAtY(WeightedObservedPoint[] points,
281                                        int startIdx,
282                                        int idxStep,
283                                        double y)
284             throws OutOfRangeException {
285             if (idxStep == 0) {
286                 throw new ZeroException();
287             }
288             final WeightedObservedPoint[] twoPoints
289                 = getInterpolationPointsForY(points, startIdx, idxStep, y);
290             final WeightedObservedPoint p1 = twoPoints[0];
291             final WeightedObservedPoint p2 = twoPoints[1];
292             if (p1.getY() == y) {
293                 return p1.getX();
294             }
295             if (p2.getY() == y) {
296                 return p2.getX();
297             }
298             return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
299                                 (p2.getY() - p1.getY()));
300         }
301 
302         /**
303          * Gets the two bounding interpolation points from the specified points
304          * suitable for determining X at the specified Y.
305          *
306          * @param points Points to use for interpolation.
307          * @param startIdx Index within points from which to start search for
308          * interpolation bounds points.
309          * @param idxStep Index step for search for interpolation bounds points.
310          * @param y Y value for which X should be determined.
311          * @return the array containing two points suitable for determining X at
312          * the specified Y.
313          * @throws ZeroException if {@code idxStep} is 0.
314          * @throws OutOfRangeException if specified {@code y} is not within the
315          * range of the specified {@code points}.
316          */
317         private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
318                                                                    int startIdx,
319                                                                    int idxStep,
320                                                                    double y)
321             throws OutOfRangeException {
322             if (idxStep == 0) {
323                 throw new ZeroException();
324             }
325             for (int i = startIdx;
326                  idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
327                  i += idxStep) {
328                 final WeightedObservedPoint p1 = points[i];
329                 final WeightedObservedPoint p2 = points[i + idxStep];
330                 if (isBetween(y, p1.getY(), p2.getY())) {
331                     if (idxStep < 0) {
332                         return new WeightedObservedPoint[] { p2, p1 };
333                     } else {
334                         return new WeightedObservedPoint[] { p1, p2 };
335                     }
336                 }
337             }
338 
339             // Boundaries are replaced by dummy values because the raised
340             // exception is caught and the message never displayed.
341             // TODO: Exceptions should not be used for flow control.
342             throw new OutOfRangeException(y,
343                                           Double.NEGATIVE_INFINITY,
344                                           Double.POSITIVE_INFINITY);
345         }
346 
347         /**
348          * Determines whether a value is between two other values.
349          *
350          * @param value Value to test whether it is between {@code boundary1}
351          * and {@code boundary2}.
352          * @param boundary1 One end of the range.
353          * @param boundary2 Other end of the range.
354          * @return {@code true} if {@code value} is between {@code boundary1} and
355          * {@code boundary2} (inclusive), {@code false} otherwise.
356          */
357         private boolean isBetween(double value,
358                                   double boundary1,
359                                   double boundary2) {
360             return (value >= boundary1 && value <= boundary2) ||
361                 (value >= boundary2 && value <= boundary1);
362         }
363     }
364 }