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.math4.legacy.fitting;
18
19 import java.util.Collection;
20 import org.apache.commons.math4.legacy.analysis.function.HarmonicOscillator;
21 import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
22 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
23 import org.apache.commons.math4.legacy.exception.ZeroException;
24 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
25 import org.apache.commons.math4.core.jdkmath.JdkMath;
26
27 /**
28 * Fits points to a {@link
29 * org.apache.commons.math4.legacy.analysis.function.HarmonicOscillator.Parametric harmonic oscillator}
30 * function.
31 * <br>
32 * The {@link #withStartPoint(double[]) initial guess values} must be passed
33 * in the following order:
34 * <ul>
35 * <li>Amplitude</li>
36 * <li>Angular frequency</li>
37 * <li>phase</li>
38 * </ul>
39 * The optimal values will be returned in the same order.
40 *
41 * @since 3.3
42 */
43 public final class HarmonicCurveFitter extends SimpleCurveFitter {
44 /** Parametric function to be fitted. */
45 private static final HarmonicOscillator.Parametric FUNCTION = new HarmonicOscillator.Parametric();
46
47 /**
48 * Constructor used by the factory methods.
49 *
50 * @param initialGuess Initial guess. If set to {@code null}, the initial guess
51 * will be estimated using the {@link ParameterGuesser}.
52 * @param maxIter Maximum number of iterations of the optimization algorithm.
53 */
54 private HarmonicCurveFitter(double[] initialGuess,
55 int maxIter) {
56 super(FUNCTION, initialGuess, new ParameterGuesser(), maxIter);
57 }
58
59 /**
60 * Creates a default curve fitter.
61 * The initial guess for the parameters will be {@link ParameterGuesser}
62 * computed automatically, and the maximum number of iterations of the
63 * optimization algorithm is set to {@link Integer#MAX_VALUE}.
64 *
65 * @return a curve fitter.
66 *
67 * @see #withStartPoint(double[])
68 * @see #withMaxIterations(int)
69 */
70 public static HarmonicCurveFitter create() {
71 return new HarmonicCurveFitter(null, Integer.MAX_VALUE);
72 }
73
74 /**
75 * This class guesses harmonic coefficients from a sample.
76 * <p>The algorithm used to guess the coefficients is as follows:</p>
77 *
78 * <p>We know \( f(t) \) at some sampling points \( t_i \) and want
79 * to find \( a \), \( \omega \) and \( \phi \) such that
80 * \( f(t) = a \cos (\omega t + \phi) \).
81 * </p>
82 *
83 * <p>From the analytical expression, we can compute two primitives :
84 * \[
85 * If2(t) = \int f^2 dt = a^2 (t + S(t)) / 2
86 * \]
87 * \[
88 * If'2(t) = \int f'^2 dt = a^2 \omega^2 (t - S(t)) / 2
89 * \]
90 * where \(S(t) = \frac{\sin(2 (\omega t + \phi))}{2\omega}\)
91 * </p>
92 *
93 * <p>We can remove \(S\) between these expressions :
94 * \[
95 * If'2(t) = a^2 \omega^2 t - \omega^2 If2(t)
96 * \]
97 * </p>
98 *
99 * <p>The preceding expression shows that \(If'2 (t)\) is a linear
100 * combination of both \(t\) and \(If2(t)\):
101 * \[
102 * If'2(t) = A t + B If2(t)
103 * \]
104 * </p>
105 *
106 * <p>From the primitive, we can deduce the same form for definite
107 * integrals between \(t_1\) and \(t_i\) for each \(t_i\) :
108 * \[
109 * If2(t_i) - If2(t_1) = A (t_i - t_1) + B (If2 (t_i) - If2(t_1))
110 * \]
111 * </p>
112 *
113 * <p>We can find the coefficients \(A\) and \(B\) that best fit the sample
114 * to this linear expression by computing the definite integrals for
115 * each sample points.
116 * </p>
117 *
118 * <p>For a bilinear expression \(z(x_i, y_i) = A x_i + B y_i\), the
119 * coefficients \(A\) and \(B\) that minimize a least-squares criterion
120 * \(\sum (z_i - z(x_i, y_i))^2\) are given by these expressions:</p>
121 * \[
122 * A = \frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i}
123 * {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
124 * \]
125 * \[
126 * B = \frac{\sum x_i x_i \sum y_i z_i - \sum x_i y_i \sum x_i z_i}
127 * {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
128 *
129 * \]
130 *
131 * <p>In fact, we can assume that both \(a\) and \(\omega\) are positive and
132 * compute them directly, knowing that \(A = a^2 \omega^2\) and that
133 * \(B = -\omega^2\). The complete algorithm is therefore:</p>
134 *
135 * For each \(t_i\) from \(t_1\) to \(t_{n-1}\), compute:
136 * \[ f(t_i) \]
137 * \[ f'(t_i) = \frac{f (t_{i+1}) - f(t_{i-1})}{t_{i+1} - t_{i-1}} \]
138 * \[ x_i = t_i - t_1 \]
139 * \[ y_i = \int_{t_1}^{t_i} f^2(t) dt \]
140 * \[ z_i = \int_{t_1}^{t_i} f'^2(t) dt \]
141 * and update the sums:
142 * \[ \sum x_i x_i, \sum y_i y_i, \sum x_i y_i, \sum x_i z_i, \sum y_i z_i \]
143 *
144 * Then:
145 * \[
146 * a = \sqrt{\frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i }
147 * {\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i }}
148 * \]
149 * \[
150 * \omega = \sqrt{\frac{\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i}
151 * {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}}
152 * \]
153 *
154 * <p>Once we know \(\omega\) we can compute:
155 * \[
156 * fc = \omega f(t) \cos(\omega t) - f'(t) \sin(\omega t)
157 * \]
158 * \[
159 * fs = \omega f(t) \sin(\omega t) + f'(t) \cos(\omega t)
160 * \]
161 * </p>
162 *
163 * <p>It appears that \(fc = a \omega \cos(\phi)\) and
164 * \(fs = -a \omega \sin(\phi)\), so we can use these
165 * expressions to compute \(\phi\). The best estimate over the sample is
166 * given by averaging these expressions.
167 * </p>
168 *
169 * <p>Since integrals and means are involved in the preceding
170 * estimations, these operations run in \(O(n)\) time, where \(n\) is the
171 * number of measurements.</p>
172 */
173 public static class ParameterGuesser extends SimpleCurveFitter.ParameterGuesser {
174 /**
175 * {@inheritDoc}
176 *
177 * @return the guessed parameters, in the following order:
178 * <ul>
179 * <li>Amplitude</li>
180 * <li>Angular frequency</li>
181 * <li>Phase</li>
182 * </ul>
183 * @throws NumberIsTooSmallException if the sample is too short.
184 * @throws ZeroException if the abscissa range is zero.
185 * @throws MathIllegalStateException when the guessing procedure cannot
186 * produce sensible results.
187 */
188 @Override
189 public double[] guess(Collection<WeightedObservedPoint> observations) {
190 if (observations.size() < 4) {
191 throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
192 observations.size(), 4, true);
193 }
194
195 final WeightedObservedPoint[] sorted
196 = sortObservations(observations).toArray(new WeightedObservedPoint[0]);
197
198 final double[] aOmega = guessAOmega(sorted);
199 final double a = aOmega[0];
200 final double omega = aOmega[1];
201
202 final double phi = guessPhi(sorted, omega);
203
204 return new double[] { a, omega, phi };
205 }
206
207 /**
208 * Estimate a first guess of the amplitude and angular frequency.
209 *
210 * @param observations Observations, sorted w.r.t. abscissa.
211 * @throws ZeroException if the abscissa range is zero.
212 * @throws MathIllegalStateException when the guessing procedure cannot
213 * produce sensible results.
214 * @return the guessed amplitude (at index 0) and circular frequency
215 * (at index 1).
216 */
217 private double[] guessAOmega(WeightedObservedPoint[] observations) {
218 final double[] aOmega = new double[2];
219
220 // initialize the sums for the linear model between the two integrals
221 double sx2 = 0;
222 double sy2 = 0;
223 double sxy = 0;
224 double sxz = 0;
225 double syz = 0;
226
227 double currentX = observations[0].getX();
228 double currentY = observations[0].getY();
229 double f2Integral = 0;
230 double fPrime2Integral = 0;
231 final double startX = currentX;
232 for (int i = 1; i < observations.length; ++i) {
233 // one step forward
234 final double previousX = currentX;
235 final double previousY = currentY;
236 currentX = observations[i].getX();
237 currentY = observations[i].getY();
238
239 // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
240 // considering a linear model for f (and therefore constant f')
241 final double dx = currentX - previousX;
242 final double dy = currentY - previousY;
243 final double f2StepIntegral =
244 dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
245 final double fPrime2StepIntegral = dy * dy / dx;
246
247 final double x = currentX - startX;
248 f2Integral += f2StepIntegral;
249 fPrime2Integral += fPrime2StepIntegral;
250
251 sx2 += x * x;
252 sy2 += f2Integral * f2Integral;
253 sxy += x * f2Integral;
254 sxz += x * fPrime2Integral;
255 syz += f2Integral * fPrime2Integral;
256 }
257
258 // compute the amplitude and pulsation coefficients
259 double c1 = sy2 * sxz - sxy * syz;
260 double c2 = sxy * sxz - sx2 * syz;
261 double c3 = sx2 * sy2 - sxy * sxy;
262 if (c1 / c2 < 0 || c2 / c3 < 0) {
263 final int last = observations.length - 1;
264 // Range of the observations, assuming that the
265 // observations are sorted.
266 final double xRange = observations[last].getX() - observations[0].getX();
267 if (xRange == 0) {
268 throw new ZeroException();
269 }
270 aOmega[1] = 2 * Math.PI / xRange;
271
272 double yMin = Double.POSITIVE_INFINITY;
273 double yMax = Double.NEGATIVE_INFINITY;
274 for (int i = 1; i < observations.length; ++i) {
275 final double y = observations[i].getY();
276 if (y < yMin) {
277 yMin = y;
278 }
279 if (y > yMax) {
280 yMax = y;
281 }
282 }
283 aOmega[0] = 0.5 * (yMax - yMin);
284 } else {
285 if (c2 == 0) {
286 // In some ill-conditioned cases (cf. MATH-844), the guesser
287 // procedure cannot produce sensible results.
288 throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
289 }
290
291 aOmega[0] = JdkMath.sqrt(c1 / c2);
292 aOmega[1] = JdkMath.sqrt(c2 / c3);
293 }
294
295 return aOmega;
296 }
297
298 /**
299 * Estimate a first guess of the phase.
300 *
301 * @param observations Observations, sorted w.r.t. abscissa.
302 * @param omega Angular frequency.
303 * @return the guessed phase.
304 */
305 private double guessPhi(WeightedObservedPoint[] observations,
306 double omega) {
307 // initialize the means
308 double fcMean = 0;
309 double fsMean = 0;
310
311 double currentX = observations[0].getX();
312 double currentY = observations[0].getY();
313 for (int i = 1; i < observations.length; ++i) {
314 // one step forward
315 final double previousX = currentX;
316 final double previousY = currentY;
317 currentX = observations[i].getX();
318 currentY = observations[i].getY();
319 final double currentYPrime = (currentY - previousY) / (currentX - previousX);
320
321 double omegaX = omega * currentX;
322 double cosine = JdkMath.cos(omegaX);
323 double sine = JdkMath.sin(omegaX);
324 fcMean += omega * currentY * cosine - currentYPrime * sine;
325 fsMean += omega * currentY * sine + currentYPrime * cosine;
326 }
327
328 return JdkMath.atan2(-fsMean, fcMean);
329 }
330 }
331 }