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

import java.util.Comparator;

import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.analysis.MultivariateFunction;

/**
 * This class implements the Nelder-Mead simplex algorithm.
 *
 * @since 3.0
 */
public class NelderMeadSimplex extends AbstractSimplex {
    /** Default value for {@link #rho}: {@value}. */
    private static final double DEFAULT_RHO = 1;
    /** Default value for {@link #khi}: {@value}. */
    private static final double DEFAULT_KHI = 2;
    /** Default value for {@link #gamma}: {@value}. */
    private static final double DEFAULT_GAMMA = 0.5;
    /** Default value for {@link #sigma}: {@value}. */
    private static final double DEFAULT_SIGMA = 0.5;
    /** Reflection coefficient. */
    private final double rho;
    /** Expansion coefficient. */
    private final double khi;
    /** Contraction coefficient. */
    private final double gamma;
    /** Shrinkage coefficient. */
    private final double sigma;

    /**
     * Build a Nelder-Mead simplex with default coefficients.
     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
     * for both gamma and sigma.
     *
     * @param n Dimension of the simplex.
     */
    public NelderMeadSimplex(final int n) {
        this(n, 1d);
    }

    /**
     * Build a Nelder-Mead simplex with default coefficients.
     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
     * for both gamma and sigma.
     *
     * @param n Dimension of the simplex.
     * @param sideLength Length of the sides of the default (hypercube)
     * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
     */
    public NelderMeadSimplex(final int n, double sideLength) {
        this(n, sideLength,
             DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
    }

    /**
     * Build a Nelder-Mead simplex with specified coefficients.
     *
     * @param n Dimension of the simplex. See
     * {@link AbstractSimplex#AbstractSimplex(int,double)}.
     * @param sideLength Length of the sides of the default (hypercube)
     * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
     * @param rho Reflection coefficient.
     * @param khi Expansion coefficient.
     * @param gamma Contraction coefficient.
     * @param sigma Shrinkage coefficient.
     */
    public NelderMeadSimplex(final int n, double sideLength,
                             final double rho, final double khi,
                             final double gamma, final double sigma) {
        super(n, sideLength);

        this.rho = rho;
        this.khi = khi;
        this.gamma = gamma;
        this.sigma = sigma;
    }

    /**
     * Build a Nelder-Mead simplex with specified coefficients.
     *
     * @param n Dimension of the simplex. See
     * {@link AbstractSimplex#AbstractSimplex(int)}.
     * @param rho Reflection coefficient.
     * @param khi Expansion coefficient.
     * @param gamma Contraction coefficient.
     * @param sigma Shrinkage coefficient.
     */
    public NelderMeadSimplex(final int n,
                             final double rho, final double khi,
                             final double gamma, final double sigma) {
        this(n, 1d, rho, khi, gamma, sigma);
    }

    /**
     * Build a Nelder-Mead simplex with default coefficients.
     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
     * for both gamma and sigma.
     *
     * @param steps Steps along the canonical axes representing box edges.
     * They may be negative but not zero. See
     */
    public NelderMeadSimplex(final double[] steps) {
        this(steps, DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
    }

    /**
     * Build a Nelder-Mead simplex with specified coefficients.
     *
     * @param steps Steps along the canonical axes representing box edges.
     * They may be negative but not zero. See
     * {@link AbstractSimplex#AbstractSimplex(double[])}.
     * @param rho Reflection coefficient.
     * @param khi Expansion coefficient.
     * @param gamma Contraction coefficient.
     * @param sigma Shrinkage coefficient.
     * @throws IllegalArgumentException if one of the steps is zero.
     */
    public NelderMeadSimplex(final double[] steps,
                             final double rho, final double khi,
                             final double gamma, final double sigma) {
        super(steps);

        this.rho = rho;
        this.khi = khi;
        this.gamma = gamma;
        this.sigma = sigma;
    }

    /**
     * Build a Nelder-Mead simplex with default coefficients.
     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
     * for both gamma and sigma.
     *
     * @param referenceSimplex Reference simplex. See
     * {@link AbstractSimplex#AbstractSimplex(double[][])}.
     */
    public NelderMeadSimplex(final double[][] referenceSimplex) {
        this(referenceSimplex, DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
    }

    /**
     * Build a Nelder-Mead simplex with specified coefficients.
     *
     * @param referenceSimplex Reference simplex. See
     * {@link AbstractSimplex#AbstractSimplex(double[][])}.
     * @param rho Reflection coefficient.
     * @param khi Expansion coefficient.
     * @param gamma Contraction coefficient.
     * @param sigma Shrinkage coefficient.
     * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
     * if the reference simplex does not contain at least one point.
     * @throws org.apache.commons.math3.exception.DimensionMismatchException
     * if there is a dimension mismatch in the reference simplex.
     */
    public NelderMeadSimplex(final double[][] referenceSimplex,
                             final double rho, final double khi,
                             final double gamma, final double sigma) {
        super(referenceSimplex);

        this.rho = rho;
        this.khi = khi;
        this.gamma = gamma;
        this.sigma = sigma;
    }

    /** {@inheritDoc} */
    @Override
    public void iterate(final MultivariateFunction evaluationFunction,
                        final Comparator<PointValuePair> comparator) {
        // The simplex has n + 1 points if dimension is n.
        final int n = getDimension();

        // Interesting values.
        final PointValuePair best = getPoint(0);
        final PointValuePair secondBest = getPoint(n - 1);
        final PointValuePair worst = getPoint(n);
        final double[] xWorst = worst.getPointRef();

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

        // compute the reflection point
        final double[] xR = new double[n];
        for (int j = 0; j < n; j++) {
            xR[j] = centroid[j] + rho * (centroid[j] - xWorst[j]);
        }
        final PointValuePair reflected
            = new PointValuePair(xR, evaluationFunction.value(xR), false);

        if (comparator.compare(best, reflected) <= 0 &&
            comparator.compare(reflected, secondBest) < 0) {
            // Accept the reflected point.
            replaceWorstPoint(reflected, comparator);
        } else if (comparator.compare(reflected, best) < 0) {
            // Compute the expansion point.
            final double[] xE = new double[n];
            for (int j = 0; j < n; j++) {
                xE[j] = centroid[j] + khi * (xR[j] - centroid[j]);
            }
            final PointValuePair expanded
                = new PointValuePair(xE, evaluationFunction.value(xE), false);

            if (comparator.compare(expanded, reflected) < 0) {
                // Accept the expansion point.
                replaceWorstPoint(expanded, comparator);
            } else {
                // Accept the reflected point.
                replaceWorstPoint(reflected, comparator);
            }
        } else {
            if (comparator.compare(reflected, worst) < 0) {
                // Perform an outside contraction.
                final double[] xC = new double[n];
                for (int j = 0; j < n; j++) {
                    xC[j] = centroid[j] + gamma * (xR[j] - centroid[j]);
                }
                final PointValuePair outContracted
                    = new PointValuePair(xC, evaluationFunction.value(xC), false);
                if (comparator.compare(outContracted, reflected) <= 0) {
                    // Accept the contraction point.
                    replaceWorstPoint(outContracted, comparator);
                    return;
                }
            } else {
                // Perform an inside contraction.
                final double[] xC = new double[n];
                for (int j = 0; j < n; j++) {
                    xC[j] = centroid[j] - gamma * (centroid[j] - xWorst[j]);
                }
                final PointValuePair inContracted
                    = new PointValuePair(xC, evaluationFunction.value(xC), false);

                if (comparator.compare(inContracted, worst) < 0) {
                    // Accept the contraction point.
                    replaceWorstPoint(inContracted, comparator);
                    return;
                }
            }

            // Perform a shrink.
            final double[] xSmallest = getPoint(0).getPointRef();
            for (int i = 1; i <= n; i++) {
                final double[] x = getPoint(i).getPoint();
                for (int j = 0; j < n; j++) {
                    x[j] = xSmallest[j] + sigma * (x[j] - xSmallest[j]);
                }
                setPoint(i, new PointValuePair(x, Double.NaN, false));
            }
            evaluate(evaluationFunction, comparator);
        }
    }
}