diff --git a/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java b/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java index 097c17c2c..47409a8e8 100644 --- a/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java @@ -18,19 +18,20 @@ 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.analysis.UnivariateRealFunction; import org.apache.commons.math.optimization.GoalType; /** * Implements Richard Brent's algorithm (from his book "Algorithms for * Minimization without Derivatives", p. 79) for finding minima of real - * univariate functions. + * univariate functions. This implementation is an adaptation partly + * based on the Python code from SciPy (module "optimize.py" v0.5). * * @version $Revision$ $Date$ * @since 2.0 */ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { - /** * Golden section. */ @@ -47,15 +48,16 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { public double optimize(final UnivariateRealFunction f, final GoalType goalType, final double min, final double max, final double startValue) throws MaxIterationsExceededException, FunctionEvaluationException { - return optimize(f, goalType, min, max); + clearResult(); + return localMin(f, goalType, min, startValue, max, + getRelativeAccuracy(), getAbsoluteAccuracy()); } /** {@inheritDoc} */ public double optimize(final UnivariateRealFunction f, final GoalType goalType, final double min, final double max) throws MaxIterationsExceededException, FunctionEvaluationException { - clearResult(); - return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy); + return optimize(f, goalType, min, max, min + GOLDEN_SECTION * (max - min)); } /** @@ -69,23 +71,41 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { * {@code eps} should be no smaller than 2 macheps and preferable not * much less than sqrt(macheps), where macheps is the relative * machine precision. {@code t} should be positive. - * @param f the function to solve + * @param f the function to solve. * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} - * or {@link GoalType#MINIMIZE} - * @param a Lower bound of the interval - * @param b Higher bound of the interval - * @param eps Relative accuracy - * @param t Absolute accuracy - * @return the point at which the function is minimal. + * or {@link GoalType#MINIMIZE}. + * @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(final UnivariateRealFunction f, final GoalType goalType, - double a, double b, final double eps, final double t) + private double localMin(UnivariateRealFunction f, + GoalType goalType, + double lo, double mid, double hi, + double eps, double t) throws MaxIterationsExceededException, FunctionEvaluationException { - double x = a + GOLDEN_SECTION * (b - a); + if (eps <= 0) { + throw new NotStrictlyPositiveException(eps); + } + if (t <= 0) { + throw new NotStrictlyPositiveException(t); + } + double a, b; + if (lo < hi) { + a = lo; + b = hi; + } else { + a = hi; + b = lo; + } + + double x = mid; double v = x; double w = x; double e = 0; @@ -99,18 +119,18 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { int count = 0; while (count < maximalIterationCount) { double m = 0.5 * (a + b); - double tol = eps * Math.abs(x) + t; - double t2 = 2 * tol; + final double tol1 = eps * Math.abs(x) + t; + final double tol2 = 2 * tol1; // Check stopping criterion. - if (Math.abs(x - m) > t2 - 0.5 * (b - a)) { + if (Math.abs(x - m) > tol2 - 0.5 * (b - a)) { double p = 0; double q = 0; double r = 0; double d = 0; double u = 0; - if (Math.abs(e) > tol) { // Fit parabola. + if (Math.abs(e) > tol1) { // Fit parabola. r = (x - w) * (fx - fv); q = (x - v) * (fx - fw); p = (x - v) * q - (x - w) * r; @@ -124,24 +144,53 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { r = e; e = d; - } - if (Math.abs(p) < Math.abs(0.5 * q * r) && - (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step. - d = p / q; - u = x + d; + if (p > q * (a - x) + && p < q * (b - x) + && Math.abs(p) < Math.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) < t2) || ((b - u) < t2)) { - d = (x < m) ? tol : -tol; + // 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; } - } else { // Golden section step. - e = ((x < m) ? b : a) - x; d = GOLDEN_SECTION * e; } - // f must not be evaluated too close to a or b. - u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol)); + // Update by at least "tol1". + if (Math.abs(d) < tol1) { + if (d >= 0) { + u = x + tol1; + } else { + u = x - tol1; + } + } else { + u = x + d; + } + double fu = computeObjectiveValue(f, u); if (goalType == GoalType.MAXIMIZE) { fu = -fu; @@ -166,12 +215,15 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { } else { b = u; } - if ((fu <= fw) || (w == x)) { + if (fu <= fw + || w == x) { v = w; fv = fw; w = u; fw = fu; - } else if ((fu <= fv) || (v == x) || (v == w)) { + } else if (fu <= fv + || v == x + || v == w) { v = u; fv = fu; } @@ -180,12 +232,8 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer { setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count); return x; } - ++count; } - throw new MaxIterationsExceededException(maximalIterationCount); - } - } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 07be4a457..d3390e02b 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,9 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Fixed several bugs in "BrentOptimizer". + Fixed inconsistency in return values in "MultiStartUnivariateRealOptimizer". diff --git a/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java b/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java index f29330228..73dd9dced 100644 --- a/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java +++ b/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java @@ -19,7 +19,7 @@ package org.apache.commons.math.analysis; import org.apache.commons.math.FunctionEvaluationException; /** - * Auxillary class for testing solvers. + * Auxiliary class for testing solvers. * * @version $Revision$ $Date$ */ diff --git a/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java b/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java index e066703b6..159588d1e 100644 --- a/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java +++ b/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java @@ -48,8 +48,8 @@ public class MultiStartUnivariateRealOptimizerTest { assertEquals(-1.0, f.value(optima[i]), 1.0e-10); assertEquals(f.value(optima[i]), optimaValues[i], 1.0e-10); } - assertTrue(minimizer.getEvaluations() > 2900); - assertTrue(minimizer.getEvaluations() < 3100); + assertTrue(minimizer.getEvaluations() > 1500); + assertTrue(minimizer.getEvaluations() < 1700); } @Test @@ -58,8 +58,9 @@ public class MultiStartUnivariateRealOptimizerTest { // The function has extrema (first derivative is zero) at 0.27195613 and 0.82221643, UnivariateRealFunction f = new QuinticFunction(); UnivariateRealOptimizer underlying = new BrentOptimizer(); + underlying.setRelativeAccuracy(1e-15); JDKRandomGenerator g = new JDKRandomGenerator(); - g.setSeed(4312000053l); + g.setSeed(4312000053L); MultiStartUnivariateRealOptimizer minimizer = new MultiStartUnivariateRealOptimizer(underlying, 5, g); minimizer.setAbsoluteAccuracy(10 * minimizer.getAbsoluteAccuracy()); @@ -82,9 +83,10 @@ public class MultiStartUnivariateRealOptimizerTest { fail("wrong exception caught"); } - assertEquals(-0.27195612846834, minimizer.optimize(f, GoalType.MINIMIZE, -0.3, -0.2), 1.0e-13); - assertEquals(-0.27195612846834, minimizer.getResult(), 1.0e-13); - assertEquals(-0.04433426954946, minimizer.getFunctionValue(), 1.0e-13); + double result = minimizer.optimize(f, GoalType.MINIMIZE, -0.3, -0.2); + assertEquals(-0.27195612525275803, result, 1.0e-13); + assertEquals(-0.27195612525275803, minimizer.getResult(), 1.0e-13); + assertEquals(-0.04433426954946637, minimizer.getFunctionValue(), 1.0e-13); double[] optima = minimizer.getOptima(); double[] optimaValues = minimizer.getOptimaValues(); @@ -92,11 +94,9 @@ public class MultiStartUnivariateRealOptimizerTest { assertEquals(f.value(optima[i]), optimaValues[i], 1.0e-10); } - assertTrue(minimizer.getEvaluations() >= 510); - assertTrue(minimizer.getEvaluations() <= 530); - assertTrue(minimizer.getIterationCount() >= 150); - assertTrue(minimizer.getIterationCount() <= 170); - + assertTrue(minimizer.getEvaluations() >= 300); + assertTrue(minimizer.getEvaluations() <= 420); + assertTrue(minimizer.getIterationCount() >= 100); + assertTrue(minimizer.getIterationCount() <= 140); } - } diff --git a/src/test/java/org/apache/commons/math/optimization/univariate/BrentMinimizerTest.java b/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java similarity index 71% rename from src/test/java/org/apache/commons/math/optimization/univariate/BrentMinimizerTest.java rename to src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java index 0358de428..90091fbb0 100644 --- a/src/test/java/org/apache/commons/math/optimization/univariate/BrentMinimizerTest.java +++ b/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java @@ -25,15 +25,16 @@ import org.apache.commons.math.MathException; import org.apache.commons.math.MaxIterationsExceededException; import org.apache.commons.math.analysis.QuinticFunction; import org.apache.commons.math.analysis.SinFunction; +import org.apache.commons.math.analysis.SincFunction; import org.apache.commons.math.analysis.UnivariateRealFunction; import org.apache.commons.math.optimization.GoalType; import org.apache.commons.math.optimization.UnivariateRealOptimizer; import org.junit.Test; /** - * @version $Revision$ $Date$ + * @version $Revision: 811685 $ $Date: 2009-09-05 19:36:48 +0200 (Sat, 05 Sep 2009) $ */ -public final class BrentMinimizerTest { +public final class BrentOptimizerTest { @Test public void testSinMin() throws MathException { @@ -54,7 +55,7 @@ public final class BrentMinimizerTest { assertEquals(3 * Math.PI / 2, minimizer.optimize(f, GoalType.MINIMIZE, 1, 5), 70 * minimizer.getAbsoluteAccuracy()); assertTrue(minimizer.getIterationCount() <= 50); assertTrue(minimizer.getEvaluations() <= 100); - assertTrue(minimizer.getEvaluations() >= 90); + assertTrue(minimizer.getEvaluations() >= 30); minimizer.setMaxEvaluations(50); try { minimizer.optimize(f, GoalType.MINIMIZE, 4, 5); @@ -68,8 +69,7 @@ public final class BrentMinimizerTest { @Test public void testQuinticMin() throws MathException { - // The quintic function has zeros at 0, +-0.5 and +-1. - // The function has extrema (first derivative is zero) at 0.27195613 and 0.82221643, + // The function has local minima at -0.27195613 and 0.82221643. UnivariateRealFunction f = new QuinticFunction(); UnivariateRealOptimizer minimizer = new BrentOptimizer(); assertEquals(-0.27195613, minimizer.optimize(f, GoalType.MINIMIZE, -0.3, -0.2), 1.0e-8); @@ -79,17 +79,48 @@ public final class BrentMinimizerTest { // search in a large interval assertEquals(-0.27195613, minimizer.optimize(f, GoalType.MINIMIZE, -1.0, 0.2), 1.0e-8); assertTrue(minimizer.getIterationCount() <= 50); + } + @Test + public void testQuinticMinPythonComparison() throws MathException { + // The function has local minima at -0.27195613 and 0.82221643. + UnivariateRealFunction f = new QuinticFunction(); + UnivariateRealOptimizer minimizer = new BrentOptimizer(); + minimizer.setRelativeAccuracy(1e-12); + minimizer.setAbsoluteAccuracy(1e-11); + + double result; + int nIter, nEval; + + result = minimizer.optimize(f, GoalType.MINIMIZE, -0.3, -0.2, -0.25); + nIter = minimizer.getIterationCount(); + nEval = minimizer.getEvaluations(); + // XXX Python: -0.27195612805911351 (instead of -0.2719561279558559). + assertEquals(-0.2719561279558559, result, 1e-12); + // XXX Python: 15 (instead of 18). + assertEquals(18, nEval); + // XXX Python: 11 (instead of 17). + assertEquals(17, nIter); + + result = minimizer.optimize(f, GoalType.MINIMIZE, 0.7, 0.9, 0.8); + nIter = minimizer.getIterationCount(); + nEval = minimizer.getEvaluations(); + // XXX Python: 0.82221643488363705 (instead of 0.8222164326561908). + assertEquals(0.8222164326561908, result, 1e-12); + // XXX Python: 25 (instead of 43). + assertEquals(43, nEval); + // XXX Python: 21 (instead of 24). + assertEquals(24, nIter); } @Test public void testQuinticMax() throws MathException { // The quintic function has zeros at 0, +-0.5 and +-1. - // The function has extrema (first derivative is zero) at 0.27195613 and 0.82221643, + // The function has a local maximum at 0.27195613. UnivariateRealFunction f = new QuinticFunction(); UnivariateRealOptimizer minimizer = new BrentOptimizer(); assertEquals(0.27195613, minimizer.optimize(f, GoalType.MAXIMIZE, 0.2, 0.3), 1.0e-8); - minimizer.setMaximalIterationCount(30); + minimizer.setMaximalIterationCount(20); try { minimizer.optimize(f, GoalType.MAXIMIZE, 0.2, 0.3); fail("an exception should have been thrown"); @@ -110,8 +141,6 @@ public final class BrentMinimizerTest { assertEquals(3 * Math.PI / 2, result, 70 * solver.getAbsoluteAccuracy()); result = solver.optimize(f, GoalType.MINIMIZE, 4, 3 * Math.PI / 2); - assertEquals(3 * Math.PI / 2, result, 70 * solver.getAbsoluteAccuracy()); - + assertEquals(3 * Math.PI / 2, result, 80 * solver.getAbsoluteAccuracy()); } - }