1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.optimization;
17
18 import com.irurueta.numerical.EvaluationException;
19 import com.irurueta.numerical.InvalidBracketRangeException;
20 import com.irurueta.numerical.LockedException;
21 import com.irurueta.numerical.NotReadyException;
22 import com.irurueta.numerical.SingleDimensionFunctionEvaluatorListener;
23
24 /**
25 * This class uses Brent algorithm to determine a local function minimum for
26 * single dimension functions.
27 * Brent's algorithm will search for a local minimum inside the provided or
28 * computed bracket of values.
29 * Brent's algorithm is not the fastest among all optimization algorithms, but
30 * it is usually one that provides good convergence for most continuous
31 * functions.
32 * It 's recommended to always set or compute a bracket of values, as the search
33 * range is reduced and results usually become more accurate.
34 * The implementation of this class is based on Numerical Recipes 3rd ed.
35 * Section 10.3. Page 496.
36 */
37 public class BrentSingleOptimizer extends BracketedSingleOptimizer {
38
39 /**
40 * Is the maximum allowed number of iterations.
41 */
42 public static final int ITMAX = 100;
43
44 /**
45 * Is the golden ratio.
46 */
47 public static final double CGOLD = 0.3819660;
48
49 /**
50 * Small number that protects against trying to achieve fractional accuracy
51 * for a minimum that happens to be exactly zero.
52 */
53 public static final double ZEPS = 1e-10;
54
55 /**
56 * Constant defining the default accuracy of the estimated minimum.
57 */
58 public static final double DEFAULT_TOLERANCE = 3e-8;
59
60 /**
61 * Minimum allowed tolerance.
62 */
63 public static final double MIN_TOLERANCE = 0.0;
64
65 /**
66 * Tolerance value. The algorithm will iterate until the result converges
67 * below this value of accuracy or until the maximum number of iterations is
68 * achieved (and in such case, convergence will be assumed to have failed).
69 */
70 private double tolerance;
71
72 /**
73 * Empty constructor.
74 */
75 public BrentSingleOptimizer() {
76 super();
77 tolerance = DEFAULT_TOLERANCE;
78 }
79
80 /**
81 * Constructor. Creates an instance with provided bracket of values.
82 *
83 * @param minEvalPoint Minimum bracket evaluation point.
84 * @param middleEvalPoint Middle bracket evaluation point.
85 * @param maxEvalPoint Maximum bracket evaluation point.
86 * @param tolerance Tolerance or accuracy to be obtained in estimated
87 * minimum.
88 * @throws InvalidBracketRangeException Raised if the following condition is
89 * not met: minEvalPoint <= middleEvalPoint <= maxEvalPoint.
90 * @throws IllegalArgumentException Raised if tolerance is negative.
91 */
92 public BrentSingleOptimizer(
93 final double minEvalPoint, final double middleEvalPoint, final double maxEvalPoint, final double tolerance)
94 throws InvalidBracketRangeException {
95 super(minEvalPoint, middleEvalPoint, maxEvalPoint);
96 internalSetTolerance(tolerance);
97 }
98
99 /**
100 * Constructor. Creates an instance with provided bracket of values and a
101 * listener to get single dimension function evaluations.
102 *
103 * @param listener Listener to evaluate a function.
104 * @param minEvalPoint Minimum bracket evaluation point.
105 * @param middleEvalPoint Middle bracket evaluation point.
106 * @param maxEvalPoint Maximum bracket evaluation point.
107 * @param tolerance Tolerance or accuracy to be obtained in estimated
108 * minimum.
109 * @throws InvalidBracketRangeException Raised if the following condition is
110 * not met: minEvalPoint <= middleEvalPoint <= maxEvalPoint.
111 * @throws IllegalArgumentException Raised if tolerance is negative.
112 */
113 public BrentSingleOptimizer(
114 final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
115 final double middleEvalPoint, final double maxEvalPoint, final double tolerance)
116 throws InvalidBracketRangeException {
117 super(listener, minEvalPoint, middleEvalPoint, maxEvalPoint);
118 internalSetTolerance(tolerance);
119 }
120
121 /**
122 * Returns tolerance value, which is the accuracy to be obtained when a
123 * minimum is estimated.
124 * The algorithm will iterate until the result converges below this value of
125 * accuracy or until the maximum number of iterations is achieved (and in
126 * such case, convergence will be assumed to have failed).
127 *
128 * @return Tolerance value.
129 */
130 public double getTolerance() {
131 return tolerance;
132 }
133
134 /**
135 * Sets algorithm's tolerance.
136 * The algorithm will iterate until the result converges below this value of
137 * accuracy or until the maximum number of iterations is achieved (an in
138 * such case, convergence will be assumed to have failed).
139 *
140 * @param tolerance Tolerance or accuracy to be obtained in estimated
141 * minimum.
142 * @throws LockedException Raised if this instance is locked. This instance
143 * will be locked while doing some operations. Attempting to change any
144 * parameter while being locked will raise this exception.
145 * @throws IllegalArgumentException Raised if tolerance is negative.
146 */
147 public void setTolerance(final double tolerance) throws LockedException {
148 if (isLocked()) {
149 throw new LockedException();
150 }
151 internalSetTolerance(tolerance);
152 }
153
154 /**
155 * This function estimates a function minimum within provided or computed
156 * bracket of values.
157 * Given a function f, and given a bracketing triplet of abscissas "ax", "bx",
158 * "cx" (such that bx is between ax and cx, and f(bx) is less than both f(ax)
159 * and f(cx), this routine isolates the minimum to a fractional prevision of
160 * about tolerance using Brent's method. The abscissa of the minimum is
161 * returned as "xmin", and the function value of the minimum is returned as
162 * "fmin", the returned function value.
163 *
164 * @throws LockedException Raised if this instance is locked, because
165 * estimation is being computed.
166 * @throws NotReadyException Raised if this instance is not ready because
167 * either a listener or a bracket has not yet been provided or computed.
168 * @throws OptimizationException Raised if the algorithm failed because of
169 * lack of convergence or because function couldn't be evaluated.
170 */
171 @Override
172 @SuppressWarnings("Duplicates")
173 public void minimize() throws LockedException, NotReadyException, OptimizationException {
174 if (isLocked()) {
175 throw new LockedException();
176 }
177 if (!isReady()) {
178 throw new NotReadyException();
179 }
180
181 locked = true;
182 final var v1 = new double[1];
183 final var v2 = new double[2];
184 final var v3 = new double[3];
185
186 try {
187 double a;
188 double b;
189 var d = 0.0;
190 double etemp;
191 double fu;
192 double fv;
193 double fw;
194 double fx;
195 double p;
196 double q;
197 double r;
198 double tol1;
199 double tol2;
200 double u;
201 double v;
202 double w;
203 double x;
204 double xm;
205 //This will be the distance moved on the step before last.
206 var e = 0.0;
207
208 //a and b must be in ascending order, but input abscissas need not
209 //be.
210 a = Math.min(ax, cx);
211 b = (bx > cx ? ax : cx);
212 //Initializations...
213 x = w = v = bx;
214 fw = fv = fx = listener.evaluate(x);
215
216 for (var iter = 0; iter < ITMAX; iter++) {
217 //Main program loop
218 xm = 0.5 * (a + b);
219 tol1 = tolerance * Math.abs(x) + ZEPS;
220 tol2 = 2.0 * tol1;
221
222 if (Math.abs(x - xm) <= (tol2 - 0.5 * (b - a))) {
223 //Test for done here.
224 fmin = fx;
225 xmin = x;
226
227 resultAvailable = true;
228 locked = false;
229 return;
230 }
231 if (Math.abs(e) > tol1) {
232 //Construct a trial parabolic fit.
233 r = (x - w) * (fx - fv);
234 q = (x - v) * (fx - fw);
235 p = (x - v) * q - (x - w) * r;
236 q = 2.0 * (q - r);
237
238 if (q > 0.0) {
239 p = -p;
240 }
241 q = Math.abs(q);
242 etemp = e;
243 e = d;
244
245 if (Math.abs(p) >= Math.abs(0.5 * q * etemp) ||
246 p <= q * (a - x) || p >= q * (b - x)) {
247 //noinspection all
248 e = x >= xm ? a - x : b - x;
249 d = CGOLD * (e);
250 //The above conditions determine the acceptability of
251 //the parabolic fit. Here we take the golden section
252 //step into the larger of the two segments.
253 } else {
254 //Take the parabolic step
255 d = p / q;
256 u = x + d;
257
258 if (u - a < tol2 || b - u < tol2) {
259 d = sign(tol1, xm - x);
260 }
261 }
262 } else {
263 //noinspection all
264 e = x >= xm ? a - x : b - x;
265 d = CGOLD * (e);
266 }
267
268 u = Math.abs(d) >= tol1 ? x + d : x + sign(tol1, d);
269 fu = listener.evaluate(u);
270
271 //This is the one function evaluation per iteration
272 if (fu <= fx) {
273 //No decide what to do with function evaluation
274 if (u >= x) {
275 a = x;
276 } else {
277 b = x;
278 }
279 //Housekeeping follows
280 v1[0] = v;
281 v2[0] = w;
282 v3[0] = x;
283 shift3(v1, v2, v3, u);
284 v = v1[0];
285 w = v2[0];
286 x = v3[0];
287
288 v1[0] = fv;
289 v2[0] = fw;
290 v3[0] = fx;
291 shift3(v1, v2, v3, fu);
292 fv = v1[0];
293 fw = v2[0];
294 fx = v3[0];
295 } else {
296 if (u < x) {
297 a = u;
298 } else {
299 b = u;
300 }
301 if (fu <= fw || w == x) {
302 v = w;
303 w = u;
304 fv = fw;
305 fw = fu;
306 } else if (fu <= fv || v == x || v == w) {
307 v = u;
308 fv = fu;
309 }
310 }
311
312 //Done with housekeeping. Back for another iteration.
313 if (iterationCompletedListener != null) {
314 iterationCompletedListener.onIterationCompleted(this, iter, ITMAX);
315 }
316 }
317 } catch (final EvaluationException e) {
318 throw new OptimizationException(e);
319 } finally {
320 locked = false;
321 }
322
323 //Too many iterations in Brent!
324 throw new OptimizationException();
325 }
326
327 /**
328 * Returns boolean indicating whether this instance is ready to start the
329 * estimation of a minimum or not.
330 * The instance is ready when both the listener and the bracket are
331 * available.
332 *
333 * @return True if this instance is ready, false otherwise.
334 */
335 @Override
336 public boolean isReady() {
337 return isListenerAvailable() && isBracketAvailable();
338 }
339
340 /**
341 * Internal method to set algorithm tolerance. This method does not check
342 * whether this instance is locked or not.
343 * The algorithm will iterate until the result converges below this value of
344 * accuracy or until the maximum number of iterations is achieved (and in
345 * such case, convergence will be assumed to have failed).
346 *
347 * @param tolerance Tolerance or accuracy to be obtained in estimated
348 * minimum.
349 * @throws IllegalArgumentException Raised if tolerance is negative.
350 */
351 private void internalSetTolerance(final double tolerance) {
352 if (tolerance < MIN_TOLERANCE) {
353 throw new IllegalArgumentException();
354 }
355 this.tolerance = tolerance;
356 }
357 }