User-defined functions in PyCont

[Maurizio De Pitta']

Dear all,
I am attempting to reconstruct an f-I curve of an IF neuron. I am using the PyCont_hybrid_osc.py as a basis, but I cannot clearly understand how to specify the PCargs.userfunc argument.

I was also looking at the not reader-friendly PyCont_van_der_Pool.py example, as suggested by the tutorial. Now if I got it correctly so far, any user-defined function should take three input arguments:

 def user_funct(C, pt, pars):
 # C    : which correspond to the PyCont object such as C.model retrieves the DS object used in the integrator
 # pt   : one (or could also be more than one?!?) of the parameters of the model... which parameter however?!?
 # pars : parameters of the model

 ...
 ...
      return <array>

So in my case I would like to define a user function cost_func that detects existence of limit cycle (i.e. periodic firing) for some value of the Iapp and then automatically retrieves the frequency for that limit cycle. It seems that a possible way to do this is save the period information in C._userdata which however seems undocumented so far.

The idea would be to modify get_cycle in PyCont_hybrid_osc.py as:

def get_cycle(DS):
    DS.compute(trajname='cont_traj', force=True)
    thresh_ts = DS.getTrajEventTimes('cont_traj')['threshold']
    vthresh = DS.getTrajEvents('cont_traj')['threshold']['v']
    return array([thresh_ts[-1]-thresh_ts[-2]],float),array([abs(vthresh[-1]-vthresh[-2]],float)

And then have a cost_func like:

def cont_func(C, pt, pars):
    DS = C.model
    DS.set(pars={'Iapp': pt['Iapp'])  # Is this correct? Assuming that I want to continue for Iapp?
    try:
        F = get_cycle(DS)[1]
    except PyDSTool_ExistError:
        print 'Problem computing orbit'
        C._userdata.problem = True
        # Where can I save Period of oscillations???
        return array([0], float)
    else:
        return F

Any help would be sincerely appreciated.
Thanks for your time.

M

[Maurizio De Pitta']

Ok, let me specify a bit more this problem, with the hope that someone could help me.
First I cannot fully understand whether I can apply PyCont as it is to my HybridModel without the use of a userfunc as in the examples. If so then, I could just simply run a continuation of my model neuron, say IF_squarepulse_model and plot the frequency of oscillations of v against the bifurcation parameter I_app. This does not seem to me possible however... correct me if I am wrong. So necessarily I need to define some userfunc (or cost functions) that deals with events in my HybridModel.

Accordingly, suppose that I want to compute the firing characteristic of the IF_squarepulse_model in the examples, starting from

par_args_linear = {..., 'I_app': 0}
IFmodel = makeIFneuron('IF_fit', par_args_linear, par_args_release)
icdict = {...}
# Compute a first trajectory to seek an equilibrium point
IFmodel.compute(trajname='eqb', tdata=[0, 100], ics=icdict)

def get_cycle(DS):
    ...

def cost_func(C,pt,pars):
    ...

 # Initialize PyCont object
 IFmodel.set(ics=IFmodel.getEndPoint('eqb'))
PyCont = ContClass(IFmodel)    
PCargs.description = """User-defined continuation of period"""
PCargs.FuncTol = 1e-6
PCargs.VarTol = 1e-6
PCargs.freepars = ['fs']
PCargs.StepSize = 3e-2
PCargs.MaxStepSize = 5e-2
PCargs.MaxNumPoints = 4 # Make larger for a more serious run (but slow)
PCargs.SaveJacobian = True
PCargs.verbosity = 4

...
<will discuss below what else to put here>
...

PyCont.newCurve(PCargs)
#PyCont['UD1'].forward()
#PyCont.display(('fs','gs'))

Now for the missing part (...), looking at PyCont_hybrid_osc.py example, I understand that ideally I should use the following arguments of PCargs:

PCargs.uservars = ... # If I understood correctly this should be a list of variables that enter my 'pt' argument in cost_func
PCargs.userpars = PyCont.gensys.query('pars') # This should be the pars argument in cost_func
PCargs.userfunc = const_func # This should be my cost_func - the one that is checked every step of the continuation to seek the tolerance and move to the next step.    
PCargs.initpoint = {'fs': PyCont.gensys.query('pars')['fs']}

There are two problems, at least, that I cannot solve:
1. First I need to understand whether I am on the right track for my cost_func and the above few lines of code. Can I define I_app both as my uservars and my freepars?
2. Second, I would like not only to detect the existence of an oscillation, but also to make available the information on the period of this oscillation (i.e. firing period) so that I can plot it after, when continuation has ended.

The only example suggesting some sort of manipulation of variables is the not-easy-to-read-through PyCont_vanDerPol.py. In particular there appears to be an aregument in the PyCont class such as _userdata where apparently I can store this information on the firing period such as retrieved by my get_cycle(DS)[0] (see my previous post). But again, it is not clear what is the structure of this _userdata. Several attributes seems to be present such as .problem, .cycle etc. but it is not obvious to me how I an save the period information (or where it is saved... maybe it is in .cycle??).

Any help would be sincerely appreciated.

Thanks.

M

[Rob Clewley]

Hi,

This is Drew's code, but I think I can answer some of this, at least. Yes, you're on the right track, to the best of my understanding of this code.

_userdata is just the name for your own, completely custom, storage object for use inside your user-defined function. The special attributes of that seem to be sgn (sign for the direction of continuation) and the problem Boolean flag. In the VdP example, there are events that detect each cycle, so the period is calculated on lines 248-256 and stored in _userdata.cycle.

All attributes of _userdata are then copied into the label information for the solution Pointset stored for that PyCont[<kind>] object (let's call it sol) created by PyCont, to become sol.labels['UD']['data']. You'll find it there to be able to plot it, as the VdP tutorial example does in the function plot_cycles.

If you get some actual errors, please mention them here and email me your code so I can try to work through it further.

[Maurizio De Pitta']

Hi Rob,
thanks for the hints. Sorry I turned back with such a delay but I eventually got stucked with some other code issues (unrelated to this problem).

Following your explanations, indeed, I am able to correctly save the period in my _userdata attribute. Nonethelss I keep on failing to correctly continue the IF_squarepulse_model in order to retrieve the f-I curve.

Basically, taking as a template PyCont_hybrid_osc, I am changing a bit the code in get_cycle and cont_func in order to just detect oscillations (i.e. firing). With this aim my code reads:

from PyDSTool import *
from IF_squarespike_model import makeIFneuron

# (Modified) Starts from small I so as not to have spikes
par_args_linear = {'Iapp': 0.1, 'gl': 0.1, 'vl': -67, 'threshval': -65, 'C': 1}
par_args_spike = {'splen': 0.75}

# Build DS
IFmodel = makeIFneuron('IF_bif', par_args_linear, par_args_spike, evtol=1e-8)
# set IC right before a spike
icdict = {'v': -66, 'excited': 0}
IFmodel.set(ics=icdict, tdata=[0,100])

# (Modified) Detect oscillations and period otherwise retrieve [0,1] 
def get_cycle(DS):
    DS.compute(trajname='cont_traj', force=True)
    # Detect oscillations (i.e. at least a termination event must have occurred)
    if bool(DS.getTrajEventTimes('cont_traj')):
        thresh_ts = DS.getTrajEventTimes('cont_traj')['threshold']
        traj = DS.getTrajEvents('cont_traj')['threshold']
        return array(r_[fabs([traj['v'][-1]-traj['v'][-2]]),thresh_ts[-1]-thresh_ts[-2]], float)
    else:
        # No firing
        return array([1,0],float)

# Set up continuation class
PyCont = ContClass(IFmodel)

# (Modified) Now return difference from v peaks before threhsold crossings and save to userdata firing period 
def cont_func(C, pt, pars):
    DS = C.model
    DS.set(pars={'Iapp': pt['Iapp']})
    try:
        cycle = get_cycle(DS) 
        F = cycle[0]
        C._userdata.cycle = cycle[1]
        print pars['Iapp'],cycle
    except PyDSTool_ExistError:
        print('Problem computing orbit')
        C._userdata.problem = True
        return array([0], float)
    else:
    return F

# Effective continuation
PCargs = args(name='UD1', type='UD-C')
PCargs.description = """User-defined continuation of period"""
PCargs.uservars = ['Iapp'] #(Modified)
PCargs.userpars = PyCont.gensys.query('pars')
PCargs.userfunc = cont_func
PCargs.FuncTol = 1e-6
PCargs.VarTol = 1e-6
PCargs.freepars = ['Iapp']
# PCargs.StepSize = 3e-2
# PCargs.MaxStepSize = 5e-2
PCargs.StepSize = 5e-2
PCargs.MaxStepSize = 0.1
PCargs.MaxNumPoints = 10 # Make larger for a more serious run (but slow)
PCargs.SaveJacobian = True
PCargs.verbosity = 4
# (Modified)
PCargs.initpoint = {'Iapp': PyCont.gensys.query('pars')['Iapp']}
PyCont.newCurve(PCargs)

print('Computing curve...')
PyCont['UD1'].forward()
PyCont['UD1'].backward()

As mentioned, when I run the above, the continuation stops right after 4 steps providing the following errors:

---------------------------------------------------------------------------
PyDSTool_ExistError                       Traceback (most recent call last)
<ipython-input-102-e469728147cf> in <module>()
      3 
      4 print('Computing curve...')
----> 5 PyCont['UD1'].forward()
      6 # PyCont['UD1'].backward()
      7

/usr/local/pydstool/PyDSTool/PyCont/Continuation.pyc in forward(self)
   1244         self.CurveInfo = PointInfo()
   1245         if self.sol is None:
-> 1246             self._compute(v0=self.initdirec)
   1247             self.sol = self._curveToPointset()
   1248

/usr/local/pydstool/PyDSTool/PyCont/Continuation.pyc in _compute(self, x0, v0, direc)
   1048 
   1049             if singular:
-> 1050                 raise PyDSTool_ExistError("Problem in _compute: Failed to compute tangent vector.")
   1051 #             v0 = zeros(self.dim, float)
   1052 #             v0 = linalg.solve(r_[c_[J_coords, J_params],

PyDSTool_ExistError: 'Problem in _compute: Failed to compute tangent vector.'

Which I suppose is related to the if--else statement inside get_cycle(...) which makes impossible to start the continuation. Nonetheless I cannot avoid that conditional statement otherwise my get_cycle works only if the neuron is firing...

I am sure computing the LIF f-I curve in PyDSTool should be much more straight forward once I get the trick.

Looking forward your suggestion, or Drew's ones.

Thanks again for your time.
Cheers,

M