Changeset 337

Show
Ignore:
Timestamp:
11/13/08 15:41:12 (2 months ago)
Author:
emuller
Message:

shotnoise_fromspikes fixed
AnalogSignal?.time_slice was changed to not assume t_start=0.0

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/src/signals.py

    r321 r337  
    22882288            t_stop  - end of the new SpikeTrain, in ms. 
    22892289        """ 
    2290         signal = self.signal[numpy.floor(t_start/self.dt):numpy.floor(t_stop/self.dt)] 
    2291         return AnalogSignal(signal, self.dt, t_start, t_stop) 
     2290        assert t_start>self.t_start 
     2291        assert t_stop<=self.t_stop 
     2292        assert t_stop>t_start 
     2293 
     2294        t = self.time_axis() 
     2295        i_start = numpy.searchsorted(t,t_start,'right')-1 
     2296        i_stop = numpy.searchsorted(t,t_stop,'right')-1 
     2297 
     2298        signal = self.signal[i_start:i_stop] 
     2299        result = AnalogSignal(signal, self.dt, 0.0, t_stop-t_start) 
     2300        result.time_offset(t_start) 
     2301        return result 
    22922302 
    22932303    def threshold_detection(self, threshold=None, format=None): 

  • trunk/src/stgen.py

    r330 r337  
    3838        a   - dimensionless 
    3939        b   - in units of seconds 
     40 
     41    See also: 
     42        inh_gamma_generator 
    4043    """ 
    4144 
    …… 
    7174        a   - dimensionless 
    7275        b   - in units of seconds 
     76 
     77    See also: 
     78        inh_gamma_generator 
     79 
    7380    """ 
    7481 
    …… 
    133140        inh_poisson_generator - inhomogeneous Poisson process (time varying rate) 
    134141        inh_gamma_generator - inhomogeneous Gamma process (time varying a,b) 
     142        inh_adaptingmarkov_generator - inhomogeneous adapting markov process (time varying) 
     143        inh_2Dadaptingmarkov_generator - inhomogeneous adapting and  
     144                                         refractory markov process (time varying) 
    135145 
    136146        Continuous time processes: 
    …… 
    178188          
    179189        See also: 
    180             inh_poisson_generator, inh_gamma_generator 
     190            inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator 
    181191        """ 
    182192 
    …… 
    261271 
    262272        See also: 
    263             poisson_generator, inh_gamma_generator 
     273            poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator 
    264274        """ 
    265275 
    …… 
    323333 
    324334        Examples: 
    325          
     335            See source:trunk/examples/stgen/inh_gamma_psth.py 
     336 
    326337        See also: 
    327338            inh_poisson_generator, gamma_hazard 
    …… 
    434445 
    435446        Examples: 
     447            See source:trunk/examples/stgen/inh_2Dmarkov_psth.py 
     448 
    436449         
    437450        See also: 
    …… 
    558571 
    559572        Examples: 
     573            See source:trunk/examples/stgen/inh_2Dmarkov_psth.py 
    560574         
    561575        See also: 
    …… 
    642656        """  
    643657        Generates an Orstein Ulbeck process using the forward euler method. The function returns 
    644         an AnalogSignal object 
     658        an AnalogSignal object. 
    645659         
    646660        Inputs: 
    …… 
    651665            t_start - start time in milliseconds 
    652666            t_stop  - end time in milliseconds 
    653             array   - if True, a numpy array of the OU signal is returned, 
    654                       rather than an AnalogSignal object. 
     667            array   - if True, the functions returns the tuple (y,t)  
     668                      where y and t are the OU signal and the time bins, respectively, 
     669                      and are both numpy arrays. 
    655670         
    656671        Examples: 
    657672            >> stgen.OU_generator(0.1, 2, 3, 0, 0, 10000) 
     673 
     674        See also: 
     675            OU_generator_weave1 
    658676        """ 
    659677 
    …… 
    693711        """  
    694712        Generates an Orstein Ulbeck process using the forward euler method. The function returns 
    695         an AnalogSignal object 
     713        an AnalogSignal object. 
    696714         
    697715        Inputs: 
    …… 
    702720            t_start - start time in milliseconds 
    703721            t_stop  - end time in milliseconds 
    704             array   - if True, a numpy array of the OU signal is returned, 
    705                       rather than an AnalogSignal object. 
     722            array   - if True, the functions returns the tuple (y,t)  
     723                      where y and t are the OU signal and the time bins, respectively, 
     724                      and are both numpy arrays. 
    706725         
    707726        Examples: 
    708727            >> stgen.OU_generator(0.1, 2, 3, 0, 0, 10000) 
     728 
     729        See also: 
     730            OU_generator_weave1 
    709731        """ 
    710732 
    …… 
    742764 
    743765    def OU_generator_weave1(self, dt,tau,sigma,y0,t_start=0.0,t_stop=1000.0,time_it=False): 
    744         """ generates an OU process using forward euler 
    745         dt = resolution 
    746         tau = correlation time 
    747         sigma = std dev of process 
    748         y0 = mean/initial value 
    749         tsim = how long to simulate 
     766        """  
     767        Generates an Orstein Ulbeck process using the forward euler method. The function returns 
     768        an AnalogSignal object. 
     769 
     770        OU_generator_weave1, as opposed to OU_generator, uses scipy.weave 
     771        and is thus much faster. 
     772         
     773        Inputs: 
     774            dt      - the time resolution in milliseconds of th signal 
     775            tau     - the correlation time in milliseconds 
     776            sigma   - std dev of the process 
     777            y0      - initial value of the process, at t_start 
     778            t_start - start time in milliseconds 
     779            t_stop  - end time in milliseconds 
     780            array   - if True, the functions returns the tuple (y,t)  
     781                      where y and t are the OU signal and the time bins, respectively, 
     782                      and are both numpy arrays. 
     783         
     784        Examples: 
     785            >> stgen.OU_generator_weave1(0.1, 2, 3, 0, 0, 10000) 
     786 
     787        See also: 
     788            OU_generator 
    750789        """ 
    751790        import scipy.weave 
    …… 
    805844# TODO fix shotnoise stuff below  ... and write tests 
    806845 
    807 def shotnoise_fromspikes(self,rate,tau,q,num,t): 
    808     """  
    809     Generate shotnoise 
    810  
    811     Inputs: 
    812         rate - the rate (in Hz) of the shotnoise 
    813         tau  - the exponential decay of the synapse 
    814         g    - the quantal increase for a spike 
    815  
    816     quantal-increase-"q"-exponential-decay-"tau" synapse 
    817     and a poisson spike train of "rate" and "num" """ 
    818  
    819     g = numpy.zeros(numpy.shape(t),float) 
    820     for i in xrange(num): 
    821  
    822         spikes = poisson_generator(rate,t[-1]) 
    823         dg=exp_conv(spikes,t,q,tau) 
    824         tmp = add(g,dg,g) 
    825  
    826     return g 
    827  
    828  
    829  
    830  
    831846# Operations on spike trains 
    832847 
    833848 
    834 def shotnoise_fromspikes(spike_train,q,tau,dt,array=False, eps = 1.0e-8): 
     849def shotnoise_fromspikes(spike_train,q,tau,dt=0.1,t_start=None, t_stop=None,array=False, eps = 1.0e-8): 
    835850    """  
    836851    Convolves the provided spike train with shot decaying exponentials 
    …… 
    842857      q - the shot jump for each spike 
    843858      tau - the shot decay time constant in milliseconds 
     859      t_start - start time of the resulting AnalogSignal 
     860                If unspecified, t_start of spike_train is used 
     861      t_stop  - stop time of the resulting AnalogSignal 
     862                If unspecified, t_stop of spike_train is used 
    844863      dt - the resolution of the resulting shotnoise in milliseconds 
    845864      array - if True, returns (shotnoise,t) as numpy arrays, otherwise an AnalogSignal. 
    …… 
    847866      the shot kernal the tail is cut.  The default is usually fine. 
    848867 
     868   Note: 
     869      Spikes in spike_train before t_start are taken into account in the convolution. 
     870 
    849871   Examples: 
    850872      >> stg = stgen.StGen() 
    851873      >> st = stg.poisson_generator(10.0,0.0,1000.0) 
    852       >> g_e = shotnoise_fromspikes(st,2.0,10.0) 
     874      >> g_e = shotnoise_fromspikes(st,2.0,10.0,dt=0.1) 
    853875 
    854876 
    855877   See also: 
    856       poisson_generator, inh_gamma_generator, OU_generator ... 
     878      poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator, OU_generator ... 
    857879   """ 
    858880 
    859881    st = spike_train 
    860882 
    861     t = numpy.arange(st.t_start,st.t_stop,dt) 
     883    assert t_stop>t_start 
    862884 
    863885    # time of vanishing significance 
    864886    vs_t = -tau*numpy.log(eps/q) 
    865887 
     888 
     889    if t_stop == None: 
     890        t_stop = st.t_stop 
     891 
     892    # need to be clever with start time 
     893    # because we want to take spikes into 
     894    # account which occured in spikes_times 
     895    # before t_start 
     896    if t_start == None: 
     897        t_start = st.t_start 
     898        window_start = st.t_start 
     899    else: 
     900        window_start = t_start 
     901        if t_start>st.t_start: 
     902            t_start = st.t_start 
     903             
     904 
     905    t = numpy.arange(t_start,t_stop,dt) 
     906 
     907 
    866908    kern = q*numpy.exp(-numpy.arange(0.0,vs_t,dt)/tau) 
    867909 
    868     idx = numpy.clip(numpy.searchsorted(t,poisson_train,'right')-1,0,len(t)-1) 
    869  
    870     a = numpy.zeros(shape(t),float) 
     910    idx = numpy.clip(numpy.searchsorted(t,st.spike_times,'right')-1,0,len(t)-1) 
     911 
     912    a = numpy.zeros(numpy.shape(t),float) 
    871913 
    872914    a[idx] = 1.0 
    873915 
    874     y = convolve(a,kern)[0:len(t)] 
    875  
    876     result = AnalogSignal(y,dt,t_start=0,t_stop=st.t_stop-st.t_start) 
    877     result.time_offset(st.t_start) 
     916    y = numpy.convolve(a,kern)[0:len(t)] 
     917 
     918    result = AnalogSignal(y,dt,t_start=0.0,t_stop=t_stop-t_start) 
     919    result.time_offset(t_start) 
     920    if window_start>t_start: 
     921        result = result.time_slice(window_start,t_stop) 
    878922    return result 
    879923 
    …… 
    882926 
    883927 
    884 def gen_g_add(spikes,tau,q,t,eps = 1.0e-8): 
     928def _gen_g_add(spikes,tau,q,t,eps = 1.0e-8): 
    885929 
    886930    #spikes = poisson_generator(rate,t[-1]) 

  • trunk/test/test_stgen.py

    r323 r337  
    294294 
    295295 
     296    def testShotNoiseFromSpikes(self): 
     297 
     298 
     299        stg = stgen.StGen() 
     300 
     301        from numpy import array 
     302 
     303        st = stg.poisson_generator(10.0,0.0,1000.0) 
     304         
     305        ge = shotnoise_fromspikes(st,2.0,10.0,dt=0.1) 
     306 
     307        assert ge.t_start==0.0 
     308        assert ge.t_stop==1000.0 
     309 
     310        st = stg.poisson_generator(10.0,0.0,1000.0) 
     311        ge = shotnoise_fromspikes(st,2.0,10.0,dt=0.1,t_start=500.0,t_stop=1500.0) 
     312         
     313        assert ge.t_start==500.0 
     314        assert ge.t_stop==1500.0 
     315 
    296316 
    297317