View Javadoc
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 org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
21  import org.apache.commons.math3.analysis.function.HarmonicOscillator;
22  import org.apache.commons.math3.exception.ZeroException;
23  import org.apache.commons.math3.exception.NumberIsTooSmallException;
24  import org.apache.commons.math3.exception.MathIllegalStateException;
25  import org.apache.commons.math3.exception.util.LocalizedFormats;
26  import org.apache.commons.math3.util.FastMath;
27  
28  /**
29   * Class that implements a curve fitting specialized for sinusoids.
30   *
31   * Harmonic fitting is a very simple case of curve fitting. The
32   * estimated coefficients are the amplitude a, the pulsation ω and
33   * the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
34   * searched by a least square estimator initialized with a rough guess
35   * based on integrals.
36   *
37   * @deprecated As of 3.1 (to be removed in 4.0).
38   * @since 2.0
39   */
40  @Deprecated
41  public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> {
42      /**
43       * Simple constructor.
44       * @param optimizer Optimizer to use for the fitting.
45       */
46      public HarmonicFitter(final DifferentiableMultivariateVectorOptimizer 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 }