[SciPy-User] bug in signal.lsim2

josef.pktd@gmai... josef.pktd@gmai...
Thu Feb 11 08:38:46 CST 2010


On Tue, Feb 9, 2010 at 1:17 PM, Warren Weckesser
<warren.weckesser@enthought.com> wrote:
> josef.pktd@gmail.com wrote:
>> On Wed, Feb 3, 2010 at 9:00 AM, Ryan Krauss <ryanlists@gmail.com> wrote:
>>
>>> FYI, I am quite happy with passing in an hmax value.  I basically
>>> copied and pasted lsim2 from signal.ltisys and adapted it just a
>>> little to make it a method of my derived class.  Then I added the hmas
>>> kwarg that gets passed to odeint.
>>>
>>> Is there any reason not to allow the user to pass in a kwargs to lsim2
>>> that gets passed to odeint?
>>>
>>
>> I don't see a reason why we cannot add a **kwargs, it should be
>> completely backwards compatible.
>> Can you file a ticket and add your adjusted version or a patch? And
>> even better, add your original example as a test case?
>>
>>
>
> Josef,
>
> I just created ticket #1112 for this.  Unless Ryan wants to adapt his
> change to lsim2, I can make a patch this week to implement the enhancement.

Thanks, Stefan (I think) or I will look at it.

For statsmodels (maximum likelihood) I was looking at a second pattern
instead of **kw that keeps the options for the optimizer, or ode in
this case, separate.

It's just a taste difference. **kw is the more usual pattern in scipy.
Using a separate keyword for the options instead, keeps the signature
a bit cleaner (and I have been reading more gauss and matlab lately),
but it requires slightly more typing in the call.

Josef

def fun1a(a, opt1=3, **kw):
    #options to second function as keywords
    print '\nfun1a a',a
    print 'fun1a opt1', opt1
    print 'fun1a kw', kw
    fun2(a, **kw)
    return 'end fun1a'

def fun1b(a, opt1=3, optfun={}):
    #options to second function as separate dictionary
    print '\nfun1b a',a
    print 'fun1 opt1', opt1
    print 'fun1b optfun', optfun
    fun2(a, **optfun)
    return 'end fun1b'


def fun2(a, tol=0, maxiter=1):
    print 'fun2 a',a
    print 'fun2 kw: tol, maxiter', tol, maxiter
    return 'end fun2'

print fun1a('1a')
print fun1a('1a', tol=5, maxiter=10)
print fun1a('1a', opt1=4, tol=5, maxiter=10)
my_options = dict(tol=5, maxiter=10)
print fun1a('1a', opt1=4, **my_options)

print fun1a('1b')
print fun1b('1b', optfun=dict(tol=5, maxiter=10))
print fun1b('1b', opt1=4, optfun={ 'tol':5, 'maxiter':10})
my_options = dict(tol=5, maxiter=10)
print fun1b('1b', opt1=4, optfun=my_options)



