1 files changed, 318 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java b/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java
new file mode 100644
index 0000000..7a52b55
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java
@@ -0,0 +1,318 @@
+/*
+ * 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.math.estimation;
+
+import java.util.Arrays;
+
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.linear.InvalidMatrixException;
+import org.apache.commons.math.linear.LUDecompositionImpl;
+import org.apache.commons.math.linear.MatrixUtils;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * Base class for implementing estimators.
+ * <p>This base class handles the boilerplates methods associated to thresholds
+ * settings, jacobian and error estimation.</p>
+ * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
+ * @since 1.2
+ * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
+ * been deprecated and replaced by package org.apache.commons.math.optimization.general
+ *
+ */
+@Deprecated
+public abstract class AbstractEstimator implements Estimator {
+
+    /** Default maximal number of cost evaluations allowed. */
+    public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
+
+    /** Array of measurements. */
+    protected WeightedMeasurement[] measurements;
+
+    /** Array of parameters. */
+    protected EstimatedParameter[] parameters;
+
+    /**
+     * Jacobian matrix.
+     * <p>This matrix is in canonical form just after the calls to
+     * {@link #updateJacobian()}, but may be modified by the solver
+     * in the derived class (the {@link LevenbergMarquardtEstimator
+     * Levenberg-Marquardt estimator} does this).</p>
+     */
+    protected double[] jacobian;
+
+    /** Number of columns of the jacobian matrix. */
+    protected int cols;
+
+    /** Number of rows of the jacobian matrix. */
+    protected int rows;
+
+    /** Residuals array.
+     * <p>This array is in canonical form just after the calls to
+     * {@link #updateJacobian()}, but may be modified by the solver
+     * in the derived class (the {@link LevenbergMarquardtEstimator
+     * Levenberg-Marquardt estimator} does this).</p>
+     */
+    protected double[] residuals;
+
+    /** Cost value (square root of the sum of the residuals). */
+    protected double cost;
+
+    /** Maximal allowed number of cost evaluations. */
+    private int maxCostEval;
+
+    /** Number of cost evaluations. */
+    private int costEvaluations;
+
+    /** Number of jacobian evaluations. */
+    private int jacobianEvaluations;
+
+    /**
+     * Build an abstract estimator for least squares problems.
+     * <p>The maximal number of cost evaluations allowed is set
+     * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p>
+     */
+    protected AbstractEstimator() {
+        setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
+    }
+
+    /**
+     * Set the maximal number of cost evaluations allowed.
+     *
+     * @param maxCostEval maximal number of cost evaluations allowed
+     * @see #estimate
+     */
+    public final void setMaxCostEval(int maxCostEval) {
+        this.maxCostEval = maxCostEval;
+    }
+
+    /**
+     * Get the number of cost evaluations.
+     *
+     * @return number of cost evaluations
+     * */
+    public final int getCostEvaluations() {
+        return costEvaluations;
+    }
+
+    /**
+     * Get the number of jacobian evaluations.
+     *
+     * @return number of jacobian evaluations
+     * */
+    public final int getJacobianEvaluations() {
+        return jacobianEvaluations;
+    }
+
+    /**
+     * Update the jacobian matrix.
+     */
+    protected void updateJacobian() {
+        incrementJacobianEvaluationsCounter();
+        Arrays.fill(jacobian, 0);
+        int index = 0;
+        for (int i = 0; i < rows; i++) {
+            WeightedMeasurement wm = measurements[i];
+            double factor = -FastMath.sqrt(wm.getWeight());
+            for (int j = 0; j < cols; ++j) {
+                jacobian[index++] = factor * wm.getPartial(parameters[j]);
+            }
+        }
+    }
+
+    /**
+     * Increment the jacobian evaluations counter.
+     */
+    protected final void incrementJacobianEvaluationsCounter() {
+      ++jacobianEvaluations;
+    }
+
+    /**
+     * Update the residuals array and cost function value.
+     * @exception EstimationException if the number of cost evaluations
+     * exceeds the maximum allowed
+     */
+    protected void updateResidualsAndCost()
+    throws EstimationException {
+
+        if (++costEvaluations > maxCostEval) {
+            throw new EstimationException(LocalizedFormats.MAX_EVALUATIONS_EXCEEDED,
+                                          maxCostEval);
+        }
+
+        cost = 0;
+        int index = 0;
+        for (int i = 0; i < rows; i++, index += cols) {
+            WeightedMeasurement wm = measurements[i];
+            double residual = wm.getResidual();
+            residuals[i] = FastMath.sqrt(wm.getWeight()) * residual;
+            cost += wm.getWeight() * residual * residual;
+        }
+        cost = FastMath.sqrt(cost);
+
+    }
+
+    /**
+     * Get the Root Mean Square value.
+     * Get the Root Mean Square value, i.e. the root of the arithmetic
+     * mean of the square of all weighted residuals. This is related to the
+     * criterion that is minimized by the estimator as follows: if
+     * <em>c</em> if the criterion, and <em>n</em> is the number of
+     * measurements, then the RMS is <em>sqrt (c/n)</em>.
+     *
+     * @param problem estimation problem
+     * @return RMS value
+     */
+    public double getRMS(EstimationProblem problem) {
+        WeightedMeasurement[] wm = problem.getMeasurements();
+        double criterion = 0;
+        for (int i = 0; i < wm.length; ++i) {
+            double residual = wm[i].getResidual();
+            criterion += wm[i].getWeight() * residual * residual;
+        }
+        return FastMath.sqrt(criterion / wm.length);
+    }
+
+    /**
+     * Get the Chi-Square value.
+     * @param problem estimation problem
+     * @return chi-square value
+     */
+    public double getChiSquare(EstimationProblem problem) {
+        WeightedMeasurement[] wm = problem.getMeasurements();
+        double chiSquare = 0;
+        for (int i = 0; i < wm.length; ++i) {
+            double residual = wm[i].getResidual();
+            chiSquare += residual * residual / wm[i].getWeight();
+        }
+        return chiSquare;
+    }
+
+    /**
+     * Get the covariance matrix of unbound estimated parameters.
+     * @param problem estimation problem
+     * @return covariance matrix
+     * @exception EstimationException if the covariance matrix
+     * cannot be computed (singular problem)
+     */
+    public double[][] getCovariances(EstimationProblem problem)
+      throws EstimationException {
+
+        // set up the jacobian
+        updateJacobian();
+
+        // compute transpose(J).J, avoiding building big intermediate matrices
+        final int n = problem.getMeasurements().length;
+        final int m = problem.getUnboundParameters().length;
+        final int max  = m * n;
+        double[][] jTj = new double[m][m];
+        for (int i = 0; i < m; ++i) {
+            for (int j = i; j < m; ++j) {
+                double sum = 0;
+                for (int k = 0; k < max; k += m) {
+                    sum += jacobian[k + i] * jacobian[k + j];
+                }
+                jTj[i][j] = sum;
+                jTj[j][i] = sum;
+            }
+        }
+
+        try {
+            // compute the covariances matrix
+            RealMatrix inverse =
+                new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
+            return inverse.getData();
+        } catch (InvalidMatrixException ime) {
+            throw new EstimationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM);
+        }
+
+    }
+
+    /**
+     * Guess the errors in unbound estimated parameters.
+     * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
+     * @param problem estimation problem
+     * @return errors in estimated parameters
+     * @exception EstimationException if the covariances matrix cannot be computed
+     * or the number of degrees of freedom is not positive (number of measurements
+     * lesser or equal to number of parameters)
+     */
+    public double[] guessParametersErrors(EstimationProblem problem)
+      throws EstimationException {
+        int m = problem.getMeasurements().length;
+        int p = problem.getUnboundParameters().length;
+        if (m <= p) {
+            throw new EstimationException(
+                    LocalizedFormats.NO_DEGREES_OF_FREEDOM,
+                    m, p);
+        }
+        double[] errors = new double[problem.getUnboundParameters().length];
+        final double c = FastMath.sqrt(getChiSquare(problem) / (m - p));
+        double[][] covar = getCovariances(problem);
+        for (int i = 0; i < errors.length; ++i) {
+            errors[i] = FastMath.sqrt(covar[i][i]) * c;
+        }
+        return errors;
+    }
+
+    /**
+     * Initialization of the common parts of the estimation.
+     * <p>This method <em>must</em> be called at the start
+     * of the {@link #estimate(EstimationProblem) estimate}
+     * method.</p>
+     * @param problem estimation problem to solve
+     */
+    protected void initializeEstimate(EstimationProblem problem) {
+
+        // reset counters
+        costEvaluations     = 0;
+        jacobianEvaluations = 0;
+
+        // retrieve the equations and the parameters
+        measurements = problem.getMeasurements();
+        parameters   = problem.getUnboundParameters();
+
+        // arrays shared with the other private methods
+        rows      = measurements.length;
+        cols      = parameters.length;
+        jacobian  = new double[rows * cols];
+        residuals = new double[rows];
+
+        cost = Double.POSITIVE_INFINITY;
+
+    }
+
+    /**
+     * Solve an estimation problem.
+     *
+     * <p>The method should set the parameters of the problem to several
+     * trial values until it reaches convergence. If this method returns
+     * normally (i.e. without throwing an exception), then the best
+     * estimate of the parameters can be retrieved from the problem
+     * itself, through the {@link EstimationProblem#getAllParameters
+     * EstimationProblem.getAllParameters} method.</p>
+     *
+     * @param problem estimation problem to solve
+     * @exception EstimationException if the problem cannot be solved
+     *
+     */
+    public abstract void estimate(EstimationProblem problem)
+    throws EstimationException;
+
+}