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.analysis.UnivariateRealFunction; 022 import org.apache.commons.math.optimization.GoalType; 023 024 /** 025 * Implements Richard Brent's algorithm (from his book "Algorithms for 026 * Minimization without Derivatives", p. 79) for finding minima of real 027 * univariate functions. 028 * 029 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 030 * @since 2.0 031 */ 032 public class BrentOptimizer extends AbstractUnivariateRealOptimizer { 033 034 /** 035 * Golden section. 036 */ 037 private static final double c = 0.5 * (3 - Math.sqrt(5)); 038 039 /** 040 * Construct a solver. 041 */ 042 public BrentOptimizer() { 043 super(100, 1E-10); 044 } 045 046 /** {@inheritDoc} */ 047 public double optimize(final UnivariateRealFunction f, final GoalType goalType, 048 final double min, final double max, final double startValue) 049 throws MaxIterationsExceededException, FunctionEvaluationException { 050 return optimize(f, goalType, min, max); 051 } 052 053 /** {@inheritDoc} */ 054 public double optimize(final UnivariateRealFunction f, final GoalType goalType, 055 final double min, final double max) 056 throws MaxIterationsExceededException, FunctionEvaluationException { 057 clearResult(); 058 return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy); 059 } 060 061 /** 062 * Find the minimum of the function {@code f} within the interval {@code (a, b)}. 063 * 064 * If the function {@code f} is defined on the interval {@code (a, b)}, then 065 * this method finds an approximation {@code x} to the point at which {@code f} 066 * attains its minimum.<br/> 067 * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and 068 * {@code f} is never evaluated at two points closer together than {@code tol}. 069 * {@code eps} should be no smaller than <em>2 macheps</em> and preferable not 070 * much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative 071 * machine precision. {@code t} should be positive. 072 * @param f the function to solve 073 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} 074 * or {@link GoalType#MINIMIZE} 075 * @param a Lower bound of the interval 076 * @param b Higher bound of the interval 077 * @param eps Relative accuracy 078 * @param t Absolute accuracy 079 * @return the point at which the function is minimal. 080 * @throws MaxIterationsExceededException if the maximum iteration count 081 * is exceeded. 082 * @throws FunctionEvaluationException if an error occurs evaluating 083 * the function. 084 */ 085 private double localMin(final UnivariateRealFunction f, final GoalType goalType, 086 double a, double b, final double eps, final double t) 087 throws MaxIterationsExceededException, FunctionEvaluationException { 088 double x = a + c * (b - a); 089 double v = x; 090 double w = x; 091 double e = 0; 092 double fx = computeObjectiveValue(f, x); 093 if (goalType == GoalType.MAXIMIZE) { 094 fx = -fx; 095 } 096 double fv = fx; 097 double fw = fx; 098 099 int count = 0; 100 while (count < maximalIterationCount) { 101 double m = 0.5 * (a + b); 102 double tol = eps * Math.abs(x) + t; 103 double t2 = 2 * tol; 104 105 // Check stopping criterion. 106 if (Math.abs(x - m) > t2 - 0.5 * (b - a)) { 107 double p = 0; 108 double q = 0; 109 double r = 0; 110 double d = 0; 111 double u = 0; 112 113 if (Math.abs(e) > tol) { // Fit parabola. 114 r = (x - w) * (fx - fv); 115 q = (x - v) * (fx - fw); 116 p = (x - v) * q - (x - w) * r; 117 q = 2 * (q - r); 118 119 if (q > 0) { 120 p = -p; 121 } else { 122 q = -q; 123 } 124 125 r = e; 126 e = d; 127 } 128 129 if (Math.abs(p) < Math.abs(0.5 * q * r) && 130 (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step. 131 d = p / q; 132 u = x + d; 133 134 // f must not be evaluated too close to a or b. 135 if (((u - a) < t2) || ((b - u) < t2)) { 136 d = (x < m) ? tol : -tol; 137 } 138 } else { // Golden section step. 139 e = ((x < m) ? b : a) - x; 140 d = c * e; 141 } 142 143 // f must not be evaluated too close to a or b. 144 u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol)); 145 double fu = computeObjectiveValue(f, u); 146 if (goalType == GoalType.MAXIMIZE) { 147 fu = -fu; 148 } 149 150 // Update a, b, v, w and x. 151 if (fu <= fx) { 152 if (u < x) { 153 b = x; 154 } else { 155 a = x; 156 } 157 v = w; 158 fv = fw; 159 w = x; 160 fw = fx; 161 x = u; 162 fx = fu; 163 } else { 164 if (u < x) { 165 a = u; 166 } else { 167 b = u; 168 } 169 if ((fu <= fw) || (w == x)) { 170 v = w; 171 fv = fw; 172 w = u; 173 fw = fu; 174 } else if ((fu <= fv) || (v == x) || (v == w)) { 175 v = u; 176 fv = fu; 177 } 178 } 179 } else { // termination 180 setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count); 181 return x; 182 } 183 184 ++count; 185 } 186 187 throw new MaxIterationsExceededException(maximalIterationCount); 188 189 } 190 191 }