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