http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.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 deleted file mode 100644 index e5e7be4..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/HarmonicCurveFitter.java +++ /dev/null @@ -1,445 +0,0 @@ -/* - * 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.ArrayList; -import java.util.Collection; -import java.util.List; - -import org.apache.commons.math3.analysis.function.HarmonicOscillator; -import org.apache.commons.math3.exception.MathIllegalStateException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.ZeroException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; -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. - * - * @since 3.3 - */ -public class HarmonicCurveFitter extends AbstractCurveFitter { - /** 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); - } - - /** - * Configure the start point (initial guess). - * @param newStart new start point (initial guess) - * @return a new instance. - */ - public HarmonicCurveFitter withStartPoint(double[] newStart) { - return new HarmonicCurveFitter(newStart.clone(), - maxIter); - } - - /** - * Configure the maximum number of iterations. - * @param newMaxIter maximum number of iterations - * @return a new instance. - */ - public HarmonicCurveFitter withMaxIterations(int newMaxIter) { - return new HarmonicCurveFitter(initialGuess, - newMaxIter); - } - - /** {@inheritDoc} */ - @Override - protected LeastSquaresProblem getProblem(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 new LeastSquaresBuilder(). - maxEvaluations(Integer.MAX_VALUE). - maxIterations(maxIter). - start(startPoint). - target(target). - weight(new DiagonalMatrix(weights)). - model(model.getModelFunction(), model.getModelFunctionJacobian()). - build(); - - } - - /** - * 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_i \) and want - * to find \( a \), \( \omega \) and \( \phi \) such that - * \( f(t) = a \cos (\omega t + \phi) \). - * </p> - * - * <p>From the analytical expression, we can compute two primitives : - * \[ - * If2(t) = \int f^2 dt = a^2 (t + S(t)) / 2 - * \] - * \[ - * If'2(t) = \int f'^2 dt = a^2 \omega^2 (t - S(t)) / 2 - * \] - * where \(S(t) = \frac{\sin(2 (\omega t + \phi))}{2\omega}\) - * </p> - * - * <p>We can remove \(S\) between these expressions : - * \[ - * If'2(t) = a^2 \omega^2 t - \omega^2 If2(t) - * \] - * </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_1\) and \(t_i\) for each \(t_i\) : - * \[ - * If2(t_i) - If2(t_1) = A (t_i - t_1) + B (If2 (t_i) - If2(t_1)) - * \] - * </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_i, y_i) = A x_i + B y_i\), the - * coefficients \(A\) and \(B\) that minimize a least-squares criterion - * \(\sum (z_i - z(x_i, y_i))^2\) are given by these expressions:</p> - * \[ - * A = \frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i} - * {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i} - * \] - * \[ - * B = \frac{\sum x_i x_i \sum y_i z_i - \sum x_i y_i \sum x_i z_i} - * {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i} - * - * \] - * - * <p>In fact, we can assume that both \(a\) and \(\omega\) are positive and - * compute them directly, knowing that \(A = a^2 \omega^2\) and that - * \(B = -\omega^2\). The complete algorithm is therefore:</p> - * - * For each \(t_i\) from \(t_1\) to \(t_{n-1}\), compute: - * \[ f(t_i) \] - * \[ f'(t_i) = \frac{f (t_{i+1}) - f(t_{i-1})}{t_{i+1} - t_{i-1}} \] - * \[ x_i = t_i - t_1 \] - * \[ y_i = \int_{t_1}^{t_i} f^2(t) dt \] - * \[ z_i = \int_{t_1}^{t_i} f'^2(t) dt \] - * and update the sums: - * \[ \sum x_i x_i, \sum y_i y_i, \sum x_i y_i, \sum x_i z_i, \sum y_i z_i \] - * - * Then: - * \[ - * a = \sqrt{\frac{\sum y_i y_i \sum x_i z_i - \sum x_i y_i \sum y_i z_i } - * {\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i }} - * \] - * \[ - * \omega = \sqrt{\frac{\sum x_i y_i \sum x_i z_i - \sum x_i x_i \sum y_i z_i} - * {\sum x_i x_i \sum y_i y_i - \sum x_i y_i \sum x_i y_i}} - * \] - * - * <p>Once we know \(\omega\) we can compute: - * \[ - * fc = \omega f(t) \cos(\omega t) - f'(t) \sin(\omega t) - * \] - * \[ - * fs = \omega f(t) \sin(\omega t) + f'(t) \cos(\omega t) - * \] - * </p> - * - * <p>It appears that \(fc = a \omega \cos(\phi)\) and - * \(fs = -a \omega \sin(\phi)\), so we can use these - * expressions to compute \(\phi\). 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); - } - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java b/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java deleted file mode 100644 index 24952a6..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java +++ /dev/null @@ -1,384 +0,0 @@ -/* - * 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 org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer; -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.util.FastMath; - -/** - * Class that implements a curve fitting specialized for sinusoids. - * - * Harmonic fitting is a very simple case of curve fitting. The - * estimated coefficients are the amplitude a, the pulsation ω and - * the phase φ: <code>f (t) = a cos (ω t + φ)</code>. They are - * searched by a least square estimator initialized with a rough guess - * based on integrals. - * - * @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. - * @param optimizer Optimizer to use for the fitting. - */ - public HarmonicFitter(final MultivariateVectorOptimizer optimizer) { - super(optimizer); - } - - /** - * Fit an harmonic function to the observed points. - * - * @param initialGuess First guess values in the following order: - * <ul> - * <li>Amplitude</li> - * <li>Angular frequency</li> - * <li>Phase</li> - * </ul> - * @return the parameters of the harmonic function that best fits the - * observed points (in the same order as above). - */ - public double[] fit(double[] initialGuess) { - return fit(new HarmonicOscillator.Parametric(), initialGuess); - } - - /** - * Fit an harmonic function to the observed points. - * An initial guess will be automatically computed. - * - * @return the parameters of the harmonic function that best fits the - * observed points (see the other {@link #fit(double[]) fit} method. - * @throws NumberIsTooSmallException if the sample is too short for the - * the first guess to be computed. - * @throws ZeroException if the first guess cannot be computed because - * the abscissa range is zero. - */ - public double[] fit() { - return fit((new ParameterGuesser(getObservations())).guess()); - } - - /** - * 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(WeightedObservedPoint[] observations) { - if (observations.length < 4) { - throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, - observations.length, 4, true); - } - - final WeightedObservedPoint[] sorted = sortObservations(observations); - - 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 WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) { - final WeightedObservedPoint[] observations = unsorted.clone(); - - // 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[0]; - for (int j = 1; j < observations.length; ++j) { - WeightedObservedPoint prec = curr; - curr = observations[j]; - if (curr.getX() < prec.getX()) { - // the current element should be inserted closer to the beginning - int i = j - 1; - WeightedObservedPoint mI = observations[i]; - while ((i >= 0) && (curr.getX() < mI.getX())) { - observations[i + 1] = mI; - if (i-- != 0) { - mI = observations[i]; - } - } - observations[i + 1] = curr; - curr = observations[j]; - } - } - - return observations; - } - - /** - * Estimate a first guess of the amplitude and angular frequency. - * This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method - * has been called previously. - * - * @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); - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/PolynomialCurveFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/PolynomialCurveFitter.java b/src/main/java/org/apache/commons/math3/fitting/PolynomialCurveFitter.java deleted file mode 100644 index 9a6f3a9..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/PolynomialCurveFitter.java +++ /dev/null @@ -1,131 +0,0 @@ -/* - * 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 org.apache.commons.math3.analysis.polynomials.PolynomialFunction; -import org.apache.commons.math3.exception.MathInternalError; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; -import org.apache.commons.math3.linear.DiagonalMatrix; - -/** - * Fits points to a {@link - * org.apache.commons.math3.analysis.polynomials.PolynomialFunction.Parametric polynomial} - * function. - * <br/> - * The size of the {@link #withStartPoint(double[]) initial guess} array defines the - * degree of the polynomial to be fitted. - * They must be sorted in increasing order of the polynomial's degree. - * The optimal values of the coefficients will be returned in the same order. - * - * @since 3.3 - */ -public class PolynomialCurveFitter extends AbstractCurveFitter { - /** Parametric function to be fitted. */ - private static final PolynomialFunction.Parametric FUNCTION = new PolynomialFunction.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. - * @param maxIter Maximum number of iterations of the optimization algorithm. - * @throws MathInternalError if {@code initialGuess} is {@code null}. - */ - private PolynomialCurveFitter(double[] initialGuess, - int maxIter) { - this.initialGuess = initialGuess; - this.maxIter = maxIter; - } - - /** - * Creates a default curve fitter. - * Zero will be used as initial guess for the coefficients, and the maximum - * number of iterations of the optimization algorithm is set to - * {@link Integer#MAX_VALUE}. - * - * @param degree Degree of the polynomial to be fitted. - * @return a curve fitter. - * - * @see #withStartPoint(double[]) - * @see #withMaxIterations(int) - */ - public static PolynomialCurveFitter create(int degree) { - return new PolynomialCurveFitter(new double[degree + 1], Integer.MAX_VALUE); - } - - /** - * Configure the start point (initial guess). - * @param newStart new start point (initial guess) - * @return a new instance. - */ - public PolynomialCurveFitter withStartPoint(double[] newStart) { - return new PolynomialCurveFitter(newStart.clone(), - maxIter); - } - - /** - * Configure the maximum number of iterations. - * @param newMaxIter maximum number of iterations - * @return a new instance. - */ - public PolynomialCurveFitter withMaxIterations(int newMaxIter) { - return new PolynomialCurveFitter(initialGuess, - newMaxIter); - } - - /** {@inheritDoc} */ - @Override - protected LeastSquaresProblem getProblem(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); - - if (initialGuess == null) { - throw new MathInternalError(); - } - - // Return a new least squares problem set up to fit a polynomial curve to the - // observed points. - return new LeastSquaresBuilder(). - maxEvaluations(Integer.MAX_VALUE). - maxIterations(maxIter). - start(initialGuess). - target(target). - weight(new DiagonalMatrix(weights)). - model(model.getModelFunction(), model.getModelFunctionJacobian()). - build(); - - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/PolynomialFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/PolynomialFitter.java b/src/main/java/org/apache/commons/math3/fitting/PolynomialFitter.java deleted file mode 100644 index ec242c1..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/PolynomialFitter.java +++ /dev/null @@ -1,72 +0,0 @@ -/* - * 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 org.apache.commons.math3.analysis.polynomials.PolynomialFunction; -import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer; - -/** - * Polynomial fitting is a very simple case of {@link CurveFitter curve fitting}. - * The estimated coefficients are the polynomial coefficients (see the - * {@link #fit(double[]) fit} method). - * - * @since 2.0 - * @deprecated As of 3.3. Please use {@link PolynomialCurveFitter} and - * {@link WeightedObservedPoints} instead. - */ -@Deprecated -public class PolynomialFitter extends CurveFitter<PolynomialFunction.Parametric> { - /** - * Simple constructor. - * - * @param optimizer Optimizer to use for the fitting. - */ - public PolynomialFitter(MultivariateVectorOptimizer optimizer) { - super(optimizer); - } - - /** - * Get the coefficients of the polynomial fitting the weighted data points. - * The degree of the fitting polynomial is {@code guess.length - 1}. - * - * @param guess First guess for the coefficients. They must be sorted in - * increasing order of the polynomial's degree. - * @param maxEval Maximum number of evaluations of the polynomial. - * @return the coefficients of the polynomial that best fits the observed points. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if - * the number of evaluations exceeds {@code maxEval}. - * @throws org.apache.commons.math3.exception.ConvergenceException - * if the algorithm failed to converge. - */ - public double[] fit(int maxEval, double[] guess) { - return fit(maxEval, new PolynomialFunction.Parametric(), guess); - } - - /** - * Get the coefficients of the polynomial fitting the weighted data points. - * The degree of the fitting polynomial is {@code guess.length - 1}. - * - * @param guess First guess for the coefficients. They must be sorted in - * increasing order of the polynomial's degree. - * @return the coefficients of the polynomial that best fits the observed points. - * @throws org.apache.commons.math3.exception.ConvergenceException - * if the algorithm failed to converge. - */ - public double[] fit(double[] guess) { - return fit(new PolynomialFunction.Parametric(), guess); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/SimpleCurveFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/SimpleCurveFitter.java b/src/main/java/org/apache/commons/math3/fitting/SimpleCurveFitter.java deleted file mode 100644 index 3f6b7d6..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/SimpleCurveFitter.java +++ /dev/null @@ -1,125 +0,0 @@ -/* - * 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 org.apache.commons.math3.analysis.ParametricUnivariateFunction; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresBuilder; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; -import org.apache.commons.math3.linear.DiagonalMatrix; - -/** - * Fits points to a user-defined {@link ParametricUnivariateFunction function}. - * - * @since 3.4 - */ -public class SimpleCurveFitter extends AbstractCurveFitter { - /** Function to fit. */ - private final ParametricUnivariateFunction function; - /** Initial guess for the parameters. */ - private final double[] initialGuess; - /** Maximum number of iterations of the optimization algorithm. */ - private final int maxIter; - - /** - * Contructor used by the factory methods. - * - * @param function Function to fit. - * @param initialGuess Initial guess. Cannot be {@code null}. Its length must - * be consistent with the number of parameters of the {@code function} to fit. - * @param maxIter Maximum number of iterations of the optimization algorithm. - */ - private SimpleCurveFitter(ParametricUnivariateFunction function, - double[] initialGuess, - int maxIter) { - this.function = function; - this.initialGuess = initialGuess; - this.maxIter = maxIter; - } - - /** - * Creates a curve fitter. - * The maximum number of iterations of the optimization algorithm is set - * to {@link Integer#MAX_VALUE}. - * - * @param f Function to fit. - * @param start Initial guess for the parameters. Cannot be {@code null}. - * Its length must be consistent with the number of parameters of the - * function to fit. - * @return a curve fitter. - * - * @see #withStartPoint(double[]) - * @see #withMaxIterations(int) - */ - public static SimpleCurveFitter create(ParametricUnivariateFunction f, - double[] start) { - return new SimpleCurveFitter(f, start, Integer.MAX_VALUE); - } - - /** - * Configure the start point (initial guess). - * @param newStart new start point (initial guess) - * @return a new instance. - */ - public SimpleCurveFitter withStartPoint(double[] newStart) { - return new SimpleCurveFitter(function, - newStart.clone(), - maxIter); - } - - /** - * Configure the maximum number of iterations. - * @param newMaxIter maximum number of iterations - * @return a new instance. - */ - public SimpleCurveFitter withMaxIterations(int newMaxIter) { - return new SimpleCurveFitter(function, - initialGuess, - newMaxIter); - } - - /** {@inheritDoc} */ - @Override - protected LeastSquaresProblem getProblem(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 count = 0; - for (WeightedObservedPoint obs : observations) { - target[count] = obs.getY(); - weights[count] = obs.getWeight(); - ++count; - } - - final AbstractCurveFitter.TheoreticalValuesFunction model - = new AbstractCurveFitter.TheoreticalValuesFunction(function, - observations); - - // Create an optimizer for fitting the curve to the observed points. - return new LeastSquaresBuilder(). - maxEvaluations(Integer.MAX_VALUE). - maxIterations(maxIter). - start(initialGuess). - target(target). - weight(new DiagonalMatrix(weights)). - model(model.getModelFunction(), model.getModelFunctionJacobian()). - build(); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoint.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoint.java b/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoint.java deleted file mode 100644 index 00292c9..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoint.java +++ /dev/null @@ -1,78 +0,0 @@ -/* - * 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.io.Serializable; - -/** - * This class is a simple container for weighted observed point in - * {@link CurveFitter curve fitting}. - * <p>Instances of this class are guaranteed to be immutable.</p> - * @since 2.0 - */ -public class WeightedObservedPoint implements Serializable { - /** Serializable version id. */ - private static final long serialVersionUID = 5306874947404636157L; - /** Weight of the measurement in the fitting process. */ - private final double weight; - /** Abscissa of the point. */ - private final double x; - /** Observed value of the function at x. */ - private final double y; - - /** - * Simple constructor. - * - * @param weight Weight of the measurement in the fitting process. - * @param x Abscissa of the measurement. - * @param y Ordinate of the measurement. - */ - public WeightedObservedPoint(final double weight, final double x, final double y) { - this.weight = weight; - this.x = x; - this.y = y; - } - - /** - * Gets the weight of the measurement in the fitting process. - * - * @return the weight of the measurement in the fitting process. - */ - public double getWeight() { - return weight; - } - - /** - * Gets the abscissa of the point. - * - * @return the abscissa of the point. - */ - public double getX() { - return x; - } - - /** - * Gets the observed value of the function at x. - * - * @return the observed value of the function at x. - */ - public double getY() { - return y; - } - -} - http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoints.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoints.java b/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoints.java deleted file mode 100644 index 9eba337..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/WeightedObservedPoints.java +++ /dev/null @@ -1,112 +0,0 @@ -/* - * 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.List; -import java.util.ArrayList; -import java.io.Serializable; - -/** - * Simple container for weighted observed points used - * in {@link AbstractCurveFitter curve fitting} algorithms. - * - * @since 3.3 - */ -public class WeightedObservedPoints implements Serializable { - /** Serializable version id. */ - private static final long serialVersionUID = 20130813L; - - /** Observed points. */ - private final List<WeightedObservedPoint> observations - = new ArrayList<WeightedObservedPoint>(); - - /** - * Adds a point to the sample. - * Calling this method is equivalent to calling - * {@code add(1.0, x, y)}. - * - * @param x Abscissa of the point. - * @param y Observed value at {@code x}. After fitting we should - * have {@code f(x)} as close as possible to this value. - * - * @see #add(double, double, double) - * @see #add(WeightedObservedPoint) - * @see #toList() - */ - public void add(double x, double y) { - add(1d, x, y); - } - - /** - * Adds a point to the sample. - * - * @param weight Weight of the observed point. - * @param x Abscissa of the point. - * @param y Observed value at {@code x}. After fitting we should - * have {@code f(x)} as close as possible to this value. - * - * @see #add(double, double) - * @see #add(WeightedObservedPoint) - * @see #toList() - */ - public void add(double weight, double x, double y) { - observations.add(new WeightedObservedPoint(weight, x, y)); - } - - /** - * Adds a point to the sample. - * - * @param observed Observed point to add. - * - * @see #add(double, double) - * @see #add(double, double, double) - * @see #toList() - */ - public void add(WeightedObservedPoint observed) { - observations.add(observed); - } - - /** - * Gets a <em>snapshot</em> of the observed points. - * The list of stored points is copied in order to ensure that - * modification of the returned instance does not affect this - * container. - * Conversely, further modification of this container (through - * the {@code add} or {@code clear} methods) will not affect the - * returned list. - * - * @return the observed points, in the order they were added to this - * container. - * - * @see #add(double, double) - * @see #add(double, double, double) - * @see #add(WeightedObservedPoint) - */ - public List<WeightedObservedPoint> toList() { - // The copy is necessary to ensure thread-safety because of the - // "clear" method (which otherwise would be able to empty the - // list of points while it is being used by another thread). - return new ArrayList<WeightedObservedPoint>(observations); - } - - /** - * Removes all observations from this container. - */ - public void clear() { - observations.clear(); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java deleted file mode 100644 index b164380..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractEvaluation.java +++ /dev/null @@ -1,87 +0,0 @@ -/* - * 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.leastsquares; - -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.DecompositionSolver; -import org.apache.commons.math3.linear.QRDecomposition; -import org.apache.commons.math3.linear.RealMatrix; -import org.apache.commons.math3.linear.RealVector; -import org.apache.commons.math3.util.FastMath; - -/** - * An implementation of {@link Evaluation} that is designed for extension. All of the - * methods implemented here use the methods that are left unimplemented. - * <p/> - * TODO cache results? - * - * @since 3.3 - */ -public abstract class AbstractEvaluation implements Evaluation { - - /** number of observations */ - private final int observationSize; - - /** - * Constructor. - * - * @param observationSize the number of observation. Needed for {@link - * #getRMS()}. - */ - AbstractEvaluation(final int observationSize) { - this.observationSize = observationSize; - } - - /** {@inheritDoc} */ - public RealMatrix getCovariances(double threshold) { - // Set up the Jacobian. - final RealMatrix j = this.getJacobian(); - - // Compute transpose(J)J. - final RealMatrix jTj = j.transpose().multiply(j); - - // Compute the covariances matrix. - final DecompositionSolver solver - = new QRDecomposition(jTj, threshold).getSolver(); - return solver.getInverse(); - } - - /** {@inheritDoc} */ - public RealVector getSigma(double covarianceSingularityThreshold) { - final RealMatrix cov = this.getCovariances(covarianceSingularityThreshold); - final int nC = cov.getColumnDimension(); - final RealVector sig = new ArrayRealVector(nC); - for (int i = 0; i < nC; ++i) { - sig.setEntry(i, FastMath.sqrt(cov.getEntry(i,i))); - } - return sig; - } - - /** {@inheritDoc} */ - public double getRMS() { - final double cost = this.getCost(); - return FastMath.sqrt(cost * cost / this.observationSize); - } - - /** {@inheritDoc} */ - public double getCost() { - final ArrayRealVector r = new ArrayRealVector(this.getResiduals()); - return FastMath.sqrt(r.dotProduct(r)); - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java deleted file mode 100644 index 89f5f1f..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/DenseWeightedEvaluation.java +++ /dev/null @@ -1,68 +0,0 @@ -/* - * 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.leastsquares; - -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; -import org.apache.commons.math3.linear.RealMatrix; -import org.apache.commons.math3.linear.RealVector; - -/** - * Applies a dense weight matrix to an evaluation. - * - * @since 3.3 - */ -class DenseWeightedEvaluation extends AbstractEvaluation { - - /** the unweighted evaluation */ - private final Evaluation unweighted; - /** reference to the weight square root matrix */ - private final RealMatrix weightSqrt; - - /** - * Create a weighted evaluation from an unweighted one. - * - * @param unweighted the evalutation before weights are applied - * @param weightSqrt the matrix square root of the weight matrix - */ - DenseWeightedEvaluation(final Evaluation unweighted, - final RealMatrix weightSqrt) { - // weight square root is square, nR=nC=number of observations - super(weightSqrt.getColumnDimension()); - this.unweighted = unweighted; - this.weightSqrt = weightSqrt; - } - - /* apply weights */ - - /** {@inheritDoc} */ - public RealMatrix getJacobian() { - return weightSqrt.multiply(this.unweighted.getJacobian()); - } - - /** {@inheritDoc} */ - public RealVector getResiduals() { - return this.weightSqrt.operate(this.unweighted.getResiduals()); - } - - /* delegate */ - - /** {@inheritDoc} */ - public RealVector getPoint() { - return unweighted.getPoint(); - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java deleted file mode 100644 index ceb5988..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/EvaluationRmsChecker.java +++ /dev/null @@ -1,75 +0,0 @@ -/* - * 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.leastsquares; - -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; -import org.apache.commons.math3.optim.ConvergenceChecker; -import org.apache.commons.math3.util.Precision; - -/** - * Check if an optimization has converged based on the change in computed RMS. - * - * @since 3.4 - */ -public class EvaluationRmsChecker implements ConvergenceChecker<Evaluation> { - - /** relative tolerance for comparisons. */ - private final double relTol; - /** absolute tolerance for comparisons. */ - private final double absTol; - - /** - * Create a convergence checker for the RMS with the same relative and absolute - * tolerance. - * - * <p>Convenience constructor for when the relative and absolute tolerances are the - * same. Same as {@code new EvaluationRmsChecker(tol, tol)}. - * - * @param tol the relative and absolute tolerance. - * @see #EvaluationRmsChecker(double, double) - */ - public EvaluationRmsChecker(final double tol) { - this(tol, tol); - } - - /** - * Create a convergence checker for the RMS with a relative and absolute tolerance. - * - * <p>The optimization has converged when the RMS of consecutive evaluations are equal - * to within the given relative tolerance or absolute tolerance. - * - * @param relTol the relative tolerance. - * @param absTol the absolute tolerance. - * @see Precision#equals(double, double, double) - * @see Precision#equalsWithRelativeTolerance(double, double, double) - */ - public EvaluationRmsChecker(final double relTol, final double absTol) { - this.relTol = relTol; - this.absTol = absTol; - } - - /** {@inheritDoc} */ - public boolean converged(final int iteration, - final Evaluation previous, - final Evaluation current) { - final double prevRms = previous.getRMS(); - final double currRms = current.getRMS(); - return Precision.equals(prevRms, currRms, this.absTol) || - Precision.equalsWithRelativeTolerance(prevRms, currRms, this.relTol); - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java deleted file mode 100644 index bc4503c..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java +++ /dev/null @@ -1,300 +0,0 @@ -/* - * 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.leastsquares; - -import org.apache.commons.math3.exception.ConvergenceException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation; -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.CholeskyDecomposition; -import org.apache.commons.math3.linear.LUDecomposition; -import org.apache.commons.math3.linear.MatrixUtils; -import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException; -import org.apache.commons.math3.linear.QRDecomposition; -import org.apache.commons.math3.linear.RealMatrix; -import org.apache.commons.math3.linear.RealVector; -import org.apache.commons.math3.linear.SingularMatrixException; -import org.apache.commons.math3.linear.SingularValueDecomposition; -import org.apache.commons.math3.optim.ConvergenceChecker; -import org.apache.commons.math3.util.Incrementor; -import org.apache.commons.math3.util.Pair; - -/** - * Gauss-Newton least-squares solver. - * <p> This class solve a least-square problem by - * solving the normal equations of the linearized problem at each iteration. Either LU - * decomposition or Cholesky decomposition can be used to solve the normal equations, - * or QR decomposition or SVD decomposition can be used to solve the linear system. LU - * decomposition is faster but QR decomposition is more robust for difficult problems, - * and SVD can compute a solution for rank-deficient problems. - * </p> - * - * @since 3.3 - */ -public class GaussNewtonOptimizer implements LeastSquaresOptimizer { - - /** The decomposition algorithm to use to solve the normal equations. */ - //TODO move to linear package and expand options? - public static enum Decomposition { - /** - * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and - * using the {@link LUDecomposition}. - * - * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the - * normal matrix and n<sup>3</sup>/3 operations (m > n) to solve the system using - * the LU decomposition. </p> - */ - LU { - @Override - protected RealVector solve(final RealMatrix jacobian, - final RealVector residuals) { - try { - final Pair<RealMatrix, RealVector> normalEquation = - computeNormalMatrix(jacobian, residuals); - final RealMatrix normal = normalEquation.getFirst(); - final RealVector jTr = normalEquation.getSecond(); - return new LUDecomposition(normal, SINGULARITY_THRESHOLD) - .getSolver() - .solve(jTr); - } catch (SingularMatrixException e) { - throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); - } - } - }, - /** - * Solve the linear least squares problem (Jx=r) using the {@link - * QRDecomposition}. - * - * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations - * (m > n) and has better numerical accuracy than any method that forms the normal - * equations. </p> - */ - QR { - @Override - protected RealVector solve(final RealMatrix jacobian, - final RealVector residuals) { - try { - return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD) - .getSolver() - .solve(residuals); - } catch (SingularMatrixException e) { - throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); - } - } - }, - /** - * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and - * using the {@link CholeskyDecomposition}. - * - * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the - * normal matrix and n<sup>3</sup>/6 operations (m > n) to solve the system using - * the Cholesky decomposition. </p> - */ - CHOLESKY { - @Override - protected RealVector solve(final RealMatrix jacobian, - final RealVector residuals) { - try { - final Pair<RealMatrix, RealVector> normalEquation = - computeNormalMatrix(jacobian, residuals); - final RealMatrix normal = normalEquation.getFirst(); - final RealVector jTr = normalEquation.getSecond(); - return new CholeskyDecomposition( - normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD) - .getSolver() - .solve(jTr); - } catch (NonPositiveDefiniteMatrixException e) { - throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e); - } - } - }, - /** - * Solve the linear least squares problem using the {@link - * SingularValueDecomposition}. - * - * <p> This method is slower, but can provide a solution for rank deficient and - * nearly singular systems. - */ - SVD { - @Override - protected RealVector solve(final RealMatrix jacobian, - final RealVector residuals) { - return new SingularValueDecomposition(jacobian) - .getSolver() - .solve(residuals); - } - }; - - /** - * Solve the linear least squares problem Jx=r. - * - * @param jacobian the Jacobian matrix, J. the number of rows >= the number or - * columns. - * @param residuals the computed residuals, r. - * @return the solution x, to the linear least squares problem Jx=r. - * @throws ConvergenceException if the matrix properties (e.g. singular) do not - * permit a solution. - */ - protected abstract RealVector solve(RealMatrix jacobian, - RealVector residuals); - } - - /** - * The singularity threshold for matrix decompositions. Determines when a {@link - * ConvergenceException} is thrown. The current value was the default value for {@link - * LUDecomposition}. - */ - private static final double SINGULARITY_THRESHOLD = 1e-11; - - /** Indicator for using LU decomposition. */ - private final Decomposition decomposition; - - /** - * Creates a Gauss Newton optimizer. - * <p/> - * The default for the algorithm is to solve the normal equations using QR - * decomposition. - */ - public GaussNewtonOptimizer() { - this(Decomposition.QR); - } - - /** - * Create a Gauss Newton optimizer that uses the given decomposition algorithm to - * solve the normal equations. - * - * @param decomposition the {@link Decomposition} algorithm. - */ - public GaussNewtonOptimizer(final Decomposition decomposition) { - this.decomposition = decomposition; - } - - /** - * Get the matrix decomposition algorithm used to solve the normal equations. - * - * @return the matrix {@link Decomposition} algoritm. - */ - public Decomposition getDecomposition() { - return this.decomposition; - } - - /** - * Configure the decomposition algorithm. - * - * @param newDecomposition the {@link Decomposition} algorithm to use. - * @return a new instance. - */ - public GaussNewtonOptimizer withDecomposition(final Decomposition newDecomposition) { - return new GaussNewtonOptimizer(newDecomposition); - } - - /** {@inheritDoc} */ - public Optimum optimize(final LeastSquaresProblem lsp) { - //create local evaluation and iteration counts - final Incrementor evaluationCounter = lsp.getEvaluationCounter(); - final Incrementor iterationCounter = lsp.getIterationCounter(); - final ConvergenceChecker<Evaluation> checker - = lsp.getConvergenceChecker(); - - // Computation will be useless without a checker (see "for-loop"). - if (checker == null) { - throw new NullArgumentException(); - } - - RealVector currentPoint = lsp.getStart(); - - // iterate until convergence is reached - Evaluation current = null; - while (true) { - iterationCounter.incrementCount(); - - // evaluate the objective function and its jacobian - Evaluation previous = current; - // Value of the objective function at "currentPoint". - evaluationCounter.incrementCount(); - current = lsp.evaluate(currentPoint); - final RealVector currentResiduals = current.getResiduals(); - final RealMatrix weightedJacobian = current.getJacobian(); - currentPoint = current.getPoint(); - - // Check convergence. - if (previous != null) { - if (checker.converged(iterationCounter.getCount(), previous, current)) { - return new OptimumImpl( - current, - evaluationCounter.getCount(), - iterationCounter.getCount()); - } - } - - // solve the linearized least squares problem - final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals); - // update the estimated parameters - currentPoint = currentPoint.add(dX); - } - } - - @Override - public String toString() { - return "GaussNewtonOptimizer{" + - "decomposition=" + decomposition + - '}'; - } - - /** - * Compute the normal matrix, J<sup>T</sup>J. - * - * @param jacobian the m by n jacobian matrix, J. Input. - * @param residuals the m by 1 residual vector, r. Input. - * @return the n by n normal matrix and the n by 1 J<sup>Tr vector. - */ - private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian, - final RealVector residuals) { - //since the normal matrix is symmetric, we only need to compute half of it. - final int nR = jacobian.getRowDimension(); - final int nC = jacobian.getColumnDimension(); - //allocate space for return values - final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC); - final RealVector jTr = new ArrayRealVector(nC); - //for each measurement - for (int i = 0; i < nR; ++i) { - //compute JTr for measurement i - for (int j = 0; j < nC; j++) { - jTr.setEntry(j, jTr.getEntry(j) + - residuals.getEntry(i) * jacobian.getEntry(i, j)); - } - - // add the the contribution to the normal matrix for measurement i - for (int k = 0; k < nC; ++k) { - //only compute the upper triangular part - for (int l = k; l < nC; ++l) { - normal.setEntry(k, l, normal.getEntry(k, l) + - jacobian.getEntry(i, k) * jacobian.getEntry(i, l)); - } - } - } - //copy the upper triangular part to the lower triangular part. - for (int i = 0; i < nC; i++) { - for (int j = 0; j < i; j++) { - normal.setEntry(i, j, normal.getEntry(j, i)); - } - } - return new Pair<RealMatrix, RealVector>(normal, jTr); - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java deleted file mode 100644 index 1c09874..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LeastSquaresAdapter.java +++ /dev/null @@ -1,77 +0,0 @@ -/* - * 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.leastsquares; - -import org.apache.commons.math3.linear.RealVector; -import org.apache.commons.math3.optim.ConvergenceChecker; -import org.apache.commons.math3.util.Incrementor; - -/** - * An adapter that delegates to another implementation of {@link LeastSquaresProblem}. - * - * @since 3.3 - */ -public class LeastSquaresAdapter implements LeastSquaresProblem { - - /** the delegate problem */ - private final LeastSquaresProblem problem; - - /** - * Delegate the {@link LeastSquaresProblem} interface to the given implementation. - * - * @param problem the delegate - */ - public LeastSquaresAdapter(final LeastSquaresProblem problem) { - this.problem = problem; - } - - /** {@inheritDoc} */ - public RealVector getStart() { - return problem.getStart(); - } - - /** {@inheritDoc} */ - public int getObservationSize() { - return problem.getObservationSize(); - } - - /** {@inheritDoc} */ - public int getParameterSize() { - return problem.getParameterSize(); - } - - /** {@inheritDoc} - * @param point*/ - public Evaluation evaluate(final RealVector point) { - return problem.evaluate(point); - } - - /** {@inheritDoc} */ - public Incrementor getEvaluationCounter() { - return problem.getEvaluationCounter(); - } - - /** {@inheritDoc} */ - public Incrementor getIterationCounter() { - return problem.getIterationCounter(); - } - - /** {@inheritDoc} */ - public ConvergenceChecker<Evaluation> getConvergenceChecker() { - return problem.getConvergenceChecker(); - } -}
