/*
 * 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.math3.optim.nonlinear.vector.jacobian;

import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.TooManyEvaluationsException;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.optim.OptimizationData;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.optim.PointVectorValuePair;
import org.apache.commons.math3.optim.nonlinear.vector.Weight;
import org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer;
import org.apache.commons.math3.util.FastMath;

Base class for implementing least-squares optimizers.
It provides methods for error estimation.
Since:
3.1
Deprecated:
All classes and interfaces in this package are deprecated. The optimizers that were provided here were moved to the leastsquares package (cf. MATH-1008)./**
 * Base class for implementing least-squares optimizers.
 * It provides methods for error estimation.
 *
 * @since 3.1
 * @deprecated All classes and interfaces in this package are deprecated.
 * The optimizers that were provided here were moved to the
 * {@link org.apache.commons.math3.fitting.leastsquares} package
 * (cf. MATH-1008).
 */
@Deprecated
public abstract class AbstractLeastSquaresOptimizer
    extends JacobianMultivariateVectorOptimizer {
    
Square-root of the weight matrix. /** Square-root of the weight matrix. */
    private RealMatrix weightMatrixSqrt;
    
Cost value (square root of the sum of the residuals). /** Cost value (square root of the sum of the residuals). */
    private double cost;

    
Params:
checker – Convergence checker./**
     * @param checker Convergence checker.
     */
    protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
        super(checker);
    }

    
Computes the weighted Jacobian matrix.
Params:
params – Model parameters at which to compute the Jacobian.
Throws:
DimensionMismatchException – if the Jacobian dimension does not
match problem dimension.
Returns:
the weighted Jacobian: W1/2 J./**
     * Computes the weighted Jacobian matrix.
     *
     * @param params Model parameters at which to compute the Jacobian.
     * @return the weighted Jacobian: W<sup>1/2</sup> J.
     * @throws DimensionMismatchException if the Jacobian dimension does not
     * match problem dimension.
     */
    protected RealMatrix computeWeightedJacobian(double[] params) {
        return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params)));
    }

    
Computes the cost.
Params:
residuals – Residuals.
See Also:
computeResiduals(double[])
Returns:
the cost./**
     * Computes the cost.
     *
     * @param residuals Residuals.
     * @return the cost.
     * @see #computeResiduals(double[])
     */
    protected double computeCost(double[] residuals) {
        final ArrayRealVector r = new ArrayRealVector(residuals);
        return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
    }

    
Gets the root-mean-square (RMS) value.
The RMS 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 c if the criterion, and n is the
number of measurements, then the RMS is sqrt (c/n).
Returns:
the RMS value./**
     * Gets the root-mean-square (RMS) value.
     *
     * The RMS 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 the RMS value.
     */
    public double getRMS() {
        return FastMath.sqrt(getChiSquare() / getTargetSize());
    }

    
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.
Returns:
chi-square value/**
     * 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;
    }

    
Gets the square-root of the weight matrix.
Returns:
the square-root of the weight matrix./**
     * Gets the square-root of the weight matrix.
     *
     * @return the square-root of the weight matrix.
     */
    public RealMatrix getWeightSquareRoot() {
        return weightMatrixSqrt.copy();
    }

    
Sets the cost.
Params:
cost – Cost value./**
     * Sets the cost.
     *
     * @param cost Cost value.
     */
    protected void setCost(double cost) {
        this.cost = cost;
    }

    
Get the covariance matrix of the optimized parameters.


Note that this operation involves the inversion of the
JTJ matrix, where J is the Jacobian matrix. The threshold parameter is a way for the caller to specify that the result of this computation should be considered meaningless, and thus trigger an exception. 
Params:
params – Model parameters.
threshold – Singularity threshold.
Throws:
SingularMatrixException – 
if the covariance matrix cannot be computed (singular problem).
Returns:
the covariance matrix./**
     * Get the covariance matrix of the optimized parameters.
     * <br/>
     * Note that this operation involves the inversion of the
     * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
     * Jacobian matrix.
     * The {@code threshold} parameter is a way for the caller to specify
     * that the result of this computation should be considered meaningless,
     * and thus trigger an exception.
     *
     * @param params Model parameters.
     * @param threshold Singularity threshold.
     * @return the covariance matrix.
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed (singular problem).
     */
    public double[][] computeCovariances(double[] params,
                                         double threshold) {
        // Set up the Jacobian.
        final RealMatrix j = computeWeightedJacobian(params);

        // Compute transpose(J)J.
        final RealMatrix jTj = j.transpose().multiply(j);

        // Compute the covariances matrix.
        final DecompositionSolver solver
            = new QRDecomposition(jTj, threshold).getSolver();
        return solver.getInverse().getData();
    }

    
