/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 * 
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.minimization;

import static com.opengamma.analytics.math.FunctionUtils.square;
import static com.opengamma.analytics.math.matrix.MatrixAlgebraFactory.OG_ALGEBRA;

import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.MathException;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;

/**
 * Standard version of Powell's method. It is intended to be used when an analytic function for the gradient is not available.
 * This implementation is taken from <i>"An Introduction to the Conjugate Gradient Method Without the Agonizing Pain", Shewchuk</i>.
 * 
 */
public class ConjugateDirectionVectorMinimizer implements Minimizer<Function1D<DoubleMatrix1D, Double>, DoubleMatrix1D> {

  private static final double SMALL = Double.MIN_NORMAL;
  private final double _eps;
  private final int _maxIterations;
  private final LineSearch _lineSearch;

  /**
   * Constructs the object with default values for tolerance (10<sup>-8</sup>) and number of iterations (100).
   * @param minimizer The minimizer, not null
   */
  public ConjugateDirectionVectorMinimizer(final ScalarMinimizer minimizer) {
    this(minimizer, 1e-8, 100);
  }

  /**
   * @param minimizer The minimizer, not null
   * @param tolerance The tolerance, must be greater than the minimum value of a double and less than one
   * @param maxIterations The maximum number of iterations
   */
  public ConjugateDirectionVectorMinimizer(final ScalarMinimizer minimizer, final double tolerance, final int maxIterations) {
    Validate.notNull(minimizer, "minimizer");
    Validate.isTrue(tolerance > SMALL && tolerance < 1, "Tolerance must be greater than " + SMALL + " and less than 1.0");
    Validate.isTrue(maxIterations >= 1, "Need at least one iteration");
    _lineSearch = new LineSearch(minimizer);
    _eps = tolerance;
    _maxIterations = maxIterations;
  }

  /**
   * {@inheritDoc}
   */
  @Override
  public DoubleMatrix1D minimize(final Function1D<DoubleMatrix1D, Double> function, final DoubleMatrix1D startPosition) {
    Validate.notNull(function, "function");
    Validate.notNull(startPosition, "start position");
    final int n = startPosition.getNumberOfElements();
    final DoubleMatrix1D[] directionSet = getDefaultDirectionSet(n);

    DoubleMatrix1D x0 = startPosition;
    for (int count = 0; count < _maxIterations; count++) {
      double delta = 0.0;
      int indexDelta = 0;
      final double startValue = function.evaluate(x0);
      double f1 = startValue;
      double f2 = 0;
      double lambda = 0.0;

      DoubleMatrix1D x = x0;
      for (int i = 0; i < n; i++) {
        final DoubleMatrix1D direction = directionSet[i];
        lambda = _lineSearch.minimise(function, direction, x);
        x = (DoubleMatrix1D) OG_ALGEBRA.add(x, OG_ALGEBRA.scale(direction, lambda));
        f2 = function.evaluate(x);
        final double temp = (f1 - f2); //TODO LineSearch should return this
        if (temp > delta) {
          delta = temp;
          indexDelta = i;
        }
        f1 = f2;
      }

      if ((startValue - f2) < _eps * (Math.abs(startValue) + Math.abs(f2)) / 2.0 + SMALL) {
        return x;
      }

      final DoubleMatrix1D deltaX = (DoubleMatrix1D) OG_ALGEBRA.subtract(x, x0);
      final DoubleMatrix1D extrapolatedPoint = (DoubleMatrix1D) OG_ALGEBRA.add(x, deltaX);

      final double extrapValue = function.evaluate(extrapolatedPoint);
      // Powell's condition for updating the direction set
      if (extrapValue < startValue && (2 * (startValue - 2 * f2 * extrapValue) * square(startValue - f2 - delta)) < (square(startValue - extrapValue) * delta)) {
        lambda = _lineSearch.minimise(function, deltaX, x);
        x = (DoubleMatrix1D) OG_ALGEBRA.add(x, OG_ALGEBRA.scale(deltaX, lambda));
        directionSet[indexDelta] = directionSet[n - 1];
        directionSet[n - 1] = deltaX;
      }

      x0 = x;
    }
    throw new MathException("ConjugateDirection failed to converge after " + _maxIterations + " iterations, with a tolerance of " + _eps + ". Final position reached was " + x0.toString());
  }

  DoubleMatrix1D[] getDefaultDirectionSet(final int dim) {
    final DoubleMatrix1D[] res = new DoubleMatrix1D[dim];
    for (int i = 0; i < dim; i++) {
      final double[] temp = new double[dim];
      temp[i] = 1.0;
      res[i] = new DoubleMatrix1D(temp);
    }
    return res;
  }

}