Changeset 77

Show
Ignore:
Timestamp:
01/04/08 20:08:42 (1 year ago)
Author:
LaurentPerrinet
Message:

Major changes to NeuroTools:
- simplifying the structure of the package
- using the SpikeList object for analysis
- examples for the single_neuron and the simple retina applied with nest1 (to be validated with other simulators)

the previous working copy is available as a branch (branches/working_copy) on the SVN

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/facets/fkbtools.py

    r75 r77  
    55 
    66try: 
    7     import NeuroTools.hdf5.FileExtension as file_extension 
    8     import NeuroTools.hdf5.Movie as movie 
     7    import NeuroTools.facets.hdf5.FileExtension as file_extension 
     8    import NeuroTools.facets.hdf5.Movie as movie 
    99    use_hdf5 = True 
    1010except ImportError: 
     11    print "Failed to load HDF5 tools" 
    1112    use_hdf5 = False 
    1213import srblib 
    …… 
    8283    parameters = _get_parameters(png_dir) 
    8384    print "Creating movie file with %d frames, saving to %s" % (len(frames),output_url) 
    84     h5file = file_extension.openHDF5File(output_url, "w", title=title)     
     85    h5file = file_extension.openHDF5File(output_url, "w", title=title) 
    8586    h5file.createCMovie(h5file.root, "Movie", frames, coding=8, complevel=9, 
    8687                       frame_duration=parameters['frame_duration']) 
    …… 
    165166    Lines beginning with '#' are taken to be comments. 
    166167    Column headings can be specified as a comment on the first line of the file. 
    167      
     168 
    168169    Example: 
    169170    # time voltage 
    …… 
    176177    contents = f.readlines() 
    177178    f.close() 
    178      
     179 
    179180    # extract header information 
    180181    firstline = contents[0] 
    …… 
    198199    if len(headers) != nfields: 
    199200        headers = ['col%d' % i for i in range(1,nfields+1)] 
    200      
     201 
    201202    # create row format 
    202     class_str = "class DataRow(tables.IsDescription):\n"   
     203    class_str = "class DataRow(tables.IsDescription):\n" 
    203204    for pos,colname in enumerate(headers): 
    204205        class_str += "  %s = tables.Float32Col(pos=%d)\n" % (colname,pos) 
    205206    exec(class_str) 
    206      
     207 
    207208    # create HDF5 file 
    208209    protocol, path = output_url.split('://') 
    …… 
    211212    assert protocol == 'file', "For now, we only support writing to local files." 
    212213    h5file = tables.openFile('%s.h5' % path, mode='w', title=filename) 
    213      
     214 
    214215    # create and fill table 
    215216    table = h5file.createTable(h5file.root, filename, DataRow, "Converted from %s" % input_url) 
    …… 
    224225        row.append() 
    225226    table.flush() 
    226      
     227 
    227228    h5file.close() 
    228229 
    …… 
    239240    assert len(h5file.listNodes(h5file.root,classname="Table")) == 1, "File contains more than one table." 
    240241    table = h5file.listNodes(h5file.root,classname="Table")[0] 
    241      
     242 
    242243    # Open and write to text file 
    243244    protocol, path = output_url.split('://') 
    …… 
    253254            values[i] = str(row[colname]) 
    254255        txtfile.write(spacer.join(values)+'\n') 
    255      
     256 
    256257    txtfile.close() 
    257258    h5file.close() 
    258      
     259 
    259260def getData(url): 
    260261    """ 
    …… 
    278279        raise Exception("Not a valid file.") 
    279280    return arr 
    280      

  • trunk/spikes.py

    r73 r77  
    11import numpy 
    2 from numpy import array 
    3 from analysis import nstats 
     2#from numpy import array 
     3#from analysis import nstats 
    44 
    55class SpikeTrain(object): 
    …… 
    2424    # TODO in the definition, should a spike train be ordered? 
    2525 
     26    # TODO : should we handle different population shapes differently? 
     27 
    2628    def __init__(self, spike_times, dt=None, t_start=None, t_stop=None): 
    2729        """ 
    28         `spike_times` is a list/numpy array of spike times 
     30        `spike_times` is a list/numpy array of spike times (in seconds) 
     31 
     32        TODO: We proposed innitially that : 
    2933        If `dt`is specified, the time values should be ints, 
    3034        and will be multiplied by `dt`, otherwise time values should be floats. 
     35 
     36        this is confusing. All times should be in seconds. period. 
     37 
    3138        If `t_start` and `t_stop` are not specified, they are inferred from the data. 
    3239 
    …… 
    3643        self.spike_times = numpy.array([numpy.float(time) for time in spike_times]) 
    3744 
    38         if not(dt==None): 
    39             if dt<=0: 
    40                 raise("dt must be greater than zero") 
    41             self.spike_times *= dt 
     45##        if not(dt==None): 
     46##            if dt<=0: 
     47##                raise("dt must be greater than zero") 
     48##            self.spike_times *= dt 
     49##        else: 
     50##            # knowing dt may be useful (for the spike matrix) 
     51##            dt = 1e4 # 0.1 ms per default 
     52##            # TODO : compute as the best descriptor for a list of floats 
    4253 
    4354        self.dt = dt 
    4455        self.t_start = t_start 
    4556        self.t_stop = t_stop 
    46          
     57 
    4758        if len(spike_times)>0: # spike list may be empty 
    48             # TODO raise an error if some data is outside [t_start, t_stop] ? 
    49             #TODO return an exception if self.t_stop==self.t_start 
    50              
    51             if self.t_start == None: 
     59            if self.t_start is None: 
    5260                self.t_start = numpy.min(self.spike_times) 
    5361 
    54             if self.t_stop == None: 
     62            if self.t_stop is None: 
    5563                self.t_stop = numpy.max(self.spike_times) 
     64 
     65        # TODO raise an error if some data is outside [t_start, t_stop] ? 
     66        #TODO return an exception if self.t_stop < self.t_start (when not empty) 
     67            if self.t_start >= self.t_stop : 
     68                raise("Incompatible time interval for the creation of the SpikeTrain") 
     69 
    5670 
    5771    def __str__(self): 
    …… 
    90104        """ 
    91105        cv_isi is the ratio between the standard deviation and the mean of the ISI 
     106 
     107          The irregularity of individual spike trains is measured by the squared 
     108        coefficient of variation of the corresponding inter-spike interval (ISI) 
     109        distribution normalized by the square of its mean. 
     110          In point processes, low values reflect more regular spiking, a 
     111        clock-like pattern yields CV2= 0. On the other hand, CV2 = 1 indicates 
     112        Poisson-type behavior. As a measure for irregularity in the network one 
     113        can use the average irregularity across all neurons. 
     114 
     115        TODO: is it useful to get the std of CV? 
    92116        """ 
    93117        isi = self.isi() 
    …… 
    95119        return numpy.std(isi)/numpy.mean(isi) 
    96120 
    97     def spike_histogram(self, binwidth, normalized=False): 
     121    def time_axis(self, binwidth): 
     122        return numpy.arange(self.t_start, self.t_stop, binwidth) 
     123 
     124    def time_histogram(self, binwidth , normalized=False): 
    98125        """ 
    99126        Bin the spikes with the specified bin width. The first and last bins 
    …… 
    102129        firing rates in spikes/second, otherwise. 
    103130        """ 
    104         bins = numpy.arange(self.t_start, self.t_stop, binwidth) 
    105         hist = nstats.histc(self.spike_times, bins) 
     131        bins = self.time_axis(binwidth) 
     132        (hist, lower_edges ) = numpy.histogram(self.spike_times, bins) 
     133 
    106134        if normalized: 
    107             hist *= 1000.0/binwidth 
     135            hist *= 1000.0/binwidth # TODO why 1000.0 ?? we are using seconds everywhere 
     136 
    108137        return hist 
    109138 
    …… 
    119148        so that each bin is associated with a particular value of the variable 
    120149        in `var_array`. 
    121          
     150 
    122151        The return value is a dictionary whose keys are the distinct values in 
    123152        `val_array`. The values in the dict depend on the arguments `method` and 
    124153        `normalized`. 
    125          
     154 
    126155        If `normalized` is False, the responses (bin values) are spike counts, 
    127156        if True, they are firing rates. 
    …… 
    130159        If `method` == "sum", each value is the summed response... 
    131160        If `method` == "mean", ... you get the idea. 
    132          
     161 
    133162        (If someone can rewrite this more clearly, please do so!) 
    134163        """ 
    135164        binwidth = (self.t_stop - self.t_start)/len(var_array) 
    136         spike_histogram = self.spike_histogram(binwidth, normalized=normalized) 
    137         assert len(spike_histogram) == len(var_array) 
     165        time_histogram = self.time_histogram(binwidth, normalized=normalized) 
     166        assert len(time_histogram) == len(var_array) 
    138167        #def _orientation_tuning_curve(time_bins, orientations, responses, method='sum'): 
    139168        tuning_curve = {} 
    140169        counts = {} 
    141         for k, x in zip(var_array, spike_histogram): 
     170        for k, x in zip(var_array, time_histogram): 
    142171            if not tuning_curve.has_key(k): 
    143172                tuning_curve[k] = 0 
    …… 
    186215        """ 
    187216        `spikes` is a list/tuple of (id,t) tuples (id in range(0,N)) 
    188         `N` is the total number of units TODO should be possibly a tuple in accordance with the dimension of a pyNN.Population 
     217        `N` is the total number of units TODO should be possibly a tuple in accordance 
     218        with the dimension of a pyNN.Population 
    189219        If `dt`is specified, the time values should be ints, 
    190220        and will be multiplied by `dt`, otherwise time values should be floats. 
    …… 
    194224 
    195225        """ 
    196  
     226        self.N = N 
     227        # taking care of some formatting 
     228        # if not(dt==None): 
     229        #    raise("not implemented, don't assign dt") 
     230 
     231        self.t_start = t_start 
     232        self.t_stop = t_stop 
     233        self.dt = dt 
    197234        # transform spikes in a spike array 
    198235        spike_array = list([ [] for i in range(N)]) 
    199236        times =[] 
    200  
    201237        for spike in spikes: 
    202238            times.append(spike[1]) 
    203239            spike_array[spike[0]].append(spike[1]) 
    204240 
    205         # taking care of some formatting 
    206         if not(dt==None): 
    207             raise("not implemented, don't assign dt") 
    208         # TODO raise an error if some data is outside [t_start, t_stop] ? 
    209         self.t_start = t_start 
    210         self.t_stop = t_stop 
    211  
    212         self.N = N 
    213         self.dt = dt 
    214  
    215  
    216         # writing spike train list 
     241        # writing as a list of SpikeTrains 
    217242        self.spiketrain_list = [] 
    218         for i_neuron in range(N): 
    219             #self.spiketrain_list.append( SpikeTrain(spike_array[i_neuron], dt=dt, t_start=t_start, t_stop=t_stop) ) 
    220             self.spiketrain_list.append( SpikeTrain(spike_array[i_neuron], dt=dt )) 
     243        for i_neuron in range(N): #TODO: handle missing fields 
     244            self.spiketrain_list.append( SpikeTrain(spike_array[i_neuron], dt=dt, t_start=t_start, t_stop=t_stop) ) 
     245            #self.spiketrain_list.append( SpikeTrain(spike_array[i_neuron], dt=dt )) 
    221246        self.__calc_startstop() 
    222247 
    223248    def __calc_startstop(self): 
     249        """ 
     250        t_start and t_stop are shared for all neurons, so we take min and max values respectively. 
     251 
     252        """ 
    224253        if len(self) > 0: 
    225254            start_times = numpy.array([spiketrain.t_start for spiketrain in self.spiketrain_list]) 
    226255            stop_times = numpy.array([spiketrain.t_stop for spiketrain in self.spiketrain_list]) 
    227256            self.t_stop = max( self.t_stop, stop_times.max() ) 
    228             if self.t_start == None: 
     257            if self.t_start is None: 
    229258                self.t_start = start_times.min() 
    230259            else: 
    …… 
    248277            self.N = id+1 
    249278        self.__calc_startstop() 
    250          
     279 
     280    def time_axis(self, binwidth): 
     281        return self[0].time_axis(binwidth) 
     282 
     283 
    251284 
    252285    # methods for different output formats 
    …… 
    289322 
    290323    def as_spikematrix(self): 
    291         pass 
     324        """ Returns a list giving for each neuron how many spikes per 
     325        time bin, usually zero everywhere, one if a spike. 
     326        Useful for plots etc. 
     327 
     328        """ 
     329        return self.firing_rate(self.dt) 
    292330 
    293331 
    …… 
    316354        pass 
    317355 
     356 
     357    def firing_rate(self, binwidth): 
     358        """ 
     359        Calculate firing rate traces (in Hz) from arrays of spike times. 
     360 
     361        Spike times and binwidth should are in seconds. 
     362        Returns a tuple (list_of_firing_rate_traces,bins) 
     363 
     364        >>> pylab.plot(output[0].time_axis(0.02),sum(output.firing_rate(0.02))) 
     365        """ 
     366 
     367        #time_size = self.time_axis(binwidth).size 
     368        time_size = int((self.t_stop-self.t_start)/binwidth) 
     369        firing_rate = numpy.zeros((self.N,time_size)) 
     370 
     371        for i_neuron in range(self.N): 
     372            firing_rate[i_neuron,:] = self.spiketrain_list[i_neuron].time_histogram(binwidth) 
     373 
     374        return firing_rate 
     375 
    318376    def fano_factor(self, binwidth): 
    319377        pass 
    …… 
    322380        pass 
    323381 
     382 
     383    def pw_corr_pearson(self,edge,bins,number_of_neurons): 
     384        """ 
     385        TODO: document/test this function 
     386 
     387        """ 
     388        #bins = edge[1]-edge[0] 
     389        cor = numpy.array([]) 
     390        [neuron_ids_array, spike_times_array] = self.as_list_id_list_time() 
     391 
     392        # Pairwise correlation 
     393        neuron_ids_unique = numpy.unique(neuron_ids_array) 
     394        for count in range(number_of_neurons): 
     395            # draw two neuron radomly out of the neuron_id_unique pool 
     396            neuron_1_tmp = neuron_ids_unique[numpy.floor(numpy.random.uniform()*len(neuron_ids_unique))] 
     397            neuron_2_tmp = neuron_ids_unique[numpy.floor(numpy.random.uniform()*len(neuron_ids_unique))] 
     398            # get the spike_times of the neurons 
     399            spike_times_1_tmp = spike_times_array[neuron_ids_array == neuron_1_tmp] 
     400            spike_times_2_tmp = spike_times_array[neuron_ids_array == neuron_2_tmp] 
     401            # hist of the spiketrains 
     402            n1_hist = numpy.histogram(spike_times_1_tmp,bins=bins,range=edge) 
     403            n2_hist = numpy.histogram(spike_times_2_tmp,bins=bins,range=edge) 
     404            # correlation 
     405            # TODO: normalize the cor, look in 1.6 in Kumar et al. 
     406            cov = numpy.corrcoef(n1_hist[0],n2_hist[0])[1][0] 
     407            #print n1_hist[0] 
     408            #print n2_hist[0] 
     409            cor = cov/numpy.sqrt(n1_hist[0].var()*n2_hist[0].var()) 
     410            #cor = numpy.append(cor,numpy.corrcoef(n1_hist[0],n2_hist[0])[1][0]) 
     411 
     412        cor_coef_mean = cor.mean() 
     413        cor_coef_std = cor.std() 
     414        return cor_coef_mean, cor_coef_std 
    324415 
    325416""" 
    …… 
    327418 
    328419""" 
    329 def tmpfile2spikelist(filename, N, t_start=None, t_stop=None): 
     420 
     421def population2spikelist(output,N, dt , t_start, t_stop): 
     422    """ 
     423    WARNING: all times in ms in PyNN, in seconds here. 
     424 
     425    """ 
     426    import os, tempfile 
     427 
     428    tmpdir = tempfile.mkdtemp() 
     429    output_filename=os.path.join(tmpdir,'output.gdf') 
     430    output.printSpikes(output_filename)# 
     431    output_DATA = tmpfile2spikelist(output_filename, N= N, dt = dt, t_start=t_start, t_stop=t_stop) 
     432    os.remove(output_filename) 
     433    os.rmdir(tmpdir) 
     434    return output_DATA 
     435 
     436def tmpfile2spikelist(filename, N, dt = None, t_start=None, t_stop=None): 
    330437    """ 
    331438    Returns a SpikeList object from the tmp file saved by PyNN-NEST. 
    …… 
    339446    * a line commented, 
    340447    * then one line per event:  absolute time in ms, GID 
    341     (see function printSpikes in nest.py) 
    342  
    343     """ 
    344     #import os, numpy 
     448    (see function printSpikes in nest1.py) 
     449 
     450    """ 
     451 
    345452    f = open(filename,'r') 
    346453    lines = f.readlines() 
    …… 
    351458    for i_line in range(1,len(lines)): 
    352459        [spike_time , neuron_id] = lines[i_line].split("\t", 1) 
    353         spikes.append( (int(neuron_id),float(spike_time)) ) 
    354  
    355     sl = SpikeList(N, spikes=spikes, dt=None, t_start=t_start, t_stop=t_stop) 
    356     print sl 
     460        #TODO: all times in ms in PyNN, in seconds here.Unify? 
     461        spikes.append( (int(neuron_id),float(spike_time)/1000.) ) 
     462 
     463    sl = SpikeList(N, spikes=spikes, dt=dt, t_start=t_start, t_stop=t_stop) 
     464    #print sl 
    357465    return sl 
    358466 

  • trunk/stgen.py

    r75 r77  
    11# -*- coding: utf8 -*- 
    22# TODO : needs refactoring 
    3 # needs python gsl wrapper 
     3# TODO: needs python gsl wrapper / convert to numpy 
     4# TODO make it generate spiketrains? 
    45 
    56import pygsl 
    6 from analysis import nstats 
     7#from analysis import nstats 
    78from numpy import array 
    89import numpy 
    …… 
    4243        If neither is specified, a gsl rng is created. 
    4344        If both are specified, the gslrng is used.""" 
    44          
     45 
    4546        if not gslrng: 
    4647            if not numpyrng: 
    …… 
    8788        rmax = numpy.max(rate) 
    8889        ps = self.poisson_generator(rmax, tsim) 
    89          
     90 
    9091        if len(ps) > 0: 
    9192            # gen uniform rand on 0,1 for each spike 
    9293            rn = array(self.rng.uniform(len(ps))) 
    93      
     94 
    9495            # instantaneous rate for each spike 
    95      
     96 
    9697            #idx = len(tbins)-1-numpy.searchsorted(-tbins[::-1],-ps) 
    9798            ## for all times before first tbin, use first tbin 
    9899            #idx = numpy.where(idx>=0,idx,0) 
    99      
     100 
    100101            # TODO : this simplifies the above code but "there must a reason why this was more complex..." 
    101102            idx=numpy.searchsorted(tbins,ps)-1 
    102103            spike_rate = numpy.choose(idx,rate) 
    103      
     104 
    104105            # thin and return spikes 
    105106            keep = numpy.nonzero(rn<spike_rate/rmax) 
    …… 
    108109        else: 
    109110            spike_train = ps 
    110              
     111 
    111112        return spike_train 
    112113