    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
    018    package org.apache.commons.math.optimization.fitting;
    019
    020    import org.apache.commons.math.exception.util.LocalizedFormats;
    021    import org.apache.commons.math.optimization.OptimizationException;
    022    import org.apache.commons.math.util.FastMath;
    023
    024    /** This class guesses harmonic coefficients from a sample.
    025
    026     * <p>The algorithm used to guess the coefficients is as follows:</p>
    027
    028     * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
    029     * &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
    030     * </p>
    031     *
    032     * <p>From the analytical expression, we can compute two primitives :
    033     * <pre>
    034     *     If2  (t) = &int; f<sup>2</sup>  = a<sup>2</sup> &times; [t + S (t)] / 2
    035     *     If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
    036     *     where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
    037     * </pre>
    038     * </p>
    039     *
    040     * <p>We can remove S between these expressions :
    041     * <pre>
    042     *     If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
    043     * </pre>
    044     * </p>
    045     *
    046     * <p>The preceding expression shows that If'2 (t) is a linear
    047     * combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
    048     * </p>
    049     *
    050     * <p>From the primitive, we can deduce the same form for definite
    051     * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
    052     * <pre>
    053     *   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>))
    054     * </pre>
    055     * </p>
    056     *
    057     * <p>We can find the coefficients A and B that best fit the sample
    058     * to this linear expression by computing the definite integrals for
    059     * each sample points.
    060     * </p>
    061     *
    062     * <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
    063     * coefficients A and B that minimize a least square criterion
    064     * &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
    065     * <pre>
    066     *
    067     *         &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>
    068     *     A = ------------------------
    069     *         &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>
    070     *
    071     *         &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>
    072     *     B = ------------------------
    073     *         &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>
    074     * </pre>
    075     * </p>
    076     *
    077     *
    078     * <p>In fact, we can assume both a and &omega; are positive and
    079     * compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
    080     * B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
    081     * <pre>
    082     *
    083     * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
    084     *   f  (t<sub>i</sub>)
    085     *   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>)
    086     *   x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
    087     *   y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
    088     *   z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
    089     *   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>
    090     * end for
    091     *
    092     *            |--------------------------
    093     *         \  | &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>
    094     * a     =  \ | ------------------------
    095     *           \| &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>
    096     *
    097     *
    098     *            |--------------------------
    099     *         \  | &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>
    100     * &omega;     =  \ | ------------------------
    101     *           \| &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>
    102     *
    103     * </pre>
    104     * </p>
    105
    106     * <p>Once we know &omega;, we can compute:
    107     * <pre>
    108     *    fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
    109     *    fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
    110     * </pre>
    111     * </p>
    112
    113     * <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
    114     * <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
    115     * expressions to compute &phi;. The best estimate over the sample is
    116     * given by averaging these expressions.
    117     * </p>
    118
    119     * <p>Since integrals and means are involved in the preceding
    120     * estimations, these operations run in O(n) time, where n is the
    121     * number of measurements.</p>
    122
    123     * @version $Revision: 1056034 $ $Date: 2011-01-06 20:41:43 +0100 (jeu. 06 janv. 2011) $
    124     * @since 2.0
    125
    126     */
    127    public class HarmonicCoefficientsGuesser {
    128
    129        /** Sampled observations. */
    130        private final WeightedObservedPoint[] observations;
    131
    132        /** Guessed amplitude. */
    133        private double a;
    134
    135        /** Guessed pulsation &omega;. */
    136        private double omega;
    137
    138        /** Guessed phase &phi;. */
    139        private double phi;
    140
    141        /** Simple constructor.
    142         * @param observations sampled observations
    143         */
    144        public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) {
    145            this.observations = observations.clone();
    146            a                 = Double.NaN;
    147            omega             = Double.NaN;
    148        }
    149
    150        /** Estimate a first guess of the coefficients.
    151         * @exception OptimizationException if the sample is too short or if
    152         * the first guess cannot be computed (when the elements under the
    153         * square roots are negative).
    154         * */
    155        public void guess() throws OptimizationException {
    156            sortObservations();
    157            guessAOmega();
    158            guessPhi();
    159        }
    160
    161        /** Sort the observations with respect to the abscissa.
    162         */
    163        private void sortObservations() {
    164
    165            // Since the samples are almost always already sorted, this
    166            // method is implemented as an insertion sort that reorders the
    167            // elements in place. Insertion sort is very efficient in this case.
    168            WeightedObservedPoint curr = observations[0];
    169            for (int j = 1; j < observations.length; ++j) {
    170                WeightedObservedPoint prec = curr;
    171                curr = observations[j];
    172                if (curr.getX() < prec.getX()) {
    173                    // the current element should be inserted closer to the beginning
    174                    int i = j - 1;
    175                    WeightedObservedPoint mI = observations[i];
    176                    while ((i >= 0) && (curr.getX() < mI.getX())) {
    177                        observations[i + 1] = mI;
    178                        if (i-- != 0) {
    179                            mI = observations[i];
    180                        }
    181                    }
    182                    observations[i + 1] = curr;
    183                    curr = observations[j];
    184                }
    185            }
    186
    187        }
    188
    189        /** Estimate a first guess of the a and &omega; coefficients.
    190         * @exception OptimizationException if the sample is too short or if
    191         * the first guess cannot be computed (when the elements under the
    192         * square roots are negative).
    193         */
    194        private void guessAOmega() throws OptimizationException {
    195
    196            // initialize the sums for the linear model between the two integrals
    197            double sx2 = 0.0;
    198            double sy2 = 0.0;
    199            double sxy = 0.0;
    200            double sxz = 0.0;
    201            double syz = 0.0;
    202
    203            double currentX        = observations[0].getX();
    204            double currentY        = observations[0].getY();
    205            double f2Integral      = 0;
    206            double fPrime2Integral = 0;
    207            final double startX = currentX;
    208            for (int i = 1; i < observations.length; ++i) {
    209
    210                // one step forward
    211                final double previousX = currentX;
    212                final double previousY = currentY;
    213                currentX = observations[i].getX();
    214                currentY = observations[i].getY();
    215
    216                // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
    217                // considering a linear model for f (and therefore constant f')
    218                final double dx = currentX - previousX;
    219                final double dy = currentY - previousY;
    220                final double f2StepIntegral =
    221                    dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
    222                final double fPrime2StepIntegral = dy * dy / dx;
    223
    224                final double x   = currentX - startX;
    225                f2Integral      += f2StepIntegral;
    226                fPrime2Integral += fPrime2StepIntegral;
    227
    228                sx2 += x * x;
    229                sy2 += f2Integral * f2Integral;
    230                sxy += x * f2Integral;
    231                sxz += x * fPrime2Integral;
    232                syz += f2Integral * fPrime2Integral;
    233
    234            }
    235
    236            // compute the amplitude and pulsation coefficients
    237            double c1 = sy2 * sxz - sxy * syz;
    238            double c2 = sxy * sxz - sx2 * syz;
    239            double c3 = sx2 * sy2 - sxy * sxy;
    240            if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) {
    241                throw new OptimizationException(LocalizedFormats.UNABLE_TO_FIRST_GUESS_HARMONIC_COEFFICIENTS);
    242            }
    243            a     = FastMath.sqrt(c1 / c2);
    244            omega = FastMath.sqrt(c2 / c3);
    245
    246        }
    247
    248        /** Estimate a first guess of the &phi; coefficient.
    249         */
    250        private void guessPhi() {
    251
    252            // initialize the means
    253            double fcMean = 0.0;
    254            double fsMean = 0.0;
    255
    256            double currentX = observations[0].getX();
    257            double currentY = observations[0].getY();
    258            for (int i = 1; i < observations.length; ++i) {
    259
    260                // one step forward
    261                final double previousX = currentX;
    262                final double previousY = currentY;
    263                currentX = observations[i].getX();
    264                currentY = observations[i].getY();
    265                final double currentYPrime = (currentY - previousY) / (currentX - previousX);
    266
    267                double   omegaX = omega * currentX;
    268                double   cosine = FastMath.cos(omegaX);
    269                double   sine   = FastMath.sin(omegaX);
    270                fcMean += omega * currentY * cosine - currentYPrime *   sine;
    271                fsMean += omega * currentY *   sine + currentYPrime * cosine;
    272
    273            }
    274
    275            phi = FastMath.atan2(-fsMean, fcMean);
    276
    277        }
    278
    279        /** Get the guessed amplitude a.
    280         * @return guessed amplitude a;
    281         */
    282        public double getGuessedAmplitude() {
    283            return a;
    284        }
    285
    286        /** Get the guessed pulsation &omega;.
    287         * @return guessed pulsation &omega;
    288         */
    289        public double getGuessedPulsation() {
    290            return omega;
    291        }
    292
    293        /** Get the guessed phase &phi;.
    294         * @return guessed phase &phi;
    295         */
    296        public double getGuessedPhase() {
    297            return phi;
    298        }
    299
    300    }