/*
 * 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.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;

/**
 * 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$ $Date$
 * @since 1.2
 *
 */
public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {

	/** Default maximal number of iterations allowed. */
	public static final int DEFAULT_MAX_ITERATIONS = 100;

	/** 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;

	/** 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;

	/** Objective function. */
	private DifferentiableMultivariateVectorialFunction f;

	/** Objective function derivatives. */
	private MultivariateMatrixFunction jF;

	/** Target value for the objective functions at optimum. */
	protected double[] target;

	/** Weight for the least squares cost computation. */
	protected double[] weights;

	/** Current point. */
	protected double[] point;

	/** Current objective function value. */
	protected double[] objective;

	/** Current residuals. */
	protected double[] residuals;

	/** Cost value (square root of the sum of the residuals). */
	protected double cost;

	/**
	 * 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 checker) {
		this.checker = checker;
	}

	/** {@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, "dimension mismatch {0} != {1}", jacobian.length, rows);
		}
		for (int i = 0; i < rows; i++) {
			final double[] ji = jacobian[i];
			final double factor = -Math.sqrt(weights[i]);
			for (int j = 0; j < cols; ++j) {
				ji[j] *= factor;
			}
		}
	}

	/**
	 * 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 = f.value(point);
		if (objective.length != rows) {
			throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}", objective.length, rows);
		}
		cost = 0;
		for (int i = 0, index = 0; i < rows; i++, index += cols) {
			final double residual = target[i] - objective[i];
			residuals[i] = residual;
			cost += weights[i] * residual * residual;
		}
		cost = Math.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() {
		double criterion = 0;
		for (int i = 0; i < rows; ++i) {
			final double residual = residuals[i];
			criterion += weights[i] * residual * residual;
		}
		return Math.sqrt(criterion / rows);
	}

	/**
	 * Get the Chi-Square value.
	 * 
	 * @return chi-square value
	 */
	public double getChiSquare() {
		double chiSquare = 0;
		for (int i = 0; i < rows; ++i) {
			final double residual = residuals[i];
			chiSquare += residual * residual / weights[i];
		}
		return chiSquare;
	}

	/**
	 * 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 += jacobian[k][i] * jacobian[k][j];
				}
				jTj[i][j] = sum;
				jTj[j][i] = sum;
			}
		}

		try {
			// compute the covariances matrix
			RealMatrix inverse = new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
			return inverse.getData();
		} catch (InvalidMatrixException ime) {
			throw new OptimizationException("unable to compute covariances: 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("no degrees of freedom ({0} measurements, {1} parameters)", rows, cols);
		}
		double[] errors = new double[cols];
		final double c = Math.sqrt(getChiSquare() / (rows - cols));
		double[][] covar = getCovariances();
		for (int i = 0; i < errors.length; ++i) {
			errors[i] = Math.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("dimension mismatch {0} != {1}", target.length, weights.length);
		}

		// reset counters
		iterations = 0;
		objectiveEvaluations = 0;
		jacobianEvaluations = 0;

		// store least squares problem characteristics
		this.f = f;
		jF = f.jacobian();
		this.target = target.clone();
		this.weights = 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];

		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
	 */
	abstract protected VectorialPointValuePair doOptimize()
			throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;

}