1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.optimization.univariate; 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.FunctionEvaluationException; 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MaxIterationsExceededException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.NotStrictlyPositiveException; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.optimization.GoalType; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Implements Richard Brent's algorithm (from his book "Algorithms for 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Minimization without Derivatives", p. 79) for finding minima of real 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * univariate functions. This implementation is an adaptation partly 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * based on the Python code from SciPy (module "optimize.py" v0.5). 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class BrentOptimizer extends AbstractUnivariateRealOptimizer { 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Golden section. 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5)); 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Construct a solver. 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BrentOptimizer() { 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMaxEvaluations(1000); 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setMaximalIterationCount(100); 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setAbsoluteAccuracy(1e-11); 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setRelativeAccuracy(1e-9); 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double doOptimize() 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws MaxIterationsExceededException, FunctionEvaluationException { 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return localMin(getGoalType() == GoalType.MINIMIZE, 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond getMin(), getStartValue(), getMax(), 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond getRelativeAccuracy(), getAbsoluteAccuracy()); 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Find the minimum of the function within the interval {@code (lo, hi)}. 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * If the function is defined on the interval {@code (lo, hi)}, then 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this method finds an approximation {@code x} to the point at which 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the function attains its minimum.<br/> 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * and the function is never evaluated at two points closer together than 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@code tol}. {@code eps} should be no smaller than <em>2 macheps</em> and 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * preferable not much less than <em>sqrt(macheps)</em>, where 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <em>macheps</em> is the relative machine precision. {@code t} should be 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * positive. 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param isMinim {@code true} when minimizing the function. 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param lo Lower bound of the interval. 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param mid Point inside the interval {@code [lo, hi]}. 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param hi Higher bound of the interval. 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param eps Relative accuracy. 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param t Absolute accuracy. 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the optimum point. 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MaxIterationsExceededException if the maximum iteration count 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * is exceeded. 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws FunctionEvaluationException if an error occurs evaluating the function. 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double localMin(boolean isMinim, 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double lo, double mid, double hi, 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double eps, double t) 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws MaxIterationsExceededException, FunctionEvaluationException { 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (eps <= 0) { 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new NotStrictlyPositiveException(eps); 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (t <= 0) { 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new NotStrictlyPositiveException(t); 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double a; 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double b; 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (lo < hi) { 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a = lo; 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b = hi; 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a = hi; 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b = lo; 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double x = mid; 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double v = x; 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double w = x; 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double d = 0; 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double e = 0; 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fx = computeObjectiveValue(x); 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!isMinim) { 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fx = -fx; 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fv = fx; 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fw = fx; 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while (true) { 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double m = 0.5 * (a + b); 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double tol1 = eps * FastMath.abs(x) + t; 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double tol2 = 2 * tol1; 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Check stopping criterion. 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.abs(x - m) > tol2 - 0.5 * (b - a)) { 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double p = 0; 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double q = 0; 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double r = 0; 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double u = 0; 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.abs(e) > tol1) { // Fit parabola. 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond r = (x - w) * (fx - fv); 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond q = (x - v) * (fx - fw); 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond p = (x - v) * q - (x - w) * r; 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond q = 2 * (q - r); 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (q > 0) { 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond p = -p; 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond q = -q; 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond r = e; 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond e = d; 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (p > q * (a - x) && 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond p < q * (b - x) && 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.abs(p) < FastMath.abs(0.5 * q * r)) { 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Parabolic interpolation step. 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond d = p / q; 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond u = x + d; 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // f must not be evaluated too close to a or b. 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (u - a < tol2 || b - u < tol2) { 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (x <= m) { 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond d = tol1; 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond d = -tol1; 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Golden section step. 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (x < m) { 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond e = b - x; 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond e = a - x; 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond d = GOLDEN_SECTION * e; 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Golden section step. 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (x < m) { 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond e = b - x; 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond e = a - x; 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond d = GOLDEN_SECTION * e; 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Update by at least "tol1". 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.abs(d) < tol1) { 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (d >= 0) { 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond u = x + tol1; 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond u = x - tol1; 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond u = x + d; 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fu = computeObjectiveValue(u); 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!isMinim) { 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fu = -fu; 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Update a, b, v, w and x. 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (fu <= fx) { 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (u < x) { 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b = x; 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a = x; 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v = w; 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fv = fw; 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond w = x; 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fw = fx; 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond x = u; 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fx = fu; 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (u < x) { 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a = u; 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b = u; 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (fu <= fw || w == x) { 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v = w; 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fv = fw; 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond w = u; 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fw = fu; 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else if (fu <= fv || v == x || v == w) { 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v = u; 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fv = fu; 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { // termination 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setFunctionValue(isMinim ? fx : -fx); 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return x; 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond incrementIterationsCounter(); 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 228