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 */ 017package org.apache.commons.math3.optim.univariate; 018 019import org.apache.commons.math3.util.FastMath; 020import org.apache.commons.math3.util.IntegerSequence; 021import org.apache.commons.math3.exception.NotStrictlyPositiveException; 022import org.apache.commons.math3.exception.TooManyEvaluationsException; 023import org.apache.commons.math3.exception.MaxCountExceededException; 024import org.apache.commons.math3.analysis.UnivariateFunction; 025import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; 026 027/** 028 * Provide an interval that brackets a local optimum of a function. 029 * This code is based on a Python implementation (from <em>SciPy</em>, 030 * module {@code optimize.py} v0.5). 031 * 032 * @since 2.2 033 */ 034public 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 IntegerSequence.Incrementor evaluations; 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, 500} (see the 076 * {@link #BracketFinder(double,int) other constructor}). 077 */ 078 public BracketFinder() { 079 this(100, 500); 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 = IntegerSequence.Incrementor.create().withMaximalCount(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, 113 GoalType goal, 114 double xA, 115 double xB) { 116 evaluations = evaluations.withStart(0); 117 final boolean isMinim = goal == GoalType.MINIMIZE; 118 119 double fA = eval(func, xA); 120 double fB = eval(func, xB); 121 if (isMinim ? 122 fA < fB : 123 fA > fB) { 124 125 double tmp = xA; 126 xA = xB; 127 xB = tmp; 128 129 tmp = fA; 130 fA = fB; 131 fB = tmp; 132 } 133 134 double xC = xB + GOLD * (xB - xA); 135 double fC = eval(func, xC); 136 137 while (isMinim ? fC < fB : fC > fB) { 138 double tmp1 = (xB - xA) * (fB - fC); 139 double tmp2 = (xB - xC) * (fB - fA); 140 141 double val = tmp2 - tmp1; 142 double denom = FastMath.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val; 143 144 double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom; 145 double wLim = xB + growLimit * (xC - xB); 146 147 double fW; 148 if ((w - xC) * (xB - w) > 0) { 149 fW = eval(func, w); 150 if (isMinim ? 151 fW < fC : 152 fW > fC) { 153 xA = xB; 154 xB = w; 155 fA = fB; 156 fB = fW; 157 break; 158 } else if (isMinim ? 159 fW > fB : 160 fW < fB) { 161 xC = w; 162 fC = fW; 163 break; 164 } 165 w = xC + GOLD * (xC - xB); 166 fW = eval(func, w); 167 } else if ((w - wLim) * (wLim - xC) >= 0) { 168 w = wLim; 169 fW = eval(func, w); 170 } else if ((w - wLim) * (xC - w) > 0) { 171 fW = eval(func, w); 172 if (isMinim ? 173 fW < fC : 174 fW > fC) { 175 xB = xC; 176 xC = w; 177 w = xC + GOLD * (xC - xB); 178 fB = fC; 179 fC =fW; 180 fW = eval(func, w); 181 } 182 } else { 183 w = xC + GOLD * (xC - xB); 184 fW = eval(func, w); 185 } 186 187 xA = xB; 188 fA = fB; 189 xB = xC; 190 fB = fC; 191 xC = w; 192 fC = fW; 193 } 194 195 lo = xA; 196 fLo = fA; 197 mid = xB; 198 fMid = fB; 199 hi = xC; 200 fHi = fC; 201 202 if (lo > hi) { 203 double tmp = lo; 204 lo = hi; 205 hi = tmp; 206 207 tmp = fLo; 208 fLo = fHi; 209 fHi = tmp; 210 } 211 } 212 213 /** 214 * @return the number of evalutations. 215 */ 216 public int getMaxEvaluations() { 217 return evaluations.getMaximalCount(); 218 } 219 220 /** 221 * @return the number of evalutations. 222 */ 223 public int getEvaluations() { 224 return evaluations.getCount(); 225 } 226 227 /** 228 * @return the lower bound of the bracket. 229 * @see #getFLo() 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 getFLo() { 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 TooManyEvaluationsException if the maximal number of evaluations is 280 * exceeded. 281 */ 282 private double eval(UnivariateFunction f, double x) { 283 try { 284 evaluations.increment(); 285 } catch (MaxCountExceededException e) { 286 throw new TooManyEvaluationsException(e.getMax()); 287 } 288 return f.value(x); 289 } 290}