HarmonicFitter.java

  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.  */

  18. package org.apache.commons.math3.optimization.fitting;

  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;

  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.     }

  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.     }

  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.     }

  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;

  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.             }

  203.             final WeightedObservedPoint[] sorted = sortObservations(observations);

  205.             final double aOmega[] = guessAOmega(sorted);
  206.             a = aOmega[0];
  207.             omega = aOmega[1];

  209.             phi = guessPhi(sorted);
  210.         }

  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.         }

  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();

  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.             }

  257.             return observations;
  258.         }

  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];

  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;

  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();

  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;

  302.                 final double x = currentX - startX;
  303.                 f2Integral += f2StepIntegral;
  304.                 fPrime2Integral += fPrime2StepIntegral;

  306.                 sx2 += x * x;
  307.                 sy2 += f2Integral * f2Integral;
  308.                 sxy += x * f2Integral;
  309.                 sxz += x * fPrime2Integral;
  310.                 syz += f2Integral * fPrime2Integral;
  311.             }

  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;

  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.                 }

  346.                 aOmega[0] = FastMath.sqrt(c1 / c2);
  347.                 aOmega[1] = FastMath.sqrt(c2 / c3);
  348.             }

  350.             return aOmega;
  351.         }

  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;

  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);

  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.             }

  381.             return FastMath.atan2(-fsMean, fcMean);
  382.         }
  383.     }
  384. }
Created with JaCoCo 0.7.2.201409121644