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

import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.optimization.GoalType;

/**
 * Implements Richard Brent's algorithm (from his book "Algorithms for
 * Minimization without Derivatives", p. 79) for finding minima of real
 * univariate functions.
 *  
 * @version $Revision$ $Date$
 * @since 2.0
 */
public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
    
    /**
     * Golden section.
     */
    private static final double c = 0.5 * (3 - Math.sqrt(5));

    /**
     * Construct a solver.
     */
    public BrentOptimizer() {
        super(100, 1E-10);
    }

    /** {@inheritDoc} */
    public double optimize(final UnivariateRealFunction f, final GoalType goalType,
                           final double min, final double max, final double startValue)
        throws MaxIterationsExceededException, FunctionEvaluationException {
        return optimize(f, goalType, min, max);
    }
    
    /** {@inheritDoc} */
    public double optimize(final UnivariateRealFunction f, final GoalType goalType,
                           final double min, final double max)
        throws MaxIterationsExceededException, FunctionEvaluationException {
        clearResult();
        return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy);
    }
    
    /**
     * Find the minimum of the function {@code f} within the interval {@code (a, b)}.
     *
     * If the function {@code f} is defined on the interval {@code (a, b)}, then
     * this method finds an approximation {@code x} to the point at which {@code f}
     * attains its minimum.<br/>
     * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and
     * {@code f} is never evaluated at two points closer together than {@code tol}.
     * {@code eps} should be no smaller than <em>2 macheps</em> and preferable not
     * much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative
     * machine precision. {@code t} should be positive.
     * @param f the function to solve
     * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
     * or {@link GoalType#MINIMIZE}
     * @param a Lower bound of the interval
     * @param b Higher bound of the interval
     * @param eps Relative accuracy
     * @param t Absolute accuracy
     * @return the point at which the function is minimal.
     * @throws MaxIterationsExceededException if the maximum iteration count
     * is exceeded.
     * @throws FunctionEvaluationException if an error occurs evaluating
     * the function. 
     */
    private double localMin(final UnivariateRealFunction f, final GoalType goalType,
                            double a, double b, final double eps, final double t)
        throws MaxIterationsExceededException, FunctionEvaluationException {
        double x = a + c * (b - a);
        double v = x;
        double w = x;
        double e = 0;
        double fx = computeObjectiveValue(f, x);
        if (goalType == GoalType.MAXIMIZE) {
            fx = -fx;
        }
        double fv = fx;
        double fw = fx;

        int count = 0;
        while (count < maximalIterationCount) {
            double m = 0.5 * (a + b);
            double tol = eps * Math.abs(x) + t;
            double t2 = 2 * tol;

            // Check stopping criterion.
            if (Math.abs(x - m) > t2 - 0.5 * (b - a)) {
                double p = 0;
                double q = 0;
                double r = 0;
                double d = 0;
                double u = 0;

                if (Math.abs(e) > tol) { // Fit parabola.
                    r = (x - w) * (fx - fv);
                    q = (x - v) * (fx - fw);
                    p = (x - v) * q - (x - w) * r;
                    q = 2 * (q - r);

                    if (q > 0) {
                        p = -p;
                    } else {
                        q = -q;
                    }

                    r = e;
                    e = d;
                }

                if (Math.abs(p) < Math.abs(0.5 * q * r) &&
                    (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step.
                    d = p / q;
                    u = x + d;

                    // f must not be evaluated too close to a or b.
                    if (((u - a) < t2) || ((b - u) < t2)) {
                        d = (x < m) ? tol : -tol;
                    }
                } else { // Golden section step.
                    e = ((x < m) ? b : a) - x;
                    d = c * e;
                }

                // f must not be evaluated too close to a or b.
                u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol));
                double fu = computeObjectiveValue(f, u);
                if (goalType == GoalType.MAXIMIZE) {
                    fu = -fu;
                }

                // Update a, b, v, w and x.
                if (fu <= fx) {
                    if (u < x) {
                        b = x;
                    } else {
                        a = x;
                    }
                    v = w;
                    fv = fw;
                    w = x;
                    fw = fx;
                    x = u;
                    fx = fu;
                } else {
                    if (u < x) {
                        a = u;
                    } else {
                        b = u;
                    }
                    if ((fu <= fw) || (w == x)) {
                        v = w;
                        fv = fw;
                        w = u;
                        fw = fu;
                    } else if ((fu <= fv) || (v == x) || (v == w)) {
                        v = u;
                        fv = fu;
                    }
                }
            } else { // termination
                setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count);
                return x;
            }

            ++count;
        }

        throw new MaxIterationsExceededException(maximalIterationCount);

    }

}
