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

import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.RealPointValuePair;
import org.apache.commons.math.optimization.general.AbstractScalarDifferentiableOptimizer;
import org.apache.commons.math.optimization.univariate.AbstractUnivariateRealOptimizer;
import org.apache.commons.math.optimization.univariate.BracketFinder;
import org.apache.commons.math.optimization.univariate.BrentOptimizer;

/**
 * Powell algorithm.
 * This code is translated and adapted from the Python version of this
 * algorithm (as implemented in module {@code optimize.py} v0.5 of
 * <em>SciPy</em>).
 *
 * @version $Revision$ $Date$
 * @since 2.2
 */
public class PowellOptimizer
    extends AbstractScalarDifferentiableOptimizer {
    /**
     * Default relative tolerance for line search ({@value}).
     */
    public static final double DEFAULT_LS_RELATIVE_TOLERANCE = 1e-7;
    /**
     * Default absolute tolerance for line search ({@value}).
     */
    public static final double DEFAULT_LS_ABSOLUTE_TOLERANCE = 1e-11;
    /**
     * Line search.
     */
    private final LineSearch line;

    /**
     * Constructor with default line search tolerances (see the
     * {@link #PowellOptimizer(double,double) other constructor}).
     */
    public PowellOptimizer() {
        this(DEFAULT_LS_RELATIVE_TOLERANCE,
             DEFAULT_LS_ABSOLUTE_TOLERANCE);
    }

    /**
     * Constructor with default absolute line search tolerances (see
     * the {@link #PowellOptimizer(double,double) other constructor}).
     *
     * @param lsRelativeTolerance Relative error tolerance for
     * the line search algorithm ({@link BrentOptimizer}).
     */
    public PowellOptimizer(double lsRelativeTolerance) {
        this(lsRelativeTolerance,
             DEFAULT_LS_ABSOLUTE_TOLERANCE);
    }

    /**
     * @param lsRelativeTolerance Relative error tolerance for
     * the line search algorithm ({@link BrentOptimizer}).
     * @param lsAbsoluteTolerance Relative error tolerance for
     * the line search algorithm ({@link BrentOptimizer}).
     */
    public PowellOptimizer(double lsRelativeTolerance,
                           double lsAbsoluteTolerance) {
        line = new LineSearch(lsRelativeTolerance,
                              lsAbsoluteTolerance);
    }

    /** {@inheritDoc} */
    @Override
    protected RealPointValuePair doOptimize()
        throws FunctionEvaluationException,
               OptimizationException {
        final double[] guess = point.clone();
        final int n = guess.length;

        final double[][] direc = new double[n][n];
        for (int i = 0; i < n; i++) {
            direc[i][i] = 1;
        }

        double[] x = guess;
        double fVal = computeObjectiveValue(x);
        double[] x1 = x.clone();
        while (true) {
            incrementIterationsCounter();

            double fX = fVal;
            double fX2 = 0;
            double delta = 0;
            int bigInd = 0;
            double alphaMin = 0;

            for (int i = 0; i < n; i++) {
                final double[] d = /* Arrays.*/ copyOf(direc[i], n); // Java 1.5 does not support Arrays.copyOf()

                fX2 = fVal;

                line.search(x, d);
                fVal = line.getValueAtOptimum();
                alphaMin = line.getOptimum();
                final double[][] result = newPointAndDirection(x, d, alphaMin);
                x = result[0];

                if ((fX2 - fVal) > delta) {
                    delta = fX2 - fVal;
                    bigInd = i;
                }
            }

            final RealPointValuePair previous = new RealPointValuePair(x1, fX);
            final RealPointValuePair current = new RealPointValuePair(x, fVal);
            if (getConvergenceChecker().converged(getIterations(), previous, current)) {
                if (goal == GoalType.MINIMIZE) {
                    return (fVal < fX) ? current : previous;
                } else {
                    return (fVal > fX) ? current : previous;
                }
            }

            final double[] d = new double[n];
            final double[] x2 = new double[n];
            for (int i = 0; i < n; i++) {
                d[i] = x[i] - x1[i];
                x2[i] = 2 * x[i] - x1[i];
            }

            x1 = x.clone();
            fX2 = computeObjectiveValue(x2);

            if (fX > fX2) {
                double t = 2 * (fX + fX2 - 2 * fVal);
                double temp = fX - fVal - delta;
                t *= temp * temp;
                temp = fX - fX2;
                t -= delta * temp * temp;

                if (t < 0.0) {
                    line.search(x, d);
                    fVal = line.getValueAtOptimum();
                    alphaMin = line.getOptimum();
                    final double[][] result = newPointAndDirection(x, d, alphaMin);
                    x = result[0];

                    final int lastInd = n - 1;
                    direc[bigInd] = direc[lastInd];
                    direc[lastInd] = result[1];
                }
            }
        }
    }

    /**
     * Compute a new point (in the original space) and a new direction
     * vector, resulting from the line search.
     * The parameters {@code p} and {@code d} will be changed in-place.
     *
     * @param p Point used in the line search.
     * @param d Direction used in the line search.
     * @param optimum Optimum found by the line search.
     * @return a 2-element array containing the new point (at index 0) and
     * the new direction (at index 1).
     */
    private double[][] newPointAndDirection(double[] p,
                                            double[] d,
                                            double optimum) {
        final int n = p.length;
        final double[][] result = new double[2][n];
        final double[] nP = result[0];
        final double[] nD = result[1];
        for (int i = 0; i < n; i++) {
            nD[i] = d[i] * optimum;
            nP[i] = p[i] + nD[i];
        }
        return result;
    }

    /**
     * Class for finding the minimum of the objective function along a given
     * direction.
     */
    private class LineSearch {
        /**
         * Optimizer.
         */
        private final AbstractUnivariateRealOptimizer optim = new BrentOptimizer();
        /**
         * Automatic bracketing.
         */
        private final BracketFinder bracket = new BracketFinder();
        /**
         * Value of the optimum.
         */
        private double optimum = Double.NaN;
        /**
         * Value of the objective function at the optimum.
         */
        private double valueAtOptimum = Double.NaN;

        /**
         * @param relativeTolerance Relative tolerance.
         * @param absoluteTolerance Absolute tolerance.
         */
        public LineSearch(double relativeTolerance,
                          double absoluteTolerance) {
            optim.setRelativeAccuracy(relativeTolerance);
            optim.setAbsoluteAccuracy(absoluteTolerance);
        }

        /**
         * Find the minimum of the function {@code f(p + alpha * d)}.
         *
         * @param p Starting point.
         * @param d Search direction.
         * @throws FunctionEvaluationException if function cannot be evaluated at some test point
         * @throws OptimizationException if algorithm fails to converge
         */
        public void search(final double[] p, final double[] d)
            throws OptimizationException, FunctionEvaluationException {

            // Reset.
            optimum = Double.NaN;
            valueAtOptimum = Double.NaN;

            try {
                final int n = p.length;
                final UnivariateRealFunction f = new UnivariateRealFunction() {
                        public double value(double alpha)
                            throws FunctionEvaluationException {

                            final double[] x = new double[n];
                            for (int i = 0; i < n; i++) {
                                x[i] = p[i] + alpha * d[i];
                            }
                            final double obj;
                            obj = computeObjectiveValue(x);
                            return obj;
                        }
                    };

                bracket.search(f, goal, 0, 1);
                optimum = optim.optimize(f, goal,
                                         bracket.getLo(),
                                         bracket.getHi(),
                                         bracket.getMid());
                valueAtOptimum = optim.getFunctionValue();
            } catch (MaxIterationsExceededException e) {
                throw new OptimizationException(e);
            }
        }

        /**
         * @return the optimum.
         */
        public double getOptimum() {
            return optimum;
        }
        /**
         * @return the value of the function at the optimum.
         */
        public double getValueAtOptimum() {
            return valueAtOptimum;
        }
    }

    /**
     * Java 1.5 does not support Arrays.copyOf()
     *
     * @param source the array to be copied
     * @param newLen the length of the copy to be returned
     * @return the copied array, truncated or padded as necessary.
     */
     private double[] copyOf(double[] source, int newLen) {
         double[] output = new double[newLen];
         System.arraycopy(source, 0, output, 0, Math.min(source.length, newLen));
         return output;
     }

}