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 */
017package org.apache.commons.math3.fitting;
018
019import java.util.ArrayList;
020import java.util.Collection;
021import java.util.List;
022
023import org.apache.commons.math3.analysis.function.HarmonicOscillator;
024import org.apache.commons.math3.exception.MathIllegalStateException;
025import org.apache.commons.math3.exception.NumberIsTooSmallException;
026import org.apache.commons.math3.exception.ZeroException;
027import org.apache.commons.math3.exception.util.LocalizedFormats;
028import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder;
029import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem;
030import org.apache.commons.math3.linear.DiagonalMatrix;
031import org.apache.commons.math3.util.FastMath;
032
033/**
034 * Fits points to a {@link
035 * org.apache.commons.math3.analysis.function.HarmonicOscillator.Parametric harmonic oscillator}
036 * function.
037 * <br/>
038 * The {@link #withStartPoint(double[]) initial guess values} must be passed
039 * in the following order:
040 * <ul>
041 *  <li>Amplitude</li>
042 *  <li>Angular frequency</li>
043 *  <li>phase</li>
044 * </ul>
045 * The optimal values will be returned in the same order.
046 *
047 * @version $Id: HarmonicCurveFitter.html 908881 2014-05-15 07:10:28Z luc $
048 * @since 3.3
049 */
050public class HarmonicCurveFitter extends AbstractCurveFitter {
051    /** Parametric function to be fitted. */
052    private static final HarmonicOscillator.Parametric FUNCTION = new HarmonicOscillator.Parametric();
053    /** Initial guess. */
054    private final double[] initialGuess;
055    /** Maximum number of iterations of the optimization algorithm. */
056    private final int maxIter;
057
058    /**
059     * Contructor used by the factory methods.
060     *
061     * @param initialGuess Initial guess. If set to {@code null}, the initial guess
062     * will be estimated using the {@link ParameterGuesser}.
063     * @param maxIter Maximum number of iterations of the optimization algorithm.
064     */
065    private HarmonicCurveFitter(double[] initialGuess,
066                                int maxIter) {
067        this.initialGuess = initialGuess;
068        this.maxIter = maxIter;
069    }
070
071    /**
072     * Creates a default curve fitter.
073     * The initial guess for the parameters will be {@link ParameterGuesser}
074     * computed automatically, and the maximum number of iterations of the
075     * optimization algorithm is set to {@link Integer#MAX_VALUE}.
076     *
077     * @return a curve fitter.
078     *
079     * @see #withStartPoint(double[])
080     * @see #withMaxIterations(int)
081     */
082    public static HarmonicCurveFitter create() {
083        return new HarmonicCurveFitter(null, Integer.MAX_VALUE);
084    }
085
086    /**
087     * Configure the start point (initial guess).
088     * @param newStart new start point (initial guess)
089     * @return a new instance.
090     */
091    public HarmonicCurveFitter withStartPoint(double[] newStart) {
092        return new HarmonicCurveFitter(newStart.clone(),
093                                       maxIter);
094    }
095
096    /**
097     * Configure the maximum number of iterations.
098     * @param newMaxIter maximum number of iterations
099     * @return a new instance.
100     */
101    public HarmonicCurveFitter withMaxIterations(int newMaxIter) {
102        return new HarmonicCurveFitter(initialGuess,
103                                       newMaxIter);
104    }
105
106    /** {@inheritDoc} */
107    @Override
108    protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> observations) {
109        // Prepare least-squares problem.
110        final int len = observations.size();
111        final double[] target  = new double[len];
112        final double[] weights = new double[len];
113
114        int i = 0;
115        for (WeightedObservedPoint obs : observations) {
116            target[i]  = obs.getY();
117            weights[i] = obs.getWeight();
118            ++i;
119        }
120
121        final AbstractCurveFitter.TheoreticalValuesFunction model
122            = new AbstractCurveFitter.TheoreticalValuesFunction(FUNCTION,
123                                                                observations);
124
125        final double[] startPoint = initialGuess != null ?
126            initialGuess :
127            // Compute estimation.
128            new ParameterGuesser(observations).guess();
129
130        // Return a new optimizer set up to fit a Gaussian curve to the
131        // observed points.
132        return new LeastSquaresBuilder().
133                maxEvaluations(Integer.MAX_VALUE).
134                maxIterations(maxIter).
135                start(startPoint).
136                target(target).
137                weight(new DiagonalMatrix(weights)).
138                model(model.getModelFunction(), model.getModelFunctionJacobian()).
139                build();
140
141    }
142
143    /**
144     * This class guesses harmonic coefficients from a sample.
145     * <p>The algorithm used to guess the coefficients is as follows:</p>
146     *
147     * <p>We know \( f(t) \) at some sampling points \( t_i \) and want
148     * to find \( a \), \( \omega \) and \( \phi \) such that
149     * \( f(t) = a \cos (\omega t + \phi) \).
150     * </p>
151     *
152     * <p>From the analytical expression, we can compute two primitives :
153     * \[
154     *     If2(t) = \int f^2 dt  = a^2 (t + S(t)) / 2
155     * \]
156     * \[
157     *     If'2(t) = \int f'^2 dt = a^2 \omega^2 (t - S(t)) / 2
158     * \]
159     * where \(S(t) = \frac{\sin(2 (\omega t + \phi))}{2\omega}\)
160     * </p>
161     *
162     * <p>We can remove \(S\) between these expressions :
163     * \[
164     *     If'2(t) = a^2 \omega^2 t - \omega^2 If2(t)
165     * \]
166     * </p>
167     *
168     * <p>The preceding expression shows that \(If'2 (t)\) is a linear
169     * combination of both \(t\) and \(If2(t)\):
170     * \[
171     *   If'2(t) = A t + B If2(t)
172     * \]
173     * </p>
174     *
175     * <p>From the primitive, we can deduce the same form for definite
176     * integrals between \(t_1\) and \(t_i\) for each \(t_i\) :
177     * \[
178     *   If2(t_i) - If2(t_1) = A (t_i - t_1) + B (If2 (t_i) - If2(t_1))
179     * \]
180     * </p>
181     *
182     * <p>We can find the coefficients \(A\) and \(B\) that best fit the sample
183     * to this linear expression by computing the definite integrals for
184     * each sample points.
185     * </p>
186     *
187     * <p>For a bilinear expression \(z(x_i, y_i) = A x_i + B y_i\), the
188     * coefficients \(A\) and \(B\) that minimize a least-squares criterion
189     * \(\sum (z_i - z(x_i, y_i))^2\) are given by these expressions:</p>
190     * \[
191     *   A = \frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i}
192     *            {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
193     * \]
194     * \[
195     *   B = \frac{\sum x_i x_i \sum y_i z_i - \sum x_i y_i \sum x_i z_i}
196     *            {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}
197     *
198     * \]
199     *
200     * <p>In fact, we can assume that both \(a\) and \(\omega\) are positive and
201     * compute them directly, knowing that \(A = a^2 \omega^2\) and that
202     * \(B = -\omega^2\). The complete algorithm is therefore:</p>
203     *
204     * For each \(t_i\) from \(t_1\) to \(t_{n-1}\), compute:
205     * \[ f(t_i) \]
206     * \[ f'(t_i) = \frac{f (t_{i+1}) - f(t_{i-1})}{t_{i+1} - t_{i-1}} \]
207     * \[ x_i = t_i  - t_1 \]
208     * \[ y_i = \int_{t_1}^{t_i} f^2(t) dt \]
209     * \[ z_i = \int_{t_1}^{t_i} f'^2(t) dt \]
210     * and update the sums:
211     * \[ \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 \]
212     *
213     * Then:
214     * \[
215     *  a = \sqrt{\frac{\sum y_i y_i  \sum x_i z_i - \sum x_i y_i \sum y_i z_i }
216     *                 {\sum x_i y_i  \sum x_i z_i - \sum x_i x_i \sum y_i z_i }}
217     * \]
218     * \[
219     *  \omega = \sqrt{\frac{\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i}
220     *                      {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}}
221     * \]
222     *
223     * <p>Once we know \(\omega\) we can compute:
224     * \[
225     *    fc = \omega f(t) \cos(\omega t) - f'(t) \sin(\omega t)
226     * \]
227     * \[
228     *    fs = \omega f(t) \sin(\omega t) + f'(t) \cos(\omega t)
229     * \]
230     * </p>
231     *
232     * <p>It appears that \(fc = a \omega \cos(\phi)\) and
233     * \(fs = -a \omega \sin(\phi)\), so we can use these
234     * expressions to compute \(\phi\). The best estimate over the sample is
235     * given by averaging these expressions.
236     * </p>
237     *
238     * <p>Since integrals and means are involved in the preceding
239     * estimations, these operations run in \(O(n)\) time, where \(n\) is the
240     * number of measurements.</p>
241     */
242    public static class ParameterGuesser {
243        /** Amplitude. */
244        private final double a;
245        /** Angular frequency. */
246        private final double omega;
247        /** Phase. */
248        private final double phi;
249
250        /**
251         * Simple constructor.
252         *
253         * @param observations Sampled observations.
254         * @throws NumberIsTooSmallException if the sample is too short.
255         * @throws ZeroException if the abscissa range is zero.
256         * @throws MathIllegalStateException when the guessing procedure cannot
257         * produce sensible results.
258         */
259        public ParameterGuesser(Collection<WeightedObservedPoint> observations) {
260            if (observations.size() < 4) {
261                throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
262                                                    observations.size(), 4, true);
263            }
264
265            final WeightedObservedPoint[] sorted
266                = sortObservations(observations).toArray(new WeightedObservedPoint[0]);
267
268            final double aOmega[] = guessAOmega(sorted);
269            a = aOmega[0];
270            omega = aOmega[1];
271
272            phi = guessPhi(sorted);
273        }
274
275        /**
276         * Gets an estimation of the parameters.
277         *
278         * @return the guessed parameters, in the following order:
279         * <ul>
280         *  <li>Amplitude</li>
281         *  <li>Angular frequency</li>
282         *  <li>Phase</li>
283         * </ul>
284         */
285        public double[] guess() {
286            return new double[] { a, omega, phi };
287        }
288
289        /**
290         * Sort the observations with respect to the abscissa.
291         *
292         * @param unsorted Input observations.
293         * @return the input observations, sorted.
294         */
295        private List<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) {
296            final List<WeightedObservedPoint> observations = new ArrayList<WeightedObservedPoint>(unsorted);
297
298            // Since the samples are almost always already sorted, this
299            // method is implemented as an insertion sort that reorders the
300            // elements in place. Insertion sort is very efficient in this case.
301            WeightedObservedPoint curr = observations.get(0);
302            final int len = observations.size();
303            for (int j = 1; j < len; j++) {
304                WeightedObservedPoint prec = curr;
305                curr = observations.get(j);
306                if (curr.getX() < prec.getX()) {
307                    // the current element should be inserted closer to the beginning
308                    int i = j - 1;
309                    WeightedObservedPoint mI = observations.get(i);
310                    while ((i >= 0) && (curr.getX() < mI.getX())) {
311                        observations.set(i + 1, mI);
312                        if (i-- != 0) {
313                            mI = observations.get(i);
314                        }
315                    }
316                    observations.set(i + 1, curr);
317                    curr = observations.get(j);
318                }
319            }
320
321            return observations;
322        }
323
324        /**
325         * Estimate a first guess of the amplitude and angular frequency.
326         *
327         * @param observations Observations, sorted w.r.t. abscissa.
328         * @throws ZeroException if the abscissa range is zero.
329         * @throws MathIllegalStateException when the guessing procedure cannot
330         * produce sensible results.
331         * @return the guessed amplitude (at index 0) and circular frequency
332         * (at index 1).
333         */
334        private double[] guessAOmega(WeightedObservedPoint[] observations) {
335            final double[] aOmega = new double[2];
336
337            // initialize the sums for the linear model between the two integrals
338            double sx2 = 0;
339            double sy2 = 0;
340            double sxy = 0;
341            double sxz = 0;
342            double syz = 0;
343
344            double currentX = observations[0].getX();
345            double currentY = observations[0].getY();
346            double f2Integral = 0;
347            double fPrime2Integral = 0;
348            final double startX = currentX;
349            for (int i = 1; i < observations.length; ++i) {
350                // one step forward
351                final double previousX = currentX;
352                final double previousY = currentY;
353                currentX = observations[i].getX();
354                currentY = observations[i].getY();
355
356                // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
357                // considering a linear model for f (and therefore constant f')
358                final double dx = currentX - previousX;
359                final double dy = currentY - previousY;
360                final double f2StepIntegral =
361                    dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
362                final double fPrime2StepIntegral = dy * dy / dx;
363
364                final double x = currentX - startX;
365                f2Integral += f2StepIntegral;
366                fPrime2Integral += fPrime2StepIntegral;
367
368                sx2 += x * x;
369                sy2 += f2Integral * f2Integral;
370                sxy += x * f2Integral;
371                sxz += x * fPrime2Integral;
372                syz += f2Integral * fPrime2Integral;
373            }
374
375            // compute the amplitude and pulsation coefficients
376            double c1 = sy2 * sxz - sxy * syz;
377            double c2 = sxy * sxz - sx2 * syz;
378            double c3 = sx2 * sy2 - sxy * sxy;
379            if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
380                final int last = observations.length - 1;
381                // Range of the observations, assuming that the
382                // observations are sorted.
383                final double xRange = observations[last].getX() - observations[0].getX();
384                if (xRange == 0) {
385                    throw new ZeroException();
386                }
387                aOmega[1] = 2 * Math.PI / xRange;
388
389                double yMin = Double.POSITIVE_INFINITY;
390                double yMax = Double.NEGATIVE_INFINITY;
391                for (int i = 1; i < observations.length; ++i) {
392                    final double y = observations[i].getY();
393                    if (y < yMin) {
394                        yMin = y;
395                    }
396                    if (y > yMax) {
397                        yMax = y;
398                    }
399                }
400                aOmega[0] = 0.5 * (yMax - yMin);
401            } else {
402                if (c2 == 0) {
403                    // In some ill-conditioned cases (cf. MATH-844), the guesser
404                    // procedure cannot produce sensible results.
405                    throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
406                }
407
408                aOmega[0] = FastMath.sqrt(c1 / c2);
409                aOmega[1] = FastMath.sqrt(c2 / c3);
410            }
411
412            return aOmega;
413        }
414
415        /**
416         * Estimate a first guess of the phase.
417         *
418         * @param observations Observations, sorted w.r.t. abscissa.
419         * @return the guessed phase.
420         */
421        private double guessPhi(WeightedObservedPoint[] observations) {
422            // initialize the means
423            double fcMean = 0;
424            double fsMean = 0;
425
426            double currentX = observations[0].getX();
427            double currentY = observations[0].getY();
428            for (int i = 1; i < observations.length; ++i) {
429                // one step forward
430                final double previousX = currentX;
431                final double previousY = currentY;
432                currentX = observations[i].getX();
433                currentY = observations[i].getY();
434                final double currentYPrime = (currentY - previousY) / (currentX - previousX);
435
436                double omegaX = omega * currentX;
437                double cosine = FastMath.cos(omegaX);
438                double sine = FastMath.sin(omegaX);
439                fcMean += omega * currentY * cosine - currentYPrime * sine;
440                fsMean += omega * currentY * sine + currentYPrime * cosine;
441            }
442
443            return FastMath.atan2(-fsMean, fcMean);
444        }
445    }
446}