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 */ 17 package org.apache.commons.math.optimization.univariate; 18 19 import org.apache.commons.math.FunctionEvaluationException; 20 import org.apache.commons.math.MaxIterationsExceededException; 21 import org.apache.commons.math.exception.NotStrictlyPositiveException; 22 import org.apache.commons.math.optimization.GoalType; 23 import 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 fvr. 2011) $ 32 * @since 2.0 33 */ 34 public 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
