| /* |
| * 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.math.optimization.univariate; |
| |
| import org.apache.commons.math.FunctionEvaluationException; |
| import org.apache.commons.math.MaxIterationsExceededException; |
| import org.apache.commons.math.exception.NotStrictlyPositiveException; |
| import org.apache.commons.math.optimization.GoalType; |
| import org.apache.commons.math.util.FastMath; |
| |
| /** |
| * Implements Richard Brent's algorithm (from his book "Algorithms for |
| * Minimization without Derivatives", p. 79) for finding minima of real |
| * univariate functions. This implementation is an adaptation partly |
| * based on the Python code from SciPy (module "optimize.py" v0.5). |
| * |
| * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ |
| * @since 2.0 |
| */ |
| public class BrentOptimizer extends AbstractUnivariateRealOptimizer { |
| /** |
| * Golden section. |
| */ |
| private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5)); |
| |
| /** |
| * Construct a solver. |
| */ |
| public BrentOptimizer() { |
| setMaxEvaluations(1000); |
| setMaximalIterationCount(100); |
| setAbsoluteAccuracy(1e-11); |
| setRelativeAccuracy(1e-9); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| protected double doOptimize() |
| throws MaxIterationsExceededException, FunctionEvaluationException { |
| return localMin(getGoalType() == GoalType.MINIMIZE, |
| getMin(), getStartValue(), getMax(), |
| getRelativeAccuracy(), getAbsoluteAccuracy()); |
| } |
| |
| /** |
| * Find the minimum of the function within the interval {@code (lo, hi)}. |
| * |
| * If the function is defined on the interval {@code (lo, hi)}, then |
| * this method finds an approximation {@code x} to the point at which |
| * the function attains its minimum.<br/> |
| * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} |
| * and the function is never evaluated at two points closer together than |
| * {@code tol}. {@code eps} should be no smaller than <em>2 macheps</em> and |
| * preferable not much less than <em>sqrt(macheps)</em>, where |
| * <em>macheps</em> is the relative machine precision. {@code t} should be |
| * positive. |
| * @param isMinim {@code true} when minimizing the function. |
| * @param lo Lower bound of the interval. |
| * @param mid Point inside the interval {@code [lo, hi]}. |
| * @param hi Higher bound of the interval. |
| * @param eps Relative accuracy. |
| * @param t Absolute accuracy. |
| * @return the optimum point. |
| * @throws MaxIterationsExceededException if the maximum iteration count |
| * is exceeded. |
| * @throws FunctionEvaluationException if an error occurs evaluating the function. |
| */ |
| private double localMin(boolean isMinim, |
| double lo, double mid, double hi, |
| double eps, double t) |
| throws MaxIterationsExceededException, FunctionEvaluationException { |
| if (eps <= 0) { |
| throw new NotStrictlyPositiveException(eps); |
| } |
| if (t <= 0) { |
| throw new NotStrictlyPositiveException(t); |
| } |
| 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; |
| |
| while (true) { |
| double m = 0.5 * (a + b); |
| final double tol1 = eps * FastMath.abs(x) + t; |
| final double tol2 = 2 * tol1; |
| |
| // Check stopping criterion. |
| if (FastMath.abs(x - m) > tol2 - 0.5 * (b - a)) { |
| double p = 0; |
| double q = 0; |
| double r = 0; |
| double u = 0; |
| |
| if (FastMath.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) && |
| FastMath.abs(p) < FastMath.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 (FastMath.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; |
| } |
| |
| // 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 || w == x) { |
| v = w; |
| fv = fw; |
| w = u; |
| fw = fu; |
| } else if (fu <= fv || v == x || v == w) { |
| v = u; |
| fv = fu; |
| } |
| } |
| } else { // termination |
| setFunctionValue(isMinim ? fx : -fx); |
| return x; |
| } |
| incrementIterationsCounter(); |
| } |
| } |
| } |