1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.optimization.univariate;
18
19import org.apache.commons.math.FunctionEvaluationException;
20import org.apache.commons.math.MaxIterationsExceededException;
21import org.apache.commons.math.exception.NotStrictlyPositiveException;
22import org.apache.commons.math.optimization.GoalType;
23import org.apache.commons.math.util.FastMath;
24
25/**
26 * Implements Richard Brent's algorithm (from his book "Algorithms for
27 * Minimization without Derivatives", p. 79) for finding minima of real
28 * univariate functions. This implementation is an adaptation partly
29 * based on the Python code from SciPy (module "optimize.py" v0.5).
30 *
31 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
32 * @since 2.0
33 */
34public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
35    /**
36     * Golden section.
37     */
38    private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
39
40    /**
41     * Construct a solver.
42     */
43    public BrentOptimizer() {
44        setMaxEvaluations(1000);
45        setMaximalIterationCount(100);
46        setAbsoluteAccuracy(1e-11);
47        setRelativeAccuracy(1e-9);
48    }
49
50    /** {@inheritDoc} */
51    @Override
52    protected double doOptimize()
53        throws MaxIterationsExceededException, FunctionEvaluationException {
54        return localMin(getGoalType() == GoalType.MINIMIZE,
55                        getMin(), getStartValue(), getMax(),
56                        getRelativeAccuracy(), getAbsoluteAccuracy());
57    }
58
59    /**
60     * Find the minimum of the function within the interval {@code (lo, hi)}.
61     *
62     * If the function is defined on the interval {@code (lo, hi)}, then
63     * this method finds an approximation {@code x} to the point at which
64     * the function attains its minimum.<br/>
65     * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t}
66     * and the function is never evaluated at two points closer together than
67     * {@code tol}. {@code eps} should be no smaller than <em>2 macheps</em> and
68     * preferable not much less than <em>sqrt(macheps)</em>, where
69     * <em>macheps</em> is the relative machine precision. {@code t} should be
70     * positive.
71     * @param isMinim {@code true} when minimizing the function.
72     * @param lo Lower bound of the interval.
73     * @param mid Point inside the interval {@code [lo, hi]}.
74     * @param hi Higher bound of the interval.
75     * @param eps Relative accuracy.
76     * @param t Absolute accuracy.
77     * @return the optimum point.
78     * @throws MaxIterationsExceededException if the maximum iteration count
79     * is exceeded.
80     * @throws FunctionEvaluationException if an error occurs evaluating the function.
81     */
82    private double localMin(boolean isMinim,
83                            double lo, double mid, double hi,
84                            double eps, double t)
85        throws MaxIterationsExceededException, FunctionEvaluationException {
86        if (eps <= 0) {
87            throw new NotStrictlyPositiveException(eps);
88        }
89        if (t <= 0) {
90            throw new NotStrictlyPositiveException(t);
91        }
92        double a;
93        double b;
94        if (lo < hi) {
95            a = lo;
96            b = hi;
97        } else {
98            a = hi;
99            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}
228