You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2015/02/16 23:39:58 UTC
[28/82] [partial] [math] Update for next development iteration:
commons-math4
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);
- }
- }
-}