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.NotAvailableException;
22 import com.irurueta.numerical.NotReadyException;
23 import com.irurueta.numerical.SingleDimensionFunctionEvaluatorListener;
24
25 /**
26 * Class to compute local minimum on single dimension functions using a
27 * modification of Brent's algorithm that takes into account the function's
28 * derivative.
29 * This class will search for a local minimum within a bracket of values.
30 * A bracket is a set of points: "a" a minimum evaluation point,
31 * "b" a middle evaluation point and "c" a maximum evaluation where a <= b
32 * <= c, and where f(b) <= f(a) and f(b) <= f(c).
33 * This class is based on the implementation of Numerical Recipes 3rd ed.
34 * Section 10.4. Page 500.
35 */
36 public class DerivativeBrentSingleOptimizer extends BracketedSingleOptimizer {
37
38 /**
39 * Maximum number of iterations to perform. If convergence is not found
40 * within this number of iterations, the minimum search will be considered
41 * as failed.
42 */
43 public static final int ITMAX = 100;
44
45 /**
46 * Constant defining machine precision.
47 */
48 public static final double ZEPS = 1e-8;
49
50 /**
51 * Default tolerance. Estimated result will be found with an accuracy below
52 * or equal to provided tolerance value.
53 */
54 public static final double DEFAULT_TOLERANCE = 3e-8;
55
56 /**
57 * Minimum allowed tolerance value.
58 */
59 public static final double MIN_TOLERANCE = 0.0;
60
61 /**
62 * Listener to evaluate the functions derivative. If the function's
63 * derivative is not know (e.g. does not have a closed expression), then
64 * a DerivativeEstimator might be used inside the listener implementation.
65 */
66 private SingleDimensionFunctionEvaluatorListener derivativeListener;
67
68 /**
69 * Tolerance. Estimated result will be found with an accuracy below or equal
70 * to provided tolerance value.
71 */
72 private double tolerance;
73
74 /**
75 * Empty constructor.
76 */
77 protected DerivativeBrentSingleOptimizer() {
78 super();
79 tolerance = DEFAULT_TOLERANCE;
80 }
81
82 /**
83 * Constructor. Creates an instance with provided bracket of values and a
84 * listener to get single dimension function evaluations.
85 *
86 * @param listener Listener to evaluate a function.
87 * @param derivativeListener Listener to get function derivative.
88 * @param minEvalPoint Minimum bracket evaluation point.
89 * @param middleEvalPoint Middle bracket evaluation point.
90 * @param maxEvalPoint Maximum bracket evaluation point.
91 * @param tolerance tolerance to find result with. Estimated result will be
92 * found with an accuracy below or equal to provided tolerance value.
93 * @throws InvalidBracketRangeException Raised if the following condition is
94 * not met: minEvalPoint <= middleEvalPoint <= maxEvalPoint.
95 * @throws IllegalArgumentException Raised if tolerance is negative.
96 */
97 protected DerivativeBrentSingleOptimizer(
98 final SingleDimensionFunctionEvaluatorListener listener,
99 final SingleDimensionFunctionEvaluatorListener derivativeListener, final double minEvalPoint,
100 final double middleEvalPoint, final double maxEvalPoint, final double tolerance)
101 throws InvalidBracketRangeException {
102 super(listener, minEvalPoint, middleEvalPoint, maxEvalPoint);
103 this.derivativeListener = derivativeListener;
104 internalSetTolerance(tolerance);
105 }
106
107 /**
108 * Returns derivative listener to get function derivative.
109 *
110 * @return Derivative listener.
111 * @throws NotAvailableException Raised if derivative listener is not
112 * available for retrieval.
113 */
114 public SingleDimensionFunctionEvaluatorListener getDerivativeListener() throws NotAvailableException {
115 if (!isDerivativeListenerAvailable()) {
116 throw new NotAvailableException();
117 }
118 return derivativeListener;
119 }
120
121 /**
122 * Sets derivative listener that gets function derivative.
123 *
124 * @param derivativeListener Sets derivative listener.
125 * @throws LockedException Raised if this instance is locked.
126 */
127 public void setDerivativeListener(final SingleDimensionFunctionEvaluatorListener derivativeListener)
128 throws LockedException {
129 if (isLocked()) {
130 throw new LockedException();
131 }
132 this.derivativeListener = derivativeListener;
133 }
134
135 /**
136 * Returns boolean indicating whether derivative listener has been provided
137 * and is available for retrieval.
138 *
139 * @return Boolean indicating whether derivative listener is available.
140 */
141 public boolean isDerivativeListenerAvailable() {
142 return derivativeListener != null;
143 }
144
145 /**
146 * Returns tolerance value. Estimated result will be found with an accuracy
147 * below or equal to provided tolerance value.
148 *
149 * @return Tolerance value.
150 */
151 public double getTolerance() {
152 return tolerance;
153 }
154
155 /**
156 * Sets tolerance value. Estimated result will be found with an accuracy
157 * below or equal to provided tolerance value.
158 *
159 * @param tolerance Tolerance value.
160 * @throws LockedException Raised if this instance is locked.
161 * @throws IllegalArgumentException Raised if tolerance is negative.
162 */
163 public void setTolerance(final double tolerance) throws LockedException {
164 if (isLocked()) {
165 throw new LockedException();
166 }
167 internalSetTolerance(tolerance);
168 }
169
170 /**
171 * This function estimates a function minimum within provided or computed
172 * bracket of values.
173 * Given a function f that computes a function and also its derivative
174 * function df, and given a bracketing triplet of abscissas "ax", "bx", "cx" (such
175 * that bx is between ax and cx, and f(bx) is less than both f(ax) and
176 * f(cx), this routine isolates the minimum to a fractional precision of
177 * about tolerance using a modification of Brent's method that uses
178 * derivatives. The abscissa of the minimum is returned as "xmin" and the
179 * minimum function value is returned as "fmin".
180 *
181 * @throws LockedException Raised if this instance is locked, because
182 * estimation is being computed.
183 * @throws NotReadyException Raised if this instance is not ready because
184 * either a listener or a bracket has not yet been provided or computed.
185 * @throws OptimizationException Raised if the algorithm failed because of
186 * lack of convergence or because function couldn't be evaluated.
187 */
188 @SuppressWarnings("DuplicatedCode")
189 @Override
190 public void minimize() throws LockedException, NotReadyException, OptimizationException {
191 if (isLocked()) {
192 throw new LockedException();
193 }
194 if (!isReady()) {
195 throw new NotReadyException();
196 }
197
198 locked = true;
199
200 final var v1 = new double[1];
201 final var v2 = new double[2];
202 final var v3 = new double[3];
203
204 try {
205 // Will be used as flags for whether proposed steps are acceptable or
206 // not
207 boolean ok1;
208 boolean ok2;
209 double a;
210 double b;
211 var d = 0.0;
212 double d1;
213 double d2;
214 double du;
215 double dv;
216 double dw;
217 double dx;
218 var e = 0.0;
219 double fu;
220 double fv;
221 double fw;
222 double fx;
223 double olde;
224 double tol1;
225 double tol2;
226 double u;
227 double u1;
228 double u2;
229 double v;
230 double w;
231 double x;
232 double xm;
233
234 // Comments following will point out only differences from the Brent
235 // single optimizer. Read that routine first.
236 a = Math.min(ax, cx);
237 b = Math.max(ax, cx);
238 x = w = v = bx;
239 fw = fv = fx = listener.evaluate(x);
240 dw = dv = dx = derivativeListener.evaluate(x);
241
242 // All out housekeeping chores are doubled by the necessity of moving
243 // around derivative values as well as function values
244 for (var iter = 0; iter < ITMAX; iter++) {
245 xm = 0.5 * (a + b);
246 tol1 = tolerance * Math.abs(x) + ZEPS;
247 tol2 = 2.0 * tol1;
248 if (Math.abs(x - xm) <= (tol2 - 0.5 * (b - a))) {
249 fmin = fx;
250 xmin = x;
251
252 resultAvailable = true;
253 locked = false;
254 return;
255 }
256
257 final var tmp = dx >= 0.0 ? a - x : b - x;
258 if (Math.abs(e) > tol1) {
259 // Initialize these d's to an out-of-bracket value
260 d1 = 2.0 * (b - a);
261 d2 = d1;
262 // Secant method with one point
263 if (dw != dx) {
264 d1 = (w - x) * dx / (dx - dw);
265 }
266 // And the other
267 if (dv != dx) {
268 d2 = (v - x) * dx / (dx - dv);
269 }
270 // Which of these two estimates of d shall we take? We will
271 // insist that they be within the bracket, and on the side
272 // pointed to by the derivative at x
273 u1 = x + d1;
274 u2 = x + d2;
275 ok1 = (a - u1) * (u1 - b) > 0.0 && dx * d1 <= 0.0;
276 ok2 = (a - u2) * (u2 - b) > 0.0 && dx * d2 <= 0.0;
277 // Movement on the step before last
278 olde = e;
279 e = d;
280 if (ok1 || ok2) {
281 // Take only an acceptable d, and if both are acceptable,
282 // then take the smallest one.
283 if (ok1 && ok2) {
284 d = Math.abs(d1) < Math.abs(d2) ? d1 : d2;
285 } else if (ok1) {
286 d = d1;
287 } else {
288 d = d2;
289 }
290
291 if (Math.abs(d) <= Math.abs(0.5 * olde)) {
292 u = x + d;
293 if (u - a < tol2 || b - u < tol2) {
294 d = sign(tol1, xm - x);
295 }
296 } else {
297 // Bisect, not golden section.
298 e = tmp;
299 d = 0.5 * (e);
300 // Decide which segment by the sign of the derivative
301 }
302 } else {
303 e = tmp;
304 d = 0.5 * e;
305 }
306 } else {
307 e = tmp;
308 d = 0.5 * e;
309 }
310
311 if (Math.abs(d) >= tol1) {
312 u = x + d;
313 fu = listener.evaluate(u);
314 } else {
315 u = x + sign(tol1, d);
316 fu = listener.evaluate(u);
317 if (fu > fx) {
318 // If the minimum step in the downhill direction takes us
319 // uphill, then we are done
320 fmin = fx;
321 xmin = x;
322
323 resultAvailable = true;
324 locked = false;
325 return;
326 }
327 }
328
329 // Now all the housekeeping, sigh
330 du = derivativeListener.evaluate(u);
331 if (fu <= fx) {
332 if (u >= x) {
333 a = x;
334 } else {
335 b = x;
336 }
337 v1[0] = v;
338 v2[0] = fv;
339 v3[0] = dv;
340 mov3(v1, v2, v3, w, fw, dw);
341 v = v1[0];
342 fv = v2[0];
343 dv = v3[0];
344
345
346 v1[0] = w;
347 v2[0] = fw;
348 v3[0] = dw;
349 mov3(v1, v2, v3, x, fx, dx);
350 w = v1[0];
351 fw = v2[0];
352 dw = v3[0];
353
354 v1[0] = x;
355 v2[0] = fx;
356 v3[0] = dx;
357 mov3(v1, v2, v3, u, fu, du);
358 x = v1[0];
359 fx = v2[0];
360 dx = v3[0];
361 } else {
362 if (u < x) {
363 a = u;
364 } else {
365 b = u;
366 }
367 if (fu <= fw || w == x) {
368 v1[0] = v;
369 v2[0] = fv;
370 v3[0] = dv;
371 mov3(v1, v2, v3, w, fw, dw);
372 v = v1[0];
373 fv = v2[0];
374 dv = v3[0];
375
376 v1[0] = w;
377 v2[0] = fw;
378 v3[0] = dw;
379 mov3(v1, v2, v3, u, fu, du);
380 w = v1[0];
381 fw = v2[0];
382 dw = v3[0];
383 } else if (fu < fv || v == x || v == w) {
384 v1[0] = v;
385 v2[0] = fv;
386 v3[0] = dv;
387 mov3(v1, v2, v3, u, fu, du);
388 v = v1[0];
389 fv = v2[0];
390 dv = v3[0];
391 }
392 }
393
394 if (iterationCompletedListener != null) {
395 iterationCompletedListener.onIterationCompleted(this, iter, ITMAX);
396 }
397 }
398
399 } catch (final EvaluationException e) {
400 throw new OptimizationException(e);
401 } finally {
402 locked = false;
403 }
404
405 // Too many iterations in Derivative Brent
406 throw new OptimizationException();
407 }
408
409 /**
410 * Returns boolean indicating whether this instance is ready to start the
411 * estimation of a local minimum.
412 * This instance will be ready once a listener, derivative listener and
413 * bracket are available.
414 *
415 * @return True if ready, false otherwise.
416 */
417 @Override
418 public boolean isReady() {
419 return isListenerAvailable() && isDerivativeListenerAvailable() && isBracketAvailable();
420 }
421
422 /**
423 * Internal method to set tolerance. Estimated result will be found with an
424 * accuracy below or equal to provided tolerance value.
425 *
426 * @param tolerance Tolerance value.
427 * @throws IllegalArgumentException Raised if tolerance is negative.
428 */
429 private void internalSetTolerance(final double tolerance) {
430 if (tolerance < MIN_TOLERANCE) {
431 throw new IllegalArgumentException();
432 }
433 this.tolerance = tolerance;
434 }
435 }