source: src/main/java/agents/anac/y2019/harddealer/math3/fitting/GaussianFitter.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 14.6 KB
Line 
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 */
17package agents.anac.y2019.harddealer.math3.fitting;
18
19import java.util.Arrays;
20import java.util.Comparator;
21import agents.anac.y2019.harddealer.math3.analysis.function.Gaussian;
22import agents.anac.y2019.harddealer.math3.exception.NullArgumentException;
23import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
24import agents.anac.y2019.harddealer.math3.exception.OutOfRangeException;
25import agents.anac.y2019.harddealer.math3.exception.ZeroException;
26import agents.anac.y2019.harddealer.math3.exception.NotStrictlyPositiveException;
27import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
28import agents.anac.y2019.harddealer.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
29import agents.anac.y2019.harddealer.math3.util.FastMath;
30
31/**
32 * Fits points to a {@link
33 * agents.anac.y2019.harddealer.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
59public 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 /** {@inheritDoc} */
85 @Override
86 public double value(double x, double ... p) {
87 double v = Double.POSITIVE_INFINITY;
88 try {
89 v = super.value(x, p);
90 } catch (NotStrictlyPositiveException e) { // NOPMD
91 // Do nothing.
92 }
93 return v;
94 }
95
96 /** {@inheritDoc} */
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 agents.anac.y2019.harddealer.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 /** {@inheritDoc} */
189 public int compare(WeightedObservedPoint p1,
190 WeightedObservedPoint p2) {
191 if (p1 == null && p2 == null) {
192 return 0;
193 }
194 if (p1 == null) {
195 return -1;
196 }
197 if (p2 == null) {
198 return 1;
199 }
200 final int cmpX = Double.compare(p1.getX(), p2.getX());
201 if (cmpX < 0) {
202 return -1;
203 }
204 if (cmpX > 0) {
205 return 1;
206 }
207 final int cmpY = Double.compare(p1.getY(), p2.getY());
208 if (cmpY < 0) {
209 return -1;
210 }
211 if (cmpY > 0) {
212 return 1;
213 }
214 final int cmpW = Double.compare(p1.getWeight(), p2.getWeight());
215 if (cmpW < 0) {
216 return -1;
217 }
218 if (cmpW > 0) {
219 return 1;
220 }
221 return 0;
222 }
223 };
224
225 Arrays.sort(observations, cmp);
226 return observations;
227 }
228
229 /**
230 * Guesses the parameters based on the specified observed points.
231 *
232 * @param points Observed points, sorted.
233 * @return the guessed parameters (normalization factor, mean and
234 * sigma).
235 */
236 private double[] basicGuess(WeightedObservedPoint[] points) {
237 final int maxYIdx = findMaxY(points);
238 final double n = points[maxYIdx].getY();
239 final double m = points[maxYIdx].getX();
240
241 double fwhmApprox;
242 try {
243 final double halfY = n + ((m - n) / 2);
244 final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
245 final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
246 fwhmApprox = fwhmX2 - fwhmX1;
247 } catch (OutOfRangeException e) {
248 // TODO: Exceptions should not be used for flow control.
249 fwhmApprox = points[points.length - 1].getX() - points[0].getX();
250 }
251 final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
252
253 return new double[] { n, m, s };
254 }
255
256 /**
257 * Finds index of point in specified points with the largest Y.
258 *
259 * @param points Points to search.
260 * @return the index in specified points array.
261 */
262 private int findMaxY(WeightedObservedPoint[] points) {
263 int maxYIdx = 0;
264 for (int i = 1; i < points.length; i++) {
265 if (points[i].getY() > points[maxYIdx].getY()) {
266 maxYIdx = i;
267 }
268 }
269 return maxYIdx;
270 }
271
272 /**
273 * Interpolates using the specified points to determine X at the
274 * specified Y.
275 *
276 * @param points Points to use for interpolation.
277 * @param startIdx Index within points from which to start the search for
278 * interpolation bounds points.
279 * @param idxStep Index step for searching interpolation bounds points.
280 * @param y Y value for which X should be determined.
281 * @return the value of X for the specified Y.
282 * @throws ZeroException if {@code idxStep} is 0.
283 * @throws OutOfRangeException if specified {@code y} is not within the
284 * range of the specified {@code points}.
285 */
286 private double interpolateXAtY(WeightedObservedPoint[] points,
287 int startIdx,
288 int idxStep,
289 double y)
290 throws OutOfRangeException {
291 if (idxStep == 0) {
292 throw new ZeroException();
293 }
294 final WeightedObservedPoint[] twoPoints
295 = getInterpolationPointsForY(points, startIdx, idxStep, y);
296 final WeightedObservedPoint p1 = twoPoints[0];
297 final WeightedObservedPoint p2 = twoPoints[1];
298 if (p1.getY() == y) {
299 return p1.getX();
300 }
301 if (p2.getY() == y) {
302 return p2.getX();
303 }
304 return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
305 (p2.getY() - p1.getY()));
306 }
307
308 /**
309 * Gets the two bounding interpolation points from the specified points
310 * suitable for determining X at the specified Y.
311 *
312 * @param points Points to use for interpolation.
313 * @param startIdx Index within points from which to start search for
314 * interpolation bounds points.
315 * @param idxStep Index step for search for interpolation bounds points.
316 * @param y Y value for which X should be determined.
317 * @return the array containing two points suitable for determining X at
318 * the specified Y.
319 * @throws ZeroException if {@code idxStep} is 0.
320 * @throws OutOfRangeException if specified {@code y} is not within the
321 * range of the specified {@code points}.
322 */
323 private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
324 int startIdx,
325 int idxStep,
326 double y)
327 throws OutOfRangeException {
328 if (idxStep == 0) {
329 throw new ZeroException();
330 }
331 for (int i = startIdx;
332 idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
333 i += idxStep) {
334 final WeightedObservedPoint p1 = points[i];
335 final WeightedObservedPoint p2 = points[i + idxStep];
336 if (isBetween(y, p1.getY(), p2.getY())) {
337 if (idxStep < 0) {
338 return new WeightedObservedPoint[] { p2, p1 };
339 } else {
340 return new WeightedObservedPoint[] { p1, p2 };
341 }
342 }
343 }
344
345 // Boundaries are replaced by dummy values because the raised
346 // exception is caught and the message never displayed.
347 // TODO: Exceptions should not be used for flow control.
348 throw new OutOfRangeException(y,
349 Double.NEGATIVE_INFINITY,
350 Double.POSITIVE_INFINITY);
351 }
352
353 /**
354 * Determines whether a value is between two other values.
355 *
356 * @param value Value to test whether it is between {@code boundary1}
357 * and {@code boundary2}.
358 * @param boundary1 One end of the range.
359 * @param boundary2 Other end of the range.
360 * @return {@code true} if {@code value} is between {@code boundary1} and
361 * {@code boundary2} (inclusive), {@code false} otherwise.
362 */
363 private boolean isBetween(double value,
364 double boundary1,
365 double boundary2) {
366 return (value >= boundary1 && value <= boundary2) ||
367 (value >= boundary2 && value <= boundary1);
368 }
369 }
370}
Note: See TracBrowser for help on using the repository browser.