source: src/main/java/agents/anac/y2019/harddealer/math3/optimization/univariate/BrentOptimizer.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 10.8 KB
Line 
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 */
17package agents.anac.y2019.harddealer.math3.optimization.univariate;
18
19import agents.anac.y2019.harddealer.math3.util.Precision;
20import agents.anac.y2019.harddealer.math3.util.FastMath;
21import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
22import agents.anac.y2019.harddealer.math3.exception.NotStrictlyPositiveException;
23import agents.anac.y2019.harddealer.math3.optimization.ConvergenceChecker;
24import agents.anac.y2019.harddealer.math3.optimization.GoalType;
25
26/**
27 * For a function defined on some interval {@code (lo, hi)}, this class
28 * finds an approximation {@code x} to the point at which the function
29 * attains its minimum.
30 * It implements Richard Brent's algorithm (from his book "Algorithms for
31 * Minimization without Derivatives", p. 79) for finding minima of real
32 * univariate functions.
33 * <br/>
34 * This code is an adaptation, partly based on the Python code from SciPy
35 * (module "optimize.py" v0.5); the original algorithm is also modified
36 * <ul>
37 * <li>to use an initial guess provided by the user,</li>
38 * <li>to ensure that the best point encountered is the one returned.</li>
39 * </ul>
40 *
41 * @deprecated As of 3.1 (to be removed in 4.0).
42 * @since 2.0
43 */
44@Deprecated
45public class BrentOptimizer extends BaseAbstractUnivariateOptimizer {
46 /**
47 * Golden section.
48 */
49 private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
50 /**
51 * Minimum relative tolerance.
52 */
53 private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
54 /**
55 * Relative threshold.
56 */
57 private final double relativeThreshold;
58 /**
59 * Absolute threshold.
60 */
61 private final double absoluteThreshold;
62
63 /**
64 * The arguments are used implement the original stopping criterion
65 * of Brent's algorithm.
66 * {@code abs} and {@code rel} define a tolerance
67 * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
68 * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
69 * where <em>macheps</em> is the relative machine precision. {@code abs} must
70 * be positive.
71 *
72 * @param rel Relative threshold.
73 * @param abs Absolute threshold.
74 * @param checker Additional, user-defined, convergence checking
75 * procedure.
76 * @throws NotStrictlyPositiveException if {@code abs <= 0}.
77 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
78 */
79 public BrentOptimizer(double rel,
80 double abs,
81 ConvergenceChecker<UnivariatePointValuePair> checker) {
82 super(checker);
83
84 if (rel < MIN_RELATIVE_TOLERANCE) {
85 throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
86 }
87 if (abs <= 0) {
88 throw new NotStrictlyPositiveException(abs);
89 }
90
91 relativeThreshold = rel;
92 absoluteThreshold = abs;
93 }
94
95 /**
96 * The arguments are used for implementing the original stopping criterion
97 * of Brent's algorithm.
98 * {@code abs} and {@code rel} define a tolerance
99 * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
100 * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
101 * where <em>macheps</em> is the relative machine precision. {@code abs} must
102 * be positive.
103 *
104 * @param rel Relative threshold.
105 * @param abs Absolute threshold.
106 * @throws NotStrictlyPositiveException if {@code abs <= 0}.
107 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
108 */
109 public BrentOptimizer(double rel,
110 double abs) {
111 this(rel, abs, null);
112 }
113
114 /** {@inheritDoc} */
115 @Override
116 protected UnivariatePointValuePair doOptimize() {
117 final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
118 final double lo = getMin();
119 final double mid = getStartValue();
120 final double hi = getMax();
121
122 // Optional additional convergence criteria.
123 final ConvergenceChecker<UnivariatePointValuePair> checker
124 = getConvergenceChecker();
125
126 double a;
127 double b;
128 if (lo < hi) {
129 a = lo;
130 b = hi;
131 } else {
132 a = hi;
133 b = lo;
134 }
135
136 double x = mid;
137 double v = x;
138 double w = x;
139 double d = 0;
140 double e = 0;
141 double fx = computeObjectiveValue(x);
142 if (!isMinim) {
143 fx = -fx;
144 }
145 double fv = fx;
146 double fw = fx;
147
148 UnivariatePointValuePair previous = null;
149 UnivariatePointValuePair current
150 = new UnivariatePointValuePair(x, isMinim ? fx : -fx);
151 // Best point encountered so far (which is the initial guess).
152 UnivariatePointValuePair best = current;
153
154 int iter = 0;
155 while (true) {
156 final double m = 0.5 * (a + b);
157 final double tol1 = relativeThreshold * FastMath.abs(x) + absoluteThreshold;
158 final double tol2 = 2 * tol1;
159
160 // Default stopping criterion.
161 final boolean stop = FastMath.abs(x - m) <= tol2 - 0.5 * (b - a);
162 if (!stop) {
163 double p = 0;
164 double q = 0;
165 double r = 0;
166 double u = 0;
167
168 if (FastMath.abs(e) > tol1) { // Fit parabola.
169 r = (x - w) * (fx - fv);
170 q = (x - v) * (fx - fw);
171 p = (x - v) * q - (x - w) * r;
172 q = 2 * (q - r);
173
174 if (q > 0) {
175 p = -p;
176 } else {
177 q = -q;
178 }
179
180 r = e;
181 e = d;
182
183 if (p > q * (a - x) &&
184 p < q * (b - x) &&
185 FastMath.abs(p) < FastMath.abs(0.5 * q * r)) {
186 // Parabolic interpolation step.
187 d = p / q;
188 u = x + d;
189
190 // f must not be evaluated too close to a or b.
191 if (u - a < tol2 || b - u < tol2) {
192 if (x <= m) {
193 d = tol1;
194 } else {
195 d = -tol1;
196 }
197 }
198 } else {
199 // Golden section step.
200 if (x < m) {
201 e = b - x;
202 } else {
203 e = a - x;
204 }
205 d = GOLDEN_SECTION * e;
206 }
207 } else {
208 // Golden section step.
209 if (x < m) {
210 e = b - x;
211 } else {
212 e = a - x;
213 }
214 d = GOLDEN_SECTION * e;
215 }
216
217 // Update by at least "tol1".
218 if (FastMath.abs(d) < tol1) {
219 if (d >= 0) {
220 u = x + tol1;
221 } else {
222 u = x - tol1;
223 }
224 } else {
225 u = x + d;
226 }
227
228 double fu = computeObjectiveValue(u);
229 if (!isMinim) {
230 fu = -fu;
231 }
232
233 // User-defined convergence checker.
234 previous = current;
235 current = new UnivariatePointValuePair(u, isMinim ? fu : -fu);
236 best = best(best,
237 best(previous,
238 current,
239 isMinim),
240 isMinim);
241
242 if (checker != null && checker.converged(iter, previous, current)) {
243 return best;
244 }
245
246 // Update a, b, v, w and x.
247 if (fu <= fx) {
248 if (u < x) {
249 b = x;
250 } else {
251 a = x;
252 }
253 v = w;
254 fv = fw;
255 w = x;
256 fw = fx;
257 x = u;
258 fx = fu;
259 } else {
260 if (u < x) {
261 a = u;
262 } else {
263 b = u;
264 }
265 if (fu <= fw ||
266 Precision.equals(w, x)) {
267 v = w;
268 fv = fw;
269 w = u;
270 fw = fu;
271 } else if (fu <= fv ||
272 Precision.equals(v, x) ||
273 Precision.equals(v, w)) {
274 v = u;
275 fv = fu;
276 }
277 }
278 } else { // Default termination (Brent's criterion).
279 return best(best,
280 best(previous,
281 current,
282 isMinim),
283 isMinim);
284 }
285 ++iter;
286 }
287 }
288
289 /**
290 * Selects the best of two points.
291 *
292 * @param a Point and value.
293 * @param b Point and value.
294 * @param isMinim {@code true} if the selected point must be the one with
295 * the lowest value.
296 * @return the best point, or {@code null} if {@code a} and {@code b} are
297 * both {@code null}. When {@code a} and {@code b} have the same function
298 * value, {@code a} is returned.
299 */
300 private UnivariatePointValuePair best(UnivariatePointValuePair a,
301 UnivariatePointValuePair b,
302 boolean isMinim) {
303 if (a == null) {
304 return b;
305 }
306 if (b == null) {
307 return a;
308 }
309
310 if (isMinim) {
311 return a.getValue() <= b.getValue() ? a : b;
312 } else {
313 return a.getValue() >= b.getValue() ? a : b;
314 }
315 }
316}
Note: See TracBrowser for help on using the repository browser.