You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2013/09/07 22:37:49 UTC
svn commit: r1520807 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math3/fitting/
test/java/org/apache/commons/math3/fitting/
Author: erans
Date: Sat Sep 7 20:37:49 2013
New Revision: 1520807
URL: http://svn.apache.org/r1520807
Log:
MATH-1014
Created "HarmonicCurveFitter" as a replacement for "HarmonicFitter".
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java (with props)
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java?rev=1520807&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java Sat Sep 7 20:37:49 2013
@@ -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.
+ * <br/>
+ * The {@link #withStartPoint(double[]) initial guess values} must be passed
+ * in the following order:
+ * <ul>
+ * <li>Amplitude</li>
+ * <li>Angular frequency</li>
+ * <li>phase</li>
+ * </ul>
+ * The optimal values will be returned in the same order.
+ *
+ * @version $Id$
+ * @since 3.3
+ */
+public class HarmonicCurveFitter extends AbstractCurveFitter<LevenbergMarquardtOptimizer>
+ implements WithStartPoint<HarmonicCurveFitter>,
+ WithMaxIterations<HarmonicCurveFitter> {
+ /** 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<WeightedObservedPoint> 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.
+ * <p>The algorithm used to guess the coefficients is as follows:</p>
+ *
+ * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
+ * ω and φ such that f (t) = a cos (ω t + φ).
+ * </p>
+ *
+ * <p>From the analytical expression, we can compute two primitives :
+ * <pre>
+ * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2
+ * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2
+ * where S (t) = sin (2 (ω t + φ)) / (2 ω)
+ * </pre>
+ * </p>
+ *
+ * <p>We can remove S between these expressions :
+ * <pre>
+ * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t)
+ * </pre>
+ * </p>
+ *
+ * <p>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)
+ * </p>
+ *
+ * <p>From the primitive, we can deduce the same form for definite
+ * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
+ * <pre>
+ * 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>))
+ * </pre>
+ * </p>
+ *
+ * <p>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.
+ * </p>
+ *
+ * <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
+ * coefficients A and B that minimize a least square criterion
+ * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
+ * <pre>
+ *
+ * ∑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>
+ * A = ------------------------
+ * ∑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>
+ *
+ * ∑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>
+ * B = ------------------------
+ * ∑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>
+ * </pre>
+ * </p>
+ *
+ *
+ * <p>In fact, we can assume both a and ω are positive and
+ * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that
+ * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p>
+ * <pre>
+ *
+ * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
+ * f (t<sub>i</sub>)
+ * 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>)
+ * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
+ * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
+ * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
+ * 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>
+ * end for
+ *
+ * |--------------------------
+ * \ | ∑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>
+ * a = \ | ------------------------
+ * \| ∑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>
+ *
+ *
+ * |--------------------------
+ * \ | ∑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>
+ * ω = \ | ------------------------
+ * \| ∑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>
+ *
+ * </pre>
+ * </p>
+ *
+ * <p>Once we know ω, we can compute:
+ * <pre>
+ * fc = ω f (t) cos (ω t) - f' (t) sin (ω t)
+ * fs = ω f (t) sin (ω t) + f' (t) cos (ω t)
+ * </pre>
+ * </p>
+ *
+ * <p>It appears that <code>fc = a ω cos (φ)</code> and
+ * <code>fs = -a ω sin (φ)</code>, so we can use these
+ * expressions to compute φ. The best estimate over the sample is
+ * given by averaging these expressions.
+ * </p>
+ *
+ * <p>Since integrals and means are involved in the preceding
+ * estimations, these operations run in O(n) time, where n is the
+ * number of measurements.</p>
+ */
+ 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<WeightedObservedPoint> 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:
+ * <ul>
+ * <li>Amplitude</li>
+ * <li>Angular frequency</li>
+ * <li>Phase</li>
+ * </ul>
+ */
+ 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<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) {
+ final List<WeightedObservedPoint> observations = new ArrayList<WeightedObservedPoint>(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 f<sup>2</sup> and f'<sup>2</sup>
+ // 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);
+ }
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java
------------------------------------------------------------------------------
svn:keywords = Id Revision
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java?rev=1520807&r1=1520806&r2=1520807&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java Sat Sep 7 20:37:49 2013
@@ -35,7 +35,10 @@ import org.apache.commons.math3.util.Fas
*
* @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<HarmonicOscillator.Parametric> {
/**
* Simple constructor.
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java?rev=1520807&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java Sat Sep 7 20:37:49 2013
@@ -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<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>();
+ 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);
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/HarmonicCurveFitterTest.java
------------------------------------------------------------------------------
svn:keywords = Id Revision