    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.optimization.univariate;
    018
    019    import org.apache.commons.math.FunctionEvaluationException;
    020    import org.apache.commons.math.MaxIterationsExceededException;
    021    import org.apache.commons.math.exception.NotStrictlyPositiveException;
    022    import org.apache.commons.math.optimization.GoalType;
    023    import org.apache.commons.math.util.FastMath;
    024
    025    /**
    026     * Implements Richard Brent's algorithm (from his book "Algorithms for
    027     * Minimization without Derivatives", p. 79) for finding minima of real
    028     * univariate functions. This implementation is an adaptation partly
    029     * based on the Python code from SciPy (module "optimize.py" v0.5).
    030     *
    031     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    032     * @since 2.0
    033     */
    034    public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
    035        /**
    036         * Golden section.
    037         */
    038        private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
    039
    040        /**
    041         * Construct a solver.
    042         */
    043        public BrentOptimizer() {
    044            setMaxEvaluations(1000);
    045            setMaximalIterationCount(100);
    046            setAbsoluteAccuracy(1e-11);
    047            setRelativeAccuracy(1e-9);
    048        }
    049
    050        /** {@inheritDoc} */
    051        @Override
    052        protected double doOptimize()
    053            throws MaxIterationsExceededException, FunctionEvaluationException {
    054            return localMin(getGoalType() == GoalType.MINIMIZE,
    055                            getMin(), getStartValue(), getMax(),
    056                            getRelativeAccuracy(), getAbsoluteAccuracy());
    057        }
    058
    059        /**
    060         * Find the minimum of the function within the interval {@code (lo, hi)}.
    061         *
    062         * If the function is defined on the interval {@code (lo, hi)}, then
    063         * this method finds an approximation {@code x} to the point at which
    064         * the function attains its minimum.<br/>
    065         * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t}
    066         * and the function is never evaluated at two points closer together than
    067         * {@code tol}. {@code eps} should be no smaller than <em>2 macheps</em> and
    068         * preferable not much less than <em>sqrt(macheps)</em>, where
    069         * <em>macheps</em> is the relative machine precision. {@code t} should be
    070         * positive.
    071         * @param isMinim {@code true} when minimizing the function.
    072         * @param lo Lower bound of the interval.
    073         * @param mid Point inside the interval {@code [lo, hi]}.
    074         * @param hi Higher bound of the interval.
    075         * @param eps Relative accuracy.
    076         * @param t Absolute accuracy.
    077         * @return the optimum point.
    078         * @throws MaxIterationsExceededException if the maximum iteration count
    079         * is exceeded.
    080         * @throws FunctionEvaluationException if an error occurs evaluating the function.
    081         */
    082        private double localMin(boolean isMinim,
    083                                double lo, double mid, double hi,
    084                                double eps, double t)
    085            throws MaxIterationsExceededException, FunctionEvaluationException {
    086            if (eps <= 0) {
    087                throw new NotStrictlyPositiveException(eps);
    088            }
    089            if (t <= 0) {
    090                throw new NotStrictlyPositiveException(t);
    091            }
    092            double a;
    093            double b;
    094            if (lo < hi) {
    095                a = lo;
    096                b = hi;
    097            } else {
    098                a = hi;
    099                b = lo;
    100            }
    101
    102            double x = mid;
    103            double v = x;
    104            double w = x;
    105            double d = 0;
    106            double e = 0;
    107            double fx = computeObjectiveValue(x);
    108            if (!isMinim) {
    109                fx = -fx;
    110            }
    111            double fv = fx;
    112            double fw = fx;
    113
    114            while (true) {
    115                double m = 0.5 * (a + b);
    116                final double tol1 = eps * FastMath.abs(x) + t;
    117                final double tol2 = 2 * tol1;
    118
    119                // Check stopping criterion.
    120                if (FastMath.abs(x - m) > tol2 - 0.5 * (b - a)) {
    121                    double p = 0;
    122                    double q = 0;
    123                    double r = 0;
    124                    double u = 0;
    125
    126                    if (FastMath.abs(e) > tol1) { // Fit parabola.
    127                        r = (x - w) * (fx - fv);
    128                        q = (x - v) * (fx - fw);
    129                        p = (x - v) * q - (x - w) * r;
    130                        q = 2 * (q - r);
    131
    132                        if (q > 0) {
    133                            p = -p;
    134                        } else {
    135                            q = -q;
    136                        }
    137
    138                        r = e;
    139                        e = d;
    140
    141                        if (p > q * (a - x) &&
    142                            p < q * (b - x) &&
    143                            FastMath.abs(p) < FastMath.abs(0.5 * q * r)) {
    144                            // Parabolic interpolation step.
    145                            d = p / q;
    146                            u = x + d;
    147
    148                            // f must not be evaluated too close to a or b.
    149                            if (u - a < tol2 || b - u < tol2) {
    150                                if (x <= m) {
    151                                    d = tol1;
    152                                } else {
    153                                    d = -tol1;
    154                                }
    155                            }
    156                        } else {
    157                            // Golden section step.
    158                            if (x < m) {
    159                                e = b - x;
    160                            } else {
    161                                e = a - x;
    162                            }
    163                            d = GOLDEN_SECTION * e;
    164                        }
    165                    } else {
    166                        // Golden section step.
    167                        if (x < m) {
    168                            e = b - x;
    169                        } else {
    170                            e = a - x;
    171                        }
    172                        d = GOLDEN_SECTION * e;
    173                    }
    174
    175                    // Update by at least "tol1".
    176                    if (FastMath.abs(d) < tol1) {
    177                        if (d >= 0) {
    178                            u = x + tol1;
    179                        } else {
    180                            u = x - tol1;
    181                        }
    182                    } else {
    183                        u = x + d;
    184                    }
    185
    186                    double fu = computeObjectiveValue(u);
    187                    if (!isMinim) {
    188                        fu = -fu;
    189                    }
    190
    191                    // Update a, b, v, w and x.
    192                    if (fu <= fx) {
    193                        if (u < x) {
    194                            b = x;
    195                        } else {
    196                            a = x;
    197                        }
    198                        v = w;
    199                        fv = fw;
    200                        w = x;
    201                        fw = fx;
    202                        x = u;
    203                        fx = fu;
    204                    } else {
    205                        if (u < x) {
    206                            a = u;
    207                        } else {
    208                            b = u;
    209                        }
    210                        if (fu <= fw || w == x) {
    211                            v = w;
    212                            fv = fw;
    213                            w = u;
    214                            fw = fu;
    215                        } else if (fu <= fv || v == x || v == w) {
    216                            v = u;
    217                            fv = fu;
    218                        }
    219                    }
    220                } else { // termination
    221                    setFunctionValue(isMinim ? fx : -fx);
    222                    return x;
    223                }
    224                incrementIterationsCounter();
    225            }
    226        }
    227    }