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    }