001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math3.optimization.fitting;
019
020import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
021import org.apache.commons.math3.analysis.function.HarmonicOscillator;
022import org.apache.commons.math3.exception.ZeroException;
023import org.apache.commons.math3.exception.NumberIsTooSmallException;
024import org.apache.commons.math3.exception.MathIllegalStateException;
025import org.apache.commons.math3.exception.util.LocalizedFormats;
026import org.apache.commons.math3.util.FastMath;
027
028/**
029 * Class that implements a curve fitting specialized for sinusoids.
030 *
031 * Harmonic fitting is a very simple case of curve fitting. The
032 * estimated coefficients are the amplitude a, the pulsation ω and
033 * the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
034 * searched by a least square estimator initialized with a rough guess
035 * based on integrals.
036 *
037 * @deprecated As of 3.1 (to be removed in 4.0).
038 * @since 2.0
039 */
040@Deprecated
041public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> {
042    /**
043     * Simple constructor.
044     * @param optimizer Optimizer to use for the fitting.
045     */
046    public HarmonicFitter(final DifferentiableMultivariateVectorOptimizer optimizer) {
047        super(optimizer);
048    }
049
050    /**
051     * Fit an harmonic function to the observed points.
052     *
053     * @param initialGuess First guess values in the following order:
054     * <ul>
055     *  <li>Amplitude</li>
056     *  <li>Angular frequency</li>
057     *  <li>Phase</li>
058     * </ul>
059     * @return the parameters of the harmonic function that best fits the
060     * observed points (in the same order as above).
061     */
062    public double[] fit(double[] initialGuess) {
063        return fit(new HarmonicOscillator.Parametric(), initialGuess);
064    }
065
066    /**
067     * Fit an harmonic function to the observed points.
068     * An initial guess will be automatically computed.
069     *
070     * @return the parameters of the harmonic function that best fits the
071     * observed points (see the other {@link #fit(double[]) fit} method.
072     * @throws NumberIsTooSmallException if the sample is too short for the
073     * the first guess to be computed.
074     * @throws ZeroException if the first guess cannot be computed because
075     * the abscissa range is zero.
076     */
077    public double[] fit() {
078        return fit((new ParameterGuesser(getObservations())).guess());
079    }
080
081    /**
082     * This class guesses harmonic coefficients from a sample.
083     * <p>The algorithm used to guess the coefficients is as follows:</p>
084     *
085     * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
086     * &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
087     * </p>
088     *
089     * <p>From the analytical expression, we can compute two primitives :
090     * <pre>
091     *     If2  (t) = &int; f<sup>2</sup>  = a<sup>2</sup> &times; [t + S (t)] / 2
092     *     If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
093     *     where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
094     * </pre>
095     * </p>
096     *
097     * <p>We can remove S between these expressions :
098     * <pre>
099     *     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}