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.optimization.univariate;
018    
019    import org.apache.commons.math3.util.Precision;
020    import org.apache.commons.math3.util.FastMath;
021    import org.apache.commons.math3.exception.NumberIsTooSmallException;
022    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
023    import org.apache.commons.math3.optimization.ConvergenceChecker;
024    import org.apache.commons.math3.optimization.GoalType;
025    
026    /**
027     * For a function defined on some interval {@code (lo, hi)}, this class
028     * finds an approximation {@code x} to the point at which the function
029     * attains its minimum.
030     * It implements Richard Brent's algorithm (from his book "Algorithms for
031     * Minimization without Derivatives", p. 79) for finding minima of real
032     * univariate functions.
033     * <br/>
034     * This code is an adaptation, partly based on the Python code from SciPy
035     * (module "optimize.py" v0.5); the original algorithm is also modified
036     * <ul>
037     *  <li>to use an initial guess provided by the user,</li>
038     *  <li>to ensure that the best point encountered is the one returned.</li>
039     * </ul>
040     *
041     * @version $Id: BrentOptimizer.java 1422230 2012-12-15 12:11:13Z erans $
042     * @deprecated As of 3.1 (to be removed in 4.0).
043     * @since 2.0
044     */
045    @Deprecated
046    public class BrentOptimizer extends BaseAbstractUnivariateOptimizer {
047        /**
048         * Golden section.
049         */
050        private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
051        /**
052         * Minimum relative tolerance.
053         */
054        private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
055        /**
056         * Relative threshold.
057         */
058        private final double relativeThreshold;
059        /**
060         * Absolute threshold.
061         */
062        private final double absoluteThreshold;
063    
064        /**
065         * The arguments are used implement the original stopping criterion
066         * of Brent's algorithm.
067         * {@code abs} and {@code rel} define a tolerance
068         * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
069         * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
070         * where <em>macheps</em> is the relative machine precision. {@code abs} must
071         * be positive.
072         *
073         * @param rel Relative threshold.
074         * @param abs Absolute threshold.
075         * @param checker Additional, user-defined, convergence checking
076         * procedure.
077         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
078         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
079         */
080        public BrentOptimizer(double rel,
081                              double abs,
082                              ConvergenceChecker<UnivariatePointValuePair> checker) {
083            super(checker);
084    
085            if (rel < MIN_RELATIVE_TOLERANCE) {
086                throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
087            }
088            if (abs <= 0) {
089                throw new NotStrictlyPositiveException(abs);
090            }
091    
092            relativeThreshold = rel;
093            absoluteThreshold = abs;
094        }
095    
096        /**
097         * The arguments are used for implementing the original stopping criterion
098         * of Brent's algorithm.
099         * {@code abs} and {@code rel} define a tolerance
100         * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
101         * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
102         * where <em>macheps</em> is the relative machine precision. {@code abs} must
103         * be positive.
104         *
105         * @param rel Relative threshold.
106         * @param abs Absolute threshold.
107         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
108         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
109         */
110        public BrentOptimizer(double rel,
111                              double abs) {
112            this(rel, abs, null);
113        }
114    
115        /** {@inheritDoc} */
116        @Override
117        protected UnivariatePointValuePair doOptimize() {
118            final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
119            final double lo = getMin();
120            final double mid = getStartValue();
121            final double hi = getMax();
122    
123            // Optional additional convergence criteria.
124            final ConvergenceChecker<UnivariatePointValuePair> checker
125                = getConvergenceChecker();
126    
127            double a;
128            double b;
129            if (lo < hi) {
130                a = lo;
131                b = hi;
132            } else {
133                a = hi;
134                b = lo;
135            }
136    
137            double x = mid;
138            double v = x;
139            double w = x;
140            double d = 0;
141            double e = 0;
142            double fx = computeObjectiveValue(x);
143            if (!isMinim) {
144                fx = -fx;
145            }
146            double fv = fx;
147            double fw = fx;
148    
149            UnivariatePointValuePair previous = null;
150            UnivariatePointValuePair current
151                = new UnivariatePointValuePair(x, isMinim ? fx : -fx);
152            // Best point encountered so far (which is the initial guess).
153            UnivariatePointValuePair best = current;
154    
155            int iter = 0;
156            while (true) {
157                final double m = 0.5 * (a + b);
158                final double tol1 = relativeThreshold * FastMath.abs(x) + absoluteThreshold;
159                final double tol2 = 2 * tol1;
160    
161                // Default stopping criterion.
162                final boolean stop = FastMath.abs(x - m) <= tol2 - 0.5 * (b - a);
163                if (!stop) {
164                    double p = 0;
165                    double q = 0;
166                    double r = 0;
167                    double u = 0;
168    
169                    if (FastMath.abs(e) > tol1) { // Fit parabola.
170                        r = (x - w) * (fx - fv);
171                        q = (x - v) * (fx - fw);
172                        p = (x - v) * q - (x - w) * r;
173                        q = 2 * (q - r);
174    
175                        if (q > 0) {
176                            p = -p;
177                        } else {
178                            q = -q;
179                        }
180    
181                        r = e;
182                        e = d;
183    
184                        if (p > q * (a - x) &&
185                            p < q * (b - x) &&
186                            FastMath.abs(p) < FastMath.abs(0.5 * q * r)) {
187                            // Parabolic interpolation step.
188                            d = p / q;
189                            u = x + d;
190    
191                            // f must not be evaluated too close to a or b.
192                            if (u - a < tol2 || b - u < tol2) {
193                                if (x <= m) {
194                                    d = tol1;
195                                } else {
196                                    d = -tol1;
197                                }
198                            }
199                        } else {
200                            // Golden section step.
201                            if (x < m) {
202                                e = b - x;
203                            } else {
204                                e = a - x;
205                            }
206                            d = GOLDEN_SECTION * e;
207                        }
208                    } else {
209                        // Golden section step.
210                        if (x < m) {
211                            e = b - x;
212                        } else {
213                            e = a - x;
214                        }
215                        d = GOLDEN_SECTION * e;
216                    }
217    
218                    // Update by at least "tol1".
219                    if (FastMath.abs(d) < tol1) {
220                        if (d >= 0) {
221                            u = x + tol1;
222                        } else {
223                            u = x - tol1;
224                        }
225                    } else {
226                        u = x + d;
227                    }
228    
229                    double fu = computeObjectiveValue(u);
230                    if (!isMinim) {
231                        fu = -fu;
232                    }
233    
234                    // User-defined convergence checker.
235                    previous = current;
236                    current = new UnivariatePointValuePair(u, isMinim ? fu : -fu);
237                    best = best(best,
238                                best(previous,
239                                     current,
240                                     isMinim),
241                                isMinim);
242    
243                    if (checker != null) {
244                        if (checker.converged(iter, previous, current)) {
245                            return best;
246                        }
247                    }
248    
249                    // Update a, b, v, w and x.
250                    if (fu <= fx) {
251                        if (u < x) {
252                            b = x;
253                        } else {
254                            a = x;
255                        }
256                        v = w;
257                        fv = fw;
258                        w = x;
259                        fw = fx;
260                        x = u;
261                        fx = fu;
262                    } else {
263                        if (u < x) {
264                            a = u;
265                        } else {
266                            b = u;
267                        }
268                        if (fu <= fw ||
269                            Precision.equals(w, x)) {
270                            v = w;
271                            fv = fw;
272                            w = u;
273                            fw = fu;
274                        } else if (fu <= fv ||
275                                   Precision.equals(v, x) ||
276                                   Precision.equals(v, w)) {
277                            v = u;
278                            fv = fu;
279                        }
280                    }
281                } else { // Default termination (Brent's criterion).
282                    return best(best,
283                                best(previous,
284                                     current,
285                                     isMinim),
286                                isMinim);
287                }
288                ++iter;
289            }
290        }
291    
292        /**
293         * Selects the best of two points.
294         *
295         * @param a Point and value.
296         * @param b Point and value.
297         * @param isMinim {@code true} if the selected point must be the one with
298         * the lowest value.
299         * @return the best point, or {@code null} if {@code a} and {@code b} are
300         * both {@code null}. When {@code a} and {@code b} have the same function
301         * value, {@code a} is returned.
302         */
303        private UnivariatePointValuePair best(UnivariatePointValuePair a,
304                                              UnivariatePointValuePair b,
305                                              boolean isMinim) {
306            if (a == null) {
307                return b;
308            }
309            if (b == null) {
310                return a;
311            }
312    
313            if (isMinim) {
314                return a.getValue() <= b.getValue() ? a : b;
315            } else {
316                return a.getValue() >= b.getValue() ? a : b;
317            }
318        }
319    }