| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math4.legacy.optim.univariate; |
| |
| import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException; |
| import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; |
| import org.apache.commons.math4.legacy.optim.ConvergenceChecker; |
| import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType; |
| import org.apache.commons.math4.core.jdkmath.JdkMath; |
| import org.apache.commons.numbers.core.Precision; |
| |
| /** |
| * For a function defined on some interval {@code (lo, hi)}, this class |
| * finds an approximation {@code x} to the point at which the function |
| * attains its minimum. |
| * It implements Richard Brent's algorithm (from his book "Algorithms for |
| * Minimization without Derivatives", p. 79) for finding minima of real |
| * univariate functions. |
| * <br> |
| * This code is an adaptation, partly based on the Python code from SciPy |
| * (module "optimize.py" v0.5); the original algorithm is also modified |
| * <ul> |
| * <li>to use an initial guess provided by the user,</li> |
| * <li>to ensure that the best point encountered is the one returned.</li> |
| * </ul> |
| * |
| * @since 2.0 |
| */ |
| public class BrentOptimizer extends UnivariateOptimizer { |
| /** |
| * Golden section. |
| */ |
| private static final double GOLDEN_SECTION = 0.5 * (3 - JdkMath.sqrt(5)); |
| /** |
| * Minimum relative tolerance. |
| */ |
| private static final double MIN_RELATIVE_TOLERANCE = 2 * JdkMath.ulp(1d); |
| /** |
| * Relative threshold. |
| */ |
| private final double relativeThreshold; |
| /** |
| * Absolute threshold. |
| */ |
| private final double absoluteThreshold; |
| |
| /** |
| * The arguments are used implement the original stopping criterion |
| * of Brent's algorithm. |
| * {@code abs} and {@code rel} define a tolerance |
| * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than |
| * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>, |
| * where <em>macheps</em> is the relative machine precision. {@code abs} must |
| * be positive. |
| * |
| * @param rel Relative threshold. |
| * @param abs Absolute threshold. |
| * @param checker Additional, user-defined, convergence checking |
| * procedure. |
| * @throws NotStrictlyPositiveException if {@code abs <= 0}. |
| * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. |
| */ |
| public BrentOptimizer(double rel, |
| double abs, |
| ConvergenceChecker<UnivariatePointValuePair> checker) { |
| super(checker); |
| |
| if (rel < MIN_RELATIVE_TOLERANCE) { |
| throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true); |
| } |
| if (abs <= 0) { |
| throw new NotStrictlyPositiveException(abs); |
| } |
| |
| relativeThreshold = rel; |
| absoluteThreshold = abs; |
| } |
| |
| /** |
| * The arguments are used for implementing the original stopping criterion |
| * of Brent's algorithm. |
| * {@code abs} and {@code rel} define a tolerance |
| * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than |
| * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>, |
| * where <em>macheps</em> is the relative machine precision. {@code abs} must |
| * be positive. |
| * |
| * @param rel Relative threshold. |
| * @param abs Absolute threshold. |
| * @throws NotStrictlyPositiveException if {@code abs <= 0}. |
| * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. |
| */ |
| public BrentOptimizer(double rel, |
| double abs) { |
| this(rel, abs, null); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| protected UnivariatePointValuePair doOptimize() { |
| final boolean isMinim = getGoalType() == GoalType.MINIMIZE; |
| final double lo = getMin(); |
| final double mid = getStartValue(); |
| final double hi = getMax(); |
| |
| // Optional additional convergence criteria. |
| final ConvergenceChecker<UnivariatePointValuePair> checker |
| = getConvergenceChecker(); |
| |
| double a; |
| double b; |
| if (lo < hi) { |
| a = lo; |
| b = hi; |
| } else { |
| a = hi; |
| b = lo; |
| } |
| |
| double x = mid; |
| double v = x; |
| double w = x; |
| double d = 0; |
| double e = 0; |
| double fx = computeObjectiveValue(x); |
| if (!isMinim) { |
| fx = -fx; |
| } |
| double fv = fx; |
| double fw = fx; |
| |
| UnivariatePointValuePair previous = null; |
| UnivariatePointValuePair current |
| = new UnivariatePointValuePair(x, isMinim ? fx : -fx); |
| // Best point encountered so far (which is the initial guess). |
| UnivariatePointValuePair best = current; |
| |
| while (true) { |
| final double m = 0.5 * (a + b); |
| final double tol1 = relativeThreshold * JdkMath.abs(x) + absoluteThreshold; |
| final double tol2 = 2 * tol1; |
| |
| // Default stopping criterion. |
| final boolean stop = JdkMath.abs(x - m) <= tol2 - 0.5 * (b - a); |
| if (!stop) { |
| double p = 0; |
| double q = 0; |
| double r = 0; |
| double u = 0; |
| |
| if (JdkMath.abs(e) > tol1) { // Fit parabola. |
| r = (x - w) * (fx - fv); |
| q = (x - v) * (fx - fw); |
| p = (x - v) * q - (x - w) * r; |
| q = 2 * (q - r); |
| |
| if (q > 0) { |
| p = -p; |
| } else { |
| q = -q; |
| } |
| |
| r = e; |
| e = d; |
| |
| if (p > q * (a - x) && |
| p < q * (b - x) && |
| JdkMath.abs(p) < JdkMath.abs(0.5 * q * r)) { |
| // Parabolic interpolation step. |
| d = p / q; |
| u = x + d; |
| |
| // f must not be evaluated too close to a or b. |
| if (u - a < tol2 || b - u < tol2) { |
| if (x <= m) { |
| d = tol1; |
| } else { |
| d = -tol1; |
| } |
| } |
| } else { |
| // Golden section step. |
| if (x < m) { |
| e = b - x; |
| } else { |
| e = a - x; |
| } |
| d = GOLDEN_SECTION * e; |
| } |
| } else { |
| // Golden section step. |
| if (x < m) { |
| e = b - x; |
| } else { |
| e = a - x; |
| } |
| d = GOLDEN_SECTION * e; |
| } |
| |
| // Update by at least "tol1". |
| if (JdkMath.abs(d) < tol1) { |
| if (d >= 0) { |
| u = x + tol1; |
| } else { |
| u = x - tol1; |
| } |
| } else { |
| u = x + d; |
| } |
| |
| double fu = computeObjectiveValue(u); |
| if (!isMinim) { |
| fu = -fu; |
| } |
| |
| // User-defined convergence checker. |
| previous = current; |
| current = new UnivariatePointValuePair(u, isMinim ? fu : -fu); |
| best = best(best, |
| best(previous, |
| current, |
| isMinim), |
| isMinim); |
| |
| if (checker != null && checker.converged(getIterations(), previous, current)) { |
| return best; |
| } |
| |
| // Update a, b, v, w and x. |
| if (fu <= fx) { |
| if (u < x) { |
| b = x; |
| } else { |
| a = x; |
| } |
| v = w; |
| fv = fw; |
| w = x; |
| fw = fx; |
| x = u; |
| fx = fu; |
| } else { |
| if (u < x) { |
| a = u; |
| } else { |
| b = u; |
| } |
| if (fu <= fw || |
| Precision.equals(w, x)) { |
| v = w; |
| fv = fw; |
| w = u; |
| fw = fu; |
| } else if (fu <= fv || |
| Precision.equals(v, x) || |
| Precision.equals(v, w)) { |
| v = u; |
| fv = fu; |
| } |
| } |
| } else { // Default termination (Brent's criterion). |
| return best(best, |
| best(previous, |
| current, |
| isMinim), |
| isMinim); |
| } |
| |
| incrementIterationCount(); |
| } |
| } |
| |
| /** |
| * Selects the best of two points. |
| * |
| * @param a Point and value. |
| * @param b Point and value. |
| * @param isMinim {@code true} if the selected point must be the one with |
| * the lowest value. |
| * @return the best point, or {@code null} if {@code a} and {@code b} are |
| * both {@code null}. When {@code a} and {@code b} have the same function |
| * value, {@code a} is returned. |
| */ |
| private UnivariatePointValuePair best(UnivariatePointValuePair a, |
| UnivariatePointValuePair b, |
| boolean isMinim) { |
| if (a == null) { |
| return b; |
| } |
| if (b == null) { |
| return a; |
| } |
| |
| if (isMinim) { |
| return a.getValue() <= b.getValue() ? a : b; |
| } else { |
| return a.getValue() >= b.getValue() ? a : b; |
| } |
| } |
| } |
