1 files changed, 374 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
new file mode 100644
index 0000000..6271f53
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
@@ -0,0 +1,374 @@
+/*
+ * 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.optimization.general;
+
+import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.MaxEvaluationsExceededException;
+import org.apache.commons.math.MaxIterationsExceededException;
+import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
+import org.apache.commons.math.analysis.MultivariateMatrixFunction;
+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.optimization.OptimizationException;
+import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
+import org.apache.commons.math.optimization.VectorialConvergenceChecker;
+import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
+import org.apache.commons.math.optimization.VectorialPointValuePair;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * Base class for implementing least squares optimizers.
+ * <p>This base class handles the boilerplate methods associated to thresholds
+ * settings, jacobian and error estimation.</p>
+ * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
+ * @since 1.2
+ *
+ */
+public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
+
+    /** Default maximal number of iterations allowed. */
+    public static final int DEFAULT_MAX_ITERATIONS = 100;
+
+    /** Convergence checker. */
+    protected VectorialConvergenceChecker checker;
+
+    /**
+     * 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 LevenbergMarquardtOptimizer
+     * Levenberg-Marquardt optimizer} 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;
+
+    /**
+     * Target value for the objective functions at optimum.
+     * @since 2.1
+     */
+    protected double[] targetValues;
+
+    /**
+     * Weight for the least squares cost computation.
+     * @since 2.1
+     */
+    protected double[] residualsWeights;
+
+    /** Current point. */
+    protected double[] point;
+
+    /** Current objective function value. */
+    protected double[] objective;
+
+    /** Current residuals. */
+    protected double[] residuals;
+
+    /** Weighted Jacobian */
+    protected double[][] wjacobian;
+
+    /** Weighted residuals */
+    protected double[] wresiduals;
+
+    /** Cost value (square root of the sum of the residuals). */
+    protected double cost;
+
+    /** Maximal number of iterations allowed. */
+    private int maxIterations;
+
+    /** Number of iterations already performed. */
+    private int iterations;
+
+    /** Maximal number of evaluations allowed. */
+    private int maxEvaluations;
+
+    /** Number of evaluations already performed. */
+    private int objectiveEvaluations;
+
+    /** Number of jacobian evaluations. */
+    private int jacobianEvaluations;
+
+    /** Objective function. */
+    private DifferentiableMultivariateVectorialFunction function;
+
+    /** Objective function derivatives. */
+    private MultivariateMatrixFunction jF;
+
+    /** Simple constructor with default settings.
+     * <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
+     * and the maximal number of evaluation is set to its default value.</p>
+     */
+    protected AbstractLeastSquaresOptimizer() {
+        setConvergenceChecker(new SimpleVectorialValueChecker());
+        setMaxIterations(DEFAULT_MAX_ITERATIONS);
+        setMaxEvaluations(Integer.MAX_VALUE);
+    }
+
+    /** {@inheritDoc} */
+    public void setMaxIterations(int maxIterations) {
+        this.maxIterations = maxIterations;
+    }
+
+    /** {@inheritDoc} */
+    public int getMaxIterations() {
+        return maxIterations;
+    }
+
+    /** {@inheritDoc} */
+    public int getIterations() {
+        return iterations;
+    }
+
+    /** {@inheritDoc} */
+    public void setMaxEvaluations(int maxEvaluations) {
+        this.maxEvaluations = maxEvaluations;
+    }
+
+    /** {@inheritDoc} */
+    public int getMaxEvaluations() {
+        return maxEvaluations;
+    }
+
+    /** {@inheritDoc} */
+    public int getEvaluations() {
+        return objectiveEvaluations;
+    }
+
+    /** {@inheritDoc} */
+    public int getJacobianEvaluations() {
+        return jacobianEvaluations;
+    }
+
+    /** {@inheritDoc} */
+    public void setConvergenceChecker(VectorialConvergenceChecker convergenceChecker) {
+        this.checker = convergenceChecker;
+    }
+
+    /** {@inheritDoc} */
+    public VectorialConvergenceChecker getConvergenceChecker() {
+        return checker;
+    }
+
+    /** Increment the iterations counter by 1.
+     * @exception OptimizationException if the maximal number
+     * of iterations is exceeded
+     */
+    protected void incrementIterationsCounter()
+        throws OptimizationException {
+        if (++iterations > maxIterations) {
+            throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
+        }
+    }
+
+    /**
+     * Update the jacobian matrix.
+     * @exception FunctionEvaluationException if the function jacobian
+     * cannot be evaluated or its dimension doesn't match problem dimension
+     */
+    protected void updateJacobian() throws FunctionEvaluationException {
+        ++jacobianEvaluations;
+        jacobian = jF.value(point);
+        if (jacobian.length != rows) {
+            throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
+                                                  jacobian.length, rows);
+        }
+        for (int i = 0; i < rows; i++) {
+            final double[] ji = jacobian[i];
+            double wi = FastMath.sqrt(residualsWeights[i]);
+            for (int j = 0; j < cols; ++j) {
+                ji[j] *=  -1.0;
+                wjacobian[i][j] = ji[j]*wi;
+            }
+        }
+    }
+
+    /**
+     * Update the residuals array and cost function value.
+     * @exception FunctionEvaluationException if the function cannot be evaluated
+     * or its dimension doesn't match problem dimension or maximal number of
+     * of evaluations is exceeded
+     */
+    protected void updateResidualsAndCost()
+        throws FunctionEvaluationException {
+
+        if (++objectiveEvaluations > maxEvaluations) {
+            throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations),
+                                                  point);
+        }
+        objective = function.value(point);
+        if (objective.length != rows) {
+            throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
+                                                  objective.length, rows);
+        }
+        cost = 0;
+        int index = 0;
+        for (int i = 0; i < rows; i++) {
+            final double residual = targetValues[i] - objective[i];
+            residuals[i] = residual;
+            wresiduals[i]= residual*FastMath.sqrt(residualsWeights[i]);
+            cost += residualsWeights[i] * residual * residual;
+            index += cols;
+        }
+        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 optimizer 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>.
+     *
+     * @return RMS value
+     */
+    public double getRMS() {
+        return FastMath.sqrt(getChiSquare() / rows);
+    }
+
+    /**
+     * Get a Chi-Square-like value assuming the N residuals follow N
+     * distinct normal distributions centered on 0 and whose variances are
+     * the reciprocal of the weights.
+     * @return chi-square value
+     */
+    public double getChiSquare() {
+        return cost*cost;
+    }
+
+    /**
+     * Get the covariance matrix of optimized parameters.
+     * @return covariance matrix
+     * @exception FunctionEvaluationException if the function jacobian cannot
+     * be evaluated
+     * @exception OptimizationException if the covariance matrix
+     * cannot be computed (singular problem)
+     */
+    public double[][] getCovariances()
+        throws FunctionEvaluationException, OptimizationException {
+
+        // set up the jacobian
+        updateJacobian();
+
+        // compute transpose(J).J, avoiding building big intermediate matrices
+        double[][] jTj = new double[cols][cols];
+        for (int i = 0; i < cols; ++i) {
+            for (int j = i; j < cols; ++j) {
+                double sum = 0;
+                for (int k = 0; k < rows; ++k) {
+                    sum += wjacobian[k][i] * wjacobian[k][j];
+                }
+                jTj[i][j] = sum;
+                jTj[j][i] = sum;
+            }
+        }
+
+        try {
+            // compute the covariance matrix
+            RealMatrix inverse =
+                new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
+            return inverse.getData();
+        } catch (InvalidMatrixException ime) {
+            throw new OptimizationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM);
+        }
+
+    }
+
+    /**
+     * Guess the errors in optimized parameters.
+     * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
+     * @return errors in optimized parameters
+     * @exception FunctionEvaluationException if the function jacobian cannot b evaluated
+     * @exception OptimizationException 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()
+        throws FunctionEvaluationException, OptimizationException {
+        if (rows <= cols) {
+            throw new OptimizationException(
+                    LocalizedFormats.NO_DEGREES_OF_FREEDOM,
+                    rows, cols);
+        }
+        double[] errors = new double[cols];
+        final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
+        double[][] covar = getCovariances();
+        for (int i = 0; i < errors.length; ++i) {
+            errors[i] = FastMath.sqrt(covar[i][i]) * c;
+        }
+        return errors;
+    }
+
+    /** {@inheritDoc} */
+    public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
+                                            final double[] target, final double[] weights,
+                                            final double[] startPoint)
+        throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
+
+        if (target.length != weights.length) {
+            throw new OptimizationException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
+                                            target.length, weights.length);
+        }
+
+        // reset counters
+        iterations           = 0;
+        objectiveEvaluations = 0;
+        jacobianEvaluations  = 0;
+
+        // store least squares problem characteristics
+        function         = f;
+        jF               = f.jacobian();
+        targetValues     = target.clone();
+        residualsWeights = weights.clone();
+        this.point       = startPoint.clone();
+        this.residuals   = new double[target.length];
+
+        // arrays shared with the other private methods
+        rows      = target.length;
+        cols      = point.length;
+        jacobian  = new double[rows][cols];
+
+        wjacobian = new double[rows][cols];
+        wresiduals = new double[rows];
+
+        cost = Double.POSITIVE_INFINITY;
+
+        return doOptimize();
+
+    }
+
+    /** Perform the bulk of optimization algorithm.
+     * @return the point/value pair giving the optimal value for objective function
+     * @exception FunctionEvaluationException if the objective function throws one during
+     * the search
+     * @exception OptimizationException if the algorithm failed to converge
+     * @exception IllegalArgumentException if the start point dimension is wrong
+     */
+    protected abstract VectorialPointValuePair doOptimize()
+        throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
+
+
+}