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.math3.optimization.fitting; 019 020 import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer; 021 import org.apache.commons.math3.analysis.function.HarmonicOscillator; 022 import org.apache.commons.math3.exception.ZeroException; 023 import org.apache.commons.math3.exception.NumberIsTooSmallException; 024 import org.apache.commons.math3.exception.MathIllegalStateException; 025 import org.apache.commons.math3.exception.util.LocalizedFormats; 026 import 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 φ: <code>f (t) = a cos (ω t + φ)</code>. They are 034 * searched by a least square estimator initialized with a rough guess 035 * based on integrals. 036 * 037 * @version $Id: HarmonicFitter.java 1422230 2012-12-15 12:11:13Z erans $ 038 * @deprecated As of 3.1 (to be removed in 4.0). 039 * @since 2.0 040 */ 041 @Deprecated 042 public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> { 043 /** 044 * Simple constructor. 045 * @param optimizer Optimizer to use for the fitting. 046 */ 047 public HarmonicFitter(final DifferentiableMultivariateVectorOptimizer optimizer) { 048 super(optimizer); 049 } 050 051 /** 052 * Fit an harmonic function to the observed points. 053 * 054 * @param initialGuess First guess values in the following order: 055 * <ul> 056 * <li>Amplitude</li> 057 * <li>Angular frequency</li> 058 * <li>Phase</li> 059 * </ul> 060 * @return the parameters of the harmonic function that best fits the 061 * observed points (in the same order as above). 062 */ 063 public double[] fit(double[] initialGuess) { 064 return fit(new HarmonicOscillator.Parametric(), initialGuess); 065 } 066 067 /** 068 * Fit an harmonic function to the observed points. 069 * An initial guess will be automatically computed. 070 * 071 * @return the parameters of the harmonic function that best fits the 072 * observed points (see the other {@link #fit(double[]) fit} method. 073 * @throws NumberIsTooSmallException if the sample is too short for the 074 * the first guess to be computed. 075 * @throws ZeroException if the first guess cannot be computed because 076 * the abscissa range is zero. 077 */ 078 public double[] fit() { 079 return fit((new ParameterGuesser(getObservations())).guess()); 080 } 081 082 /** 083 * This class guesses harmonic coefficients from a sample. 084 * <p>The algorithm used to guess the coefficients is as follows:</p> 085 * 086 * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a, 087 * ω and φ such that f (t) = a cos (ω t + φ). 088 * </p> 089 * 090 * <p>From the analytical expression, we can compute two primitives : 091 * <pre> 092 * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2 093 * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2 094 * where S (t) = sin (2 (ω t + φ)) / (2 ω) 095 * </pre> 096 * </p> 097 * 098 * <p>We can remove S between these expressions : 099 * <pre> 100 * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t) 101 * </pre> 102 * </p> 103 * 104 * <p>The preceding expression shows that If'2 (t) is a linear 105 * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t) 106 * </p> 107 * 108 * <p>From the primitive, we can deduce the same form for definite 109 * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> : 110 * <pre> 111 * If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A × (t<sub>i</sub> - t<sub>1</sub>) + B × (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>)) 112 * </pre> 113 * </p> 114 * 115 * <p>We can find the coefficients A and B that best fit the sample 116 * to this linear expression by computing the definite integrals for 117 * each sample points. 118 * </p> 119 * 120 * <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A × x<sub>i</sub> + B × y<sub>i</sub>, the 121 * coefficients A and B that minimize a least square criterion 122 * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p> 123 * <pre> 124 * 125 * ∑y<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> 126 * A = ------------------------ 127 * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub> 128 * 129 * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> 130 * B = ------------------------ 131 * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub> 132 * </pre> 133 * </p> 134 * 135 * 136 * <p>In fact, we can assume both a and ω are positive and 137 * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that 138 * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p> 139 * <pre> 140 * 141 * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute: 142 * f (t<sub>i</sub>) 143 * 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>) 144 * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub> 145 * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 146 * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 147 * update the sums ∑x<sub>i</sub>x<sub>i</sub>, ∑y<sub>i</sub>y<sub>i</sub>, ∑x<sub>i</sub>y<sub>i</sub>, ∑x<sub>i</sub>z<sub>i</sub> and ∑y<sub>i</sub>z<sub>i</sub> 148 * end for 149 * 150 * |-------------------------- 151 * \ | ∑y<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> 152 * a = \ | ------------------------ 153 * \| ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> 154 * 155 * 156 * |-------------------------- 157 * \ | ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> 158 * ω = \ | ------------------------ 159 * \| ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub> 160 * 161 * </pre> 162 * </p> 163 * 164 * <p>Once we know ω, we can compute: 165 * <pre> 166 * fc = ω f (t) cos (ω t) - f' (t) sin (ω t) 167 * fs = ω f (t) sin (ω t) + f' (t) cos (ω t) 168 * </pre> 169 * </p> 170 * 171 * <p>It appears that <code>fc = a ω cos (φ)</code> and 172 * <code>fs = -a ω sin (φ)</code>, so we can use these 173 * expressions to compute φ. The best estimate over the sample is 174 * given by averaging these expressions. 175 * </p> 176 * 177 * <p>Since integrals and means are involved in the preceding 178 * estimations, these operations run in O(n) time, where n is the 179 * number of measurements.</p> 180 */ 181 public static class ParameterGuesser { 182 /** Amplitude. */ 183 private final double a; 184 /** Angular frequency. */ 185 private final double omega; 186 /** Phase. */ 187 private final double phi; 188 189 /** 190 * Simple constructor. 191 * 192 * @param observations Sampled observations. 193 * @throws NumberIsTooSmallException if the sample is too short. 194 * @throws ZeroException if the abscissa range is zero. 195 * @throws MathIllegalStateException when the guessing procedure cannot 196 * produce sensible results. 197 */ 198 public ParameterGuesser(WeightedObservedPoint[] observations) { 199 if (observations.length < 4) { 200 throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, 201 observations.length, 4, true); 202 } 203 204 final WeightedObservedPoint[] sorted = sortObservations(observations); 205 206 final double aOmega[] = guessAOmega(sorted); 207 a = aOmega[0]; 208 omega = aOmega[1]; 209 210 phi = guessPhi(sorted); 211 } 212 213 /** 214 * Gets an estimation of the parameters. 215 * 216 * @return the guessed parameters, in the following order: 217 * <ul> 218 * <li>Amplitude</li> 219 * <li>Angular frequency</li> 220 * <li>Phase</li> 221 * </ul> 222 */ 223 public double[] guess() { 224 return new double[] { a, omega, phi }; 225 } 226 227 /** 228 * Sort the observations with respect to the abscissa. 229 * 230 * @param unsorted Input observations. 231 * @return the input observations, sorted. 232 */ 233 private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) { 234 final WeightedObservedPoint[] observations = unsorted.clone(); 235 236 // Since the samples are almost always already sorted, this 237 // method is implemented as an insertion sort that reorders the 238 // elements in place. Insertion sort is very efficient in this case. 239 WeightedObservedPoint curr = observations[0]; 240 for (int j = 1; j < observations.length; ++j) { 241 WeightedObservedPoint prec = curr; 242 curr = observations[j]; 243 if (curr.getX() < prec.getX()) { 244 // the current element should be inserted closer to the beginning 245 int i = j - 1; 246 WeightedObservedPoint mI = observations[i]; 247 while ((i >= 0) && (curr.getX() < mI.getX())) { 248 observations[i + 1] = mI; 249 if (i-- != 0) { 250 mI = observations[i]; 251 } 252 } 253 observations[i + 1] = curr; 254 curr = observations[j]; 255 } 256 } 257 258 return observations; 259 } 260 261 /** 262 * Estimate a first guess of the amplitude and angular frequency. 263 * This method assumes that the {@link #sortObservations()} method 264 * has been called previously. 265 * 266 * @param observations Observations, sorted w.r.t. abscissa. 267 * @throws ZeroException if the abscissa range is zero. 268 * @throws MathIllegalStateException when the guessing procedure cannot 269 * produce sensible results. 270 * @return the guessed amplitude (at index 0) and circular frequency 271 * (at index 1). 272 */ 273 private double[] guessAOmega(WeightedObservedPoint[] observations) { 274 final double[] aOmega = new double[2]; 275 276 // initialize the sums for the linear model between the two integrals 277 double sx2 = 0; 278 double sy2 = 0; 279 double sxy = 0; 280 double sxz = 0; 281 double syz = 0; 282 283 double currentX = observations[0].getX(); 284 double currentY = observations[0].getY(); 285 double f2Integral = 0; 286 double fPrime2Integral = 0; 287 final double startX = currentX; 288 for (int i = 1; i < observations.length; ++i) { 289 // one step forward 290 final double previousX = currentX; 291 final double previousY = currentY; 292 currentX = observations[i].getX(); 293 currentY = observations[i].getY(); 294 295 // update the integrals of f<sup>2</sup> and f'<sup>2</sup> 296 // considering a linear model for f (and therefore constant f') 297 final double dx = currentX - previousX; 298 final double dy = currentY - previousY; 299 final double f2StepIntegral = 300 dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3; 301 final double fPrime2StepIntegral = dy * dy / dx; 302 303 final double x = currentX - startX; 304 f2Integral += f2StepIntegral; 305 fPrime2Integral += fPrime2StepIntegral; 306 307 sx2 += x * x; 308 sy2 += f2Integral * f2Integral; 309 sxy += x * f2Integral; 310 sxz += x * fPrime2Integral; 311 syz += f2Integral * fPrime2Integral; 312 } 313 314 // compute the amplitude and pulsation coefficients 315 double c1 = sy2 * sxz - sxy * syz; 316 double c2 = sxy * sxz - sx2 * syz; 317 double c3 = sx2 * sy2 - sxy * sxy; 318 if ((c1 / c2 < 0) || (c2 / c3 < 0)) { 319 final int last = observations.length - 1; 320 // Range of the observations, assuming that the 321 // observations are sorted. 322 final double xRange = observations[last].getX() - observations[0].getX(); 323 if (xRange == 0) { 324 throw new ZeroException(); 325 } 326 aOmega[1] = 2 * Math.PI / xRange; 327 328 double yMin = Double.POSITIVE_INFINITY; 329 double yMax = Double.NEGATIVE_INFINITY; 330 for (int i = 1; i < observations.length; ++i) { 331 final double y = observations[i].getY(); 332 if (y < yMin) { 333 yMin = y; 334 } 335 if (y > yMax) { 336 yMax = y; 337 } 338 } 339 aOmega[0] = 0.5 * (yMax - yMin); 340 } else { 341 if (c2 == 0) { 342 // In some ill-conditioned cases (cf. MATH-844), the guesser 343 // procedure cannot produce sensible results. 344 throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR); 345 } 346 347 aOmega[0] = FastMath.sqrt(c1 / c2); 348 aOmega[1] = FastMath.sqrt(c2 / c3); 349 } 350 351 return aOmega; 352 } 353 354 /** 355 * Estimate a first guess of the phase. 356 * 357 * @param observations Observations, sorted w.r.t. abscissa. 358 * @return the guessed phase. 359 */ 360 private double guessPhi(WeightedObservedPoint[] observations) { 361 // initialize the means 362 double fcMean = 0; 363 double fsMean = 0; 364 365 double currentX = observations[0].getX(); 366 double currentY = observations[0].getY(); 367 for (int i = 1; i < observations.length; ++i) { 368 // one step forward 369 final double previousX = currentX; 370 final double previousY = currentY; 371 currentX = observations[i].getX(); 372 currentY = observations[i].getY(); 373 final double currentYPrime = (currentY - previousY) / (currentX - previousX); 374 375 double omegaX = omega * currentX; 376 double cosine = FastMath.cos(omegaX); 377 double sine = FastMath.sin(omegaX); 378 fcMean += omega * currentY * cosine - currentYPrime * sine; 379 fsMean += omega * currentY * sine + currentYPrime * cosine; 380 } 381 382 return FastMath.atan2(-fsMean, fcMean); 383 } 384 } 385 }