From f2b1cf5e5dc5e432bdafdeba25f532f7636caec9 Mon Sep 17 00:00:00 2001 From: Gilles Sadowski Date: Sat, 7 Sep 2013 20:37:49 +0000 Subject: [PATCH] MATH-1014 Created "HarmonicCurveFitter" as a replacement for "HarmonicFitter". git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1520807 13f79535-47bb-0310-9956-ffa450edef68 --- .../math3/fitting/HarmonicCurveFitter.java | 439 ++++++++++++++++++ .../commons/math3/fitting/HarmonicFitter.java | 3 + .../fitting/HarmonicCurveFitterTest.java | 183 ++++++++ 3 files changed, 625 insertions(+) create mode 100644 src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java diff --git a/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java b/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java new file mode 100644 index 000000000..2ff63fdbb --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java @@ -0,0 +1,439 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting; + +import java.util.Collection; +import java.util.List; +import java.util.ArrayList; +import org.apache.commons.math3.analysis.function.HarmonicOscillator; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer; +import org.apache.commons.math3.fitting.leastsquares.WithStartPoint; +import org.apache.commons.math3.fitting.leastsquares.WithMaxIterations; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.util.FastMath; + +/** + * Fits points to a {@link + * org.apache.commons.math3.analysis.function.HarmonicOscillator.Parametric harmonic oscillator} + * function. + *
+ * The {@link #withStartPoint(double[]) initial guess values} must be passed + * in the following order: + * + * The optimal values will be returned in the same order. + * + * @version $Id$ + * @since 3.3 + */ +public class HarmonicCurveFitter extends AbstractCurveFitter + implements WithStartPoint, + WithMaxIterations { + /** Parametric function to be fitted. */ + private static final HarmonicOscillator.Parametric FUNCTION = new HarmonicOscillator.Parametric(); + /** Initial guess. */ + private final double[] initialGuess; + /** Maximum number of iterations of the optimization algorithm. */ + private final int maxIter; + + /** + * Contructor used by the factory methods. + * + * @param initialGuess Initial guess. If set to {@code null}, the initial guess + * will be estimated using the {@link ParameterGuesser}. + * @param maxIter Maximum number of iterations of the optimization algorithm. + */ + private HarmonicCurveFitter(double[] initialGuess, + int maxIter) { + this.initialGuess = initialGuess; + this.maxIter = maxIter; + } + + /** + * Creates a default curve fitter. + * The initial guess for the parameters will be {@link ParameterGuesser} + * computed automatically, and the maximum number of iterations of the + * optimization algorithm is set to {@link Integer#MAX_VALUE}. + * + * @return a curve fitter. + * + * @see #withStartPoint(double[]) + * @see #withMaxIterations(int) + */ + public static HarmonicCurveFitter create() { + return new HarmonicCurveFitter(null, Integer.MAX_VALUE); + } + + /** {@inheritDoc} */ + public HarmonicCurveFitter withStartPoint(double[] start) { + return new HarmonicCurveFitter(start.clone(), + maxIter); + } + + /** {@inheritDoc} */ + public HarmonicCurveFitter withMaxIterations(int max) { + return new HarmonicCurveFitter(initialGuess, + max); + } + + /** {@inheritDoc} */ + @Override + protected LevenbergMarquardtOptimizer getOptimizer(Collection observations) { + // Prepare least-squares problem. + final int len = observations.size(); + final double[] target = new double[len]; + final double[] weights = new double[len]; + + int i = 0; + for (WeightedObservedPoint obs : observations) { + target[i] = obs.getY(); + weights[i] = obs.getWeight(); + ++i; + } + + final AbstractCurveFitter.TheoreticalValuesFunction model + = new AbstractCurveFitter.TheoreticalValuesFunction(FUNCTION, + observations); + + final double[] startPoint = initialGuess != null ? + initialGuess : + // Compute estimation. + new ParameterGuesser(observations).guess(); + + // Return a new optimizer set up to fit a Gaussian curve to the + // observed points. + return LevenbergMarquardtOptimizer.create() + .withMaxEvaluations(Integer.MAX_VALUE) + .withMaxIterations(maxIter) + .withStartPoint(startPoint) + .withTarget(target) + .withWeight(new DiagonalMatrix(weights)) + .withModelAndJacobian(model.getModelFunction(), + model.getModelFunctionJacobian()); + } + + /** + * This class guesses harmonic coefficients from a sample. + *

The algorithm used to guess the coefficients is as follows:

+ * + *

We know f (t) at some sampling points ti and want to find a, + * ω and φ such that f (t) = a cos (ω t + φ). + *

+ * + *

From the analytical expression, we can compute two primitives : + * 

+     *     If2  (t) = ∫ f2  = a2 × [t + S (t)] / 2
+     *     If'2 (t) = ∫ f'2 = a2 ω2 × [t - S (t)] / 2
+     *     where S (t) = sin (2 (ω t + φ)) / (2 ω)
+     *
+ *

+ * + *

We can remove S between these expressions : + * 

+     *     If'2 (t) = a2 ω2 t - ω2 If2 (t)
+     *
+ *

+ * + *

The preceding expression shows that If'2 (t) is a linear + * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t) + *

+ * + *

From the primitive, we can deduce the same form for definite + * integrals between t1 and ti for each ti : + * 

+     *   If2 (ti) - If2 (t1) = A × (ti - t1) + B × (If2 (ti) - If2 (t1))
+     *
+ *

+ * + *

We can find the coefficients A and B that best fit the sample + * to this linear expression by computing the definite integrals for + * each sample points. + *

+ * + *

For a bilinear expression z (xi, yi) = A × xi + B × yi, the + * coefficients A and B that minimize a least square criterion + * ∑ (zi - z (xi, yi))2 are given by these expressions:

+ * 
+     *
+     *         ∑yiyi ∑xizi - ∑xiyi ∑yizi
+     *     A = ------------------------
+     *         ∑xixi ∑yiyi - ∑xiyi ∑xiyi
+     *
+     *         ∑xixi ∑yizi - ∑xiyi ∑xizi
+     *     B = ------------------------
+     *         ∑xixi ∑yiyi - ∑xiyi ∑xiyi
+     *
+ *

+ * + * + *

In fact, we can assume both a and ω are positive and + * compute them directly, knowing that A = a2 ω2 and that + * B = - ω2. The complete algorithm is therefore:

+ * 
+     *
+     * for each ti from t1 to tn-1, compute:
+     *   f  (ti)
+     *   f' (ti) = (f (ti+1) - f(ti-1)) / (ti+1 - ti-1)
+     *   xi = ti - t1
+     *   yi = ∫ f2 from t1 to ti
+     *   zi = ∫ f'2 from t1 to ti
+     *   update the sums ∑xixi, ∑yiyi, ∑xiyi, ∑xizi and ∑yizi
+     * end for
+     *
+     *            |--------------------------
+     *         \  | ∑yiyi ∑xizi - ∑xiyi ∑yizi
+     * a     =  \ | ------------------------
+     *           \| ∑xiyi ∑xizi - ∑xixi ∑yizi
+     *
+     *
+     *            |--------------------------
+     *         \  | ∑xiyi ∑xizi - ∑xixi ∑yizi
+     * ω     =  \ | ------------------------
+     *           \| ∑xixi ∑yiyi - ∑xiyi ∑xiyi
+     *
+     *
+ *

+ * + *

Once we know ω, we can compute: + * 

+     *    fc = ω f (t) cos (ω t) - f' (t) sin (ω t)
+     *    fs = ω f (t) sin (ω t) + f' (t) cos (ω t)
+     *
+ *

+ * + *

It appears that fc = a ω cos (φ) and + * fs = -a ω sin (φ), so we can use these + * expressions to compute φ. The best estimate over the sample is + * given by averaging these expressions. + *

+ * + *

Since integrals and means are involved in the preceding + * estimations, these operations run in O(n) time, where n is the + * number of measurements.

+ */ + public static class ParameterGuesser { + /** Amplitude. */ + private final double a; + /** Angular frequency. */ + private final double omega; + /** Phase. */ + private final double phi; + + /** + * Simple constructor. + * + * @param observations Sampled observations. + * @throws NumberIsTooSmallException if the sample is too short. + * @throws ZeroException if the abscissa range is zero. + * @throws MathIllegalStateException when the guessing procedure cannot + * produce sensible results. + */ + public ParameterGuesser(Collection observations) { + if (observations.size() < 4) { + throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, + observations.size(), 4, true); + } + + final WeightedObservedPoint[] sorted + = sortObservations(observations).toArray(new WeightedObservedPoint[0]); + + final double aOmega[] = guessAOmega(sorted); + a = aOmega[0]; + omega = aOmega[1]; + + phi = guessPhi(sorted); + } + + /** + * Gets an estimation of the parameters. + * + * @return the guessed parameters, in the following order: + *
    + *
  • Amplitude
    • + *
  • Angular frequency
    • + *
  • Phase
    • + *
+ */ + public double[] guess() { + return new double[] { a, omega, phi }; + } + + /** + * Sort the observations with respect to the abscissa. + * + * @param unsorted Input observations. + * @return the input observations, sorted. + */ + private List sortObservations(Collection unsorted) { + final List observations = new ArrayList(unsorted); + + // Since the samples are almost always already sorted, this + // method is implemented as an insertion sort that reorders the + // elements in place. Insertion sort is very efficient in this case. + WeightedObservedPoint curr = observations.get(0); + final int len = observations.size(); + for (int j = 1; j < len; j++) { + WeightedObservedPoint prec = curr; + curr = observations.get(j); + if (curr.getX() < prec.getX()) { + // the current element should be inserted closer to the beginning + int i = j - 1; + WeightedObservedPoint mI = observations.get(i); + while ((i >= 0) && (curr.getX() < mI.getX())) { + observations.set(i + 1, mI); + if (i-- != 0) { + mI = observations.get(i); + } + } + observations.set(i + 1, curr); + curr = observations.get(j); + } + } + + return observations; + } + + /** + * Estimate a first guess of the amplitude and angular frequency. + * + * @param observations Observations, sorted w.r.t. abscissa. + * @throws ZeroException if the abscissa range is zero. + * @throws MathIllegalStateException when the guessing procedure cannot + * produce sensible results. + * @return the guessed amplitude (at index 0) and circular frequency + * (at index 1). + */ + private double[] guessAOmega(WeightedObservedPoint[] observations) { + final double[] aOmega = new double[2]; + + // initialize the sums for the linear model between the two integrals + double sx2 = 0; + double sy2 = 0; + double sxy = 0; + double sxz = 0; + double syz = 0; + + double currentX = observations[0].getX(); + double currentY = observations[0].getY(); + double f2Integral = 0; + double fPrime2Integral = 0; + final double startX = currentX; + for (int i = 1; i < observations.length; ++i) { + // one step forward + final double previousX = currentX; + final double previousY = currentY; + currentX = observations[i].getX(); + currentY = observations[i].getY(); + + // update the integrals of f2 and f'2 + // considering a linear model for f (and therefore constant f') + final double dx = currentX - previousX; + final double dy = currentY - previousY; + final double f2StepIntegral = + dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3; + final double fPrime2StepIntegral = dy * dy / dx; + + final double x = currentX - startX; + f2Integral += f2StepIntegral; + fPrime2Integral += fPrime2StepIntegral; + + sx2 += x * x; + sy2 += f2Integral * f2Integral; + sxy += x * f2Integral; + sxz += x * fPrime2Integral; + syz += f2Integral * fPrime2Integral; + } + + // compute the amplitude and pulsation coefficients + double c1 = sy2 * sxz - sxy * syz; + double c2 = sxy * sxz - sx2 * syz; + double c3 = sx2 * sy2 - sxy * sxy; + if ((c1 / c2 < 0) || (c2 / c3 < 0)) { + final int last = observations.length - 1; + // Range of the observations, assuming that the + // observations are sorted. + final double xRange = observations[last].getX() - observations[0].getX(); + if (xRange == 0) { + throw new ZeroException(); + } + aOmega[1] = 2 * Math.PI / xRange; + + double yMin = Double.POSITIVE_INFINITY; + double yMax = Double.NEGATIVE_INFINITY; + for (int i = 1; i < observations.length; ++i) { + final double y = observations[i].getY(); + if (y < yMin) { + yMin = y; + } + if (y > yMax) { + yMax = y; + } + } + aOmega[0] = 0.5 * (yMax - yMin); + } else { + if (c2 == 0) { + // In some ill-conditioned cases (cf. MATH-844), the guesser + // procedure cannot produce sensible results. + throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR); + } + + aOmega[0] = FastMath.sqrt(c1 / c2); + aOmega[1] = FastMath.sqrt(c2 / c3); + } + + return aOmega; + } + + /** + * Estimate a first guess of the phase. + * + * @param observations Observations, sorted w.r.t. abscissa. + * @return the guessed phase. + */ + private double guessPhi(WeightedObservedPoint[] observations) { + // initialize the means + double fcMean = 0; + double fsMean = 0; + + double currentX = observations[0].getX(); + double currentY = observations[0].getY(); + for (int i = 1; i < observations.length; ++i) { + // one step forward + final double previousX = currentX; + final double previousY = currentY; + currentX = observations[i].getX(); + currentY = observations[i].getY(); + final double currentYPrime = (currentY - previousY) / (currentX - previousX); + + double omegaX = omega * currentX; + double cosine = FastMath.cos(omegaX); + double sine = FastMath.sin(omegaX); + fcMean += omega * currentY * cosine - currentYPrime * sine; + fsMean += omega * currentY * sine + currentYPrime * cosine; + } + + return FastMath.atan2(-fsMean, fcMean); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java b/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java index 12badd2b1..94183eef7 100644 --- a/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java +++ b/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java @@ -35,7 +35,10 @@ import org.apache.commons.math3.util.FastMath; * * @version $Id: HarmonicFitter.java 1416643 2012-12-03 19:37:14Z tn $ * @since 2.0 + * @deprecated As of 3.3. Please use {@link HarmonicCurveFitter} and + * {@link WeightedObservedPoints} instead. */ +@Deprecated public class HarmonicFitter extends CurveFitter { /** * Simple constructor. diff --git a/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java b/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java new file mode 100644 index 000000000..9b797ea67 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java @@ -0,0 +1,183 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.fitting; + +import java.util.Random; +import java.util.List; +import java.util.ArrayList; +import org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer; +import org.apache.commons.math3.analysis.function.HarmonicOscillator; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +import org.junit.Test; +import org.junit.Assert; + +public class HarmonicCurveFitterTest { + /** + * Zero points is not enough observed points. + */ + @Test(expected=NumberIsTooSmallException.class) + public void testPreconditions1() { + HarmonicCurveFitter.create().fit(new WeightedObservedPoints().toList()); + } + + @Test + public void testNoError() { + final double a = 0.2; + final double w = 3.4; + final double p = 4.1; + final HarmonicOscillator f = new HarmonicOscillator(a, w, p); + + final WeightedObservedPoints points = new WeightedObservedPoints(); + for (double x = 0.0; x < 1.3; x += 0.01) { + points.add(1, x, f.value(x)); + } + + final HarmonicCurveFitter fitter = HarmonicCurveFitter.create(); + final double[] fitted = fitter.fit(points.toList()); + Assert.assertEquals(a, fitted[0], 1.0e-13); + Assert.assertEquals(w, fitted[1], 1.0e-13); + Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1e-13); + + final HarmonicOscillator ff = new HarmonicOscillator(fitted[0], fitted[1], fitted[2]); + for (double x = -1.0; x < 1.0; x += 0.01) { + Assert.assertTrue(FastMath.abs(f.value(x) - ff.value(x)) < 1e-13); + } + } + + @Test + public void test1PercentError() { + final Random randomizer = new Random(64925784252L); + final double a = 0.2; + final double w = 3.4; + final double p = 4.1; + final HarmonicOscillator f = new HarmonicOscillator(a, w, p); + + final WeightedObservedPoints points = new WeightedObservedPoints(); + for (double x = 0.0; x < 10.0; x += 0.1) { + points.add(1, x, f.value(x) + 0.01 * randomizer.nextGaussian()); + } + + final HarmonicCurveFitter fitter = HarmonicCurveFitter.create(); + final double[] fitted = fitter.fit(points.toList()); + Assert.assertEquals(a, fitted[0], 7.6e-4); + Assert.assertEquals(w, fitted[1], 2.7e-3); + Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.3e-2); + } + + @Test + public void testTinyVariationsData() { + final Random randomizer = new Random(64925784252L); + + final WeightedObservedPoints points = new WeightedObservedPoints(); + for (double x = 0.0; x < 10.0; x += 0.1) { + points.add(1, x, 1e-7 * randomizer.nextGaussian()); + } + + final HarmonicCurveFitter fitter = HarmonicCurveFitter.create(); + final double[] fitted = fitter.fit(points.toList()); + + // This test serves to cover the part of the code of "guessAOmega" + // when the algorithm using integrals fails. + } + + @Test + public void testInitialGuess() { + final Random randomizer = new Random(45314242L); + final double a = 0.2; + final double w = 3.4; + final double p = 4.1; + final HarmonicOscillator f = new HarmonicOscillator(a, w, p); + + final WeightedObservedPoints points = new WeightedObservedPoints(); + for (double x = 0.0; x < 10.0; x += 0.1) { + points.add(1, x, f.value(x) + 0.01 * randomizer.nextGaussian()); + } + + final HarmonicCurveFitter fitter = HarmonicCurveFitter.create() + .withStartPoint(new double[] { 0.15, 3.6, 4.5 }); + final double[] fitted = fitter.fit(points.toList()); + Assert.assertEquals(a, fitted[0], 1.2e-3); + Assert.assertEquals(w, fitted[1], 3.3e-3); + Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.7e-2); + } + + @Test + public void testUnsorted() { + Random randomizer = new Random(64925784252L); + final double a = 0.2; + final double w = 3.4; + final double p = 4.1; + final HarmonicOscillator f = new HarmonicOscillator(a, w, p); + + // Build a regularly spaced array of measurements. + final int size = 100; + final double[] xTab = new double[size]; + final double[] yTab = new double[size]; + for (int i = 0; i < size; i++) { + xTab[i] = 0.1 * i; + yTab[i] = f.value(xTab[i]) + 0.01 * randomizer.nextGaussian(); + } + + // shake it + for (int i = 0; i < size; i++) { + int i1 = randomizer.nextInt(size); + int i2 = randomizer.nextInt(size); + double xTmp = xTab[i1]; + double yTmp = yTab[i1]; + xTab[i1] = xTab[i2]; + yTab[i1] = yTab[i2]; + xTab[i2] = xTmp; + yTab[i2] = yTmp; + } + + // Pass it to the fitter. + final WeightedObservedPoints points = new WeightedObservedPoints(); + for (int i = 0; i < size; ++i) { + points.add(1, xTab[i], yTab[i]); + } + + final HarmonicCurveFitter fitter = HarmonicCurveFitter.create(); + final double[] fitted = fitter.fit(points.toList()); + Assert.assertEquals(a, fitted[0], 7.6e-4); + Assert.assertEquals(w, fitted[1], 3.5e-3); + Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.5e-2); + } + + @Test(expected=MathIllegalStateException.class) + public void testMath844() { + final double[] y = { 0, 1, 2, 3, 2, 1, + 0, -1, -2, -3, -2, -1, + 0, 1, 2, 3, 2, 1, + 0, -1, -2, -3, -2, -1, + 0, 1, 2, 3, 2, 1, 0 }; + final List points = new ArrayList(); + for (int i = 0; i < y.length; i++) { + points.add(new WeightedObservedPoint(1, i, y[i])); + } + + // The guesser fails because the function is far from an harmonic + // function: It is a triangular periodic function with amplitude 3 + // and period 12, and all sample points are taken at integer abscissae + // so function values all belong to the integer subset {-3, -2, -1, 0, + // 1, 2, 3}. + new HarmonicCurveFitter.ParameterGuesser(points); + } +}