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