>
> Warren
>
>
>> Josef
>>
>>
>>
>>> On Fri, Jan 29, 2010 at 6:44 AM, Ryan Krauss <ryanlists@gmail.com> wrote:
>>>
>>>> Thanks to Warren and Josef for their time and thoughts.  I feel like I
>>>> now understand the underlying problem and have some good options to
>>>> solve my short term issues (I assigned the project last night and they
>>>> need to be able to start working on it immediately).  I actually use a
>>>> TransferFunction class that derives from ltisys.  I could override its
>>>> lsim2 method to try out some of these solutions quickly and fairly
>>>> easily.
>>>>
>>>> Ryan
>>>>
>>>> On Thu, Jan 28, 2010 at 10:00 PM,  <josef.pktd@gmail.com> wrote:
>>>>
>>>>> On Thu, Jan 28, 2010 at 10:33 PM, Warren Weckesser
>>>>> <warren.weckesser@enthought.com> wrote:
>>>>>
>>>>>> josef.pktd@gmail.com wrote:
>>>>>>
>>>>>>> On Thu, Jan 28, 2010 at 8:50 PM, Warren Weckesser
>>>>>>> <warren.weckesser@enthought.com> wrote:
>>>>>>>
>>>>>>>
>>>>>>>> Ryan,
>>>>>>>>
>>>>>>>> The problem is that the ODE solver used by lsim2 is too good. :)
>>>>>>>>
>>>>>>>> It uses scipy.integrate.odeint, which in turn uses the Fortran library
>>>>>>>> LSODA.  Like any good solver, LSODA is an adaptive solver--it adjusts its
>>>>>>>> step size to be as large as possible while keeping estimates of the error
>>>>>>>> bounded.  For the problem you are solving, with initial condition 0, the
>>>>>>>> exact solution is initially exactly 0.  This is such a nice smooth solution
>>>>>>>> that the solver's step size quickly grows--so big, in fact, that it skips
>>>>>>>> right over your pulse and never sees it.
>>>>>>>>
>>>>>>>> So how does it create all those intermediate points at the requested time
>>>>>>>> values?  It uses interpolation between the steps that it computed to create
>>>>>>>> the solution values at the times that you requested.  So using a finer grid
>>>>>>>> of time values won't help.  (If lsim2 gave you a hook into the parameters
>>>>>>>> passed to odeint, you could set odeint's 'hmax' to a value smaller than your
>>>>>>>> pulse width, which would force the solver to see the pulse.  But there is no
>>>>>>>> way to set that parameter from lsim2.)
>>>>>>>>
>>>>>>>>
>>>>>>> It's something what I suspected. I don't know much about odeint, but
>>>>>>> do you think it would be useful to let lsim2 pass through some
>>>>>>> parameters to odeint?
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>> Sounds useful to me.  A simple implementation is an optional keyword
>>>>>> argument that is a dict of odeint arguments.   But this would almost
>>>>>> certainly break if lsim2 were ever reimplemented with a different
>>>>>> solver.  So perhaps it should allow a common set of ODE solver
>>>>>> parameters (e.g. absolute and relative error tolerances, max and min
>>>>>> step sizes, others?).
>>>>>>
>>>>>> Perhaps this should wait until after the ODE solver redesign that is
>>>>>> occasionally discussed:
>>>>>>    http://projects.scipy.org/scipy/wiki/OdeintRedesign
>>>>>> Then the solver itself could be an optional argument to lsim2.
>>>>>>
>>>>> I was just thinking of adding to the argument list a **kwds argument
>>>>> that is directly passed on to whatever ODE solver is used. This should
>>>>> be pretty flexible for any changes and be backwards compatible.
>>>>>
>>>>> I've seen and used it in a similar way for calls to optimization
>>>>> routines, e.g. also optimize.curve_fit, does it. What are actually
>>>>> valid keywords would depend on which function is called.
>>>>>
>>>>> (But I'm not a user of lsim, I'm just stealing some ideas from lti and
>>>>> friends for time series analysis.)
>>>>>
>>>>> Josef
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>> Warren
>>>>>>
>>>>>>
>>>>>>> Josef
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> The basic problem is you are passing in a discontinuous function to a solver
>>>>>>>> that expects a smooth function.  A better way to solve this problem is to
>>>>>>>> explicitly account for the discontinuity. One possibility is the attached
>>>>>>>> script.
>>>>>>>>
>>>>>>>> This is an excellent "learning opportunity" for your students on the hazards
>>>>>>>> of numerical computing!
>>>>>>>>
>>>>>>>> Warren
>>>>>>>>
>>>>>>>>
>>>>>>>> Ryan Krauss wrote:
>>>>>>>>
>>>>>>>>
>>>>>>>>> I believe I have discovered a bug in signal.lsim2.  I believe the
>>>>>>>>> short attached script illustrates the problem.  I was trying to
>>>>>>>>> predict the response of a transfer function with a pure integrator:
>>>>>>>>>
>>>>>>>>>             g
>>>>>>>>> G = -------------
>>>>>>>>>         s(s+p)
>>>>>>>>>
>>>>>>>>> to a finite width pulse.  lsim2 seems to handle the step response just
>>>>>>>>> fine, but says that the pulse response is exactly 0.0 for the entire
>>>>>>>>> time of the simulation.  Obviously, this isn't the right answer.
>>>>>>>>>
>>>>>>>>> I am running scipy 0.7.0 and numpy 1.2.1 on Ubuntu 9.04, but I also
>>>>>>>>> have the same problem on Windows running 0.7.1 and 1.4.0.
>>>>>>>>>
>>>>>>>>> Thanks,
>>>>>>>>>
>>>>>>>>> Ryan
>>>>>>>>>  ------------------------------------------------------------------------
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> SciPy-User mailing list
>>>>>>>>> SciPy-User@scipy.org
>>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>>
>>>>>>>>>
>>>>>>>> from pylab import *
>>>>>>>> from scipy import signal
>>>>>>>>
>>>>>>>>
>>>>>>>> g = 100.0
>>>>>>>> p = 15.0
>>>>>>>> G = signal.ltisys.lti(g, [1,p,0])
>>>>>>>>
>>>>>>>> t = arange(0, 1.0, 0.002)
>>>>>>>> N = len(t)
>>>>>>>>
>>>>>>>> # u for the whole interval (not used in lsim2, only for plotting later).
>>>>>>>> amp = 50.0
>>>>>>>> u = zeros(N)
>>>>>>>> k1 = 50
>>>>>>>> k2 = 100
>>>>>>>> u[k1:k2] = amp
>>>>>>>>
>>>>>>>> # Create input functions for each smooth interval. (This could be simpler,
>>>>>>>> since u
>>>>>>>> # is constant on each interval.)
>>>>>>>> a = float(k1)/N
>>>>>>>> b = float(k2)/N
>>>>>>>> T1 = linspace(0, a, 201)
>>>>>>>> u1 = zeros_like(T1)
>>>>>>>> T2 = linspace(a, b, 201)
>>>>>>>> u2 = amp*ones_like(T2)
>>>>>>>> T3 = linspace(b, 1.0, 201)
>>>>>>>> u3 = zeros_like(T3)
>>>>>>>>
>>>>>>>> # Solve on each interval; use the final value of one solution as the
>>>>>>>> starting
>>>>>>>> # point of the next solution.
>>>>>>>> # (We could skip the first calculation, since we know the solution will be
>>>>>>>> 0.)
>>>>>>>> (t1, y1, x1) = signal.lsim2(G,u1,T1)
>>>>>>>> (t2, y2, x2) = signal.lsim2(G, u2, T2, X0=x1[-1])
>>>>>>>> (t3, y3, x3) = signal.lsim2(G, u3, T3, X0=x2[-1])
>>>>>>>>
>>>>>>>> figure(1)
>>>>>>>> clf()
>>>>>>>> plot(t, u, 'k', linewidth=3)
>>>>>>>> plot(t1, y1, 'y', linewidth=3)
>>>>>>>> plot(t2, y2, 'b', linewidth=3)
>>>>>>>> plot(t3, y3, 'g', linewidth=3)
>>>>>>>>
>>>>>>>> show()
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> SciPy-User mailing list
>>>>>>>> SciPy-User@scipy.org
>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> SciPy-User mailing list
>>>>>>> SciPy-User@scipy.org
>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>
>>>>>>>
>>>>>> _______________________________________________
>>>>>> SciPy-User mailing list
>>>>>> SciPy-User@scipy.org
>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>
>>>>>>
>>>>> _______________________________________________
>>>>> SciPy-User mailing list
>>>>> SciPy-User@scipy.org
>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>
>>>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list