    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math.estimation;
    019
    020    import java.util.Arrays;
    021
    022    import org.apache.commons.math.exception.util.LocalizedFormats;
    023    import org.apache.commons.math.linear.InvalidMatrixException;
    024    import org.apache.commons.math.linear.LUDecompositionImpl;
    025    import org.apache.commons.math.linear.MatrixUtils;
    026    import org.apache.commons.math.linear.RealMatrix;
    027    import org.apache.commons.math.util.FastMath;
    028
    029    /**
    030     * Base class for implementing estimators.
    031     * <p>This base class handles the boilerplates methods associated to thresholds
    032     * settings, jacobian and error estimation.</p>
    033     * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
    034     * @since 1.2
    035     * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
    036     * been deprecated and replaced by package org.apache.commons.math.optimization.general
    037     *
    038     */
    039    @Deprecated
    040    public abstract class AbstractEstimator implements Estimator {
    041
    042        /** Default maximal number of cost evaluations allowed. */
    043        public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
    044
    045        /** Array of measurements. */
    046        protected WeightedMeasurement[] measurements;
    047
    048        /** Array of parameters. */
    049        protected EstimatedParameter[] parameters;
    050
    051        /**
    052         * Jacobian matrix.
    053         * <p>This matrix is in canonical form just after the calls to
    054         * {@link #updateJacobian()}, but may be modified by the solver
    055         * in the derived class (the {@link LevenbergMarquardtEstimator
    056         * Levenberg-Marquardt estimator} does this).</p>
    057         */
    058        protected double[] jacobian;
    059
    060        /** Number of columns of the jacobian matrix. */
    061        protected int cols;
    062
    063        /** Number of rows of the jacobian matrix. */
    064        protected int rows;
    065
    066        /** Residuals array.
    067         * <p>This array is in canonical form just after the calls to
    068         * {@link #updateJacobian()}, but may be modified by the solver
    069         * in the derived class (the {@link LevenbergMarquardtEstimator
    070         * Levenberg-Marquardt estimator} does this).</p>
    071         */
    072        protected double[] residuals;
    073
    074        /** Cost value (square root of the sum of the residuals). */
    075        protected double cost;
    076
    077        /** Maximal allowed number of cost evaluations. */
    078        private int maxCostEval;
    079
    080        /** Number of cost evaluations. */
    081        private int costEvaluations;
    082
    083        /** Number of jacobian evaluations. */
    084        private int jacobianEvaluations;
    085
    086        /**
    087         * Build an abstract estimator for least squares problems.
    088         * <p>The maximal number of cost evaluations allowed is set
    089         * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p>
    090         */
    091        protected AbstractEstimator() {
    092            setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
    093        }
    094
    095        /**
    096         * Set the maximal number of cost evaluations allowed.
    097         *
    098         * @param maxCostEval maximal number of cost evaluations allowed
    099         * @see #estimate
    100         */
    101        public final void setMaxCostEval(int maxCostEval) {
    102            this.maxCostEval = maxCostEval;
    103        }
    104
    105        /**
    106         * Get the number of cost evaluations.
    107         *
    108         * @return number of cost evaluations
    109         * */
    110        public final int getCostEvaluations() {
    111            return costEvaluations;
    112        }
    113
    114        /**
    115         * Get the number of jacobian evaluations.
    116         *
    117         * @return number of jacobian evaluations
    118         * */
    119        public final int getJacobianEvaluations() {
    120            return jacobianEvaluations;
    121        }
    122
    123        /**
    124         * Update the jacobian matrix.
    125         */
    126        protected void updateJacobian() {
    127            incrementJacobianEvaluationsCounter();
    128            Arrays.fill(jacobian, 0);
    129            int index = 0;
    130            for (int i = 0; i < rows; i++) {
    131                WeightedMeasurement wm = measurements[i];
    132                double factor = -FastMath.sqrt(wm.getWeight());
    133                for (int j = 0; j < cols; ++j) {
    134                    jacobian[index++] = factor * wm.getPartial(parameters[j]);
    135                }
    136            }
    137        }
    138
    139        /**
    140         * Increment the jacobian evaluations counter.
    141         */
    142        protected final void incrementJacobianEvaluationsCounter() {
    143          ++jacobianEvaluations;
    144        }
    145
    146        /**
    147         * Update the residuals array and cost function value.
    148         * @exception EstimationException if the number of cost evaluations
    149         * exceeds the maximum allowed
    150         */
    151        protected void updateResidualsAndCost()
    152        throws EstimationException {
    153
    154            if (++costEvaluations > maxCostEval) {
    155                throw new EstimationException(LocalizedFormats.MAX_EVALUATIONS_EXCEEDED,
    156                                              maxCostEval);
    157            }
    158
    159            cost = 0;
    160            int index = 0;
    161            for (int i = 0; i < rows; i++, index += cols) {
    162                WeightedMeasurement wm = measurements[i];
    163                double residual = wm.getResidual();
    164                residuals[i] = FastMath.sqrt(wm.getWeight()) * residual;
    165                cost += wm.getWeight() * residual * residual;
    166            }
    167            cost = FastMath.sqrt(cost);
    168
    169        }
    170
    171        /**
    172         * Get the Root Mean Square value.
    173         * Get the Root Mean Square value, i.e. the root of the arithmetic
    174         * mean of the square of all weighted residuals. This is related to the
    175         * criterion that is minimized by the estimator as follows: if
    176         * <em>c</em> if the criterion, and <em>n</em> is the number of
    177         * measurements, then the RMS is <em>sqrt (c/n)</em>.
    178         *
    179         * @param problem estimation problem
    180         * @return RMS value
    181         */
    182        public double getRMS(EstimationProblem problem) {
    183            WeightedMeasurement[] wm = problem.getMeasurements();
    184            double criterion = 0;
    185            for (int i = 0; i < wm.length; ++i) {
    186                double residual = wm[i].getResidual();
    187                criterion += wm[i].getWeight() * residual * residual;
    188            }
    189            return FastMath.sqrt(criterion / wm.length);
    190        }
    191
    192        /**
    193         * Get the Chi-Square value.
    194         * @param problem estimation problem
    195         * @return chi-square value
    196         */
    197        public double getChiSquare(EstimationProblem problem) {
    198            WeightedMeasurement[] wm = problem.getMeasurements();
    199            double chiSquare = 0;
    200            for (int i = 0; i < wm.length; ++i) {
    201                double residual = wm[i].getResidual();
    202                chiSquare += residual * residual / wm[i].getWeight();
    203            }
    204            return chiSquare;
    205        }
    206
    207        /**
    208         * Get the covariance matrix of unbound estimated parameters.
    209         * @param problem estimation problem
    210         * @return covariance matrix
    211         * @exception EstimationException if the covariance matrix
    212         * cannot be computed (singular problem)
    213         */
    214        public double[][] getCovariances(EstimationProblem problem)
    215          throws EstimationException {
    216
    217            // set up the jacobian
    218            updateJacobian();
    219
    220            // compute transpose(J).J, avoiding building big intermediate matrices
    221            final int n = problem.getMeasurements().length;
    222            final int m = problem.getUnboundParameters().length;
    223            final int max  = m * n;
    224            double[][] jTj = new double[m][m];
    225            for (int i = 0; i < m; ++i) {
    226                for (int j = i; j < m; ++j) {
    227                    double sum = 0;
    228                    for (int k = 0; k < max; k += m) {
    229                        sum += jacobian[k + i] * jacobian[k + j];
    230                    }
    231                    jTj[i][j] = sum;
    232                    jTj[j][i] = sum;
    233                }
    234            }
    235
    236            try {
    237                // compute the covariances matrix
    238                RealMatrix inverse =
    239                    new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
    240                return inverse.getData();
    241            } catch (InvalidMatrixException ime) {
    242                throw new EstimationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM);
    243            }
    244
    245        }
    246
    247        /**
    248         * Guess the errors in unbound estimated parameters.
    249         * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
    250         * @param problem estimation problem
    251         * @return errors in estimated parameters
    252         * @exception EstimationException if the covariances matrix cannot be computed
    253         * or the number of degrees of freedom is not positive (number of measurements
    254         * lesser or equal to number of parameters)
    255         */
    256        public double[] guessParametersErrors(EstimationProblem problem)
    257          throws EstimationException {
    258            int m = problem.getMeasurements().length;
    259            int p = problem.getUnboundParameters().length;
    260            if (m <= p) {
    261                throw new EstimationException(
    262                        LocalizedFormats.NO_DEGREES_OF_FREEDOM,
    263                        m, p);
    264            }
    265            double[] errors = new double[problem.getUnboundParameters().length];
    266            final double c = FastMath.sqrt(getChiSquare(problem) / (m - p));
    267            double[][] covar = getCovariances(problem);
    268            for (int i = 0; i < errors.length; ++i) {
    269                errors[i] = FastMath.sqrt(covar[i][i]) * c;
    270            }
    271            return errors;
    272        }
    273
    274        /**
    275         * Initialization of the common parts of the estimation.
    276         * <p>This method <em>must</em> be called at the start
    277         * of the {@link #estimate(EstimationProblem) estimate}
    278         * method.</p>
    279         * @param problem estimation problem to solve
    280         */
    281        protected void initializeEstimate(EstimationProblem problem) {
    282
    283            // reset counters
    284            costEvaluations     = 0;
    285            jacobianEvaluations = 0;
    286
    287            // retrieve the equations and the parameters
    288            measurements = problem.getMeasurements();
    289            parameters   = problem.getUnboundParameters();
    290
    291            // arrays shared with the other private methods
    292            rows      = measurements.length;
    293            cols      = parameters.length;
    294            jacobian  = new double[rows * cols];
    295            residuals = new double[rows];
    296
    297            cost = Double.POSITIVE_INFINITY;
    298
    299        }
    300
    301        /**
    302         * Solve an estimation problem.
    303         *
    304         * <p>The method should set the parameters of the problem to several
    305         * trial values until it reaches convergence. If this method returns
    306         * normally (i.e. without throwing an exception), then the best
    307         * estimate of the parameters can be retrieved from the problem
    308         * itself, through the {@link EstimationProblem#getAllParameters
    309         * EstimationProblem.getAllParameters} method.</p>
    310         *
    311         * @param problem estimation problem to solve
    312         * @exception EstimationException if the problem cannot be solved
    313         *
    314         */
    315        public abstract void estimate(EstimationProblem problem)
    316        throws EstimationException;
    317
    318    }