Computes an estimate of the standard deviation of the parameters. The returned values are the square root of the diagonal coefficients of the covariance matrix, sd(a[i]) ~= sqrt(C[i][i]), where a[i] is the optimized value of the i-th parameter, and C is the covariance matrix. 
Params:
params – Model parameters.
covarianceSingularityThreshold – Singularity threshold (see computeCovariances).
Throws:
SingularMatrixException – 
if the covariance matrix cannot be computed.
Returns:
an estimate of the standard deviation of the optimized parameters/**
     * Computes an estimate of the standard deviation of the parameters. The
     * returned values are the square root of the diagonal coefficients of the
     * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
     * is the optimized value of the {@code i}-th parameter, and {@code C} is
     * the covariance matrix.
     *
     * @param params Model parameters.
     * @param covarianceSingularityThreshold Singularity threshold (see
     * {@link #computeCovariances(double[],double) computeCovariances}).
     * @return an estimate of the standard deviation of the optimized parameters
     * @throws org.apache.commons.math3.linear.SingularMatrixException
     * if the covariance matrix cannot be computed.
     */
    public double[] computeSigma(double[] params,
                                 double covarianceSingularityThreshold) {
        final int nC = params.length;
        final double[] sig = new double[nC];
        final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
        for (int i = 0; i < nC; ++i) {
            sig[i] = FastMath.sqrt(cov[i][i]);
        }
        return sig;
    }

    
{@inheritDoc}
Params:
optData – Optimization data. In addition to those documented in 
JacobianMultivariateVectorOptimizer, this method will register the following data: 

 
Weight
Throws:
TooManyEvaluationsException – if the maximal number of
evaluations is exceeded.
DimensionMismatchException – if the initial guess, target, and weight
arguments have inconsistent dimensions.
Returns:
{@inheritDoc}/**
     * {@inheritDoc}
     *
     * @param optData Optimization data. In addition to those documented in
     * {@link JacobianMultivariateVectorOptimizer#parseOptimizationData(OptimizationData[])
     * JacobianMultivariateVectorOptimizer}, this method will register the following data:
     * <ul>
     *  <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Weight}</li>
     * </ul>
     * @return {@inheritDoc}
     * @throws TooManyEvaluationsException if the maximal number of
     * evaluations is exceeded.
     * @throws DimensionMismatchException if the initial guess, target, and weight
     * arguments have inconsistent dimensions.
     */
    @Override
    public PointVectorValuePair optimize(OptimizationData... optData)
        throws TooManyEvaluationsException {
        // Set up base class and perform computation.
        return super.optimize(optData);
    }

    
Computes the residuals.
The residual is the difference between the observed (target)
values and the model (objective function) value.
There is one residual for each element of the vector-valued
function.
Params:
objectiveValue – Value of the the objective function. This is the value returned from a call to computeObjectiveValue (whose array argument contains the model parameters).
Throws:
DimensionMismatchException – if params has a wrong length.
Returns:
the residuals./**
     * Computes the residuals.
     * The residual is the difference between the observed (target)
     * values and the model (objective function) value.
     * There is one residual for each element of the vector-valued
     * function.
     *
     * @param objectiveValue Value of the the objective function. This is
     * the value returned from a call to
     * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
     * (whose array argument contains the model parameters).
     * @return the residuals.
     * @throws DimensionMismatchException if {@code params} has a wrong
     * length.
     */
    protected double[] computeResiduals(double[] objectiveValue) {
        final double[] target = getTarget();
        if (objectiveValue.length != target.length) {
            throw new DimensionMismatchException(target.length,
                                                 objectiveValue.length);
        }

        final double[] residuals = new double[target.length];
        for (int i = 0; i < target.length; i++) {
            residuals[i] = target[i] - objectiveValue[i];
        }

        return residuals;
    }

    
Scans the list of (required and optional) optimization data that characterize the problem. If the weight matrix is specified, the weightMatrixSqrt field is recomputed. 
Params:
optData – Optimization data. The following data will be looked for:

 
Weight
/**
     * Scans the list of (required and optional) optimization data that
     * characterize the problem.
     * If the weight matrix is specified, the {@link #weightMatrixSqrt}
     * field is recomputed.
     *
     * @param optData Optimization data. The following data will be looked for:
     * <ul>
     *  <li>{@link Weight}</li>
     * </ul>
     */
    @Override
    protected void parseOptimizationData(OptimizationData... optData) {
        // Allow base class to register its own data.
        super.parseOptimizationData(optData);

        // The existing values (as set by the previous call) are reused if
        // not provided in the argument list.
        for (OptimizationData data : optData) {
            if (data instanceof Weight) {
                weightMatrixSqrt = squareRoot(((Weight) data).getWeight());
                // If more data must be parsed, this statement _must_ be
                // changed to "continue".
                break;
            }
        }
    }

    
Computes the square-root of the weight matrix.
Params:
m – Symmetric, positive-definite (weight) matrix.
Returns:
the square-root of the weight matrix./**
     * Computes the square-root of the weight matrix.
     *
     * @param m Symmetric, positive-definite (weight) matrix.
     * @return the square-root of the weight matrix.
     */
    private RealMatrix squareRoot(RealMatrix m) {
        if (m instanceof DiagonalMatrix) {
            final int dim = m.getRowDimension();
            final RealMatrix sqrtM = new DiagonalMatrix(dim);
            for (int i = 0; i < dim; i++) {
                sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
            }
            return sqrtM;
        } else {
            final EigenDecomposition dec = new EigenDecomposition(m);
            return dec.getSquareRoot();
        }
    }
}