source: src/main/java/agents/anac/y2019/harddealer/math3/fitting/HarmonicFitter.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: 16.2 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 agents.anac.y2019.harddealer.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
20import agents.anac.y2019.harddealer.math3.analysis.function.HarmonicOscillator;
21import agents.anac.y2019.harddealer.math3.exception.ZeroException;
22import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
23import agents.anac.y2019.harddealer.math3.exception.MathIllegalStateException;
24import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
25import agents.anac.y2019.harddealer.math3.util.FastMath;
26
27/**
28 * Class that implements a curve fitting specialized for sinusoids.
29 *
30 * Harmonic fitting is a very simple case of curve fitting. The
31 * estimated coefficients are the amplitude a, the pulsation ω and
32 * the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
33 * searched by a least square estimator initialized with a rough guess
34 * based on integrals.
35 *
36 * @since 2.0
37 * @deprecated As of 3.3. Please use {@link HarmonicCurveFitter} and
38 * {@link WeightedObservedPoints} instead.
39 */
40@Deprecated
41public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> {
42 /**
43 * Simple constructor.
44 * @param optimizer Optimizer to use for the fitting.
45 */
46 public HarmonicFitter(final MultivariateVectorOptimizer optimizer) {
47 super(optimizer);
48 }
49
50 /**
51 * Fit an harmonic function to the observed points.
52 *
53 * @param initialGuess First guess values in the following order:
54 * <ul>
55 * <li>Amplitude</li>
56 * <li>Angular frequency</li>
57 * <li>Phase</li>
58 * </ul>
59 * @return the parameters of the harmonic function that best fits the
60 * observed points (in the same order as above).
61 */
62 public double[] fit(double[] initialGuess) {
63 return fit(new HarmonicOscillator.Parametric(), initialGuess);
64 }
65
66 /**
67 * Fit an harmonic function to the observed points.
68 * An initial guess will be automatically computed.
69 *
70 * @return the parameters of the harmonic function that best fits the
71 * observed points (see the other {@link #fit(double[]) fit} method.
72 * @throws NumberIsTooSmallException if the sample is too short for the
73 * the first guess to be computed.
74 * @throws ZeroException if the first guess cannot be computed because
75 * the abscissa range is zero.
76 */
77 public double[] fit() {
78 return fit((new ParameterGuesser(getObservations())).guess());
79 }
80
81 /**
82 * This class guesses harmonic coefficients from a sample.
83 * <p>The algorithm used to guess the coefficients is as follows:</p>
84 *
85 * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
86 * &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
87 * </p>
88 *
89 * <p>From the analytical expression, we can compute two primitives :
90 * <pre>
91 * If2 (t) = &int; f<sup>2</sup> = a<sup>2</sup> &times; [t + S (t)] / 2
92 * If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
93 * where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
94 * </pre>
95 * </p>
96 *
97 * <p>We can remove S between these expressions :
98 * <pre>
99 * If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
100 * </pre>
101 * </p>
102 *
103 * <p>The preceding expression shows that If'2 (t) is a linear
104 * combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
105 * </p>
106 *
107 * <p>From the primitive, we can deduce the same form for definite
108 * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
109 * <pre>
110 * If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A &times; (t<sub>i</sub> - t<sub>1</sub>) + B &times; (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
111 * </pre>
112 * </p>
113 *
114 * <p>We can find the coefficients A and B that best fit the sample
115 * to this linear expression by computing the definite integrals for
116 * each sample points.
117 * </p>
118 *
119 * <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A &times; x<sub>i</sub> + B &times; y<sub>i</sub>, the
120 * coefficients A and B that minimize a least square criterion
121 * &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
122 * <pre>
123 *
124 * &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
125 * A = ------------------------
126 * &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
127 *
128 * &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub>
129 * B = ------------------------
130 * &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
131 * </pre>
132 * </p>
133 *
134 *
135 * <p>In fact, we can assume both a and &omega; are positive and
136 * compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
137 * B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
138 * <pre>
139 *
140 * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
141 * f (t<sub>i</sub>)
142 * f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>)
143 * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
144 * y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
145 * z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
146 * update the sums &sum;x<sub>i</sub>x<sub>i</sub>, &sum;y<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>z<sub>i</sub> and &sum;y<sub>i</sub>z<sub>i</sub>
147 * end for
148 *
149 * |--------------------------
150 * \ | &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
151 * a = \ | ------------------------
152 * \| &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
153 *
154 *
155 * |--------------------------
156 * \ | &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
157 * &omega; = \ | ------------------------
158 * \| &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
159 *
160 * </pre>
161 * </p>
162 *
163 * <p>Once we know &omega;, we can compute:
164 * <pre>
165 * fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
166 * fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
167 * </pre>
168 * </p>
169 *
170 * <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
171 * <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
172 * expressions to compute &phi;. The best estimate over the sample is
173 * given by averaging these expressions.
174 * </p>
175 *
176 * <p>Since integrals and means are involved in the preceding
177 * estimations, these operations run in O(n) time, where n is the
178 * number of measurements.</p>
179 */
180 public static class ParameterGuesser {
181 /** Amplitude. */
182 private final double a;
183 /** Angular frequency. */
184 private final double omega;
185 /** Phase. */
186 private final double phi;
187
188 /**
189 * Simple constructor.
190 *
191 * @param observations Sampled observations.
192 * @throws NumberIsTooSmallException if the sample is too short.
193 * @throws ZeroException if the abscissa range is zero.
194 * @throws MathIllegalStateException when the guessing procedure cannot
195 * produce sensible results.
196 */
197 public ParameterGuesser(WeightedObservedPoint[] observations) {
198 if (observations.length < 4) {
199 throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
200 observations.length, 4, true);
201 }
202
203 final WeightedObservedPoint[] sorted = sortObservations(observations);
204
205 final double aOmega[] = guessAOmega(sorted);
206 a = aOmega[0];
207 omega = aOmega[1];
208
209 phi = guessPhi(sorted);
210 }
211
212 /**
213 * Gets an estimation of the parameters.
214 *
215 * @return the guessed parameters, in the following order:
216 * <ul>
217 * <li>Amplitude</li>
218 * <li>Angular frequency</li>
219 * <li>Phase</li>
220 * </ul>
221 */
222 public double[] guess() {
223 return new double[] { a, omega, phi };
224 }
225
226 /**
227 * Sort the observations with respect to the abscissa.
228 *
229 * @param unsorted Input observations.
230 * @return the input observations, sorted.
231 */
232 private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
233 final WeightedObservedPoint[] observations = unsorted.clone();
234
235 // Since the samples are almost always already sorted, this
236 // method is implemented as an insertion sort that reorders the
237 // elements in place. Insertion sort is very efficient in this case.
238 WeightedObservedPoint curr = observations[0];
239 for (int j = 1; j < observations.length; ++j) {
240 WeightedObservedPoint prec = curr;
241 curr = observations[j];
242 if (curr.getX() < prec.getX()) {
243 // the current element should be inserted closer to the beginning
244 int i = j - 1;
245 WeightedObservedPoint mI = observations[i];
246 while ((i >= 0) && (curr.getX() < mI.getX())) {
247 observations[i + 1] = mI;
248 if (i-- != 0) {
249 mI = observations[i];
250 }
251 }
252 observations[i + 1] = curr;
253 curr = observations[j];
254 }
255 }
256
257 return observations;
258 }
259
260 /**
261 * Estimate a first guess of the amplitude and angular frequency.
262 * This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method
263 * has been called previously.
264 *
265 * @param observations Observations, sorted w.r.t. abscissa.
266 * @throws ZeroException if the abscissa range is zero.
267 * @throws MathIllegalStateException when the guessing procedure cannot
268 * produce sensible results.
269 * @return the guessed amplitude (at index 0) and circular frequency
270 * (at index 1).
271 */
272 private double[] guessAOmega(WeightedObservedPoint[] observations) {
273 final double[] aOmega = new double[2];
274
275 // initialize the sums for the linear model between the two integrals
276 double sx2 = 0;
277 double sy2 = 0;
278 double sxy = 0;
279 double sxz = 0;
280 double syz = 0;
281
282 double currentX = observations[0].getX();
283 double currentY = observations[0].getY();
284 double f2Integral = 0;
285 double fPrime2Integral = 0;
286 final double startX = currentX;
287 for (int i = 1; i < observations.length; ++i) {
288 // one step forward
289 final double previousX = currentX;
290 final double previousY = currentY;
291 currentX = observations[i].getX();
292 currentY = observations[i].getY();
293
294 // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
295 // considering a linear model for f (and therefore constant f')
296 final double dx = currentX - previousX;
297 final double dy = currentY - previousY;
298 final double f2StepIntegral =
299 dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
300 final double fPrime2StepIntegral = dy * dy / dx;
301
302 final double x = currentX - startX;
303 f2Integral += f2StepIntegral;
304 fPrime2Integral += fPrime2StepIntegral;
305
306 sx2 += x * x;
307 sy2 += f2Integral * f2Integral;
308 sxy += x * f2Integral;
309 sxz += x * fPrime2Integral;
310 syz += f2Integral * fPrime2Integral;
311 }
312
313 // compute the amplitude and pulsation coefficients
314 double c1 = sy2 * sxz - sxy * syz;
315 double c2 = sxy * sxz - sx2 * syz;
316 double c3 = sx2 * sy2 - sxy * sxy;
317 if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
318 final int last = observations.length - 1;
319 // Range of the observations, assuming that the
320 // observations are sorted.
321 final double xRange = observations[last].getX() - observations[0].getX();
322 if (xRange == 0) {
323 throw new ZeroException();
324 }
325 aOmega[1] = 2 * Math.PI / xRange;
326
327 double yMin = Double.POSITIVE_INFINITY;
328 double yMax = Double.NEGATIVE_INFINITY;
329 for (int i = 1; i < observations.length; ++i) {
330 final double y = observations[i].getY();
331 if (y < yMin) {
332 yMin = y;
333 }
334 if (y > yMax) {
335 yMax = y;
336 }
337 }
338 aOmega[0] = 0.5 * (yMax - yMin);
339 } else {
340 if (c2 == 0) {
341 // In some ill-conditioned cases (cf. MATH-844), the guesser
342 // procedure cannot produce sensible results.
343 throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
344 }
345
346 aOmega[0] = FastMath.sqrt(c1 / c2);
347 aOmega[1] = FastMath.sqrt(c2 / c3);
348 }
349
350 return aOmega;
351 }
352
353 /**
354 * Estimate a first guess of the phase.
355 *
356 * @param observations Observations, sorted w.r.t. abscissa.
357 * @return the guessed phase.
358 */
359 private double guessPhi(WeightedObservedPoint[] observations) {
360 // initialize the means
361 double fcMean = 0;
362 double fsMean = 0;
363
364 double currentX = observations[0].getX();
365 double currentY = observations[0].getY();
366 for (int i = 1; i < observations.length; ++i) {
367 // one step forward
368 final double previousX = currentX;
369 final double previousY = currentY;
370 currentX = observations[i].getX();
371 currentY = observations[i].getY();
372 final double currentYPrime = (currentY - previousY) / (currentX - previousX);
373
374 double omegaX = omega * currentX;
375 double cosine = FastMath.cos(omegaX);
376 double sine = FastMath.sin(omegaX);
377 fcMean += omega * currentY * cosine - currentYPrime * sine;
378 fsMean += omega * currentY * sine + currentYPrime * cosine;
379 }
380
381 return FastMath.atan2(-fsMean, fcMean);
382 }
383 }
384}
Note: See TracBrowser for help on using the repository browser.