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    *      https://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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.fitting;
23  
24  import java.util.ArrayList;
25  import java.util.Collection;
26  import java.util.Collections;
27  import java.util.Comparator;
28  import java.util.List;
29  
30  import org.hipparchus.analysis.function.Gaussian;
31  import org.hipparchus.exception.LocalizedCoreFormats;
32  import org.hipparchus.exception.MathIllegalArgumentException;
33  import org.hipparchus.linear.DiagonalMatrix;
34  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresBuilder;
35  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem;
36  import org.hipparchus.util.FastMath;
37  import org.hipparchus.util.MathUtils;
38  
39  /**
40   * Fits points to a {@link
41   * org.hipparchus.analysis.function.Gaussian.Parametric Gaussian}
42   * function.
43   * <br>
44   * The {@link #withStartPoint(double[]) initial guess values} must be passed
45   * in the following order:
46   * <ul>
47   *  <li>Normalization</li>
48   *  <li>Mean</li>
49   *  <li>Sigma</li>
50   * </ul>
51   * The optimal values will be returned in the same order.
52   *
53   * <p>
54   * Usage example:
55   * <pre>
56   *   WeightedObservedPoints obs = new WeightedObservedPoints();
57   *   obs.add(4.0254623,  531026.0);
58   *   obs.add(4.03128248, 984167.0);
59   *   obs.add(4.03839603, 1887233.0);
60   *   obs.add(4.04421621, 2687152.0);
61   *   obs.add(4.05132976, 3461228.0);
62   *   obs.add(4.05326982, 3580526.0);
63   *   obs.add(4.05779662, 3439750.0);
64   *   obs.add(4.0636168,  2877648.0);
65   *   obs.add(4.06943698, 2175960.0);
66   *   obs.add(4.07525716, 1447024.0);
67   *   obs.add(4.08237071, 717104.0);
68   *   obs.add(4.08366408, 620014.0);
69   *   double[] parameters = GaussianCurveFitter.create().fit(obs.toList());
70   * </pre>
71   *
72   */
73  public class GaussianCurveFitter extends AbstractCurveFitter {
74      /** Parametric function to be fitted. */
75      private static final Gaussian.Parametric FUNCTION = new Gaussian.Parametric() {
76              /** {@inheritDoc} */
77              @Override
78              public double value(double x, double ... p) {
79                  double v = Double.POSITIVE_INFINITY;
80                  try {
81                      v = super.value(x, p);
82                  } catch (MathIllegalArgumentException e) { // NOPMD
83                      // Do nothing.
84                  }
85                  return v;
86              }
87  
88              /** {@inheritDoc} */
89              @Override
90              public double[] gradient(double x, double ... p) {
91                  double[] v = { Double.POSITIVE_INFINITY,
92                                 Double.POSITIVE_INFINITY,
93                                 Double.POSITIVE_INFINITY };
94                  try {
95                      v = super.gradient(x, p);
96                  } catch (MathIllegalArgumentException e) { // NOPMD
97                      // Do nothing.
98                  }
99                  return v;
100             }
101         };
102     /** Initial guess. */
103     private final double[] initialGuess;
104     /** Maximum number of iterations of the optimization algorithm. */
105     private final int maxIter;
106 
107     /**
108      * Constructor used by the factory methods.
109      *
110      * @param initialGuess Initial guess. If set to {@code null}, the initial guess
111      * will be estimated using the {@link ParameterGuesser}.
112      * @param maxIter Maximum number of iterations of the optimization algorithm.
113      */
114     private GaussianCurveFitter(double[] initialGuess, int maxIter) {
115         this.initialGuess = initialGuess == null ? null : initialGuess.clone();
116         this.maxIter = maxIter;
117     }
118 
119     /**
120      * Creates a default curve fitter.
121      * The initial guess for the parameters will be {@link ParameterGuesser}
122      * computed automatically, and the maximum number of iterations of the
123      * optimization algorithm is set to {@link Integer#MAX_VALUE}.
124      *
125      * @return a curve fitter.
126      *
127      * @see #withStartPoint(double[])
128      * @see #withMaxIterations(int)
129      */
130     public static GaussianCurveFitter create() {
131         return new GaussianCurveFitter(null, Integer.MAX_VALUE);
132     }
133 
134     /**
135      * Configure the start point (initial guess).
136      * @param newStart new start point (initial guess)
137      * @return a new instance.
138      */
139     public GaussianCurveFitter withStartPoint(double[] newStart) {
140         return new GaussianCurveFitter(newStart.clone(),
141                                        maxIter);
142     }
143 
144     /**
145      * Configure the maximum number of iterations.
146      * @param newMaxIter maximum number of iterations
147      * @return a new instance.
148      */
149     public GaussianCurveFitter withMaxIterations(int newMaxIter) {
150         return new GaussianCurveFitter(initialGuess,
151                                        newMaxIter);
152     }
153 
154     /** {@inheritDoc} */
155     @Override
156     protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> observations) {
157 
158         // Prepare least-squares problem.
159         final int len = observations.size();
160         final double[] target  = new double[len];
161         final double[] weights = new double[len];
162 
163         int i = 0;
164         for (WeightedObservedPoint obs : observations) {
165             target[i]  = obs.getY();
166             weights[i] = obs.getWeight();
167             ++i;
168         }
169 
170         final AbstractCurveFitter.TheoreticalValuesFunction model =
171                 new AbstractCurveFitter.TheoreticalValuesFunction(FUNCTION, observations);
172 
173         final double[] startPoint = initialGuess != null ?
174             initialGuess :
175             // Compute estimation.
176             new ParameterGuesser(observations).guess();
177 
178         // Return a new least squares problem set up to fit a Gaussian curve to the
179         // observed points.
180         return new LeastSquaresBuilder().
181                 maxEvaluations(Integer.MAX_VALUE).
182                 maxIterations(maxIter).
183                 start(startPoint).
184                 target(target).
185                 weight(new DiagonalMatrix(weights)).
186                 model(model.getModelFunction(), model.getModelFunctionJacobian()).
187                 build();
188 
189     }
190 
191     /**
192      * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
193      * of a {@link org.hipparchus.analysis.function.Gaussian.Parametric}
194      * based on the specified observed points.
195      */
196     public static class ParameterGuesser {
197         /** Normalization factor. */
198         private final double norm;
199         /** Mean. */
200         private final double mean;
201         /** Standard deviation. */
202         private final double sigma;
203 
204         /**
205          * Constructs instance with the specified observed points.
206          *
207          * @param observations Observed points from which to guess the
208          * parameters of the Gaussian.
209          * @throws org.hipparchus.exception.NullArgumentException if {@code observations} is
210          * {@code null}.
211          * @throws MathIllegalArgumentException if there are less than 3
212          * observations.
213          */
214         public ParameterGuesser(Collection<WeightedObservedPoint> observations) {
215             MathUtils.checkNotNull(observations);
216             if (observations.size() < 3) {
217                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
218                                                        observations.size(), 3);
219             }
220 
221             final List<WeightedObservedPoint> sorted = sortObservations(observations);
222             final double[] params = basicGuess(sorted.toArray(new WeightedObservedPoint[0]));
223 
224             norm = params[0];
225             mean = params[1];
226             sigma = params[2];
227         }
228 
229         /**
230          * Gets an estimation of the parameters.
231          *
232          * @return the guessed parameters, in the following order:
233          * <ul>
234          *  <li>Normalization factor</li>
235          *  <li>Mean</li>
236          *  <li>Standard deviation</li>
237          * </ul>
238          */
239         public double[] guess() {
240             return new double[] { norm, mean, sigma };
241         }
242 
243         /**
244          * Sort the observations.
245          *
246          * @param unsorted Input observations.
247          * @return the input observations, sorted.
248          */
249         private List<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) {
250             final List<WeightedObservedPoint> observations = new ArrayList<>(unsorted);
251 
252             final Comparator<WeightedObservedPoint> cmp = new Comparator<WeightedObservedPoint>() {
253                 /** {@inheritDoc} */
254                 @Override
255                 public int compare(WeightedObservedPoint p1,
256                                    WeightedObservedPoint p2) {
257                     if (p1 == null && p2 == null) {
258                         return 0;
259                     }
260                     if (p1 == null) {
261                         return -1;
262                     }
263                     if (p2 == null) {
264                         return 1;
265                     }
266                     int comp = Double.compare(p1.getX(), p2.getX());
267                     if (comp != 0) {
268                         return comp;
269                     }
270                     comp = Double.compare(p1.getY(), p2.getY());
271                     if (comp != 0) {
272                         return comp;
273                     }
274                     comp = Double.compare(p1.getWeight(), p2.getWeight());
275                     if (comp != 0) {
276                         return comp;
277                     }
278                     return 0;
279                 }
280             };
281 
282             Collections.sort(observations, cmp);
283             return observations;
284         }
285 
286         /**
287          * Guesses the parameters based on the specified observed points.
288          *
289          * @param points Observed points, sorted.
290          * @return the guessed parameters (normalization factor, mean and
291          * sigma).
292          */
293         private double[] basicGuess(WeightedObservedPoint[] points) {
294             final int maxYIdx = findMaxY(points);
295             final double n = points[maxYIdx].getY();
296             final double m = points[maxYIdx].getX();
297 
298             double fwhmApprox;
299             try {
300                 final double halfY = n + ((m - n) / 2);
301                 final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
302                 final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
303                 fwhmApprox = fwhmX2 - fwhmX1;
304             } catch (MathIllegalArgumentException e) {
305                 // TODO: Exceptions should not be used for flow control.
306                 fwhmApprox = points[points.length - 1].getX() - points[0].getX();
307             }
308             final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
309 
310             return new double[] { n, m, s };
311         }
312 
313         /**
314          * Finds index of point in specified points with the largest Y.
315          *
316          * @param points Points to search.
317          * @return the index in specified points array.
318          */
319         private int findMaxY(WeightedObservedPoint[] points) {
320             int maxYIdx = 0;
321             for (int i = 1; i < points.length; i++) {
322                 if (points[i].getY() > points[maxYIdx].getY()) {
323                     maxYIdx = i;
324                 }
325             }
326             return maxYIdx;
327         }
328 
329         /**
330          * Interpolates using the specified points to determine X at the
331          * specified Y.
332          *
333          * @param points Points to use for interpolation.
334          * @param startIdx Index within points from which to start the search for
335          * interpolation bounds points.
336          * @param idxStep Index step for searching interpolation bounds points.
337          * @param y Y value for which X should be determined.
338          * @return the value of X for the specified Y.
339          * @throws MathIllegalArgumentException if {@code idxStep} is 0.
340          * @throws MathIllegalArgumentException if specified {@code y} is not within the
341          * range of the specified {@code points}.
342          */
343         private double interpolateXAtY(WeightedObservedPoint[] points,
344                                        int startIdx,
345                                        int idxStep,
346                                        double y)
347             throws MathIllegalArgumentException {
348             if (idxStep == 0) {
349                 throw new MathIllegalArgumentException(LocalizedCoreFormats.ZERO_NOT_ALLOWED);
350             }
351             final WeightedObservedPoint[] twoPoints
352                 = getInterpolationPointsForY(points, startIdx, idxStep, y);
353             final WeightedObservedPoint p1 = twoPoints[0];
354             final WeightedObservedPoint p2 = twoPoints[1];
355             if (p1.getY() == y) {
356                 return p1.getX();
357             }
358             if (p2.getY() == y) {
359                 return p2.getX();
360             }
361             return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
362                                 (p2.getY() - p1.getY()));
363         }
364 
365         /**
366          * Gets the two bounding interpolation points from the specified points
367          * suitable for determining X at the specified Y.
368          *
369          * @param points Points to use for interpolation.
370          * @param startIdx Index within points from which to start search for
371          * interpolation bounds points.
372          * @param idxStep Index step for search for interpolation bounds points.
373          * @param y Y value for which X should be determined.
374          * @return the array containing two points suitable for determining X at
375          * the specified Y.
376          * @throws MathIllegalArgumentException if {@code idxStep} is 0.
377          * @throws MathIllegalArgumentException if specified {@code y} is not within the
378          * range of the specified {@code points}.
379          */
380         private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
381                                                                    int startIdx,
382                                                                    int idxStep,
383                                                                    double y)
384             throws MathIllegalArgumentException {
385             if (idxStep == 0) {
386                 throw new MathIllegalArgumentException(LocalizedCoreFormats.ZERO_NOT_ALLOWED);
387             }
388             for (int i = startIdx;
389                  idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
390                  i += idxStep) {
391                 final WeightedObservedPoint p1 = points[i];
392                 final WeightedObservedPoint p2 = points[i + idxStep];
393                 if (isBetween(y, p1.getY(), p2.getY())) {
394                     if (idxStep < 0) {
395                         return new WeightedObservedPoint[] { p2, p1 };
396                     } else {
397                         return new WeightedObservedPoint[] { p1, p2 };
398                     }
399                 }
400             }
401 
402             // Boundaries are replaced by dummy values because the raised
403             // exception is caught and the message never displayed.
404             // TODO: Exceptions should not be used for flow control.
405             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
406                                                    y, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
407         }
408 
409         /**
410          * Determines whether a value is between two other values.
411          *
412          * @param value Value to test whether it is between {@code boundary1}
413          * and {@code boundary2}.
414          * @param boundary1 One end of the range.
415          * @param boundary2 Other end of the range.
416          * @return {@code true} if {@code value} is between {@code boundary1} and
417          * {@code boundary2} (inclusive), {@code false} otherwise.
418          */
419         private boolean isBetween(double value,
420                                   double boundary1,
421                                   double boundary2) {
422             return (value >= boundary1 && value <= boundary2) ||
423                 (value >= boundary2 && value <= boundary1);
424         }
425     }
426 }