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.analysis.UnivariateRealFunction;
022    import org.apache.commons.math.optimization.GoalType;
023    
024    /**
025     * Implements Richard Brent's algorithm (from his book "Algorithms for
026     * Minimization without Derivatives", p. 79) for finding minima of real
027     * univariate functions.
028     *  
029     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
030     * @since 2.0
031     */
032    public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
033        
034        /**
035         * Golden section.
036         */
037        private static final double c = 0.5 * (3 - Math.sqrt(5));
038    
039        /**
040         * Construct a solver.
041         */
042        public BrentOptimizer() {
043            super(100, 1E-10);
044        }
045    
046        /** {@inheritDoc} */
047        public double optimize(final UnivariateRealFunction f, final GoalType goalType,
048                               final double min, final double max, final double startValue)
049            throws MaxIterationsExceededException, FunctionEvaluationException {
050            return optimize(f, goalType, min, max);
051        }
052        
053        /** {@inheritDoc} */
054        public double optimize(final UnivariateRealFunction f, final GoalType goalType,
055                               final double min, final double max)
056            throws MaxIterationsExceededException, FunctionEvaluationException {
057            clearResult();
058            return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy);
059        }
060        
061        /**
062         * Find the minimum of the function {@code f} within the interval {@code (a, b)}.
063         *
064         * If the function {@code f} is defined on the interval {@code (a, b)}, then
065         * this method finds an approximation {@code x} to the point at which {@code f}
066         * attains its minimum.<br/>
067         * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and
068         * {@code f} is never evaluated at two points closer together than {@code tol}.
069         * {@code eps} should be no smaller than <em>2 macheps</em> and preferable not
070         * much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative
071         * machine precision. {@code t} should be positive.
072         * @param f the function to solve
073         * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
074         * or {@link GoalType#MINIMIZE}
075         * @param a Lower bound of the interval
076         * @param b Higher bound of the interval
077         * @param eps Relative accuracy
078         * @param t Absolute accuracy
079         * @return the point at which the function is minimal.
080         * @throws MaxIterationsExceededException if the maximum iteration count
081         * is exceeded.
082         * @throws FunctionEvaluationException if an error occurs evaluating
083         * the function. 
084         */
085        private double localMin(final UnivariateRealFunction f, final GoalType goalType,
086                                double a, double b, final double eps, final double t)
087            throws MaxIterationsExceededException, FunctionEvaluationException {
088            double x = a + c * (b - a);
089            double v = x;
090            double w = x;
091            double e = 0;
092            double fx = computeObjectiveValue(f, x);
093            if (goalType == GoalType.MAXIMIZE) {
094                fx = -fx;
095            }
096            double fv = fx;
097            double fw = fx;
098    
099            int count = 0;
100            while (count < maximalIterationCount) {
101                double m = 0.5 * (a + b);
102                double tol = eps * Math.abs(x) + t;
103                double t2 = 2 * tol;
104    
105                // Check stopping criterion.
106                if (Math.abs(x - m) > t2 - 0.5 * (b - a)) {
107                    double p = 0;
108                    double q = 0;
109                    double r = 0;
110                    double d = 0;
111                    double u = 0;
112    
113                    if (Math.abs(e) > tol) { // Fit parabola.
114                        r = (x - w) * (fx - fv);
115                        q = (x - v) * (fx - fw);
116                        p = (x - v) * q - (x - w) * r;
117                        q = 2 * (q - r);
118    
119                        if (q > 0) {
120                            p = -p;
121                        } else {
122                            q = -q;
123                        }
124    
125                        r = e;
126                        e = d;
127                    }
128    
129                    if (Math.abs(p) < Math.abs(0.5 * q * r) &&
130                        (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step.
131                        d = p / q;
132                        u = x + d;
133    
134                        // f must not be evaluated too close to a or b.
135                        if (((u - a) < t2) || ((b - u) < t2)) {
136                            d = (x < m) ? tol : -tol;
137                        }
138                    } else { // Golden section step.
139                        e = ((x < m) ? b : a) - x;
140                        d = c * e;
141                    }
142    
143                    // f must not be evaluated too close to a or b.
144                    u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol));
145                    double fu = computeObjectiveValue(f, u);
146                    if (goalType == GoalType.MAXIMIZE) {
147                        fu = -fu;
148                    }
149    
150                    // Update a, b, v, w and x.
151                    if (fu <= fx) {
152                        if (u < x) {
153                            b = x;
154                        } else {
155                            a = x;
156                        }
157                        v = w;
158                        fv = fw;
159                        w = x;
160                        fw = fx;
161                        x = u;
162                        fx = fu;
163                    } else {
164                        if (u < x) {
165                            a = u;
166                        } else {
167                            b = u;
168                        }
169                        if ((fu <= fw) || (w == x)) {
170                            v = w;
171                            fv = fw;
172                            w = u;
173                            fw = fu;
174                        } else if ((fu <= fv) || (v == x) || (v == w)) {
175                            v = u;
176                            fv = fu;
177                        }
178                    }
179                } else { // termination
180                    setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count);
181                    return x;
182                }
183    
184                ++count;
185            }
186    
187            throw new MaxIterationsExceededException(maximalIterationCount);
188    
189        }
190    
191    }