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.math4.legacy.optim.univariate;
18
19 import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
20 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
21 import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
22 import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
23 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
24 import org.apache.commons.math4.core.jdkmath.JdkMath;
25 import org.apache.commons.math4.legacy.core.IntegerSequence.Incrementor;
26
27 /**
28 * Provide an interval that brackets a local optimum of a function.
29 * This code is based on a Python implementation (from <em>SciPy</em>,
30 * module {@code optimize.py} v0.5).
31 *
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 * Number of allowed function evaluations.
47 */
48 private final int maxEvaluations;
49 /**
50 * Number of function evaluations performed in the last search.
51 */
52 private int evaluations;
53 /**
54 * Lower bound of the bracket.
55 */
56 private double lo;
57 /**
58 * Higher bound of the bracket.
59 */
60 private double hi;
61 /**
62 * Point inside the bracket.
63 */
64 private double mid;
65 /**
66 * Function value at {@link #lo}.
67 */
68 private double fLo;
69 /**
70 * Function value at {@link #hi}.
71 */
72 private double fHi;
73 /**
74 * Function value at {@link #mid}.
75 */
76 private double fMid;
77
78 /**
79 * Constructor with default values {@code 100, 500} (see the
80 * {@link #BracketFinder(double,int) other constructor}).
81 */
82 public BracketFinder() {
83 this(100, 500);
84 }
85
86 /**
87 * Create a bracketing interval finder.
88 *
89 * @param growLimit Expanding factor.
90 * @param maxEvaluations Maximum number of evaluations allowed for finding
91 * a bracketing interval.
92 */
93 public BracketFinder(double growLimit,
94 int maxEvaluations) {
95 if (growLimit <= 0) {
96 throw new NotStrictlyPositiveException(growLimit);
97 }
98 if (maxEvaluations <= 0) {
99 throw new NotStrictlyPositiveException(maxEvaluations);
100 }
101
102 this.growLimit = growLimit;
103 this.maxEvaluations = maxEvaluations;
104 }
105
106 /**
107 * Search new points that bracket a local optimum of the function.
108 *
109 * @param func Function whose optimum should be bracketed.
110 * @param goal {@link GoalType Goal type}.
111 * @param xA Initial point.
112 * @param xB Initial point.
113 * @throws TooManyEvaluationsException if the maximum number of evaluations
114 * is exceeded.
115 */
116 public void search(UnivariateFunction func,
117 GoalType goal,
118 double xA,
119 double xB) {
120 final FunctionEvaluator eval = new FunctionEvaluator(func);
121 final boolean isMinim = goal == GoalType.MINIMIZE;
122
123 double fA = eval.value(xA);
124 double fB = eval.value(xB);
125 if (isMinim ?
126 fA < fB :
127 fA > fB) {
128
129 double tmp = xA;
130 xA = xB;
131 xB = tmp;
132
133 tmp = fA;
134 fA = fB;
135 fB = tmp;
136 }
137
138 double xC = xB + GOLD * (xB - xA);
139 double fC = eval.value(xC);
140
141 while (isMinim ? fC < fB : fC > fB) {
142 double tmp1 = (xB - xA) * (fB - fC);
143 double tmp2 = (xB - xC) * (fB - fA);
144
145 double val = tmp2 - tmp1;
146 double denom = JdkMath.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;
147
148 double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom;
149 double wLim = xB + growLimit * (xC - xB);
150
151 double fW;
152 if ((w - xC) * (xB - w) > 0) {
153 fW = eval.value(w);
154 if (isMinim ?
155 fW < fC :
156 fW > fC) {
157 xA = xB;
158 xB = w;
159 fA = fB;
160 fB = fW;
161 break;
162 } else if (isMinim ?
163 fW > fB :
164 fW < fB) {
165 xC = w;
166 fC = fW;
167 break;
168 }
169 w = xC + GOLD * (xC - xB);
170 fW = eval.value(w);
171 } else if ((w - wLim) * (wLim - xC) >= 0) {
172 w = wLim;
173 fW = eval.value(w);
174 } else if ((w - wLim) * (xC - w) > 0) {
175 fW = eval.value(w);
176 if (isMinim ?
177 fW < fC :
178 fW > fC) {
179 xB = xC;
180 xC = w;
181 w = xC + GOLD * (xC - xB);
182 fB = fC;
183 fC =fW;
184 fW = eval.value(w);
185 }
186 } else {
187 w = xC + GOLD * (xC - xB);
188 fW = eval.value(w);
189 }
190
191 xA = xB;
192 fA = fB;
193 xB = xC;
194 fB = fC;
195 xC = w;
196 fC = fW;
197 }
198
199 lo = xA;
200 fLo = fA;
201 mid = xB;
202 fMid = fB;
203 hi = xC;
204 fHi = fC;
205
206 if (lo > hi) {
207 double tmp = lo;
208 lo = hi;
209 hi = tmp;
210
211 tmp = fLo;
212 fLo = fHi;
213 fHi = tmp;
214 }
215 }
216
217 /**
218 * @return the number of evaluations.
219 */
220 public int getMaxEvaluations() {
221 return maxEvaluations;
222 }
223
224 /**
225 * @return the number of evaluations.
226 */
227 public int getEvaluations() {
228 return evaluations;
229 }
230
231 /**
232 * @return the lower bound of the bracket.
233 * @see #getFLo()
234 */
235 public double getLo() {
236 return lo;
237 }
238
239 /**
240 * Get function value at {@link #getLo()}.
241 * @return function value at {@link #getLo()}
242 */
243 public double getFLo() {
244 return fLo;
245 }
246
247 /**
248 * @return the higher bound of the bracket.
249 * @see #getFHi()
250 */
251 public double getHi() {
252 return hi;
253 }
254
255 /**
256 * Get function value at {@link #getHi()}.
257 * @return function value at {@link #getHi()}
258 */
259 public double getFHi() {
260 return fHi;
261 }
262
263 /**
264 * @return a point in the middle of the bracket.
265 * @see #getFMid()
266 */
267 public double getMid() {
268 return mid;
269 }
270
271 /**
272 * Get function value at {@link #getMid()}.
273 * @return function value at {@link #getMid()}
274 */
275 public double getFMid() {
276 return fMid;
277 }
278
279 /**
280 * Utility for incrementing a counter at each function evaluation.
281 */
282 private final class FunctionEvaluator {
283 /** Function. */
284 private final UnivariateFunction func;
285 /** Counter. */
286 private final Incrementor inc;
287
288 /**
289 * @param func Function.
290 */
291 FunctionEvaluator(UnivariateFunction func) {
292 this.func = func;
293 inc = Incrementor.create().withMaximalCount(maxEvaluations);
294 evaluations = 0;
295 }
296
297 /**
298 * @param x Argument.
299 * @return {@code f(x)}
300 * @throws TooManyEvaluationsException if the maximal number of evaluations is
301 * exceeded.
302 */
303 double value(double x) {
304 try {
305 inc.increment();
306 evaluations = inc.getCount();
307 } catch (MaxCountExceededException e) {
308 throw new TooManyEvaluationsException(e.getMax());
309 }
310
311 return func.value(x);
312 }
313 }
314 }