src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java

/*
 * 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.exception.NotStrictlyPositiveException;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.util.FastMath;

/**
 * Implements Richard Brent's algorithm (from his book "Algorithms for
 * Minimization without Derivatives", p. 79) for finding minima of real
 * univariate functions. This implementation is an adaptation partly
 * based on the Python code from SciPy (module "optimize.py" v0.5).
 *
 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
 * @since 2.0
 */
public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
    /**
     * Golden section.
     */
    private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));

    /**
     * Construct a solver.
     */
    public BrentOptimizer() {
        setMaxEvaluations(1000);
        setMaximalIterationCount(100);
        setAbsoluteAccuracy(1e-11);
        setRelativeAccuracy(1e-9);
    }

    /** {@inheritDoc} */
    @Override
    protected double doOptimize()
        throws MaxIterationsExceededException, FunctionEvaluationException {
        return localMin(getGoalType() == GoalType.MINIMIZE,
                        getMin(), getStartValue(), getMax(),
                        getRelativeAccuracy(), getAbsoluteAccuracy());
    }

    /**
     * Find the minimum of the function within the interval {@code (lo, hi)}.
     *
     * If the function is defined on the interval {@code (lo, hi)}, then
     * this method finds an approximation {@code x} to the point at which
     * the function attains its minimum.<br/>
     * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t}
     * and the function 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 isMinim {@code true} when minimizing the function.
     * @param lo Lower bound of the interval.
     * @param mid Point inside the interval {@code [lo, hi]}.
     * @param hi Higher bound of the interval.
     * @param eps Relative accuracy.
     * @param t Absolute accuracy.
     * @return the optimum point.
     * @throws MaxIterationsExceededException if the maximum iteration count
     * is exceeded.
     * @throws FunctionEvaluationException if an error occurs evaluating the function.
     */
    private double localMin(boolean isMinim,
                            double lo, double mid, double hi,
                            double eps, double t)
        throws MaxIterationsExceededException, FunctionEvaluationException {
        if (eps <= 0) {
            throw new NotStrictlyPositiveException(eps);
        }
        if (t <= 0) {
            throw new NotStrictlyPositiveException(t);
        }
        double a;
        double b;
        if (lo < hi) {
            a = lo;
            b = hi;
        } else {
            a = hi;
            b = lo;
        }

        double x = mid;
        double v = x;
        double w = x;
        double d = 0;
        double e = 0;
        double fx = computeObjectiveValue(x);
        if (!isMinim) {
            fx = -fx;
        }
        double fv = fx;
        double fw = fx;

        while (true) {
            double m = 0.5 * (a + b);
            final double tol1 = eps * FastMath.abs(x) + t;
            final double tol2 = 2 * tol1;

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

                if (FastMath.abs(e) > tol1) { // 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 (p > q * (a - x) &&
                        p < q * (b - x) &&
                        FastMath.abs(p) < FastMath.abs(0.5 * q * r)) {
                        // Parabolic interpolation step.
                        d = p / q;
                        u = x + d;

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

                // Update by at least "tol1".
                if (FastMath.abs(d) < tol1) {
                    if (d >= 0) {
                        u = x + tol1;
                    } else {
                        u = x - tol1;
                    }
                } else {
                    u = x + d;
                }

                double fu = computeObjectiveValue(u);
                if (!isMinim) {
                    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
                setFunctionValue(isMinim ? fx : -fx);
                return x;
            }
            incrementIterationsCounter();
        }
    }
}