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 }