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.exception.NotStrictlyPositiveException; 21 import org.apache.commons.math.MaxIterationsExceededException; 22 import org.apache.commons.math.analysis.UnivariateRealFunction; 23 import org.apache.commons.math.optimization.GoalType; 24 25 /** 26 * Provide an interval that brackets a local optimum of a function. 27 * This code is based on a Python implementation (from <em>SciPy</em>, 28 * module {@code optimize.py} v0.5). 29 * @version $Revision$ $Date$ 30 * @since 2.2 31 */ 32 public class BracketFinder { 33 /** Tolerance to avoid division by zero. */ 34 private static final double EPS_MIN = 1e-21; 35 /** 36 * Golden section. 37 */ 38 private static final double GOLD = 1.618034; 39 /** 40 * Factor for expanding the interval. 41 */ 42 private final double growLimit; 43 /** 44 * Maximum number of iterations. 45 */ 46 private final int maxIterations; 47 /** 48 * Number of iterations. 49 */ 50 private int iterations; 51 /** 52 * Number of function evaluations. 53 */ 54 private int evaluations; 55 /** 56 * Lower bound of the bracket. 57 */ 58 private double lo; 59 /** 60 * Higher bound of the bracket. 61 */ 62 private double hi; 63 /** 64 * Point inside the bracket. 65 */ 66 private double mid; 67 /** 68 * Function value at {@link #lo}. 69 */ 70 private double fLo; 71 /** 72 * Function value at {@link #hi}. 73 */ 74 private double fHi; 75 /** 76 * Function value at {@link #mid}. 77 */ 78 private double fMid; 79 80 /** 81 * Constructor with default values {@code 100, 50} (see the 82 * {@link #BracketFinder(double,int) other constructor}). 83 */ 84 public BracketFinder() { 85 this(100, 50); 86 } 87 88 /** 89 * Create a bracketing interval finder. 90 * 91 * @param growLimit Expanding factor. 92 * @param maxIterations Maximum number of iterations allowed for finding 93 * a bracketing interval. 94 */ 95 public BracketFinder(double growLimit, 96 int maxIterations) { 97 if (growLimit <= 0) { 98 throw new NotStrictlyPositiveException(growLimit); 99 } 100 if (maxIterations <= 0) { 101 throw new NotStrictlyPositiveException(maxIterations); 102 } 103 104 this.growLimit = growLimit; 105 this.maxIterations = maxIterations; 106 } 107 108 /** 109 * Search new points that bracket a local optimum of the function. 110 * 111 * @param func Function whose optimum should be bracketted. 112 * @param goal {@link GoalType Goal type}. 113 * @param xA Initial point. 114 * @param xB Initial point. 115 * @throws MaxIterationsExceededException if the maximum iteration count 116 * is exceeded. 117 * @throws FunctionEvaluationException if an error occurs evaluating the function. 118 */ 119 public void search(UnivariateRealFunction func, 120 GoalType goal, 121 double xA, 122 double xB) 123 throws MaxIterationsExceededException, FunctionEvaluationException { 124 reset(); 125 final boolean isMinim = goal == GoalType.MINIMIZE; 126 127 double fA = eval(func, xA); 128 double fB = eval(func, xB); 129 if (isMinim ? 130 fA < fB : 131 fA > fB) { 132 double tmp = xA; 133 xA = xB; 134 xB = tmp; 135 136 tmp = fA; 137 fA = fB; 138 fB = tmp; 139 } 140 141 double xC = xB + GOLD * (xB - xA); 142 double fC = eval(func, xC); 143 144 while (isMinim ? fC < fB : fC > fB) { 145 if (++iterations > maxIterations) { 146 throw new MaxIterationsExceededException(maxIterations); 147 } 148 149 double tmp1 = (xB - xA) * (fB - fC); 150 double tmp2 = (xB - xC) * (fB - fA); 151 152 double val = tmp2 - tmp1; 153 double denom = Math.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val; 154 155 double w = xB - ((xB - xC) * tmp2 - (xB -xA) * tmp1) / denom; 156 double wLim = xB + growLimit * (xC - xB); 157 158 double fW; 159 if ((w - xC) * (xB - w) > 0) { 160 fW = eval(func, w); 161 if (isMinim ? 162 fW < fC : 163 fW > fC) { 164 xA = xB; 165 xB = w; 166 fA = fB; 167 fB = fW; 168 break; 169 } else if (isMinim ? 170 fW > fB : 171 fW < fB) { 172 xC = w; 173 fC = fW; 174 break; 175 } 176 w = xC + GOLD * (xC - xB); 177 fW = eval(func, w); 178 } else if ((w - wLim) * (wLim - xC) >= 0) { 179 w = wLim; 180 fW = eval(func, w); 181 } else if ((w - wLim) * (xC - w) > 0) { 182 fW = eval(func, w); 183 if (isMinim ? 184 fW < fC : 185 fW > fC) { 186 xB = xC; 187 xC = w; 188 w = xC + GOLD * (xC -xB); 189 fB = fC; 190 fC =fW; 191 fW = eval(func, w); 192 } 193 } else { 194 w = xC + GOLD * (xC - xB); 195 fW = eval(func, w); 196 } 197 198 xA = xB; 199 xB = xC; 200 xC = w; 201 fA = fB; 202 fB = fC; 203 fC = fW; 204 } 205 206 lo = xA; 207 mid = xB; 208 hi = xC; 209 fLo = fA; 210 fMid = fB; 211 fHi = fC; 212 } 213 214 /** 215 * @return the number of iterations. 216 */ 217 public int getIterations() { 218 return iterations; 219 } 220 /** 221 * @return the number of evaluations. 222 */ 223 public int getEvaluations() { 224 return evaluations; 225 } 226 227 /** 228 * @return the lower bound of the bracket. 229 * @see #getFLow() 230 */ 231 public double getLo() { 232 return lo; 233 } 234 235 /** 236 * Get function value at {@link #getLo()}. 237 * @return function value at {@link #getLo()} 238 */ 239 public double getFLow() { 240 return fLo; 241 } 242 243 /** 244 * @return the higher bound of the bracket. 245 * @see #getFHi() 246 */ 247 public double getHi() { 248 return hi; 249 } 250 251 /** 252 * Get function value at {@link #getHi()}. 253 * @return function value at {@link #getHi()} 254 */ 255 public double getFHi() { 256 return fHi; 257 } 258 259 /** 260 * @return a point in the middle of the bracket. 261 * @see #getFMid() 262 */ 263 public double getMid() { 264 return mid; 265 } 266 267 /** 268 * Get function value at {@link #getMid()}. 269 * @return function value at {@link #getMid()} 270 */ 271 public double getFMid() { 272 return fMid; 273 } 274 275 /** 276 * @param f Function. 277 * @param x Argument. 278 * @return {@code f(x)} 279 * @throws FunctionEvaluationException if function cannot be evaluated at x 280 */ 281 private double eval(UnivariateRealFunction f, double x) 282 throws FunctionEvaluationException { 283 284 ++evaluations; 285 return f.value(x); 286 } 287 288 /** 289 * Reset internal state. 290 */ 291 private void reset() { 292 iterations = 0; 293 evaluations = 0; 294 } 295 } 296