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.exception.NotStrictlyPositiveException;
021 import org.apache.commons.math.MaxIterationsExceededException;
022 import org.apache.commons.math.analysis.UnivariateRealFunction;
023 import org.apache.commons.math.optimization.GoalType;
024
025 /**
026 * Provide an interval that brackets a local optimum of a function.
027 * This code is based on a Python implementation (from <em>SciPy</em>,
028 * module {@code optimize.py} v0.5).
029 * @version $Revision$ $Date$
030 * @since 2.2
031 */
032 public class BracketFinder {
033 /** Tolerance to avoid division by zero. */
034 private static final double EPS_MIN = 1e-21;
035 /**
036 * Golden section.
037 */
038 private static final double GOLD = 1.618034;
039 /**
040 * Factor for expanding the interval.
041 */
042 private final double growLimit;
043 /**
044 * Maximum number of iterations.
045 */
046 private final int maxIterations;
047 /**
048 * Number of iterations.
049 */
050 private int iterations;
051 /**
052 * Number of function evaluations.
053 */
054 private int evaluations;
055 /**
056 * Lower bound of the bracket.
057 */
058 private double lo;
059 /**
060 * Higher bound of the bracket.
061 */
062 private double hi;
063 /**
064 * Point inside the bracket.
065 */
066 private double mid;
067 /**
068 * Function value at {@link #lo}.
069 */
070 private double fLo;
071 /**
072 * Function value at {@link #hi}.
073 */
074 private double fHi;
075 /**
076 * Function value at {@link #mid}.
077 */
078 private double fMid;
079
080 /**
081 * Constructor with default values {@code 100, 50} (see the
082 * {@link #BracketFinder(double,int) other constructor}).
083 */
084 public BracketFinder() {
085 this(100, 50);
086 }
087
088 /**
089 * Create a bracketing interval finder.
090 *
091 * @param growLimit Expanding factor.
092 * @param maxIterations Maximum number of iterations allowed for finding
093 * a bracketing interval.
094 */
095 public BracketFinder(double growLimit,
096 int maxIterations) {
097 if (growLimit <= 0) {
098 throw new NotStrictlyPositiveException(growLimit);
099 }
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 }