commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r979257 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/optimization/univariate/ site/xdoc/ test/java/org/apache/commons/math/analysis/ test/java/org/apache/commons/math/optimization/ test/java/org/apache/commons/math...
Date Mon, 26 Jul 2010 12:17:45 GMT
Author: erans
Date: Mon Jul 26 12:17:45 2010
New Revision: 979257

URL: http://svn.apache.org/viewvc?rev=979257&view=rev
Log:
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.

Added:
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java
  (with props)
Removed:
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentMinimizerTest.java
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java?rev=979257&r1=979256&r2=979257&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/univariate/BrentOptimizer.java
Mon Jul 26 12:17:45 2010
@@ -18,19 +18,20 @@ package org.apache.commons.math.optimiza
 
 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 Abst
     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 Abst
      * {@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}
-     * @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 Abst
         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 Abst
 
                     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;
 
-                    // f must not be evaluated too close to a or b.
-                    if (((u - a) < t2) || ((b - u) < t2)) {
-                        d = (x < m) ? tol : -tol;
+                    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 < 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 Abst
                     } 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 Abst
                 setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count);
                 return x;
             }
-
             ++count;
         }
-
         throw new MaxIterationsExceededException(maximalIterationCount);
-
     }
-
 }

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=979257&r1=979256&r2=979257&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Jul 26 12:17:45 2010
@@ -52,6 +52,9 @@ The <action> type attribute can be add,u
     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>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java?rev=979257&r1=979256&r2=979257&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/QuinticFunction.java
Mon Jul 26 12:17:45 2010
@@ -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$
  */

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java?rev=979257&r1=979256&r2=979257&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/MultiStartUnivariateRealOptimizerTest.java
Mon Jul 26 12:17:45 2010
@@ -48,8 +48,8 @@ public class MultiStartUnivariateRealOpt
             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 MultiStartUnivariateRealOpt
         // 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 MultiStartUnivariateRealOpt
             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 MultiStartUnivariateRealOpt
             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);
     }
-
 }

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java?rev=979257&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java
(added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java
Mon Jul 26 12:17:45 2010
@@ -0,0 +1,146 @@
+/*
+ * 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 static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
+import static org.junit.Assert.fail;
+
+import org.apache.commons.math.FunctionEvaluationException;
+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: 811685 $ $Date: 2009-09-05 19:36:48 +0200 (Sat, 05 Sep 2009) $
+ */
+public final class BrentOptimizerTest {
+
+    @Test
+    public void testSinMin() throws MathException {
+        UnivariateRealFunction f = new SinFunction();
+        UnivariateRealOptimizer minimizer = new BrentOptimizer();
+        minimizer.setMaxEvaluations(200);
+        assertEquals(200, minimizer.getMaxEvaluations());
+        try {
+            minimizer.getResult();
+            fail("an exception should have been thrown");
+        } catch (IllegalStateException ise) {
+            // expected
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        assertEquals(3 * Math.PI / 2, minimizer.optimize(f, GoalType.MINIMIZE, 4, 5), 70
* minimizer.getAbsoluteAccuracy());
+        assertTrue(minimizer.getIterationCount() <= 50);
+        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()    >=  30);
+        minimizer.setMaxEvaluations(50);
+        try {
+            minimizer.optimize(f, GoalType.MINIMIZE, 4, 5);
+            fail("an exception should have been thrown");
+        } catch (FunctionEvaluationException fee) {
+            // expected
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
+    @Test
+    public void testQuinticMin() throws MathException {
+        // 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);
+        assertEquals( 0.82221643, minimizer.optimize(f, GoalType.MINIMIZE,  0.3,  0.9), 1.0e-8);
+        assertTrue(minimizer.getIterationCount() <= 50);
+
+        // 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 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(20);
+        try {
+            minimizer.optimize(f, GoalType.MAXIMIZE, 0.2, 0.3);
+            fail("an exception should have been thrown");
+        } catch (MaxIterationsExceededException miee) {
+            // expected
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
+    @Test
+    public void testMinEndpoints() throws Exception {
+        UnivariateRealFunction f = new SinFunction();
+        UnivariateRealOptimizer solver = new BrentOptimizer();
+
+        // endpoint is minimum
+        double result = solver.optimize(f, GoalType.MINIMIZE, 3 * Math.PI / 2, 5);
+        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, 80 * solver.getAbsoluteAccuracy());
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/univariate/BrentOptimizerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native



Mime
View raw message