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