[SciPy-Dev] optimize.fsolve endless loop with nan

[josef.pktd at gmail.com]

On Wed, Mar 13, 2013 at 1:24 PM,  <josef.pktd at gmail.com> wrote:
> On Wed, Mar 13, 2013 at 12:59 PM,  <josef.pktd at gmail.com> wrote:
>> On Wed, Mar 13, 2013 at 12:17 PM,  <josef.pktd at gmail.com> wrote:
>>> On Wed, Mar 13, 2013 at 10:51 AM,  <josef.pktd at gmail.com> wrote:
>>>> preliminary question, I didn't have time yet to look closely
>>>>
>>>>>>> scipy.__version__
>>>> '0.9.0'
>>>>
>>>> I have a problem where fsolve goes into a range where the values are
>>>> nan. After that it goes into an endless loop, as far as I can tell.
>>>>
>>>> Something like this has been fixed for optimize.fmin_bfgs. Was there a
>>>> fix for this also for fsolve, since 0.9.0?
>>>>
>>>>
>>>> (The weirder story: I rearranged some test, and made unfortunately
>>>> also some other changes, and now when I run nosetests it never
>>>> returns. Ctrl+C kills nosetests, but leaves a python process running.
>>>> I have no clue why the test sequence should matter.)
>>>
>>> I had left the starting values in a module global even after I started
>>> to adjust them in one of the cases.
>>>
>>> The starting value for fsolve was in a range where the curvature is
>>> very flat, and fsolve made huge steps into the invalid range. After
>>> getting nans, it went AWOL.
>>>
>>> If I return np.inf as soon as I get a nan, then fsolve seems to stop
>>> right away. Is there a way to induce fsolve to stay out of the nan
>>> zone, for example returning something else than inf?
>>>
>>> I don't want to find a very specific solution, because I'm throwing
>>> lot's of different cases at the same generic method.
>>
>> same result with python 2.7, scipy version 0.11.0b1
>>
>>>"C:\Programs\Python27\python.exe" fsolve_endless_nan.py
>> scipy version 0.11.0b1
>> args 0.3 [100] 0.05
>> args 0.3 [ 100.] 0.05
>> args 0.3 [ 100.] 0.05
>> args 0.3 [ 100.00000149] 0.05
>> args 0.3 [-132.75434239] 0.05
>> fsolve_endless_nan.py:36: RuntimeWarning: invalid value encountered in sqrt
>>   pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
>> fsolve_endless_nan.py:39: RuntimeWarning: invalid value encountered in sqrt
>>   pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
>> args 0.3 [ nan] 0.05
>>
>>
>> standalone test case (from my power branch)
>>
>> Don't run in an interpreter (session) that you want to keep alive!
>> And open TaskManager if you are on Windows :)
>>
>> ------------
>> # -*- coding: utf-8 -*-
>> """Warning: endless loop in runaway process, requires hard kill of process
>>
>> Created on Wed Mar 13 12:44:15 2013
>>
>> Author: Josef Perktold
>> """
>>
>> import numpy as np
>> from scipy import stats, optimize
>>
>> import scipy
>> print "scipy version", scipy.__version__
>>
>>
>> def ttest_power(effect_size, nobs, alpha, df=None, alternative='two-sided'):
>>     '''Calculate power of a ttest
>>     '''
>>     print 'args', effect_size, nobs, alpha
>>     d = effect_size
>>     if df is None:
>>         df = nobs - 1
>>
>>     if alternative in ['two-sided', '2s']:
>>         alpha_ = alpha / 2.  #no inplace changes, doesn't work
>>     elif alternative in ['smaller', 'larger']:
>>         alpha_ = alpha
>>     else:
>>         raise ValueError("alternative has to be 'two-sided', 'larger' " +
>>                          "or 'smaller'")
>>
>>     pow_ = 0
>>     if alternative in ['two-sided', '2s', 'larger']:
>>         crit_upp = stats.t.isf(alpha_, df)
>>         # use private methods, generic methods return nan with negative d
>>         pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
>>     if alternative in ['two-sided', '2s', 'smaller']:
>>         crit_low = stats.t.ppf(alpha_, df)
>>         pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
>>     return pow_
>>
>> func = lambda nobs, *args: ttest_power(args[0], nobs, args[1])
>>
>> print optimize.fsolve(func, 100, args=(0.3, 0.05))
>> ------------
>
> correction to get a solution that would make sense (last two lines)
>
> func = lambda nobs, *args: ttest_power(args[0], nobs, args[1]) - args[2]
>
> print optimize.fsolve(func, 10, args=(0.76638635, 0.1, 0.8))
> #converges to 12
> print optimize.fsolve(func, 100, args=(0.76638635, 0.1, 0.8))      #runaway
>
> Josef
> "Lost in Translation"

(continuing my monologue)

more problems with fsolve, looks like on Ubuntu

Skipper uses scipy master, but python-xy daily testing uses
python-scipy (0.10.1+dfsg1-4)
TravisCI, which also runs Ubuntu doesn't have any problems (using
python-scipy amd64 0.9.0+dfsg1-1ubuntu2 )
https://travis-ci.org/statsmodels/statsmodels/jobs/5585028

(just realized: I've tested on Windows so far only with scipy 0.9)

from Skipper's run ( https://github.com/statsmodels/statsmodels/issues/710 ):

-----------
[~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
[67]: optimize.fsolve(func, .14, full_output=True)
[67]:
(array([ 0.05]),
 {'fjac': array([[-1.]]),
  'fvec': array([ 0.]),
  'nfev': 12,
  'qtf': array([ 0.]),
  'r': array([-3.408293])},
 1,
 'The solution converged.')

[~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
[68]: optimize.fsolve(func, .15, full_output=True)
[68]:
(array([ 0.15]),
 {'fjac': array([[-1.]]),
  'fvec': array([ 0.20178728]),
  'nfev': 4,
  'qtf': array([-0.20178728]),
  'r': array([ inf])},
 1,
 'The solution converged.')

It looks like the QR factorization is failing (r == inf) and then it's
reporting convergence still.
---------

That it stops with a "solution converged" also doesn't trigger my
backup root-finding.


I can work around this since I have no idea how to debug this.
However, there might be something fishy with fsolve.

I never had problems with fsolve before, and scipy.stats.distributions
is a heavy user of it.

Josef

>
>>
>>>
>>> Josef
>>>
>>>>
>>>> Josef