    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.math3.analysis.solvers;
    018
    019
    020    import org.apache.commons.math3.exception.NoBracketingException;
    021    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    022    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    023    import org.apache.commons.math3.util.FastMath;
    024    import org.apache.commons.math3.util.Precision;
    025
    026    /**
    027     * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
    028     * Brent algorithm</a> for finding zeros of real univariate functions.
    029     * The function should be continuous but not necessarily smooth.
    030     * The {@code solve} method returns a zero {@code x} of the function {@code f}
    031     * in the given interval {@code [a, b]} to within a tolerance
    032     * {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and
    033     * {@code t} is the absolute accuracy.
    034     * The given interval must bracket the root.
    035     *
    036     * @version $Id: BrentSolver.java 1379560 2012-08-31 19:40:30Z erans $
    037     */
    038    public class BrentSolver extends AbstractUnivariateSolver {
    039
    040        /** Default absolute accuracy. */
    041        private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
    042
    043        /**
    044         * Construct a solver with default accuracy (1e-6).
    045         */
    046        public BrentSolver() {
    047            this(DEFAULT_ABSOLUTE_ACCURACY);
    048        }
    049        /**
    050         * Construct a solver.
    051         *
    052         * @param absoluteAccuracy Absolute accuracy.
    053         */
    054        public BrentSolver(double absoluteAccuracy) {
    055            super(absoluteAccuracy);
    056        }
    057        /**
    058         * Construct a solver.
    059         *
    060         * @param relativeAccuracy Relative accuracy.
    061         * @param absoluteAccuracy Absolute accuracy.
    062         */
    063        public BrentSolver(double relativeAccuracy,
    064                           double absoluteAccuracy) {
    065            super(relativeAccuracy, absoluteAccuracy);
    066        }
    067        /**
    068         * Construct a solver.
    069         *
    070         * @param relativeAccuracy Relative accuracy.
    071         * @param absoluteAccuracy Absolute accuracy.
    072         * @param functionValueAccuracy Function value accuracy.
    073         */
    074        public BrentSolver(double relativeAccuracy,
    075                           double absoluteAccuracy,
    076                           double functionValueAccuracy) {
    077            super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
    078        }
    079
    080        /**
    081         * {@inheritDoc}
    082         */
    083        @Override
    084        protected double doSolve()
    085            throws NoBracketingException,
    086                   TooManyEvaluationsException,
    087                   NumberIsTooLargeException {
    088            double min = getMin();
    089            double max = getMax();
    090            final double initial = getStartValue();
    091            final double functionValueAccuracy = getFunctionValueAccuracy();
    092
    093            verifySequence(min, initial, max);
    094
    095            // Return the initial guess if it is good enough.
    096            double yInitial = computeObjectiveValue(initial);
    097            if (FastMath.abs(yInitial) <= functionValueAccuracy) {
    098                return initial;
    099            }
    100
    101            // Return the first endpoint if it is good enough.
    102            double yMin = computeObjectiveValue(min);
    103            if (FastMath.abs(yMin) <= functionValueAccuracy) {
    104                return min;
    105            }
    106
    107            // Reduce interval if min and initial bracket the root.
    108            if (yInitial * yMin < 0) {
    109                return brent(min, initial, yMin, yInitial);
    110            }
    111
    112            // Return the second endpoint if it is good enough.
    113            double yMax = computeObjectiveValue(max);
    114            if (FastMath.abs(yMax) <= functionValueAccuracy) {
    115                return max;
    116            }
    117
    118            // Reduce interval if initial and max bracket the root.
    119            if (yInitial * yMax < 0) {
    120                return brent(initial, max, yInitial, yMax);
    121            }
    122
    123            throw new NoBracketingException(min, max, yMin, yMax);
    124        }
    125
    126        /**
    127         * Search for a zero inside the provided interval.
    128         * This implementation is based on the algorithm described at page 58 of
    129         * the book
    130         * <quote>
    131         *  <b>Algorithms for Minimization Without Derivatives</b>
    132         *  <it>Richard P. Brent</it>
    133         *  Dover 0-486-41998-3
    134         * </quote>
    135         *
    136         * @param lo Lower bound of the search interval.
    137         * @param hi Higher bound of the search interval.
    138         * @param fLo Function value at the lower bound of the search interval.
    139         * @param fHi Function value at the higher bound of the search interval.
    140         * @return the value where the function is zero.
    141         */
    142        private double brent(double lo, double hi,
    143                             double fLo, double fHi) {
    144            double a = lo;
    145            double fa = fLo;
    146            double b = hi;
    147            double fb = fHi;
    148            double c = a;
    149            double fc = fa;
    150            double d = b - a;
    151            double e = d;
    152
    153            final double t = getAbsoluteAccuracy();
    154            final double eps = getRelativeAccuracy();
    155
    156            while (true) {
    157                if (FastMath.abs(fc) < FastMath.abs(fb)) {
    158                    a = b;
    159                    b = c;
    160                    c = a;
    161                    fa = fb;
    162                    fb = fc;
    163                    fc = fa;
    164                }
    165
    166                final double tol = 2 * eps * FastMath.abs(b) + t;
    167                final double m = 0.5 * (c - b);
    168
    169                if (FastMath.abs(m) <= tol ||
    170                    Precision.equals(fb, 0))  {
    171                    return b;
    172                }
    173                if (FastMath.abs(e) < tol ||
    174                    FastMath.abs(fa) <= FastMath.abs(fb)) {
    175                    // Force bisection.
    176                    d = m;
    177                    e = d;
    178                } else {
    179                    double s = fb / fa;
    180                    double p;
    181                    double q;
    182                    // The equality test (a == c) is intentional,
    183                    // it is part of the original Brent's method and
    184                    // it should NOT be replaced by proximity test.
    185                    if (a == c) {
    186                        // Linear interpolation.
    187                        p = 2 * m * s;
    188                        q = 1 - s;
    189                    } else {
    190                        // Inverse quadratic interpolation.
    191                        q = fa / fc;
    192                        final double r = fb / fc;
    193                        p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
    194                        q = (q - 1) * (r - 1) * (s - 1);
    195                    }
    196                    if (p > 0) {
    197                        q = -q;
    198                    } else {
    199                        p = -p;
    200                    }
    201                    s = e;
    202                    e = d;
    203                    if (p >= 1.5 * m * q - FastMath.abs(tol * q) ||
    204                        p >= FastMath.abs(0.5 * s * q)) {
    205                        // Inverse quadratic interpolation gives a value
    206                        // in the wrong direction, or progress is slow.
    207                        // Fall back to bisection.
    208                        d = m;
    209                        e = d;
    210                    } else {
    211                        d = p / q;
    212                    }
    213                }
    214                a = b;
    215                fa = fb;
    216
    217                if (FastMath.abs(d) > tol) {
    218                    b += d;
    219                } else if (m > 0) {
    220                    b += tol;
    221                } else {
    222                    b -= tol;
    223                }
    224                fb = computeObjectiveValue(b);
    225                if ((fb > 0 && fc > 0) ||
    226                    (fb <= 0 && fc <= 0)) {
    227                    c = a;
    228                    fc = fa;
    229                    d = b - a;
    230                    e = d;
    231                }
    232            }
    233        }
    234    }