Hi Ondrej, As requested I have attached a modification and a test case. The modification checks a value for equality to zero, but maybe it should more appropriately handle values near zero. Also, other functions in nonlin.py may benefit from similar modifications. Alex Ondrej Certik wrote:
On Wed, Mar 19, 2008 at 9:30 PM, argriffi <argriffi@ncsu.edu> wrote:
Hi Ondrej,
Running the attached script should demonstrate the problem I mentioned. I removed dependencies on other modules I had written, but let me know if you want it reduced even further.
Thanks, that is enough. I can reproduce the problem now:
------------------
test_solver_20 (__main__.TestCodonInversion) ... ok test_solver_40 (__main__.TestCodonInversion) ... ERROR
====================================================================== ERROR: test_solver_40 (__main__.TestCodonInversion) ---------------------------------------------------------------------- Traceback (most recent call last): File "CodonInversion.py", line 181, in test_solver_40 self._test_solver_helper(40) File "CodonInversion.py", line 173, in _test_solver_helper raise RuntimeError('failed on iteration %d of %d: %s' % (len(objective_function.guesses), iterations, str(e))) RuntimeError: failed on iteration 22 of 40: detected possibly invalid input to the objective function: (nan, nan, nan)
---------------------------------------------------------------------- Ran 2 tests in 0.029s
FAILED (errors=1)
-------------------------
And the reason is simple. If you look at this code:
for n in range(iter): deltaxm=-Gm*Fxm xm=xm+deltaxm Fxm1=myF(F,xm) deltaFxm=Fxm1-Fxm Fxm=Fxm1 Gm=Gm+(deltaxm-Gm*deltaFxm)*deltaFxm.T/norm(deltaFxm)**2
And if you debug it, as I just did (for example using winpdb), you'll notice, that deltaFxm gets zero, after a lot of iterations, because you converged to a solution already. Then norm(deltaFxm) is 0 and after dividing by 0, you get a NaN in Gm, thus you get NaN in deltaxm and xm and you get NaN as an input into your function and here we are.
The solution is to check that deltaFxm is near zero and stop iterating. Do you think you could please provide such a patch? :)
Also with a test case.
Thanks, Ondrej
P.S. I CCed scipy-dev as I think you forgot to forward your email to it. I hope it's ok.
import unittest import scipy.optimize import numpy import math def easy_function(X): return tuple(x*x - 1 for x in X) class MyTest(unittest.TestCase): def helper(self, solve): guess = (0, 0, 0) iterations = 20 result = solve(easy_function, guess, iterations) for x in result: self.failUnless(x in (-1.0, 1.0), x) def test_current(self): "The current version of broyden2 fails." self.helper(scipy.optimize.nonlin.broyden2) def test_modified(self): "A modified version of broyden2 does not fail so hard." def norm(v): return math.sqrt(numpy.sum((numpy.array(v)**2).flat)) def myF(F,xm): return numpy.matrix(F(tuple(xm.flat))).T def broyden2(F, xin, iter=10, alpha=0.4, verbose = False): xm=numpy.matrix(xin).T Fxm=myF(F,xm) Gm=-alpha*numpy.matrix(numpy.identity(len(xin))) for n in range(iter): deltaxm=-Gm*Fxm xm=xm+deltaxm Fxm1=myF(F,xm) # begin modification deltaFxm=Fxm1-Fxm denominator = norm(deltaFxm)**2 if not denominator: break Fxm=Fxm1 Gm=Gm+(deltaxm-Gm*deltaFxm)*deltaFxm.T/denominator # end modification if verbose: print "%d: |F(x)|=%.3f"%(n+1, norm(Fxm)) return xm.flat self.helper(broyden2) unittest.main()