HarmonicCurveFitter.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.  */
  17. package org.apache.commons.math4.legacy.fitting;

  18. import java.util.Collection;
  19. import org.apache.commons.math4.legacy.analysis.function.HarmonicOscillator;
  20. import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
  21. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math4.legacy.exception.ZeroException;
  23. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  24. import org.apache.commons.math4.core.jdkmath.JdkMath;

  25. /**
  26.  * Fits points to a {@link
  27.  * org.apache.commons.math4.legacy.analysis.function.HarmonicOscillator.Parametric harmonic oscillator}
  28.  * function.
  29.  * <br>
  30.  * The {@link #withStartPoint(double[]) initial guess values} must be passed
  31.  * in the following order:
  32.  * <ul>
  33.  *  <li>Amplitude</li>
  34.  *  <li>Angular frequency</li>
  35.  *  <li>phase</li>
  36.  * </ul>
  37.  * The optimal values will be returned in the same order.
  38.  *
  39.  * @since 3.3
  40.  */
  41. public final class HarmonicCurveFitter extends SimpleCurveFitter {
  42.     /** Parametric function to be fitted. */
  43.     private static final HarmonicOscillator.Parametric FUNCTION = new HarmonicOscillator.Parametric();

  44.     /**
  45.      * Constructor used by the factory methods.
  46.      *
  47.      * @param initialGuess Initial guess. If set to {@code null}, the initial guess
  48.      * will be estimated using the {@link ParameterGuesser}.
  49.      * @param maxIter Maximum number of iterations of the optimization algorithm.
  50.      */
  51.     private HarmonicCurveFitter(double[] initialGuess,
  52.                                 int maxIter) {
  53.         super(FUNCTION, initialGuess, new ParameterGuesser(), maxIter);
  54.     }

  55.     /**
  56.      * Creates a default curve fitter.
  57.      * The initial guess for the parameters will be {@link ParameterGuesser}
  58.      * computed automatically, and the maximum number of iterations of the
  59.      * optimization algorithm is set to {@link Integer#MAX_VALUE}.
  60.      *
  61.      * @return a curve fitter.
  62.      *
  63.      * @see #withStartPoint(double[])
  64.      * @see #withMaxIterations(int)
  65.      */
  66.     public static HarmonicCurveFitter create() {
  67.         return new HarmonicCurveFitter(null, Integer.MAX_VALUE);
  68.     }

  69.     /**
  70.      * This class guesses harmonic coefficients from a sample.
  71.      * <p>The algorithm used to guess the coefficients is as follows:</p>
  72.      *
  73.      * <p>We know \( f(t) \) at some sampling points \( t_i \) and want
  74.      * to find \( a \), \( \omega \) and \( \phi \) such that
  75.      * \( f(t) = a \cos (\omega t + \phi) \).
  76.      * </p>
  77.      *
  78.      * <p>From the analytical expression, we can compute two primitives :
  79.      * \[
  80.      *     If2(t) = \int f^2 dt  = a^2 (t + S(t)) / 2
  81.      * \]
  82.      * \[
  83.      *     If'2(t) = \int f'^2 dt = a^2 \omega^2 (t - S(t)) / 2
  84.      * \]
  85.      * where \(S(t) = \frac{\sin(2 (\omega t + \phi))}{2\omega}\)
  86.      * </p>
  87.      *
  88.      * <p>We can remove \(S\) between these expressions :
  89.      * \[
  90.      *     If'2(t) = a^2 \omega^2 t - \omega^2 If2(t)
  91.      * \]
  92.      * </p>
  93.      *
  94.      * <p>The preceding expression shows that \(If'2 (t)\) is a linear
  95.      * combination of both \(t\) and \(If2(t)\):
  96.      * \[
  97.      *   If'2(t) = A t + B If2(t)
  98.      * \]
  99.      * </p>
  100.      *
  101.      * <p>From the primitive, we can deduce the same form for definite
  102.      * integrals between \(t_1\) and \(t_i\) for each \(t_i\) :
  103.      * \[
  104.      *   If2(t_i) - If2(t_1) = A (t_i - t_1) + B (If2 (t_i) - If2(t_1))
  105.      * \]
  106.      * </p>
  107.      *
  108.      * <p>We can find the coefficients \(A\) and \(B\) that best fit the sample
  109.      * to this linear expression by computing the definite integrals for
  110.      * each sample points.
  111.      * </p>
  112.      *
  113.      * <p>For a bilinear expression \(z(x_i, y_i) = A x_i + B y_i\), the
  114.      * coefficients \(A\) and \(B\) that minimize a least-squares criterion
  115.      * \(\sum (z_i - z(x_i, y_i))^2\) are given by these expressions:</p>
  116.      * \[
  117.      *   A = \frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i}
  118.      *            {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
  119.      * \]
  120.      * \[
  121.      *   B = \frac{\sum x_i x_i \sum y_i z_i - \sum x_i y_i \sum x_i z_i}
  122.      *            {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
  123.      *
  124.      * \]
  125.      *
  126.      * <p>In fact, we can assume that both \(a\) and \(\omega\) are positive and
  127.      * compute them directly, knowing that \(A = a^2 \omega^2\) and that
  128.      * \(B = -\omega^2\). The complete algorithm is therefore:</p>
  129.      *
  130.      * For each \(t_i\) from \(t_1\) to \(t_{n-1}\), compute:
  131.      * \[ f(t_i) \]
  132.      * \[ f'(t_i) = \frac{f (t_{i+1}) - f(t_{i-1})}{t_{i+1} - t_{i-1}} \]
  133.      * \[ x_i = t_i  - t_1 \]
  134.      * \[ y_i = \int_{t_1}^{t_i} f^2(t) dt \]
  135.      * \[ z_i = \int_{t_1}^{t_i} f'^2(t) dt \]
  136.      * and update the sums:
  137.      * \[ \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 \]
  138.      *
  139.      * Then:
  140.      * \[
  141.      *  a = \sqrt{\frac{\sum y_i y_i  \sum x_i z_i - \sum x_i y_i \sum y_i z_i }
  142.      *                 {\sum x_i y_i  \sum x_i z_i - \sum x_i x_i \sum y_i z_i }}
  143.      * \]
  144.      * \[
  145.      *  \omega = \sqrt{\frac{\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i}
  146.      *                      {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}}
  147.      * \]
  148.      *
  149.      * <p>Once we know \(\omega\) we can compute:
  150.      * \[
  151.      *    fc = \omega f(t) \cos(\omega t) - f'(t) \sin(\omega t)
  152.      * \]
  153.      * \[
  154.      *    fs = \omega f(t) \sin(\omega t) + f'(t) \cos(\omega t)
  155.      * \]
  156.      * </p>
  157.      *
  158.      * <p>It appears that \(fc = a \omega \cos(\phi)\) and
  159.      * \(fs = -a \omega \sin(\phi)\), so we can use these
  160.      * expressions to compute \(\phi\). The best estimate over the sample is
  161.      * given by averaging these expressions.
  162.      * </p>
  163.      *
  164.      * <p>Since integrals and means are involved in the preceding
  165.      * estimations, these operations run in \(O(n)\) time, where \(n\) is the
  166.      * number of measurements.</p>
  167.      */
  168.     public static class ParameterGuesser extends SimpleCurveFitter.ParameterGuesser {
  169.         /**
  170.          * {@inheritDoc}
  171.          *
  172.          * @return the guessed parameters, in the following order:
  173.          * <ul>
  174.          *  <li>Amplitude</li>
  175.          *  <li>Angular frequency</li>
  176.          *  <li>Phase</li>
  177.          * </ul>
  178.          * @throws NumberIsTooSmallException if the sample is too short.
  179.          * @throws ZeroException if the abscissa range is zero.
  180.          * @throws MathIllegalStateException when the guessing procedure cannot
  181.          * produce sensible results.
  182.          */
  183.         @Override
  184.         public double[] guess(Collection<WeightedObservedPoint> observations) {
  185.             if (observations.size() < 4) {
  186.                 throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
  187.                                                     observations.size(), 4, true);
  188.             }

  189.             final WeightedObservedPoint[] sorted
  190.                 = sortObservations(observations).toArray(new WeightedObservedPoint[0]);

  191.             final double[] aOmega = guessAOmega(sorted);
  192.             final double a = aOmega[0];
  193.             final double omega = aOmega[1];

  194.             final double phi = guessPhi(sorted, omega);

  195.             return new double[] { a, omega, phi };
  196.         }

  197.         /**
  198.          * Estimate a first guess of the amplitude and angular frequency.
  199.          *
  200.          * @param observations Observations, sorted w.r.t. abscissa.
  201.          * @throws ZeroException if the abscissa range is zero.
  202.          * @throws MathIllegalStateException when the guessing procedure cannot
  203.          * produce sensible results.
  204.          * @return the guessed amplitude (at index 0) and circular frequency
  205.          * (at index 1).
  206.          */
  207.         private double[] guessAOmega(WeightedObservedPoint[] observations) {
  208.             final double[] aOmega = new double[2];

  209.             // initialize the sums for the linear model between the two integrals
  210.             double sx2 = 0;
  211.             double sy2 = 0;
  212.             double sxy = 0;
  213.             double sxz = 0;
  214.             double syz = 0;

  215.             double currentX = observations[0].getX();
  216.             double currentY = observations[0].getY();
  217.             double f2Integral = 0;
  218.             double fPrime2Integral = 0;
  219.             final double startX = currentX;
  220.             for (int i = 1; i < observations.length; ++i) {
  221.                 // one step forward
  222.                 final double previousX = currentX;
  223.                 final double previousY = currentY;
  224.                 currentX = observations[i].getX();
  225.                 currentY = observations[i].getY();

  226.                 // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
  227.                 // considering a linear model for f (and therefore constant f')
  228.                 final double dx = currentX - previousX;
  229.                 final double dy = currentY - previousY;
  230.                 final double f2StepIntegral =
  231.                     dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
  232.                 final double fPrime2StepIntegral = dy * dy / dx;

  233.                 final double x = currentX - startX;
  234.                 f2Integral += f2StepIntegral;
  235.                 fPrime2Integral += fPrime2StepIntegral;

  236.                 sx2 += x * x;
  237.                 sy2 += f2Integral * f2Integral;
  238.                 sxy += x * f2Integral;
  239.                 sxz += x * fPrime2Integral;
  240.                 syz += f2Integral * fPrime2Integral;
  241.             }

  242.             // compute the amplitude and pulsation coefficients
  243.             double c1 = sy2 * sxz - sxy * syz;
  244.             double c2 = sxy * sxz - sx2 * syz;
  245.             double c3 = sx2 * sy2 - sxy * sxy;
  246.             if (c1 / c2 < 0 || c2 / c3 < 0) {
  247.                 final int last = observations.length - 1;
  248.                 // Range of the observations, assuming that the
  249.                 // observations are sorted.
  250.                 final double xRange = observations[last].getX() - observations[0].getX();
  251.                 if (xRange == 0) {
  252.                     throw new ZeroException();
  253.                 }
  254.                 aOmega[1] = 2 * Math.PI / xRange;

  255.                 double yMin = Double.POSITIVE_INFINITY;
  256.                 double yMax = Double.NEGATIVE_INFINITY;
  257.                 for (int i = 1; i < observations.length; ++i) {
  258.                     final double y = observations[i].getY();
  259.                     if (y < yMin) {
  260.                         yMin = y;
  261.                     }
  262.                     if (y > yMax) {
  263.                         yMax = y;
  264.                     }
  265.                 }
  266.                 aOmega[0] = 0.5 * (yMax - yMin);
  267.             } else {
  268.                 if (c2 == 0) {
  269.                     // In some ill-conditioned cases (cf. MATH-844), the guesser
  270.                     // procedure cannot produce sensible results.
  271.                     throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
  272.                 }

  273.                 aOmega[0] = JdkMath.sqrt(c1 / c2);
  274.                 aOmega[1] = JdkMath.sqrt(c2 / c3);
  275.             }

  276.             return aOmega;
  277.         }

  278.         /**
  279.          * Estimate a first guess of the phase.
  280.          *
  281.          * @param observations Observations, sorted w.r.t. abscissa.
  282.          * @param omega Angular frequency.
  283.          * @return the guessed phase.
  284.          */
  285.         private double guessPhi(WeightedObservedPoint[] observations,
  286.                                 double omega) {
  287.             // initialize the means
  288.             double fcMean = 0;
  289.             double fsMean = 0;

  290.             double currentX = observations[0].getX();
  291.             double currentY = observations[0].getY();
  292.             for (int i = 1; i < observations.length; ++i) {
  293.                 // one step forward
  294.                 final double previousX = currentX;
  295.                 final double previousY = currentY;
  296.                 currentX = observations[i].getX();
  297.                 currentY = observations[i].getY();
  298.                 final double currentYPrime = (currentY - previousY) / (currentX - previousX);

  299.                 double omegaX = omega * currentX;
  300.                 double cosine = JdkMath.cos(omegaX);
  301.                 double sine = JdkMath.sin(omegaX);
  302.                 fcMean += omega * currentY * cosine - currentYPrime * sine;
  303.                 fsMean += omega * currentY * sine + currentYPrime * cosine;
  304.             }

  305.             return JdkMath.atan2(-fsMean, fcMean);
  306.         }
  307.     }
  308. }