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.math3.optim.univariate; 018 019 import org.apache.commons.math3.util.Incrementor; 020 import org.apache.commons.math3.exception.NotStrictlyPositiveException; 021 import org.apache.commons.math3.exception.TooManyEvaluationsException; 022 import org.apache.commons.math3.exception.MaxCountExceededException; 023 import org.apache.commons.math3.analysis.UnivariateFunction; 024 import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; 025 026 /** 027 * Provide an interval that brackets a local optimum of a function. 028 * This code is based on a Python implementation (from <em>SciPy</em>, 029 * module {@code optimize.py} v0.5). 030 * 031 * @version $Id: BracketFinder.java 1413186 2012-11-24 13:47:59Z erans $ 032 * @since 2.2 033 */ 034 public class BracketFinder { 035 /** Tolerance to avoid division by zero. */ 036 private static final double EPS_MIN = 1e-21; 037 /** 038 * Golden section. 039 */ 040 private static final double GOLD = 1.618034; 041 /** 042 * Factor for expanding the interval. 043 */ 044 private final double growLimit; 045 /** 046 * Counter for function evaluations. 047 */ 048 private final Incrementor evaluations = new Incrementor(); 049 /** 050 * Lower bound of the bracket. 051 */ 052 private double lo; 053 /** 054 * Higher bound of the bracket. 055 */ 056 private double hi; 057 /** 058 * Point inside the bracket. 059 */ 060 private double mid; 061 /** 062 * Function value at {@link #lo}. 063 */ 064 private double fLo; 065 /** 066 * Function value at {@link #hi}. 067 */ 068 private double fHi; 069 /** 070 * Function value at {@link #mid}. 071 */ 072 private double fMid; 073 074 /** 075 * Constructor with default values {@code 100, 50} (see the 076 * {@link #BracketFinder(double,int) other constructor}). 077 */ 078 public BracketFinder() { 079 this(100, 50); 080 } 081 082 /** 083 * Create a bracketing interval finder. 084 * 085 * @param growLimit Expanding factor. 086 * @param maxEvaluations Maximum number of evaluations allowed for finding 087 * a bracketing interval. 088 */ 089 public BracketFinder(double growLimit, 090 int maxEvaluations) { 091 if (growLimit <= 0) { 092 throw new NotStrictlyPositiveException(growLimit); 093 } 094 if (maxEvaluations <= 0) { 095 throw new NotStrictlyPositiveException(maxEvaluations); 096 } 097 098 this.growLimit = growLimit; 099 evaluations.setMaximalCount(maxEvaluations); 100 } 101 102 /** 103 * Search new points that bracket a local optimum of the function. 104 * 105 * @param func Function whose optimum should be bracketed. 106 * @param goal {@link GoalType Goal type}. 107 * @param xA Initial point. 108 * @param xB Initial point. 109 * @throws TooManyEvaluationsException if the maximum number of evaluations 110 * is exceeded. 111 */ 112 public void search(UnivariateFunction func, GoalType goal, double xA, double xB) { 113 evaluations.resetCount(); 114 final boolean isMinim = goal == GoalType.MINIMIZE; 115 116 double fA = eval(func, xA); 117 double fB = eval(func, xB); 118 if (isMinim ? 119 fA < fB : 120 fA > fB) { 121 122 double tmp = xA; 123 xA = xB; 124 xB = tmp; 125 126 tmp = fA; 127 fA = fB; 128 fB = tmp; 129 } 130 131 double xC = xB + GOLD * (xB - xA); 132 double fC = eval(func, xC); 133 134 while (isMinim ? fC < fB : fC > fB) { 135 double tmp1 = (xB - xA) * (fB - fC); 136 double tmp2 = (xB - xC) * (fB - fA); 137 138 double val = tmp2 - tmp1; 139 double denom = Math.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val; 140 141 double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom; 142 double wLim = xB + growLimit * (xC - xB); 143 144 double fW; 145 if ((w - xC) * (xB - w) > 0) { 146 fW = eval(func, w); 147 if (isMinim ? 148 fW < fC : 149 fW > fC) { 150 xA = xB; 151 xB = w; 152 fA = fB; 153 fB = fW; 154 break; 155 } else if (isMinim ? 156 fW > fB : 157 fW < fB) { 158 xC = w; 159 fC = fW; 160 break; 161 } 162 w = xC + GOLD * (xC - xB); 163 fW = eval(func, w); 164 } else if ((w - wLim) * (wLim - xC) >= 0) { 165 w = wLim; 166 fW = eval(func, w); 167 } else if ((w - wLim) * (xC - w) > 0) { 168 fW = eval(func, w); 169 if (isMinim ? 170 fW < fC : 171 fW > fC) { 172 xB = xC; 173 xC = w; 174 w = xC + GOLD * (xC - xB); 175 fB = fC; 176 fC =fW; 177 fW = eval(func, w); 178 } 179 } else { 180 w = xC + GOLD * (xC - xB); 181 fW = eval(func, w); 182 } 183 184 xA = xB; 185 fA = fB; 186 xB = xC; 187 fB = fC; 188 xC = w; 189 fC = fW; 190 } 191 192 lo = xA; 193 fLo = fA; 194 mid = xB; 195 fMid = fB; 196 hi = xC; 197 fHi = fC; 198 199 if (lo > hi) { 200 double tmp = lo; 201 lo = hi; 202 hi = tmp; 203 204 tmp = fLo; 205 fLo = fHi; 206 fHi = tmp; 207 } 208 } 209 210 /** 211 * @return the number of evalutations. 212 */ 213 public int getMaxEvaluations() { 214 return evaluations.getMaximalCount(); 215 } 216 217 /** 218 * @return the number of evalutations. 219 */ 220 public int getEvaluations() { 221 return evaluations.getCount(); 222 } 223 224 /** 225 * @return the lower bound of the bracket. 226 * @see #getFLo() 227 */ 228 public double getLo() { 229 return lo; 230 } 231 232 /** 233 * Get function value at {@link #getLo()}. 234 * @return function value at {@link #getLo()} 235 */ 236 public double getFLo() { 237 return fLo; 238 } 239 240 /** 241 * @return the higher bound of the bracket. 242 * @see #getFHi() 243 */ 244 public double getHi() { 245 return hi; 246 } 247 248 /** 249 * Get function value at {@link #getHi()}. 250 * @return function value at {@link #getHi()} 251 */ 252 public double getFHi() { 253 return fHi; 254 } 255 256 /** 257 * @return a point in the middle of the bracket. 258 * @see #getFMid() 259 */ 260 public double getMid() { 261 return mid; 262 } 263 264 /** 265 * Get function value at {@link #getMid()}. 266 * @return function value at {@link #getMid()} 267 */ 268 public double getFMid() { 269 return fMid; 270 } 271 272 /** 273 * @param f Function. 274 * @param x Argument. 275 * @return {@code f(x)} 276 * @throws TooManyEvaluationsException if the maximal number of evaluations is 277 * exceeded. 278 */ 279 private double eval(UnivariateFunction f, double x) { 280 try { 281 evaluations.incrementCount(); 282 } catch (MaxCountExceededException e) { 283 throw new TooManyEvaluationsException(e.getMax()); 284 } 285 return f.value(x); 286 } 287 }