http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/filter/KalmanFilter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/filter/KalmanFilter.java b/src/main/java/org/apache/commons/math3/filter/KalmanFilter.java deleted file mode 100644 index 2f204bf..0000000 --- a/src/main/java/org/apache/commons/math3/filter/KalmanFilter.java +++ /dev/null @@ -1,388 +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.filter; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.CholeskyDecomposition; -import org.apache.commons.math3.linear.MatrixDimensionMismatchException; -import org.apache.commons.math3.linear.MatrixUtils; -import org.apache.commons.math3.linear.NonSquareMatrixException; -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.util.MathUtils; - -/** - * Implementation of a Kalman filter to estimate the state <i>x<sub>k</sub></i> - * of a discrete-time controlled process that is governed by the linear - * stochastic difference equation: - * - * <pre> - * <i>x<sub>k</sub></i> = <b>A</b><i>x<sub>k-1</sub></i> + <b>B</b><i>u<sub>k-1</sub></i> + <i>w<sub>k-1</sub></i> - * </pre> - * - * with a measurement <i>x<sub>k</sub></i> that is - * - * <pre> - * <i>z<sub>k</sub></i> = <b>H</b><i>x<sub>k</sub></i> + <i>v<sub>k</sub></i>. - * </pre> - * - * <p> - * The random variables <i>w<sub>k</sub></i> and <i>v<sub>k</sub></i> represent - * the process and measurement noise and are assumed to be independent of each - * other and distributed with normal probability (white noise). - * <p> - * The Kalman filter cycle involves the following steps: - * <ol> - * <li>predict: project the current state estimate ahead in time</li> - * <li>correct: adjust the projected estimate by an actual measurement</li> - * </ol> - * <p> - * The Kalman filter is initialized with a {@link ProcessModel} and a - * {@link MeasurementModel}, which contain the corresponding transformation and - * noise covariance matrices. The parameter names used in the respective models - * correspond to the following names commonly used in the mathematical - * literature: - * <ul> - * <li>A - state transition matrix</li> - * <li>B - control input matrix</li> - * <li>H - measurement matrix</li> - * <li>Q - process noise covariance matrix</li> - * <li>R - measurement noise covariance matrix</li> - * <li>P - error covariance matrix</li> - * </ul> - * - * @see <a href="http://www.cs.unc.edu/~welch/kalman/";>Kalman filter - * resources</a> - * @see <a href="http://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdf";>An - * introduction to the Kalman filter by Greg Welch and Gary Bishop</a> - * @see <a href="http://academic.csuohio.edu/simond/courses/eec644/kalman.pdf";> - * Kalman filter example by Dan Simon</a> - * @see ProcessModel - * @see MeasurementModel - * @since 3.0 - */ -public class KalmanFilter { - /** The process model used by this filter instance. */ - private final ProcessModel processModel; - /** The measurement model used by this filter instance. */ - private final MeasurementModel measurementModel; - /** The transition matrix, equivalent to A. */ - private RealMatrix transitionMatrix; - /** The transposed transition matrix. */ - private RealMatrix transitionMatrixT; - /** The control matrix, equivalent to B. */ - private RealMatrix controlMatrix; - /** The measurement matrix, equivalent to H. */ - private RealMatrix measurementMatrix; - /** The transposed measurement matrix. */ - private RealMatrix measurementMatrixT; - /** The internal state estimation vector, equivalent to x hat. */ - private RealVector stateEstimation; - /** The error covariance matrix, equivalent to P. */ - private RealMatrix errorCovariance; - - /** - * Creates a new Kalman filter with the given process and measurement models. - * - * @param process - * the model defining the underlying process dynamics - * @param measurement - * the model defining the given measurement characteristics - * @throws NullArgumentException - * if any of the given inputs is null (except for the control matrix) - * @throws NonSquareMatrixException - * if the transition matrix is non square - * @throws DimensionMismatchException - * if the column dimension of the transition matrix does not match the dimension of the - * initial state estimation vector - * @throws MatrixDimensionMismatchException - * if the matrix dimensions do not fit together - */ - public KalmanFilter(final ProcessModel process, final MeasurementModel measurement) - throws NullArgumentException, NonSquareMatrixException, DimensionMismatchException, - MatrixDimensionMismatchException { - - MathUtils.checkNotNull(process); - MathUtils.checkNotNull(measurement); - - this.processModel = process; - this.measurementModel = measurement; - - transitionMatrix = processModel.getStateTransitionMatrix(); - MathUtils.checkNotNull(transitionMatrix); - transitionMatrixT = transitionMatrix.transpose(); - - // create an empty matrix if no control matrix was given - if (processModel.getControlMatrix() == null) { - controlMatrix = new Array2DRowRealMatrix(); - } else { - controlMatrix = processModel.getControlMatrix(); - } - - measurementMatrix = measurementModel.getMeasurementMatrix(); - MathUtils.checkNotNull(measurementMatrix); - measurementMatrixT = measurementMatrix.transpose(); - - // check that the process and measurement noise matrices are not null - // they will be directly accessed from the model as they may change - // over time - RealMatrix processNoise = processModel.getProcessNoise(); - MathUtils.checkNotNull(processNoise); - RealMatrix measNoise = measurementModel.getMeasurementNoise(); - MathUtils.checkNotNull(measNoise); - - // set the initial state estimate to a zero vector if it is not - // available from the process model - if (processModel.getInitialStateEstimate() == null) { - stateEstimation = new ArrayRealVector(transitionMatrix.getColumnDimension()); - } else { - stateEstimation = processModel.getInitialStateEstimate(); - } - - if (transitionMatrix.getColumnDimension() != stateEstimation.getDimension()) { - throw new DimensionMismatchException(transitionMatrix.getColumnDimension(), - stateEstimation.getDimension()); - } - - // initialize the error covariance to the process noise if it is not - // available from the process model - if (processModel.getInitialErrorCovariance() == null) { - errorCovariance = processNoise.copy(); - } else { - errorCovariance = processModel.getInitialErrorCovariance(); - } - - // sanity checks, the control matrix B may be null - - // A must be a square matrix - if (!transitionMatrix.isSquare()) { - throw new NonSquareMatrixException( - transitionMatrix.getRowDimension(), - transitionMatrix.getColumnDimension()); - } - - // row dimension of B must be equal to A - // if no control matrix is available, the row and column dimension will be 0 - if (controlMatrix != null && - controlMatrix.getRowDimension() > 0 && - controlMatrix.getColumnDimension() > 0 && - controlMatrix.getRowDimension() != transitionMatrix.getRowDimension()) { - throw new MatrixDimensionMismatchException(controlMatrix.getRowDimension(), - controlMatrix.getColumnDimension(), - transitionMatrix.getRowDimension(), - controlMatrix.getColumnDimension()); - } - - // Q must be equal to A - MatrixUtils.checkAdditionCompatible(transitionMatrix, processNoise); - - // column dimension of H must be equal to row dimension of A - if (measurementMatrix.getColumnDimension() != transitionMatrix.getRowDimension()) { - throw new MatrixDimensionMismatchException(measurementMatrix.getRowDimension(), - measurementMatrix.getColumnDimension(), - measurementMatrix.getRowDimension(), - transitionMatrix.getRowDimension()); - } - - // row dimension of R must be equal to row dimension of H - if (measNoise.getRowDimension() != measurementMatrix.getRowDimension()) { - throw new MatrixDimensionMismatchException(measNoise.getRowDimension(), - measNoise.getColumnDimension(), - measurementMatrix.getRowDimension(), - measNoise.getColumnDimension()); - } - } - - /** - * Returns the dimension of the state estimation vector. - * - * @return the state dimension - */ - public int getStateDimension() { - return stateEstimation.getDimension(); - } - - /** - * Returns the dimension of the measurement vector. - * - * @return the measurement vector dimension - */ - public int getMeasurementDimension() { - return measurementMatrix.getRowDimension(); - } - - /** - * Returns the current state estimation vector. - * - * @return the state estimation vector - */ - public double[] getStateEstimation() { - return stateEstimation.toArray(); - } - - /** - * Returns a copy of the current state estimation vector. - * - * @return the state estimation vector - */ - public RealVector getStateEstimationVector() { - return stateEstimation.copy(); - } - - /** - * Returns the current error covariance matrix. - * - * @return the error covariance matrix - */ - public double[][] getErrorCovariance() { - return errorCovariance.getData(); - } - - /** - * Returns a copy of the current error covariance matrix. - * - * @return the error covariance matrix - */ - public RealMatrix getErrorCovarianceMatrix() { - return errorCovariance.copy(); - } - - /** - * Predict the internal state estimation one time step ahead. - */ - public void predict() { - predict((RealVector) null); - } - - /** - * Predict the internal state estimation one time step ahead. - * - * @param u - * the control vector - * @throws DimensionMismatchException - * if the dimension of the control vector does not fit - */ - public void predict(final double[] u) throws DimensionMismatchException { - predict(new ArrayRealVector(u, false)); - } - - /** - * Predict the internal state estimation one time step ahead. - * - * @param u - * the control vector - * @throws DimensionMismatchException - * if the dimension of the control vector does not match - */ - public void predict(final RealVector u) throws DimensionMismatchException { - // sanity checks - if (u != null && - u.getDimension() != controlMatrix.getColumnDimension()) { - throw new DimensionMismatchException(u.getDimension(), - controlMatrix.getColumnDimension()); - } - - // project the state estimation ahead (a priori state) - // xHat(k)- = A * xHat(k-1) + B * u(k-1) - stateEstimation = transitionMatrix.operate(stateEstimation); - - // add control input if it is available - if (u != null) { - stateEstimation = stateEstimation.add(controlMatrix.operate(u)); - } - - // project the error covariance ahead - // P(k)- = A * P(k-1) * A' + Q - errorCovariance = transitionMatrix.multiply(errorCovariance) - .multiply(transitionMatrixT) - .add(processModel.getProcessNoise()); - } - - /** - * Correct the current state estimate with an actual measurement. - * - * @param z - * the measurement vector - * @throws NullArgumentException - * if the measurement vector is {@code null} - * @throws DimensionMismatchException - * if the dimension of the measurement vector does not fit - * @throws SingularMatrixException - * if the covariance matrix could not be inverted - */ - public void correct(final double[] z) - throws NullArgumentException, DimensionMismatchException, SingularMatrixException { - correct(new ArrayRealVector(z, false)); - } - - /** - * Correct the current state estimate with an actual measurement. - * - * @param z - * the measurement vector - * @throws NullArgumentException - * if the measurement vector is {@code null} - * @throws DimensionMismatchException - * if the dimension of the measurement vector does not fit - * @throws SingularMatrixException - * if the covariance matrix could not be inverted - */ - public void correct(final RealVector z) - throws NullArgumentException, DimensionMismatchException, SingularMatrixException { - - // sanity checks - MathUtils.checkNotNull(z); - if (z.getDimension() != measurementMatrix.getRowDimension()) { - throw new DimensionMismatchException(z.getDimension(), - measurementMatrix.getRowDimension()); - } - - // S = H * P(k) * H' + R - RealMatrix s = measurementMatrix.multiply(errorCovariance) - .multiply(measurementMatrixT) - .add(measurementModel.getMeasurementNoise()); - - // Inn = z(k) - H * xHat(k)- - RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); - - // calculate gain matrix - // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 - // K(k) = P(k)- * H' * S^-1 - - // instead of calculating the inverse of S we can rearrange the formula, - // and then solve the linear equation A x X = B with A = S', X = K' and B = (H * P)' - - // K(k) * S = P(k)- * H' - // S' * K(k)' = H * P(k)-' - RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver() - .solve(measurementMatrix.multiply(errorCovariance.transpose())) - .transpose(); - - // update estimate with measurement z(k) - // xHat(k) = xHat(k)- + K * Inn - stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); - - // update covariance of prediction error - // P(k) = (I - K * H) * P(k)- - RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); - errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/filter/MeasurementModel.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/filter/MeasurementModel.java b/src/main/java/org/apache/commons/math3/filter/MeasurementModel.java deleted file mode 100644 index 2e0a379..0000000 --- a/src/main/java/org/apache/commons/math3/filter/MeasurementModel.java +++ /dev/null @@ -1,44 +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.filter; - -import org.apache.commons.math3.linear.RealMatrix; - -/** - * Defines the measurement model for the use with a {@link KalmanFilter}. - * - * @since 3.0 - */ -public interface MeasurementModel { - /** - * Returns the measurement matrix. - * - * @return the measurement matrix - */ - RealMatrix getMeasurementMatrix(); - - /** - * Returns the measurement noise matrix. This method is called by the {@link KalmanFilter} every - * correction step, so implementations of this interface may return a modified measurement noise - * depending on the current iteration step. - * - * @return the measurement noise matrix - * @see KalmanFilter#correct(double[]) - * @see KalmanFilter#correct(org.apache.commons.math3.linear.RealVector) - */ - RealMatrix getMeasurementNoise(); -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/filter/ProcessModel.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/filter/ProcessModel.java b/src/main/java/org/apache/commons/math3/filter/ProcessModel.java deleted file mode 100644 index 179ed1b..0000000 --- a/src/main/java/org/apache/commons/math3/filter/ProcessModel.java +++ /dev/null @@ -1,73 +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.filter; - -import org.apache.commons.math3.linear.RealMatrix; -import org.apache.commons.math3.linear.RealVector; - -/** - * Defines the process dynamics model for the use with a {@link KalmanFilter}. - * - * @since 3.0 - */ -public interface ProcessModel { - /** - * Returns the state transition matrix. - * - * @return the state transition matrix - */ - RealMatrix getStateTransitionMatrix(); - - /** - * Returns the control matrix. - * - * @return the control matrix - */ - RealMatrix getControlMatrix(); - - /** - * Returns the process noise matrix. This method is called by the {@link KalmanFilter} every - * prediction step, so implementations of this interface may return a modified process noise - * depending on the current iteration step. - * - * @return the process noise matrix - * @see KalmanFilter#predict() - * @see KalmanFilter#predict(double[]) - * @see KalmanFilter#predict(RealVector) - */ - RealMatrix getProcessNoise(); - - /** - * Returns the initial state estimation vector. - * <p> - * <b>Note:</b> if the return value is zero, the Kalman filter will initialize the - * state estimation with a zero vector. - * - * @return the initial state estimation vector - */ - RealVector getInitialStateEstimate(); - - /** - * Returns the initial error covariance matrix. - * <p> - * <b>Note:</b> if the return value is zero, the Kalman filter will initialize the - * error covariance with the process noise matrix. - * - * @return the initial error covariance matrix - */ - RealMatrix getInitialErrorCovariance(); -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/filter/package-info.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/filter/package-info.java b/src/main/java/org/apache/commons/math3/filter/package-info.java deleted file mode 100644 index 4e76536..0000000 --- a/src/main/java/org/apache/commons/math3/filter/package-info.java +++ /dev/null @@ -1,20 +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. - */ -/** - * Implementations of common discrete-time linear filters. - */ -package org.apache.commons.math3.filter; http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/AbstractCurveFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/AbstractCurveFitter.java b/src/main/java/org/apache/commons/math3/fitting/AbstractCurveFitter.java deleted file mode 100644 index 6f39c23..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/AbstractCurveFitter.java +++ /dev/null @@ -1,147 +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.MultivariateVectorFunction; -import org.apache.commons.math3.analysis.MultivariateMatrixFunction; -import org.apache.commons.math3.analysis.ParametricUnivariateFunction; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer; -import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; -import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer; - -/** - * Base class that contains common code for fitting parametric univariate - * real functions <code>y = f(p<sub>i</sub>;x)</code>, where {@code x} is - * the independent variable and the <code>p<sub>i</sub></code> are the - * <em>parameters</em>. - * <br/> - * A fitter will find the optimal values of the parameters by - * <em>fitting</em> the curve so it remains very close to a set of - * {@code N} observed points <code>(x<sub>k</sub>, y<sub>k</sub>)</code>, - * {@code 0 <= k < N}. - * <br/> - * An algorithm usually performs the fit by finding the parameter - * values that minimizes the objective function - * <pre><code> - * ∑y<sub>k</sub> - f(x<sub>k</sub>)<sup>2</sup>, - * </code></pre> - * which is actually a least-squares problem. - * This class contains boilerplate code for calling the - * {@link #fit(Collection)} method for obtaining the parameters. - * The problem setup, such as the choice of optimization algorithm - * for fitting a specific function is delegated to subclasses. - * - * @since 3.3 - */ -public abstract class AbstractCurveFitter { - /** - * Fits a curve. - * This method computes the coefficients of the curve that best - * fit the sample of observed points. - * - * @param points Observations. - * @return the fitted parameters. - */ - public double[] fit(Collection<WeightedObservedPoint> points) { - // Perform the fit. - return getOptimizer().optimize(getProblem(points)).getPoint().toArray(); - } - - /** - * Creates an optimizer set up to fit the appropriate curve. - * <p> - * The default implementation uses a {@link LevenbergMarquardtOptimizer - * Levenberg-Marquardt} optimizer. - * </p> - * @return the optimizer to use for fitting the curve to the - * given {@code points}. - */ - protected LeastSquaresOptimizer getOptimizer() { - return new LevenbergMarquardtOptimizer(); - } - - /** - * Creates a least squares problem corresponding to the appropriate curve. - * - * @param points Sample points. - * @return the least squares problem to use for fitting the curve to the - * given {@code points}. - */ - protected abstract LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> points); - - /** - * Vector function for computing function theoretical values. - */ - protected static class TheoreticalValuesFunction { - /** Function to fit. */ - private final ParametricUnivariateFunction f; - /** Observations. */ - private final double[] points; - - /** - * @param f function to fit. - * @param observations Observations. - */ - public TheoreticalValuesFunction(final ParametricUnivariateFunction f, - final Collection<WeightedObservedPoint> observations) { - this.f = f; - - final int len = observations.size(); - this.points = new double[len]; - int i = 0; - for (WeightedObservedPoint obs : observations) { - this.points[i++] = obs.getX(); - } - } - - /** - * @return the model function values. - */ - public MultivariateVectorFunction getModelFunction() { - return new MultivariateVectorFunction() { - /** {@inheritDoc} */ - public double[] value(double[] p) { - final int len = points.length; - final double[] values = new double[len]; - for (int i = 0; i < len; i++) { - values[i] = f.value(points[i], p); - } - - return values; - } - }; - } - - /** - * @return the model function Jacobian. - */ - public MultivariateMatrixFunction getModelFunctionJacobian() { - return new MultivariateMatrixFunction() { - public double[][] value(double[] p) { - final int len = points.length; - final double[][] jacobian = new double[len][]; - for (int i = 0; i < len; i++) { - jacobian[i] = f.gradient(points[i], p); - } - return jacobian; - } - }; - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/CurveFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/CurveFitter.java b/src/main/java/org/apache/commons/math3/fitting/CurveFitter.java deleted file mode 100644 index 0ddeacc..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/CurveFitter.java +++ /dev/null @@ -1,232 +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.List; -import org.apache.commons.math3.analysis.MultivariateVectorFunction; -import org.apache.commons.math3.analysis.MultivariateMatrixFunction; -import org.apache.commons.math3.analysis.ParametricUnivariateFunction; -import org.apache.commons.math3.optim.MaxEval; -import org.apache.commons.math3.optim.InitialGuess; -import org.apache.commons.math3.optim.PointVectorValuePair; -import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer; -import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction; -import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian; -import org.apache.commons.math3.optim.nonlinear.vector.Target; -import org.apache.commons.math3.optim.nonlinear.vector.Weight; - -/** - * Fitter for parametric univariate real functions y = f(x). - * <br/> - * When a univariate real function y = f(x) does depend on some - * unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>, - * this class can be used to find these parameters. It does this - * by <em>fitting</em> the curve so it remains very close to a set of - * observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>, - * y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting - * is done by finding the parameters values that minimizes the objective - * function ∑(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is - * really a least squares problem. - * - * @param <T> Function to use for the fit. - * - * @since 2.0 - * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and - * {@link WeightedObservedPoints} instead. - */ -@Deprecated -public class CurveFitter<T extends ParametricUnivariateFunction> { - /** Optimizer to use for the fitting. */ - private final MultivariateVectorOptimizer optimizer; - /** Observed points. */ - private final List<WeightedObservedPoint> observations; - - /** - * Simple constructor. - * - * @param optimizer Optimizer to use for the fitting. - * @since 3.1 - */ - public CurveFitter(final MultivariateVectorOptimizer optimizer) { - this.optimizer = optimizer; - observations = new ArrayList<WeightedObservedPoint>(); - } - - /** Add an observed (x,y) point to the sample with unit weight. - * <p>Calling this method is equivalent to call - * {@code addObservedPoint(1.0, x, y)}.</p> - * @param x abscissa of the point - * @param y observed value of the point at x, after fitting we should - * have f(x) as close as possible to this value - * @see #addObservedPoint(double, double, double) - * @see #addObservedPoint(WeightedObservedPoint) - * @see #getObservations() - */ - public void addObservedPoint(double x, double y) { - addObservedPoint(1.0, x, y); - } - - /** Add an observed weighted (x,y) point to the sample. - * @param weight weight of the observed point in the fit - * @param x abscissa of the point - * @param y observed value of the point at x, after fitting we should - * have f(x) as close as possible to this value - * @see #addObservedPoint(double, double) - * @see #addObservedPoint(WeightedObservedPoint) - * @see #getObservations() - */ - public void addObservedPoint(double weight, double x, double y) { - observations.add(new WeightedObservedPoint(weight, x, y)); - } - - /** Add an observed weighted (x,y) point to the sample. - * @param observed observed point to add - * @see #addObservedPoint(double, double) - * @see #addObservedPoint(double, double, double) - * @see #getObservations() - */ - public void addObservedPoint(WeightedObservedPoint observed) { - observations.add(observed); - } - - /** Get the observed points. - * @return observed points - * @see #addObservedPoint(double, double) - * @see #addObservedPoint(double, double, double) - * @see #addObservedPoint(WeightedObservedPoint) - */ - public WeightedObservedPoint[] getObservations() { - return observations.toArray(new WeightedObservedPoint[observations.size()]); - } - - /** - * Remove all observations. - */ - public void clearObservations() { - observations.clear(); - } - - /** - * Fit a curve. - * This method compute the coefficients of the curve that best - * fit the sample of observed points previously given through calls - * to the {@link #addObservedPoint(WeightedObservedPoint) - * addObservedPoint} method. - * - * @param f parametric function to fit. - * @param initialGuess first guess of the function parameters. - * @return the fitted parameters. - * @throws org.apache.commons.math3.exception.DimensionMismatchException - * if the start point dimension is wrong. - */ - public double[] fit(T f, final double[] initialGuess) { - return fit(Integer.MAX_VALUE, f, initialGuess); - } - - /** - * Fit a curve. - * This method compute the coefficients of the curve that best - * fit the sample of observed points previously given through calls - * to the {@link #addObservedPoint(WeightedObservedPoint) - * addObservedPoint} method. - * - * @param f parametric function to fit. - * @param initialGuess first guess of the function parameters. - * @param maxEval Maximum number of function evaluations. - * @return the fitted parameters. - * @throws org.apache.commons.math3.exception.TooManyEvaluationsException - * if the number of allowed evaluations is exceeded. - * @throws org.apache.commons.math3.exception.DimensionMismatchException - * if the start point dimension is wrong. - * @since 3.0 - */ - public double[] fit(int maxEval, T f, - final double[] initialGuess) { - // Prepare least squares problem. - double[] target = new double[observations.size()]; - double[] weights = new double[observations.size()]; - int i = 0; - for (WeightedObservedPoint point : observations) { - target[i] = point.getY(); - weights[i] = point.getWeight(); - ++i; - } - - // Input to the optimizer: the model and its Jacobian. - final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f); - - // Perform the fit. - final PointVectorValuePair optimum - = optimizer.optimize(new MaxEval(maxEval), - model.getModelFunction(), - model.getModelFunctionJacobian(), - new Target(target), - new Weight(weights), - new InitialGuess(initialGuess)); - // Extract the coefficients. - return optimum.getPointRef(); - } - - /** Vectorial function computing function theoretical values. */ - private class TheoreticalValuesFunction { - /** Function to fit. */ - private final ParametricUnivariateFunction f; - - /** - * @param f function to fit. - */ - public TheoreticalValuesFunction(final ParametricUnivariateFunction f) { - this.f = f; - } - - /** - * @return the model function values. - */ - public ModelFunction getModelFunction() { - return new ModelFunction(new MultivariateVectorFunction() { - /** {@inheritDoc} */ - public double[] value(double[] point) { - // compute the residuals - final double[] values = new double[observations.size()]; - int i = 0; - for (WeightedObservedPoint observed : observations) { - values[i++] = f.value(observed.getX(), point); - } - - return values; - } - }); - } - - /** - * @return the model function Jacobian. - */ - public ModelFunctionJacobian getModelFunctionJacobian() { - return new ModelFunctionJacobian(new MultivariateMatrixFunction() { - public double[][] value(double[] point) { - final double[][] jacobian = new double[observations.size()][]; - int i = 0; - for (WeightedObservedPoint observed : observations) { - jacobian[i++] = f.gradient(observed.getX(), point); - } - return jacobian; - } - }); - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/GaussianCurveFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/GaussianCurveFitter.java b/src/main/java/org/apache/commons/math3/fitting/GaussianCurveFitter.java deleted file mode 100644 index 64e0bda..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/GaussianCurveFitter.java +++ /dev/null @@ -1,430 +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.Collections; -import java.util.Comparator; -import java.util.List; - -import org.apache.commons.math3.analysis.function.Gaussian; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.OutOfRangeException; -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.Gaussian.Parametric Gaussian} - * function. - * <br/> - * The {@link #withStartPoint(double[]) initial guess values} must be passed - * in the following order: - * <ul> - * <li>Normalization</li> - * <li>Mean</li> - * <li>Sigma</li> - * </ul> - * The optimal values will be returned in the same order. - * - * <p> - * Usage example: - * <pre> - * WeightedObservedPoints obs = new WeightedObservedPoints(); - * obs.add(4.0254623, 531026.0); - * obs.add(4.03128248, 984167.0); - * obs.add(4.03839603, 1887233.0); - * obs.add(4.04421621, 2687152.0); - * obs.add(4.05132976, 3461228.0); - * obs.add(4.05326982, 3580526.0); - * obs.add(4.05779662, 3439750.0); - * obs.add(4.0636168, 2877648.0); - * obs.add(4.06943698, 2175960.0); - * obs.add(4.07525716, 1447024.0); - * obs.add(4.08237071, 717104.0); - * obs.add(4.08366408, 620014.0); - * double[] parameters = GaussianCurveFitter.create().fit(obs.toList()); - * </pre> - * - * @since 3.3 - */ -public class GaussianCurveFitter extends AbstractCurveFitter { - /** Parametric function to be fitted. */ - private static final Gaussian.Parametric FUNCTION = new Gaussian.Parametric() { - @Override - public double value(double x, double ... p) { - double v = Double.POSITIVE_INFINITY; - try { - v = super.value(x, p); - } catch (NotStrictlyPositiveException e) { // NOPMD - // Do nothing. - } - return v; - } - - @Override - public double[] gradient(double x, double ... p) { - double[] v = { Double.POSITIVE_INFINITY, - Double.POSITIVE_INFINITY, - Double.POSITIVE_INFINITY }; - try { - v = super.gradient(x, p); - } catch (NotStrictlyPositiveException e) { // NOPMD - // Do nothing. - } - return v; - } - }; - /** 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 GaussianCurveFitter(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 GaussianCurveFitter create() { - return new GaussianCurveFitter(null, Integer.MAX_VALUE); - } - - /** - * Configure the start point (initial guess). - * @param newStart new start point (initial guess) - * @return a new instance. - */ - public GaussianCurveFitter withStartPoint(double[] newStart) { - return new GaussianCurveFitter(newStart.clone(), - maxIter); - } - - /** - * Configure the maximum number of iterations. - * @param newMaxIter maximum number of iterations - * @return a new instance. - */ - public GaussianCurveFitter withMaxIterations(int newMaxIter) { - return new GaussianCurveFitter(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 least squares problem 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(); - - } - - /** - * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma} - * of a {@link org.apache.commons.math3.analysis.function.Gaussian.Parametric} - * based on the specified observed points. - */ - public static class ParameterGuesser { - /** Normalization factor. */ - private final double norm; - /** Mean. */ - private final double mean; - /** Standard deviation. */ - private final double sigma; - - /** - * Constructs instance with the specified observed points. - * - * @param observations Observed points from which to guess the - * parameters of the Gaussian. - * @throws NullArgumentException if {@code observations} is - * {@code null}. - * @throws NumberIsTooSmallException if there are less than 3 - * observations. - */ - public ParameterGuesser(Collection<WeightedObservedPoint> observations) { - if (observations == null) { - throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY); - } - if (observations.size() < 3) { - throw new NumberIsTooSmallException(observations.size(), 3, true); - } - - final List<WeightedObservedPoint> sorted = sortObservations(observations); - final double[] params = basicGuess(sorted.toArray(new WeightedObservedPoint[0])); - - norm = params[0]; - mean = params[1]; - sigma = params[2]; - } - - /** - * Gets an estimation of the parameters. - * - * @return the guessed parameters, in the following order: - * <ul> - * <li>Normalization factor</li> - * <li>Mean</li> - * <li>Standard deviation</li> - * </ul> - */ - public double[] guess() { - return new double[] { norm, mean, sigma }; - } - - /** - * Sort the observations. - * - * @param unsorted Input observations. - * @return the input observations, sorted. - */ - private List<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) { - final List<WeightedObservedPoint> observations = new ArrayList<WeightedObservedPoint>(unsorted); - - final Comparator<WeightedObservedPoint> cmp = new Comparator<WeightedObservedPoint>() { - public int compare(WeightedObservedPoint p1, - WeightedObservedPoint p2) { - if (p1 == null && p2 == null) { - return 0; - } - if (p1 == null) { - return -1; - } - if (p2 == null) { - return 1; - } - if (p1.getX() < p2.getX()) { - return -1; - } - if (p1.getX() > p2.getX()) { - return 1; - } - if (p1.getY() < p2.getY()) { - return -1; - } - if (p1.getY() > p2.getY()) { - return 1; - } - if (p1.getWeight() < p2.getWeight()) { - return -1; - } - if (p1.getWeight() > p2.getWeight()) { - return 1; - } - return 0; - } - }; - - Collections.sort(observations, cmp); - return observations; - } - - /** - * Guesses the parameters based on the specified observed points. - * - * @param points Observed points, sorted. - * @return the guessed parameters (normalization factor, mean and - * sigma). - */ - private double[] basicGuess(WeightedObservedPoint[] points) { - final int maxYIdx = findMaxY(points); - final double n = points[maxYIdx].getY(); - final double m = points[maxYIdx].getX(); - - double fwhmApprox; - try { - final double halfY = n + ((m - n) / 2); - final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY); - final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY); - fwhmApprox = fwhmX2 - fwhmX1; - } catch (OutOfRangeException e) { - // TODO: Exceptions should not be used for flow control. - fwhmApprox = points[points.length - 1].getX() - points[0].getX(); - } - final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2))); - - return new double[] { n, m, s }; - } - - /** - * Finds index of point in specified points with the largest Y. - * - * @param points Points to search. - * @return the index in specified points array. - */ - private int findMaxY(WeightedObservedPoint[] points) { - int maxYIdx = 0; - for (int i = 1; i < points.length; i++) { - if (points[i].getY() > points[maxYIdx].getY()) { - maxYIdx = i; - } - } - return maxYIdx; - } - - /** - * Interpolates using the specified points to determine X at the - * specified Y. - * - * @param points Points to use for interpolation. - * @param startIdx Index within points from which to start the search for - * interpolation bounds points. - * @param idxStep Index step for searching interpolation bounds points. - * @param y Y value for which X should be determined. - * @return the value of X for the specified Y. - * @throws ZeroException if {@code idxStep} is 0. - * @throws OutOfRangeException if specified {@code y} is not within the - * range of the specified {@code points}. - */ - private double interpolateXAtY(WeightedObservedPoint[] points, - int startIdx, - int idxStep, - double y) - throws OutOfRangeException { - if (idxStep == 0) { - throw new ZeroException(); - } - final WeightedObservedPoint[] twoPoints - = getInterpolationPointsForY(points, startIdx, idxStep, y); - final WeightedObservedPoint p1 = twoPoints[0]; - final WeightedObservedPoint p2 = twoPoints[1]; - if (p1.getY() == y) { - return p1.getX(); - } - if (p2.getY() == y) { - return p2.getX(); - } - return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) / - (p2.getY() - p1.getY())); - } - - /** - * Gets the two bounding interpolation points from the specified points - * suitable for determining X at the specified Y. - * - * @param points Points to use for interpolation. - * @param startIdx Index within points from which to start search for - * interpolation bounds points. - * @param idxStep Index step for search for interpolation bounds points. - * @param y Y value for which X should be determined. - * @return the array containing two points suitable for determining X at - * the specified Y. - * @throws ZeroException if {@code idxStep} is 0. - * @throws OutOfRangeException if specified {@code y} is not within the - * range of the specified {@code points}. - */ - private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points, - int startIdx, - int idxStep, - double y) - throws OutOfRangeException { - if (idxStep == 0) { - throw new ZeroException(); - } - for (int i = startIdx; - idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length; - i += idxStep) { - final WeightedObservedPoint p1 = points[i]; - final WeightedObservedPoint p2 = points[i + idxStep]; - if (isBetween(y, p1.getY(), p2.getY())) { - if (idxStep < 0) { - return new WeightedObservedPoint[] { p2, p1 }; - } else { - return new WeightedObservedPoint[] { p1, p2 }; - } - } - } - - // Boundaries are replaced by dummy values because the raised - // exception is caught and the message never displayed. - // TODO: Exceptions should not be used for flow control. - throw new OutOfRangeException(y, - Double.NEGATIVE_INFINITY, - Double.POSITIVE_INFINITY); - } - - /** - * Determines whether a value is between two other values. - * - * @param value Value to test whether it is between {@code boundary1} - * and {@code boundary2}. - * @param boundary1 One end of the range. - * @param boundary2 Other end of the range. - * @return {@code true} if {@code value} is between {@code boundary1} and - * {@code boundary2} (inclusive), {@code false} otherwise. - */ - private boolean isBetween(double value, - double boundary1, - double boundary2) { - return (value >= boundary1 && value <= boundary2) || - (value >= boundary2 && value <= boundary1); - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/fitting/GaussianFitter.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/fitting/GaussianFitter.java b/src/main/java/org/apache/commons/math3/fitting/GaussianFitter.java deleted file mode 100644 index 3946540..0000000 --- a/src/main/java/org/apache/commons/math3/fitting/GaussianFitter.java +++ /dev/null @@ -1,364 +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.Arrays; -import java.util.Comparator; -import org.apache.commons.math3.analysis.function.Gaussian; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.ZeroException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer; -import org.apache.commons.math3.util.FastMath; - -/** - * Fits points to a {@link - * org.apache.commons.math3.analysis.function.Gaussian.Parametric Gaussian} function. - * <p> - * Usage example: - * <pre> - * GaussianFitter fitter = new GaussianFitter( - * new LevenbergMarquardtOptimizer()); - * fitter.addObservedPoint(4.0254623, 531026.0); - * fitter.addObservedPoint(4.03128248, 984167.0); - * fitter.addObservedPoint(4.03839603, 1887233.0); - * fitter.addObservedPoint(4.04421621, 2687152.0); - * fitter.addObservedPoint(4.05132976, 3461228.0); - * fitter.addObservedPoint(4.05326982, 3580526.0); - * fitter.addObservedPoint(4.05779662, 3439750.0); - * fitter.addObservedPoint(4.0636168, 2877648.0); - * fitter.addObservedPoint(4.06943698, 2175960.0); - * fitter.addObservedPoint(4.07525716, 1447024.0); - * fitter.addObservedPoint(4.08237071, 717104.0); - * fitter.addObservedPoint(4.08366408, 620014.0); - * double[] parameters = fitter.fit(); - * </pre> - * - * @since 2.2 - * @deprecated As of 3.3. Please use {@link GaussianCurveFitter} and - * {@link WeightedObservedPoints} instead. - */ -@Deprecated -public class GaussianFitter extends CurveFitter<Gaussian.Parametric> { - /** - * Constructs an instance using the specified optimizer. - * - * @param optimizer Optimizer to use for the fitting. - */ - public GaussianFitter(MultivariateVectorOptimizer optimizer) { - super(optimizer); - } - - /** - * Fits a Gaussian function to the observed points. - * - * @param initialGuess First guess values in the following order: - * <ul> - * <li>Norm</li> - * <li>Mean</li> - * <li>Sigma</li> - * </ul> - * @return the parameters of the Gaussian function that best fits the - * observed points (in the same order as above). - * @since 3.0 - */ - public double[] fit(double[] initialGuess) { - final Gaussian.Parametric f = new Gaussian.Parametric() { - @Override - public double value(double x, double ... p) { - double v = Double.POSITIVE_INFINITY; - try { - v = super.value(x, p); - } catch (NotStrictlyPositiveException e) { // NOPMD - // Do nothing. - } - return v; - } - - @Override - public double[] gradient(double x, double ... p) { - double[] v = { Double.POSITIVE_INFINITY, - Double.POSITIVE_INFINITY, - Double.POSITIVE_INFINITY }; - try { - v = super.gradient(x, p); - } catch (NotStrictlyPositiveException e) { // NOPMD - // Do nothing. - } - return v; - } - }; - - return fit(f, initialGuess); - } - - /** - * Fits a Gaussian function to the observed points. - * - * @return the parameters of the Gaussian function that best fits the - * observed points (in the same order as above). - */ - public double[] fit() { - final double[] guess = (new ParameterGuesser(getObservations())).guess(); - return fit(guess); - } - - /** - * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma} - * of a {@link org.apache.commons.math3.analysis.function.Gaussian.Parametric} - * based on the specified observed points. - */ - public static class ParameterGuesser { - /** Normalization factor. */ - private final double norm; - /** Mean. */ - private final double mean; - /** Standard deviation. */ - private final double sigma; - - /** - * Constructs instance with the specified observed points. - * - * @param observations Observed points from which to guess the - * parameters of the Gaussian. - * @throws NullArgumentException if {@code observations} is - * {@code null}. - * @throws NumberIsTooSmallException if there are less than 3 - * observations. - */ - public ParameterGuesser(WeightedObservedPoint[] observations) { - if (observations == null) { - throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY); - } - if (observations.length < 3) { - throw new NumberIsTooSmallException(observations.length, 3, true); - } - - final WeightedObservedPoint[] sorted = sortObservations(observations); - final double[] params = basicGuess(sorted); - - norm = params[0]; - mean = params[1]; - sigma = params[2]; - } - - /** - * Gets an estimation of the parameters. - * - * @return the guessed parameters, in the following order: - * <ul> - * <li>Normalization factor</li> - * <li>Mean</li> - * <li>Standard deviation</li> - * </ul> - */ - public double[] guess() { - return new double[] { norm, mean, sigma }; - } - - /** - * Sort the observations. - * - * @param unsorted Input observations. - * @return the input observations, sorted. - */ - private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) { - final WeightedObservedPoint[] observations = unsorted.clone(); - final Comparator<WeightedObservedPoint> cmp - = new Comparator<WeightedObservedPoint>() { - public int compare(WeightedObservedPoint p1, - WeightedObservedPoint p2) { - if (p1 == null && p2 == null) { - return 0; - } - if (p1 == null) { - return -1; - } - if (p2 == null) { - return 1; - } - if (p1.getX() < p2.getX()) { - return -1; - } - if (p1.getX() > p2.getX()) { - return 1; - } - if (p1.getY() < p2.getY()) { - return -1; - } - if (p1.getY() > p2.getY()) { - return 1; - } - if (p1.getWeight() < p2.getWeight()) { - return -1; - } - if (p1.getWeight() > p2.getWeight()) { - return 1; - } - return 0; - } - }; - - Arrays.sort(observations, cmp); - return observations; - } - - /** - * Guesses the parameters based on the specified observed points. - * - * @param points Observed points, sorted. - * @return the guessed parameters (normalization factor, mean and - * sigma). - */ - private double[] basicGuess(WeightedObservedPoint[] points) { - final int maxYIdx = findMaxY(points); - final double n = points[maxYIdx].getY(); - final double m = points[maxYIdx].getX(); - - double fwhmApprox; - try { - final double halfY = n + ((m - n) / 2); - final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY); - final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY); - fwhmApprox = fwhmX2 - fwhmX1; - } catch (OutOfRangeException e) { - // TODO: Exceptions should not be used for flow control. - fwhmApprox = points[points.length - 1].getX() - points[0].getX(); - } - final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2))); - - return new double[] { n, m, s }; - } - - /** - * Finds index of point in specified points with the largest Y. - * - * @param points Points to search. - * @return the index in specified points array. - */ - private int findMaxY(WeightedObservedPoint[] points) { - int maxYIdx = 0; - for (int i = 1; i < points.length; i++) { - if (points[i].getY() > points[maxYIdx].getY()) { - maxYIdx = i; - } - } - return maxYIdx; - } - - /** - * Interpolates using the specified points to determine X at the - * specified Y. - * - * @param points Points to use for interpolation. - * @param startIdx Index within points from which to start the search for - * interpolation bounds points. - * @param idxStep Index step for searching interpolation bounds points. - * @param y Y value for which X should be determined. - * @return the value of X for the specified Y. - * @throws ZeroException if {@code idxStep} is 0. - * @throws OutOfRangeException if specified {@code y} is not within the - * range of the specified {@code points}. - */ - private double interpolateXAtY(WeightedObservedPoint[] points, - int startIdx, - int idxStep, - double y) - throws OutOfRangeException { - if (idxStep == 0) { - throw new ZeroException(); - } - final WeightedObservedPoint[] twoPoints - = getInterpolationPointsForY(points, startIdx, idxStep, y); - final WeightedObservedPoint p1 = twoPoints[0]; - final WeightedObservedPoint p2 = twoPoints[1]; - if (p1.getY() == y) { - return p1.getX(); - } - if (p2.getY() == y) { - return p2.getX(); - } - return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) / - (p2.getY() - p1.getY())); - } - - /** - * Gets the two bounding interpolation points from the specified points - * suitable for determining X at the specified Y. - * - * @param points Points to use for interpolation. - * @param startIdx Index within points from which to start search for - * interpolation bounds points. - * @param idxStep Index step for search for interpolation bounds points. - * @param y Y value for which X should be determined. - * @return the array containing two points suitable for determining X at - * the specified Y. - * @throws ZeroException if {@code idxStep} is 0. - * @throws OutOfRangeException if specified {@code y} is not within the - * range of the specified {@code points}. - */ - private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points, - int startIdx, - int idxStep, - double y) - throws OutOfRangeException { - if (idxStep == 0) { - throw new ZeroException(); - } - for (int i = startIdx; - idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length; - i += idxStep) { - final WeightedObservedPoint p1 = points[i]; - final WeightedObservedPoint p2 = points[i + idxStep]; - if (isBetween(y, p1.getY(), p2.getY())) { - if (idxStep < 0) { - return new WeightedObservedPoint[] { p2, p1 }; - } else { - return new WeightedObservedPoint[] { p1, p2 }; - } - } - } - - // Boundaries are replaced by dummy values because the raised - // exception is caught and the message never displayed. - // TODO: Exceptions should not be used for flow control. - throw new OutOfRangeException(y, - Double.NEGATIVE_INFINITY, - Double.POSITIVE_INFINITY); - } - - /** - * Determines whether a value is between two other values. - * - * @param value Value to test whether it is between {@code boundary1} - * and {@code boundary2}. - * @param boundary1 One end of the range. - * @param boundary2 Other end of the range. - * @return {@code true} if {@code value} is between {@code boundary1} and - * {@code boundary2} (inclusive), {@code false} otherwise. - */ - private boolean isBetween(double value, - double boundary1, - double boundary2) { - return (value >= boundary1 && value <= boundary2) || - (value >= boundary2 && value <= boundary1); - } - } -}
