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 }