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