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    }