1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package org.apache.commons.math.estimation;
19
20import java.util.Arrays;
21
22import org.apache.commons.math.exception.util.LocalizedFormats;
23import org.apache.commons.math.linear.InvalidMatrixException;
24import org.apache.commons.math.linear.LUDecompositionImpl;
25import org.apache.commons.math.linear.MatrixUtils;
26import org.apache.commons.math.linear.RealMatrix;
27import org.apache.commons.math.util.FastMath;
28
29/**
30 * Base class for implementing estimators.
31 * <p>This base class handles the boilerplates methods associated to thresholds
32 * settings, jacobian and error estimation.</p>
33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
34 * @since 1.2
35 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
36 * been deprecated and replaced by package org.apache.commons.math.optimization.general
37 *
38 */
39@Deprecated
40public abstract class AbstractEstimator implements Estimator {
41
42    /** Default maximal number of cost evaluations allowed. */
43    public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
44
45    /** Array of measurements. */
46    protected WeightedMeasurement[] measurements;
47
48    /** Array of parameters. */
49    protected EstimatedParameter[] parameters;
50
51    /**
52     * Jacobian matrix.
53     * <p>This matrix is in canonical form just after the calls to
54     * {@link #updateJacobian()}, but may be modified by the solver
55     * in the derived class (the {@link LevenbergMarquardtEstimator
56     * Levenberg-Marquardt estimator} does this).</p>
57     */
58    protected double[] jacobian;
59
60    /** Number of columns of the jacobian matrix. */
61    protected int cols;
62
63    /** Number of rows of the jacobian matrix. */
64    protected int rows;
65
66    /** Residuals array.
67     * <p>This array is in canonical form just after the calls to
68     * {@link #updateJacobian()}, but may be modified by the solver
69     * in the derived class (the {@link LevenbergMarquardtEstimator
70     * Levenberg-Marquardt estimator} does this).</p>
71     */
72    protected double[] residuals;
73
74    /** Cost value (square root of the sum of the residuals). */
75    protected double cost;
76
77    /** Maximal allowed number of cost evaluations. */
78    private int maxCostEval;
79
80    /** Number of cost evaluations. */
81    private int costEvaluations;
82
83    /** Number of jacobian evaluations. */
84    private int jacobianEvaluations;
85
86    /**
87     * Build an abstract estimator for least squares problems.
88     * <p>The maximal number of cost evaluations allowed is set
89     * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p>
90     */
91    protected AbstractEstimator() {
92        setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
93    }
94
95    /**
96     * Set the maximal number of cost evaluations allowed.
97     *
98     * @param maxCostEval maximal number of cost evaluations allowed
99     * @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}
319