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  /*
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.List;
27  
28  import org.hipparchus.analysis.function.HarmonicOscillator;
29  import org.hipparchus.exception.LocalizedCoreFormats;
30  import org.hipparchus.exception.MathIllegalArgumentException;
31  import org.hipparchus.exception.MathIllegalStateException;
32  import org.hipparchus.linear.DiagonalMatrix;
33  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresBuilder;
34  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem;
35  import org.hipparchus.util.FastMath;
36  import org.hipparchus.util.SinCos;
37  
38  /**
39   * Fits points to a {@link
40   * org.hipparchus.analysis.function.HarmonicOscillator.Parametric harmonic oscillator}
41   * function.
42   * <br/>
43   * The {@link #withStartPoint(double[]) initial guess values} must be passed
44   * in the following order:
45   * <ul>
46   *  <li>Amplitude</li>
47   *  <li>Angular frequency</li>
48   *  <li>phase</li>
49   * </ul>
50   * The optimal values will be returned in the same order.
51   *
52   */
53  public class HarmonicCurveFitter extends AbstractCurveFitter {
54      /** Parametric function to be fitted. */
55      private static final HarmonicOscillator.Parametric FUNCTION = new HarmonicOscillator.Parametric();
56      /** Initial guess. */
57      private final double[] initialGuess;
58      /** Maximum number of iterations of the optimization algorithm. */
59      private final int maxIter;
60  
61      /**
62       * Constructor used by the factory methods.
63       *
64       * @param initialGuess Initial guess. If set to {@code null}, the initial guess
65       * will be estimated using the {@link ParameterGuesser}.
66       * @param maxIter Maximum number of iterations of the optimization algorithm.
67       */
68      private HarmonicCurveFitter(double[] initialGuess, int maxIter) {
69          this.initialGuess = initialGuess == null ? null : initialGuess.clone();
70          this.maxIter = maxIter;
71      }
72  
73      /**
74       * Creates a default curve fitter.
75       * The initial guess for the parameters will be {@link ParameterGuesser}
76       * computed automatically, and the maximum number of iterations of the
77       * optimization algorithm is set to {@link Integer#MAX_VALUE}.
78       *
79       * @return a curve fitter.
80       *
81       * @see #withStartPoint(double[])
82       * @see #withMaxIterations(int)
83       */
84      public static HarmonicCurveFitter create() {
85          return new HarmonicCurveFitter(null, Integer.MAX_VALUE);
86      }
87  
88      /**
89       * Configure the start point (initial guess).
90       * @param newStart new start point (initial guess)
91       * @return a new instance.
92       */
93      public HarmonicCurveFitter withStartPoint(double[] newStart) {
94          return new HarmonicCurveFitter(newStart.clone(),
95                                         maxIter);
96      }
97  
98      /**
99       * Configure the maximum number of iterations.
100      * @param newMaxIter maximum number of iterations
101      * @return a new instance.
102      */
103     public HarmonicCurveFitter withMaxIterations(int newMaxIter) {
104         return new HarmonicCurveFitter(initialGuess,
105                                        newMaxIter);
106     }
107 
108     /** {@inheritDoc} */
109     @Override
110     protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> observations) {
111         // Prepare least-squares problem.
112         final int len = observations.size();
113         final double[] target  = new double[len];
114         final double[] weights = new double[len];
115 
116         int i = 0;
117         for (WeightedObservedPoint obs : observations) {
118             target[i]  = obs.getY();
119             weights[i] = obs.getWeight();
120             ++i;
121         }
122 
123         final AbstractCurveFitter.TheoreticalValuesFunction model
124             = new AbstractCurveFitter.TheoreticalValuesFunction(FUNCTION,
125                                                                 observations);
126 
127         final double[] startPoint = initialGuess != null ?
128             initialGuess :
129             // Compute estimation.
130             new ParameterGuesser(observations).guess();
131 
132         // Return a new optimizer set up to fit a Gaussian curve to the
133         // observed points.
134         return new LeastSquaresBuilder().
135                 maxEvaluations(Integer.MAX_VALUE).
136                 maxIterations(maxIter).
137                 start(startPoint).
138                 target(target).
139                 weight(new DiagonalMatrix(weights)).
140                 model(model.getModelFunction(), model.getModelFunctionJacobian()).
141                 build();
142 
143     }
144 
145     /**
146      * This class guesses harmonic coefficients from a sample.
147      * <p>The algorithm used to guess the coefficients is as follows:</p>
148      *
149      * <p>We know \( f(t) \) at some sampling points \( t_i \) and want
150      * to find \( a \), \( \omega \) and \( \phi \) such that
151      * \( f(t) = a \cos (\omega t + \phi) \).
152      * </p>
153      *
154      * <p>From the analytical expression, we can compute two primitives :
155      * \[
156      *     If2(t) = \int f^2 dt  = a^2 (t + S(t)) / 2
157      * \]
158      * \[
159      *     If'2(t) = \int f'^2 dt = a^2 \omega^2 (t - S(t)) / 2
160      * \]
161      * where \(S(t) = \frac{\sin(2 (\omega t + \phi))}{2\omega}\)
162      * </p>
163      *
164      * <p>We can remove \(S\) between these expressions :
165      * \[
166      *     If'2(t) = a^2 \omega^2 t - \omega^2 If2(t)
167      * \]
168      * </p>
169      *
170      * <p>The preceding expression shows that \(If'2 (t)\) is a linear
171      * combination of both \(t\) and \(If2(t)\):
172      * \[
173      *   If'2(t) = A t + B If2(t)
174      * \]
175      * </p>
176      *
177      * <p>From the primitive, we can deduce the same form for definite
178      * integrals between \(t_1\) and \(t_i\) for each \(t_i\) :
179      * \[
180      *   If2(t_i) - If2(t_1) = A (t_i - t_1) + B (If2 (t_i) - If2(t_1))
181      * \]
182      * </p>
183      *
184      * <p>We can find the coefficients \(A\) and \(B\) that best fit the sample
185      * to this linear expression by computing the definite integrals for
186      * each sample points.
187      * </p>
188      *
189      * <p>For a bilinear expression \(z(x_i, y_i) = A x_i + B y_i\), the
190      * coefficients \(A\) and \(B\) that minimize a least-squares criterion
191      * \(\sum (z_i - z(x_i, y_i))^2\) are given by these expressions:</p>
192      * \[
193      *   A = \frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i}
194      *            {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
195      * \]
196      * \[
197      *   B = \frac{\sum x_i x_i \sum y_i z_i - \sum x_i y_i \sum x_i z_i}
198      *            {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
199      *
200      * \]
201      *
202      * <p>In fact, we can assume that both \(a\) and \(\omega\) are positive and
203      * compute them directly, knowing that \(A = a^2 \omega^2\) and that
204      * \(B = -\omega^2\). The complete algorithm is therefore:</p>
205      *
206      * For each \(t_i\) from \(t_1\) to \(t_{n-1}\), compute:
207      * \[ f(t_i) \]
208      * \[ f'(t_i) = \frac{f (t_{i+1}) - f(t_{i-1})}{t_{i+1} - t_{i-1}} \]
209      * \[ x_i = t_i  - t_1 \]
210      * \[ y_i = \int_{t_1}^{t_i} f^2(t) dt \]
211      * \[ z_i = \int_{t_1}^{t_i} f'^2(t) dt \]
212      * and update the sums:
213      * \[ \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 \]
214      *
215      * Then:
216      * \[
217      *  a = \sqrt{\frac{\sum y_i y_i  \sum x_i z_i - \sum x_i y_i \sum y_i z_i }
218      *                 {\sum x_i y_i  \sum x_i z_i - \sum x_i x_i \sum y_i z_i }}
219      * \]
220      * \[
221      *  \omega = \sqrt{\frac{\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i}
222      *                      {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}}
223      * \]
224      *
225      * <p>Once we know \(\omega\) we can compute:
226      * \[
227      *    fc = \omega f(t) \cos(\omega t) - f'(t) \sin(\omega t)
228      * \]
229      * \[
230      *    fs = \omega f(t) \sin(\omega t) + f'(t) \cos(\omega t)
231      * \]
232      * </p>
233      *
234      * <p>It appears that \(fc = a \omega \cos(\phi)\) and
235      * \(fs = -a \omega \sin(\phi)\), so we can use these
236      * expressions to compute \(\phi\). The best estimate over the sample is
237      * given by averaging these expressions.
238      * </p>
239      *
240      * <p>Since integrals and means are involved in the preceding
241      * estimations, these operations run in \(O(n)\) time, where \(n\) is the
242      * number of measurements.</p>
243      */
244     public static class ParameterGuesser {
245         /** Amplitude. */
246         private final double a;
247         /** Angular frequency. */
248         private final double omega;
249         /** Phase. */
250         private final double phi;
251 
252         /**
253          * Simple constructor.
254          *
255          * @param observations Sampled observations.
256          * @throws MathIllegalArgumentException if the sample is too short.
257          * @throws MathIllegalArgumentException if the abscissa range is zero.
258          * @throws MathIllegalStateException when the guessing procedure cannot
259          * produce sensible results.
260          */
261         public ParameterGuesser(Collection<WeightedObservedPoint> observations) {
262             if (observations.size() < 4) {
263                 throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
264                                                     observations.size(), 4, true);
265             }
266 
267             final WeightedObservedPoint[] sorted
268                 = sortObservations(observations).toArray(new WeightedObservedPoint[0]);
269 
270             final double aOmega[] = guessAOmega(sorted);
271             a = aOmega[0];
272             omega = aOmega[1];
273 
274             phi = guessPhi(sorted);
275         }
276 
277         /**
278          * Gets an estimation of the parameters.
279          *
280          * @return the guessed parameters, in the following order:
281          * <ul>
282          *  <li>Amplitude</li>
283          *  <li>Angular frequency</li>
284          *  <li>Phase</li>
285          * </ul>
286          */
287         public double[] guess() {
288             return new double[] { a, omega, phi };
289         }
290 
291         /**
292          * Sort the observations with respect to the abscissa.
293          *
294          * @param unsorted Input observations.
295          * @return the input observations, sorted.
296          */
297         private List<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) {
298             final List<WeightedObservedPoint> observations = new ArrayList<>(unsorted);
299 
300             // Since the samples are almost always already sorted, this
301             // method is implemented as an insertion sort that reorders the
302             // elements in place. Insertion sort is very efficient in this case.
303             WeightedObservedPoint curr = observations.get(0);
304             final int len = observations.size();
305             for (int j = 1; j < len; j++) {
306                 WeightedObservedPoint prec = curr;
307                 curr = observations.get(j);
308                 if (curr.getX() < prec.getX()) {
309                     // the current element should be inserted closer to the beginning
310                     int i = j - 1;
311                     WeightedObservedPoint mI = observations.get(i);
312                     while ((i >= 0) && (curr.getX() < mI.getX())) {
313                         observations.set(i + 1, mI);
314                         if (i != 0) {
315                             mI = observations.get(i - 1);
316                         }
317                         --i;
318                     }
319                     observations.set(i + 1, curr);
320                     curr = observations.get(j);
321                 }
322             }
323 
324             return observations;
325         }
326 
327         /**
328          * Estimate a first guess of the amplitude and angular frequency.
329          *
330          * @param observations Observations, sorted w.r.t. abscissa.
331          * @throws MathIllegalArgumentException if the abscissa range is zero.
332          * @throws MathIllegalStateException when the guessing procedure cannot
333          * produce sensible results.
334          * @return the guessed amplitude (at index 0) and circular frequency
335          * (at index 1).
336          */
337         private double[] guessAOmega(WeightedObservedPoint[] observations) {
338             final double[] aOmega = new double[2];
339 
340             // initialize the sums for the linear model between the two integrals
341             double sx2 = 0;
342             double sy2 = 0;
343             double sxy = 0;
344             double sxz = 0;
345             double syz = 0;
346 
347             double currentX = observations[0].getX();
348             double currentY = observations[0].getY();
349             double f2Integral = 0;
350             double fPrime2Integral = 0;
351             final double startX = currentX;
352             for (int i = 1; i < observations.length; ++i) {
353                 // one step forward
354                 final double previousX = currentX;
355                 final double previousY = currentY;
356                 currentX = observations[i].getX();
357                 currentY = observations[i].getY();
358 
359                 // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
360                 // considering a linear model for f (and therefore constant f')
361                 final double dx = currentX - previousX;
362                 final double dy = currentY - previousY;
363                 final double f2StepIntegral =
364                     dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
365                 final double fPrime2StepIntegral = dy * dy / dx;
366 
367                 final double x = currentX - startX;
368                 f2Integral += f2StepIntegral;
369                 fPrime2Integral += fPrime2StepIntegral;
370 
371                 sx2 += x * x;
372                 sy2 += f2Integral * f2Integral;
373                 sxy += x * f2Integral;
374                 sxz += x * fPrime2Integral;
375                 syz += f2Integral * fPrime2Integral;
376             }
377 
378             // compute the amplitude and pulsation coefficients
379             double c1 = sy2 * sxz - sxy * syz;
380             double c2 = sxy * sxz - sx2 * syz;
381             double c3 = sx2 * sy2 - sxy * sxy;
382             if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
383                 final int last = observations.length - 1;
384                 // Range of the observations, assuming that the
385                 // observations are sorted.
386                 final double xRange = observations[last].getX() - observations[0].getX();
387                 if (xRange == 0) {
388                     throw new MathIllegalArgumentException(LocalizedCoreFormats.ZERO_NOT_ALLOWED);
389                 }
390                 aOmega[1] = 2 * Math.PI / xRange;
391 
392                 double yMin = Double.POSITIVE_INFINITY;
393                 double yMax = Double.NEGATIVE_INFINITY;
394                 for (int i = 1; i < observations.length; ++i) {
395                     final double y = observations[i].getY();
396                     if (y < yMin) {
397                         yMin = y;
398                     }
399                     if (y > yMax) {
400                         yMax = y;
401                     }
402                 }
403                 aOmega[0] = 0.5 * (yMax - yMin);
404             } else {
405                 if (c2 == 0) {
406                     // In some ill-conditioned cases (cf. MATH-844), the guesser
407                     // procedure cannot produce sensible results.
408                     throw new MathIllegalStateException(LocalizedCoreFormats.ZERO_DENOMINATOR);
409                 }
410 
411                 aOmega[0] = FastMath.sqrt(c1 / c2);
412                 aOmega[1] = FastMath.sqrt(c2 / c3);
413             }
414 
415             return aOmega;
416         }
417 
418         /**
419          * Estimate a first guess of the phase.
420          *
421          * @param observations Observations, sorted w.r.t. abscissa.
422          * @return the guessed phase.
423          */
424         private double guessPhi(WeightedObservedPoint[] observations) {
425             // initialize the means
426             double fcMean = 0;
427             double fsMean = 0;
428 
429             double currentX = observations[0].getX();
430             double currentY = observations[0].getY();
431             for (int i = 1; i < observations.length; ++i) {
432                 // one step forward
433                 final double previousX = currentX;
434                 final double previousY = currentY;
435                 currentX = observations[i].getX();
436                 currentY = observations[i].getY();
437                 final double currentYPrime = (currentY - previousY) / (currentX - previousX);
438 
439                 double omegaX = omega * currentX;
440                 SinCos sc     = FastMath.sinCos(omegaX);
441                 fcMean += omega * currentY * sc.cos() - currentYPrime * sc.sin();
442                 fsMean += omega * currentY * sc.sin() + currentYPrime * sc.cos();
443             }
444 
445             return FastMath.atan2(-fsMean, fcMean);
446         }
447     }
448 }