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    }