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