View Javadoc

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  
18  package org.apache.commons.math.optimization.fitting;
19  
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.MathRuntimeException;
22  import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
23  import org.apache.commons.math.optimization.OptimizationException;
24  
25  /** This class implements a curve fitting specialized for sinusoids.
26   * <p>Harmonic fitting is a very simple case of curve fitting. The
27   * estimated coefficients are the amplitude a, the pulsation &omega; and
28   * the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
29   * searched by a least square estimator initialized with a rough guess
30   * based on integrals.</p>
31   * @version $Revision: 786479 $ $Date: 2009-06-19 08:36:16 -0400 (Fri, 19 Jun 2009) $
32   * @since 2.0
33   */
34  public class HarmonicFitter {
35  
36      /** Fitter for the coefficients. */
37      private final CurveFitter fitter;
38  
39      /** Values for amplitude, pulsation &omega; and phase &phi;. */
40      private double[] parameters;
41  
42      /** Simple constructor.
43       * @param optimizer optimizer to use for the fitting
44       */
45      public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) {
46          this.fitter = new CurveFitter(optimizer);
47          parameters  = null;
48      }
49  
50      /** Simple constructor.
51       * <p>This constructor can be used when a first guess of the
52       * coefficients is already known.</p>
53       * @param optimizer optimizer to use for the fitting
54       * @param initialGuess guessed values for amplitude (index 0),
55       * pulsation &omega; (index 1) and phase &phi; (index 2)
56       */
57      public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer,
58                            final double[] initialGuess) {
59          this.fitter     = new CurveFitter(optimizer);
60          this.parameters = initialGuess.clone();
61      }
62  
63      /** Add an observed weighted (x,y) point to the sample.
64       * @param weight weight of the observed point in the fit
65       * @param x abscissa of the point
66       * @param y observed value of the point at x, after fitting we should
67       * have P(x) as close as possible to this value
68       */
69      public void addObservedPoint(double weight, double x, double y) {
70          fitter.addObservedPoint(weight, x, y);
71      }
72  
73      /** Fit an harmonic function to the observed points.
74       * @return harmonic function best fitting the observed points
75       * @throws OptimizationException if the sample is too short or if
76       * the first guess cannot be computed
77       */
78      public HarmonicFunction fit() throws OptimizationException {
79          try {
80  
81              // shall we compute the first guess of the parameters ourselves ?
82              if (parameters == null) {
83                  final WeightedObservedPoint[] observations = fitter.getObservations();
84                  if (observations.length < 4) {
85                      throw new OptimizationException("sample contains {0} observed points, at least {1} are required",
86                                                      observations.length, 4);
87                  }
88  
89                  HarmonicCoefficientsGuesser guesser = new HarmonicCoefficientsGuesser(observations);
90                  guesser.guess();
91                  parameters = new double[] {
92                                   guesser.getGuessedAmplitude(),
93                                   guesser.getGuessedPulsation(),
94                                   guesser.getGuessedPhase()
95                              };
96  
97              }
98  
99              double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters);
100             return new HarmonicFunction(fitted[0], fitted[1], fitted[2]);
101 
102         } catch (FunctionEvaluationException fee) {
103             // this should never happen
104             throw MathRuntimeException.createInternalError(fee);
105         }
106     }
107 
108     /** Parametric harmonic function. */
109     private static class ParametricHarmonicFunction implements ParametricRealFunction {
110 
111         /** {@inheritDoc} */
112         public double value(double x, double[] parameters) {
113             final double a     = parameters[0];
114             final double omega = parameters[1];
115             final double phi   = parameters[2];
116             return a * Math.cos(omega * x + phi);
117         }
118 
119         /** {@inheritDoc} */
120         public double[] gradient(double x, double[] parameters) {
121             final double a     = parameters[0];
122             final double omega = parameters[1];
123             final double phi   = parameters[2];
124             final double alpha = omega * x + phi;
125             final double cosAlpha = Math.cos(alpha);
126             final double sinAlpha = Math.sin(alpha);
127             return new double[] { cosAlpha, -a * x * sinAlpha, -a * sinAlpha };
128         }
129 
130     }
131 
132 }