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 }