You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by lu...@apache.org on 2007/05/08 20:16:53 UTC

svn commit: r536283 - in /jakarta/commons/proper/math/trunk: src/java/org/apache/commons/math/analysis/ src/test/org/apache/commons/math/analysis/ xdocs/

Author: luc
Date: Tue May  8 11:16:52 2007
New Revision: 536283

URL: http://svn.apache.org/viewvc?view=rev&rev=536283
Log:
[MATH-156] use initial guess provided by the caller to BrentSolver.solve(),
thus improving speed

Added:
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/MonitoredFunction.java   (with props)
Modified:
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/BrentSolver.java
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/BrentSolverTest.java
    jakarta/commons/proper/math/trunk/xdocs/changes.xml

Modified: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/BrentSolver.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/BrentSolver.java?view=diff&rev=536283&r1=536282&r2=536283
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/BrentSolver.java (original)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/BrentSolver.java Tue May  8 11:16:52 2007
@@ -43,24 +43,67 @@
     }
 
     /**
-     * Find a zero in the given interval.
-     * <p>
-     * Throws <code>ConvergenceException</code> if the values of the function
-     * at the endpoints of the interval have the same sign.
+     * Find a zero in the given interval with an initial guess.
+     * <p>Throws <code>IllegalArgumentException</code> if the values of the
+     * function at the three points have the same sign (note that it is
+     * allowed to have endpoints with the same signe if the initial point has
+     * opposite sign function-wise).</p>
      * 
      * @param min the lower bound for the interval.
      * @param max the upper bound for the interval.
-     * @param initial the start value to use (ignored).
+     * @param initial the start value to use (must be set to x0 if no
+     * initial point is known).
      * @return the value where the function is zero
-     * @throws MaxIterationsExceededException the maximum iteration count is exceeded 
+     * @throws MaxIterationsExceededException the maximum iteration count
+     * is exceeded 
      * @throws FunctionEvaluationException if an error occurs evaluating
      *  the function
      * @throws IllegalArgumentException if initial is not between min and max
+     * (even if it <em>is</em> a root)
      */
     public double solve(double min, double max, double initial)
         throws MaxIterationsExceededException, FunctionEvaluationException {
-            
-        return solve(min, max);
+
+        if (((initial - min) * (max -initial)) < 0) {
+            throw new IllegalArgumentException("Initial guess is not in search"
+                    + " interval." + "  Initial: " + initial
+                    +  "  Endpoints: [" + min + "," + max + "]");
+        }
+
+        // return the initial guess if it is good enough
+        double yInitial = f.value(initial);
+        if (Math.abs(yInitial) <= functionValueAccuracy) {
+            setResult(initial, 0);
+            return result;
+        }
+
+        // return the first endpoint if it is good enough
+        double yMin = f.value(min);
+        if (Math.abs(yMin) <= functionValueAccuracy) {
+            setResult(yMin, 0);
+            return result;
+        }
+
+        // reduce interval if min and initial bracket the root
+        if (yInitial * yMin < 0) {
+            return solve(min, yMin, initial, yInitial, min, yMin);
+        }
+
+        // return the second endpoint if it is good enough
+        double yMax = f.value(max);
+        if (Math.abs(yMax) <= functionValueAccuracy) {
+            setResult(yMax, 0);
+            return result;
+        }
+
+        // reduce interval if initial and max bracket the root
+        if (yInitial * yMax < 0) {
+            return solve(initial, yInitial, max, yMax, initial, yInitial);
+        }
+
+        // full Brent algorithm starting with provided initial guess
+        return solve(min, yMin, max, yMax, initial, yInitial);
+
     }
     
     /**
@@ -85,32 +128,50 @@
         clearResult();
         verifyInterval(min, max);
         
-        // Index 0 is the old approximation for the root.
-        // Index 1 is the last calculated approximation  for the root.
-        // Index 2 is a bracket for the root with respect to x1.
-        double x0 = min;
-        double x1 = max;
-        double y0;
-        double y1;
-        y0 = f.value(x0);
-        y1 = f.value(x1);
+        double yMin = f.value(min);
+        double yMax = f.value(max);
         
         // Verify bracketing
-        if (y0 * y1 >= 0) {
+        if (yMin * yMax >= 0) {
             throw new IllegalArgumentException
             ("Function values at endpoints do not have different signs." +
                     "  Endpoints: [" + min + "," + max + "]" + 
-                    "  Values: [" + y0 + "," + y1 + "]");       
+                    "  Values: [" + yMin + "," + yMax + "]");       
         }
-   
-        double x2 = x0;
-        double y2 = y0;
+
+        // solve using only the first endpoint as initial guess
+        return solve(min, yMin, max, yMax, min, yMin);
+
+    }
+        
+    /**
+     * Find a zero starting search according to the three provided points.
+     * @param x0 old approximation for the root
+     * @param y0 function value at the approximation for the root
+     * @param x1 last calculated approximation for the root
+     * @param y1 function value at the last calculated approximation
+     * for the root
+     * @param x2 bracket point (must be set to x0 if no bracket point is
+     * known, this will force starting with linear interpolation)
+     * @param y3 function value at the bracket point.
+     * @return the value where the function is zero
+     * @throws MaxIterationsExceededException if the maximum iteration count
+     * is exceeded
+     * @throws FunctionEvaluationException if an error occurs evaluating
+     * the function 
+     */
+    private double solve(double x0, double y0,
+                         double x1, double y1,
+                         double x2, double y2)
+    throws MaxIterationsExceededException, FunctionEvaluationException {
+
         double delta = x1 - x0;
         double oldDelta = delta;
 
         int i = 0;
         while (i < maximalIterationCount) {
             if (Math.abs(y2) < Math.abs(y1)) {
+                // use the bracket point if is better than last approximation
                 x0 = x1;
                 x1 = x2;
                 x2 = x0;

Modified: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/BrentSolverTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/BrentSolverTest.java?view=diff&rev=536283&r1=536282&r2=536283
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/BrentSolverTest.java (original)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/BrentSolverTest.java Tue May  8 11:16:52 2007
@@ -284,4 +284,48 @@
             // expected
         }
     }
+
+    public void testInitialGuess() throws MathException {
+
+        MonitoredFunction f = new MonitoredFunction(new QuinticFunction());
+        UnivariateRealSolver solver = new BrentSolver(f);
+        double result;
+
+        // no guess
+        result = solver.solve(0.6, 7.0);
+        assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
+        int referenceCallsCount = f.getCallsCount();
+        assertTrue(referenceCallsCount >= 13);
+ 
+        // invalid guess (it *is* a root, but outside of the range)
+        try {
+          result = solver.solve(0.6, 7.0, 0.0);
+          fail("an IllegalArgumentException was expected");
+        } catch (IllegalArgumentException iae) {
+            // expected behaviour
+        } catch (Exception e) {
+            fail("wrong exception caught: " + e.getMessage());
+        }
+ 
+        // bad guess
+        f.setCallsCount(0);
+        result = solver.solve(0.6, 7.0, 0.61);
+        assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
+        assertTrue(f.getCallsCount() > referenceCallsCount);
+ 
+        // good guess
+        f.setCallsCount(0);
+        result = solver.solve(0.6, 7.0, 0.999999);
+        assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
+        assertTrue(f.getCallsCount() < referenceCallsCount);
+
+        // perfect guess
+        f.setCallsCount(0);
+        result = solver.solve(0.6, 7.0, 1.0);
+        assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
+        assertEquals(0, solver.getIterationCount());
+        assertEquals(1, f.getCallsCount());
+ 
+    }
+    
 }

Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/MonitoredFunction.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/MonitoredFunction.java?view=auto&rev=536283
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/MonitoredFunction.java (added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/MonitoredFunction.java Tue May  8 11:16:52 2007
@@ -0,0 +1,49 @@
+/*
+ * 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.analysis;
+
+import org.apache.commons.math.FunctionEvaluationException;
+
+/**
+ * Wrapper class for counting functions calls.
+ *
+ * @version $Revision: 480442 $ $Date: 2006-11-29 08:21:22 +0100 (mer., 29 nov. 2006) $ 
+ */
+public class MonitoredFunction implements UnivariateRealFunction {
+
+    public MonitoredFunction(UnivariateRealFunction f) {
+        callsCount = 0;
+        this.f = f;
+    }
+
+    public void setCallsCount(int callsCount) {
+        this.callsCount = callsCount;
+    }
+
+    public int getCallsCount() {
+        return callsCount;
+    }
+
+    public double value(double x) throws FunctionEvaluationException {
+        ++callsCount;
+        return f.value(x);
+    }
+
+    private int callsCount;
+    private UnivariateRealFunction f;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/MonitoredFunction.java
------------------------------------------------------------------------------
    svn:eol-style = native

Modified: jakarta/commons/proper/math/trunk/xdocs/changes.xml
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/xdocs/changes.xml?view=diff&rev=536283&r1=536282&r2=536283
==============================================================================
--- jakarta/commons/proper/math/trunk/xdocs/changes.xml (original)
+++ jakarta/commons/proper/math/trunk/xdocs/changes.xml Tue May  8 11:16:52 2007
@@ -40,6 +40,10 @@
   </properties>
   <body>
     <release version="1.2-SNAPSHOT" date="TBD">
+      <action dev="luc" type="fix" issue="MATH-156" due-to="Tyler Ward">
+        Use the initial guess provided by the user in BrentSolver.solve(), thus
+        improving speed.
+      </action>
       <action dev="luc" type="update">
         Added the estimation optimization, geometry and ode package from the
         Mantissa library.



---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org