[SciPy-dev] nonlin patch

[alex]

Ondrej Certik wrote:
> [...]
> 
> Just checkout the svn, copy your working scipy files over it, do "svn
> di" and send us the patch.
> 

Here is a bug-exposing patch that adds three tests to test_nonlin.py.
One of them currently fails, and the fact that the other two pass helps
to show the nature of the problem (if you use a small enough number of
iterations or start your guess far enough from the true answer, then the
bug is not triggered).

This is probably fixed by Pauli Virtanen's work, but I haven't checked this.

Alex




Index: scipy/optimize/tests/test_nonlin.py
===================================================================
--- scipy/optimize/tests/test_nonlin.py	(revision 5597)
+++ scipy/optimize/tests/test_nonlin.py	(working copy)
@@ -3,12 +3,13 @@
 May 2007
 """

+import numpy
+
 from numpy.testing import *

 from scipy.optimize import nonlin
 from numpy import matrix, diag

-
 def F(x):
     def p3(y):
         return float(y.T*y)*y
@@ -24,6 +25,40 @@

     return tuple(f.flat)

+def G(x):
+    return [v**3 - 1 for v in x]
+
+class TestNonlinEasy(TestCase):
+    """ Test case for a stupidly easy function optimization problem.
+    """
+    def broyden_helper(self, initial_guess, expected_result,
iterations=None):
+        if iterations:
+            x = nonlin.broyden2(G, initial_guess, iter=iterations)
+        else:
+            x = nonlin.broyden2(G, initial_guess)
+        x_array = numpy.array(x)
+        xout_array = numpy.array(expected_result)
+        eps = 1e-9
+        errmsg = 'got %s but expected %s' % (str(x_array), str(xout_array))
+        difference = numpy.array(x_array - xout_array)
+        assert nonlin.norm(difference) < eps, errmsg
+        assert nonlin.norm(G(x)) < eps
+
+    def test_broyden2_near_default_iterations(self):
+        initial_guess = [1.1, 1.1, 1.1]
+        expected_result = [1, 1, 1]
+        self.broyden_helper(initial_guess, expected_result)
+
+    def test_broyden2_near_fewer_iterations(self):
+        initial_guess = [1.1, 1.1, 1.1]
+        expected_result = [1, 1, 1]
+        self.broyden_helper(initial_guess, expected_result, 8)
+
+    def test_broyden2_far(self):
+        initial_guess = [2, 2, 2]
+        expected_result = [1, 1, 1]
+        self.broyden_helper(initial_guess, expected_result)
+
 class TestNonlin(TestCase):
     """ Test case for a simple constrained entropy maximization problem
     (the machine translation example of Berger et al in