package org.spaceroots.mantissa.optimization;

/** This class implements the Nelder-Mead direct search method.

 * @author Luc Maisonobe
 * @version $Id: NelderMead.java 1709 2006-12-03 21:16:50Z luc $
 * @see MultiDirectional
 */
public class NelderMead
  extends DirectSearchOptimizer {

  /** Build a Nelder-Mead optimizer with default coefficients.
   * 
The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
   * for both gamma and sigma.
   */
  public NelderMead() {
    super();
    this.rho   = 1.0;
    this.khi   = 2.0;
    this.gamma = 0.5;
    this.sigma = 0.5;
  }

  /** Build a Nelder-Mead optimizer with specified coefficients.
   * @param rho reflection coefficient
   * @param khi expansion coefficient
   * @param gamma contraction coefficient
   * @param sigma shrinkage coefficient
   */
  public NelderMead(double rho, double khi, double gamma, double sigma) {
    super();
    this.rho   = rho;
    this.khi   = khi;
    this.gamma = gamma;
    this.sigma = sigma;
  }

  /** Compute the next simplex of the algorithm.
   */
  protected void iterateSimplex()
    throws CostException {

    // the simplex has n+1 point if dimension is n
    int n = simplex.length - 1;

    // interesting costs
    double   smallest      = simplex[0].cost;
    double   secondLargest = simplex[n-1].cost;
    double   largest       = simplex[n].cost;
    double[] xLargest      = simplex[n].point;

    // compute the centroid of the best vertices
    // (dismissing the worst point at index n)
    double[] centroid = new double[n];
    for (int i = 0; i < n; ++i) {
      double[] x = simplex[i].point;
      for (int j = 0; j < n; ++j) {
        centroid[j] += x[j];
      }
    }
    double scaling = 1.0 / n;
    for (int j = 0; j < n; ++j) {
      centroid[j] *= scaling;
    }

    // compute the reflection point
    double[] xR       = new double[n];
    for (int j = 0; j < n; ++j) {
      xR[j] = centroid[j] + rho * (centroid[j] - xLargest[j]);
    }
    double costR = evaluateCost(xR);

    if ((smallest <= costR) && (costR < secondLargest)) {

      // accept the reflected point
      replaceWorstPoint(new PointCostPair(xR, costR));

    } else if (costR < smallest) {

      // compute the expansion point
      double[] xE = new double[n];
      for (int j = 0; j < n; ++j) {
        xE[j] = centroid[j] + khi * (xR[j] - centroid[j]);
      }
      double costE = evaluateCost(xE);

      if (costE < costR) {
        // accept the expansion point
        replaceWorstPoint(new PointCostPair(xE, costE));
      } else {
        // accept the reflected point
        replaceWorstPoint(new PointCostPair(xR, costR));
      }

    } else {

      if (costR < largest) {

        // perform an outside contraction
        double[] xC = new double[n];
        for (int j = 0; j < n; ++j) {
          xC[j] = centroid[j] + gamma * (xR[j] - centroid[j]);
        }
        double costC = evaluateCost(xC);

        if (costC <= costR) {
          // accept the contraction point
          replaceWorstPoint(new PointCostPair(xC, costC));
          return;
        }

      } else {

        // perform an inside contraction
        double[] xC = new double[n];
        for (int j = 0; j < n; ++j) {
          xC[j] = centroid[j] - gamma * (centroid[j] - xLargest[j]);
        }
        double costC = evaluateCost(xC);

        if (costC < largest) {
          // accept the contraction point
          replaceWorstPoint(new PointCostPair(xC, costC));
          return;
        }

      }

      // perform a shrink
      double[] xSmallest = simplex[0].point;
      for (int i = 1; i < simplex.length; ++i) {
        double[] x = simplex[i].point;
        for (int j = 0; j < n; ++j) {
          x[j] = xSmallest[j] + sigma * (x[j] - xSmallest[j]);
        }
        simplex[i] = new PointCostPair(x, Double.NaN);
      }
      evaluateSimplex();

    }

  }

  /** Reflection coefficient. */
  private double rho;

  /** Expansion coefficient. */
  private double khi;

  /** Contraction coefficient. */
  private double gamma;

  /** Shrinkage coefficient. */
  private double sigma;

}