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