1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math3.optim.univariate;
18
19 import org.apache.commons.math3.util.Incrementor;
20 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
21 import org.apache.commons.math3.exception.TooManyEvaluationsException;
22 import org.apache.commons.math3.exception.MaxCountExceededException;
23 import org.apache.commons.math3.analysis.UnivariateFunction;
24 import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
25
26 /**
27 * Provide an interval that brackets a local optimum of a function.
28 * This code is based on a Python implementation (from <em>SciPy</em>,
29 * module {@code optimize.py} v0.5).
30 *
31 * @version $Id: BracketFinder.java 1435539 2013-01-19 13:27:24Z tn $
32 * @since 2.2
33 */
34 public class BracketFinder {
35 /** Tolerance to avoid division by zero. */
36 private static final double EPS_MIN = 1e-21;
37 /**
38 * Golden section.
39 */
40 private static final double GOLD = 1.618034;
41 /**
42 * Factor for expanding the interval.
43 */
44 private final double growLimit;
45 /**
46 * Counter for function evaluations.
47 */
48 private final Incrementor evaluations = new Incrementor();
49 /**
50 * Lower bound of the bracket.
51 */
52 private double lo;
53 /**
54 * Higher bound of the bracket.
55 */
56 private double hi;
57 /**
58 * Point inside the bracket.
59 */
60 private double mid;
61 /**
62 * Function value at {@link #lo}.
63 */
64 private double fLo;
65 /**
66 * Function value at {@link #hi}.
67 */
68 private double fHi;
69 /**
70 * Function value at {@link #mid}.
71 */
72 private double fMid;
73
74 /**
75 * Constructor with default values {@code 100, 50} (see the
76 * {@link #BracketFinder(double,int) other constructor}).
77 */
78 public BracketFinder() {
79 this(100, 50);
80 }
81
82 /**
83 * Create a bracketing interval finder.
84 *
85 * @param growLimit Expanding factor.
86 * @param maxEvaluations Maximum number of evaluations allowed for finding
87 * a bracketing interval.
88 */
89 public BracketFinder(double growLimit,
90 int maxEvaluations) {
91 if (growLimit <= 0) {
92 throw new NotStrictlyPositiveException(growLimit);
93 }
94 if (maxEvaluations <= 0) {
95 throw new NotStrictlyPositiveException(maxEvaluations);
96 }
97
98 this.growLimit = growLimit;
99 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 }