/*
    * 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.math4.legacy.optim.nonlinear.scalar.noderiv;
  
  import java.util.Arrays;
  import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
  import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
  import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
  import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
  import org.apache.commons.math4.legacy.optim.PointValuePair;
  import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
  import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
  import org.apache.commons.math4.legacy.optim.univariate.UnivariatePointValuePair;
  import org.apache.commons.math4.core.jdkmath.JdkMath;
  
  /**
   * Powell's 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>).
   * <br>
   * The default stopping criterion is based on the differences of the
   * function value between two successive iterations. It is however possible
   * to define a custom convergence checker that might terminate the algorithm
   * earlier.
   * <br>
   * Constraints are not supported: the call to
   * {@link #optimize(org.apache.commons.math4.legacy.optim.OptimizationData...)} will throw
   * {@link MathUnsupportedOperationException} if bounds are passed to it.
   * In order to impose simple constraints, the objective function must be
   * wrapped in an adapter like
   * {@link org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
   * MultivariateFunctionMappingAdapter} or
   * {@link org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
   * MultivariateFunctionPenaltyAdapter}.
   *
   * @since 2.2
   */
  public class PowellOptimizer
      extends MultivariateOptimizer {
      /**
       * Minimum relative tolerance.
       */
      private static final double MIN_RELATIVE_TOLERANCE = 2 * JdkMath.ulp(1d);
      /** Relative threshold. */
      private final double relativeThreshold;
      /** Absolute threshold. */
      private final double absoluteThreshold;
      /** Relative threshold. */
      private final double lineSearchRelativeThreshold;
      /** Absolute threshold. */
      private final double lineSearchAbsoluteThreshold;
  
      /**
       * This constructor allows to specify a user-defined convergence checker,
       * in addition to the parameters that control the default convergence
       * checking procedure.
       * <br>
       * The internal line search tolerances are set to the square-root of their
       * corresponding value in the multivariate optimizer.
       *
       * @param rel Relative threshold.
       * @param abs Absolute threshold.
       * @param checker Convergence checker.
       * @throws NotStrictlyPositiveException if {@code abs <= 0}.
       * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
       */
      public PowellOptimizer(double rel,
                             double abs,
                             ConvergenceChecker<PointValuePair> checker) {
          this(rel, abs, JdkMath.sqrt(rel), JdkMath.sqrt(abs), checker);
      }
  
      /**
       * This constructor allows to specify a user-defined convergence checker,
       * in addition to the parameters that control the default convergence
       * checking procedure and the line search tolerances.
       *
       * @param rel Relative threshold for this optimizer.
       * @param abs Absolute threshold for this optimizer.
       * @param lineRel Relative threshold for the internal line search optimizer.
       * @param lineAbs Absolute threshold for the internal line search optimizer.
       * @param checker Convergence checker.
      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
      */
     public PowellOptimizer(double rel,
                            double abs,
                            double lineRel,
                            double lineAbs,
                            ConvergenceChecker<PointValuePair> checker) {
         super(checker);
 
         if (rel < MIN_RELATIVE_TOLERANCE) {
             throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
         }
         if (abs <= 0) {
             throw new NotStrictlyPositiveException(abs);
         }
 
         relativeThreshold = rel;
         absoluteThreshold = abs;
         lineSearchRelativeThreshold = lineRel;
         lineSearchAbsoluteThreshold = lineAbs;
     }
 
     /**
      * The parameters control the default convergence checking procedure.
      * <br>
      * The internal line search tolerances are set to the square-root of their
      * corresponding value in the multivariate optimizer.
      *
      * @param rel Relative threshold.
      * @param abs Absolute threshold.
      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
      */
     public PowellOptimizer(double rel,
                            double abs) {
         this(rel, abs, null);
     }
 
     /**
      * Builds an instance with the default convergence checking procedure.
      *
      * @param rel Relative threshold.
      * @param abs Absolute threshold.
      * @param lineRel Relative threshold for the internal line search optimizer.
      * @param lineAbs Absolute threshold for the internal line search optimizer.
      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
      */
     public PowellOptimizer(double rel,
                            double abs,
                            double lineRel,
                            double lineAbs) {
         this(rel, abs, lineRel, lineAbs, null);
     }
 
     /** {@inheritDoc} */
     @Override
     protected PointValuePair doOptimize() {
         checkParameters();
 
         // Line search optimizer.
         createLineSearch();
 
         final GoalType goal = getGoalType();
         final double[] guess = getStartPoint();
         final MultivariateFunction func = getObjectiveFunction();
         final int n = guess.length;
 
         final double[][] direc = new double[n][n];
         for (int i = 0; i < n; i++) {
             direc[i][i] = 1;
         }
 
         final ConvergenceChecker<PointValuePair> checker
             = getConvergenceChecker();
 
         double[] x = guess;
         double fVal = func.value(x);
         double[] x1 = x.clone();
         while (true) {
             incrementIterationCount();
 
             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], direc[i].length);
 
                 fX2 = fVal;
 
                 final UnivariatePointValuePair optimum = lineSearch(x, d);
                 fVal = optimum.getValue();
                 alphaMin = optimum.getPoint();
                 final double[][] result = newPointAndDirection(x, d, alphaMin);
                 x = result[0];
 
                 if ((fX2 - fVal) > delta) {
                     delta = fX2 - fVal;
                     bigInd = i;
                 }
             }
 
             // Default convergence check.
             boolean stop = 2 * (fX - fVal) <=
                 (relativeThreshold * (JdkMath.abs(fX) + JdkMath.abs(fVal)) +
                  absoluteThreshold);
 
             final PointValuePair previous = new PointValuePair(x1, fX);
             final PointValuePair current = new PointValuePair(x, fVal);
             if (!stop && checker != null) { // User-defined stopping criteria.
                 stop = checker.converged(getIterations(), previous, current);
             }
             if (stop) {
                 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 = func.value(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) {
                     final UnivariatePointValuePair optimum = lineSearch(x, d);
                     fVal = optimum.getValue();
                     alphaMin = optimum.getPoint();
                     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.
      *
      * @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[] nP = new double[n];
         final double[] nD = new double[n];
         for (int i = 0; i < n; i++) {
             nD[i] = d[i] * optimum;
             nP[i] = p[i] + nD[i];
         }
 
         final double[][] result = new double[2][];
         result[0] = nP;
         result[1] = nD;
 
         return result;
     }
 
     /**
      * @throws MathUnsupportedOperationException if bounds were passed to the
      * {@link #optimize(org.apache.commons.math4.legacy.optim.OptimizationData[]) optimize} method.
      */
     private void checkParameters() {
         if (getLowerBound() != null ||
             getUpperBound() != null) {
             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
         }
     }
 }