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 * @version $Id: GaussianFitter.java 1416643 2012-12-03 19:37:14Z tn $
56 */
57 public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
58 /**
59 * Constructs an instance using the specified optimizer.
60 *
61 * @param optimizer Optimizer to use for the fitting.
62 */
63 public GaussianFitter(MultivariateVectorOptimizer optimizer) {
64 super(optimizer);
65 }
66
67 /**
68 * Fits a Gaussian function to the observed points.
69 *
70 * @param initialGuess First guess values in the following order:
71 * <ul>
72 * <li>Norm</li>
73 * <li>Mean</li>
74 * <li>Sigma</li>
75 * </ul>
76 * @return the parameters of the Gaussian function that best fits the
77 * observed points (in the same order as above).
78 * @since 3.0
79 */
80 public double[] fit(double[] initialGuess) {
81 final Gaussian.Parametric f = new Gaussian.Parametric() {
82 @Override
83 public double value(double x, double ... p) {
84 double v = Double.POSITIVE_INFINITY;
85 try {
86 v = super.value(x, p);
87 } catch (NotStrictlyPositiveException e) { // NOPMD
88 // Do nothing.
89 }
90 return v;
91 }
92
93 @Override
94 public double[] gradient(double x, double ... p) {
95 double[] v = { Double.POSITIVE_INFINITY,
96 Double.POSITIVE_INFINITY,
97 Double.POSITIVE_INFINITY };
98 try {
99 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 }