Fixed bugs in "BrentOptimizer".

Renamed "BrentMinimizerTest" to "BrentOptimizerTest".
Modified "MultiStartUnivariateRealOptimizerTest" because it uses
"BrentOptimizer" as the underlying optimizer, and so was affected
by the changes.


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@979257 13f79535-47bb-0310-9956-ffa450edef68

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 optimize(f, goalType, min, max, min + GOLDEN_SECTION * (max - min));
-        return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy);
     }
 
     /**
@@ -69,23 +71,41 @@ public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
      * {@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 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}
+     * or {@link GoalType#MINIMIZE}.
-     * @param a Lower bound of the interval
+     * @param lo Lower bound of the interval.
-     * @param b Higher bound of the interval
+     * @param mid Point inside the interval {@code [lo, hi]}.
-     * @param eps Relative accuracy
+     * @param hi Higher bound of the interval.
-     * @param t Absolute accuracy
+     * @param eps Relative accuracy.
-     * @return the point at which the function is minimal.
+     * @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,
+    private double localMin(UnivariateRealFunction f,
-                            double a, double b, final double eps, final double t)
+                            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;
+            final double tol1 = eps * Math.abs(x) + t;
-            double t2 = 2 * tol;
+            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) &&
+                    if (p > q * (a - x)
-                    (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step.
+                        && p < q * (b - x)
-                    d = p / q;
+                        && Math.abs(p) < Math.abs(0.5 * q * r)) {
-                    u = x + d;
+                        // Parabolic interpolation step.
+                        d = p / q;
+                        u = x + d;
 
-                    // f must not be evaluated too close to a or b.
+                        // f must not be evaluated too close to a or b.
-                    if (((u - a) < t2) || ((b - u) < t2)) {
+                        if (u - a < tol2
-                        d = (x < m) ? tol : -tol;
+                            || 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.
+                // Update by at least "tol1".
-                u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol));
+                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);
-
     }
-
 }

src/site/xdoc/changes.xml

@@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="2.2" date="TBD" description="TBD">
+      <action dev="erans" type="fix" issue="MATH-395">
+        Fixed several bugs in "BrentOptimizer".
+      </action>
       <action dev="erans" type="fix" issue="MATH-393">
         Fixed inconsistency in return values in "MultiStartUnivariateRealOptimizer".
       </action>

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$
  */

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() > 1500);
-        assertTrue(minimizer.getEvaluations() < 3100);
+        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);
+        double result = minimizer.optimize(f, GoalType.MINIMIZE, -0.3, -0.2);
-        assertEquals(-0.27195612846834, minimizer.getResult(), 1.0e-13);
+        assertEquals(-0.27195612525275803, result, 1.0e-13);
-        assertEquals(-0.04433426954946, minimizer.getFunctionValue(), 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()    >= 300);
-        assertTrue(minimizer.getEvaluations()    <= 530);
+        assertTrue(minimizer.getEvaluations()    <= 420);
-        assertTrue(minimizer.getIterationCount() >= 150);
+        assertTrue(minimizer.getIterationCount() >= 100);
-        assertTrue(minimizer.getIterationCount() <= 170);
+        assertTrue(minimizer.getIterationCount() <= 140);
-
     }
-
 }

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 local minima at -0.27195613 and 0.82221643.
-        // The function has extrema (first derivative is zero) 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());
-
     }
-
 }