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