root/trunk/src/signals.py

"""
NeuroTools.signals
==================

A collection of functions to create, manipulate and play with spikes and analog
signals. 

Classes
-------

SpikeTrain       - object representing a spike train, for one cell. Useful for plots, 
                   calculations such as ISI, CV, mean rate(), ...
SpikeList        - object representing the activity of a population of neurons. Functions as a
                   dictionary of SpikeTrain objects, with methods to compute firing rate,
                   ISI, CV, cross-correlations, and so on. 
AnalogSignal     - object representing an analog signal, with its data. Can be used to do 
                   threshold detection, event triggered averages, ...
AnalogSignalList - list of AnalogSignal objects, again with methods such as mean, std, plot, 
                   and so on
VmList           - AnalogSignalList object used for Vm traces
ConductanceList  - AnalogSignalList object used for conductance traces
CurrentList      - AnalogSignalList object used for current traces

Functions
---------

load_spikelist       - load a SpikeList object from a file. Expects a particular format.
                       Can also load data in a different format, but then you have
                       to write your own File object that will know how to read the data (see io.py)
load_vmlist          - function to load a VmList object (inherits from AnalogSignalList) from a file.
                       Same comments on format as previously.
load_currentlist     - function to load a CurrentList object (inherits from AnalogSignalList) from a file.
                       Same comments on format as previously.
load_conductancelist - function to load a ConductanceList object (inherits from AnalogSignalList) from a file.
                       Same comments on format as previously. load_conductancelist returns two 
                       ConductanceLists, one for the excitatory conductance and one for the inhibitory conductance
load                 - a generic loader for all the previous load methods.
"""

import os, re
from NeuroTools import check_dependency, check_numpy_version
from NeuroTools import analysis
from NeuroTools.io import *
from NeuroTools.plotting import get_display, set_axis_limits, set_labels, SimpleMultiplot

if check_dependency('psyco'):
    import psyco
    psyco.full()

import numpy
newnum = check_numpy_version()

HAVE_PYLAB = check_dependency('pylab')
HAVE_MATPLOTLIB = check_dependency('matplotlib')
if HAVE_PYLAB:
    import pylab
else:
    PYLAB_ERROR = "The pylab package was not detected"
if not HAVE_MATPLOTLIB:
    MATPLOTLIB_ERROR = "The matplotlib package was not detected"

class SpikeTrain(object):
    """
    SpikeTrain(spikes_times, t_start=None, t_stop=None)
    This class defines a spike train as a list of times events.

    Event times are given in a list (sparse representation) in milliseconds.

    Inputs:
        spike_times - a list/numpy array of spike times (in milliseconds)
        t_start     - beginning of the SpikeTrain (if not, this is infered)
        t_stop      - end of the SpikeTrain (if not, this is infered)

    Examples:
        >> s1 = SpikeTrain([0.0, 0.1, 0.2, 0.5])
        >> s1.isi()
            array([ 0.1,  0.1,  0.3])
        >> s1.mean_rate()
            8.0
        >> s1.cv_isi()
            0.565685424949
    """

    #######################################################################
    ## Constructor and key methods to manipulate the SpikeTrain objects  ##
    #######################################################################
    def __init__(self, spike_times, t_start=None, t_stop=None):
        """
        Constructor of the SpikeTrain object
        
        See also
            SpikeTrain
        """

        self.t_start     = t_start
        self.t_stop      = t_stop
        self.spike_times = numpy.array(spike_times, float)

        # If t_start is not None, we resize the spike_train keeping only
        # the spikes with t >= t_start
        if self.t_start is not None:
            self.spike_times = numpy.extract((self.spike_times >= self.t_start), self.spike_times)

        # If t_stop is not None, we resize the spike_train keeping only
        # the spikes with t <= t_stop
        if self.t_stop is not None:
            self.spike_times = numpy.extract((self.spike_times <= self.t_stop), self.spike_times)

        # We sort the spike_times. May be slower, but is necessary by the way for quite a 
        # lot of methods...
        self.spike_times = numpy.sort(self.spike_times, kind="quicksort")
        # Here we deal with the t_start and t_stop values if the SpikeTrain
        # is empty, with only one element or several elements, if we
        # need to guess t_start and t_stop
        # no element : t_start = 0, t_stop = 0.1 
        # 1 element  : t_start = time, t_stop = time + 0.1
        # several    : t_start = min(time), t_stop = max(time)
        
        size = len(self.spike_times)
        if size == 0:
            if self.t_start is None: 
                self.t_start = 0
            if self.t_stop is None:
                self.t_stop  = 0.1
        elif size == 1: # spike list may be empty
            if self.t_start is None:
                self.t_start = self.spike_times[0]
            if self.t_stop is None:
                self.t_stop = self.spike_times[0] + 0.1
        elif size > 1:
            if self.t_start is None:
                self.t_start = numpy.min(self.spike_times)
            if numpy.any(self.spike_times < self.t_start):
                raise ValueError("Spike times must not be less than t_start")
            if self.t_stop is None:
                self.t_stop = numpy.max(self.spike_times)
            if numpy.any(self.spike_times > self.t_stop):
                raise ValueError("Spike times must not be greater than t_stop")

        if self.t_start >= self.t_stop :
            raise Exception("Incompatible time interval : t_start = %s, t_stop = %s" % (self.t_start, self.t_stop))
        if self.t_start < 0:
            raise ValueError("t_start must not be negative")
        if numpy.any(self.spike_times < 0):
            raise ValueError("Spike times must not be negative")

    def __str__(self):
        return str(self.spike_times)

    def __len__(self):
        return len(self.spike_times)
    
    def __getslice__(self, i, j):
        """
        Return a sublist of the spike_times vector of the SpikeTrain
        """
        return self.spike_times[i:j]
    
    def time_parameters(self):
        """
        Return the time parameters of the SpikeTrain (t_start, t_stop)
        """
        return (self.t_start, self.t_stop)
    
    def is_equal(self, spktrain):
        """
        Return True if the SpikeTrain object is equal to one other SpikeTrain, i.e
        if they have same time parameters and same spikes_times
        
        Inputs:
            spktrain - A SpikeTrain object
        
        See also:
            time_parameters()
        """
        test = (self.time_parameters() == spktrain.time_parameters())
        return numpy.all(self.spike_times == spktrain.spike_times) and test
    
    def copy(self):
        """
        Return a copy of the SpikeTrain object
        """
        return SpikeTrain(self.spike_times, self.t_start, self.t_stop)


    def duration(self):
        """
        Return the duration of the SpikeTrain
        """
        return self.t_stop - self.t_start
    
    
    def merge(self, spiketrain):
        """
        Add the spike times from a spiketrain to the current SpikeTrain
        
        Inputs:
            spiketrain - The SpikeTrain that should be added
        
        Examples:
            >> a = SpikeTrain(range(0,100,10),0.1,0,100)
            >> b = SpikeTrain(range(400,500,10),0.1,400,500)
            >> a.merge(b)
            >> a.spike_times
                [   0.,   10.,   20.,   30.,   40.,   50.,   60.,   70.,   80.,
                90.,  400.,  410.,  420.,  430.,  440.,  450.,  460.,  470.,
                480.,  490.]
            >> a.t_stop
                500
        """
        self.spike_times = numpy.insert(self.spike_times, self.spike_times.searchsorted(spiketrain.spike_times), \
                                        spiketrain.spike_times)
        self.t_start     = min(self.t_start, spiketrain.t_start)
        self.t_stop      = max(self.t_stop, spiketrain.t_stop)

    def format(self, relative=False, quantized=False):
        """
        Return an array with a new representation of the spike times
        
        Inputs:
            relative  - if True, spike times are expressed in a relative
                       time compared to the previsous one
            quantized - a value to divide spike times with before rounding
            
        Examples:
            >> st.spikes_times=[0, 2.1, 3.1, 4.4]
            >> st.format(relative=True)
                [0, 2.1, 1, 1.3]
            >> st.format(quantized=2)
                [0, 1, 2, 2]
        """
        spike_times = self.spike_times.copy()

        if relative and len(spike_times) > 0:
            spike_times[1:] = spike_times[1:] - spike_times[:-1]

        if quantized:
            assert quantized > 0, "quantized must either be False or a positive number"
            #spike_times =  numpy.array([time/self.quantized for time in spike_times],int)
            spike_times = (spike_times/quantized).round().astype('int')

        return spike_times

    def jitter(self,jitter):
        """
        Returns a new SpikeTrain with spiketimes jittered by a normal distribution.

        Inputs:
              jitter - sigma of the normal distribution

        Examples:
              >> st_jittered = st.jitter(2.0)
        """
        
        return SpikeTrain(self.spike_times+jitter*(numpy.random.normal(loc=0.0,scale=1.0,size=self.spike_times.shape[0])),t_start=self.t_start,t_stop=self.t_stop)
        


    #######################################################################
    ## Analysis methods that can be applied to a SpikeTrain object       ##
    #######################################################################
    
    def isi(self):
        """
        Return an array with the inter-spike intervals of the SpikeTrain
        
        Examples:
            >> st.spikes_times=[0, 2.1, 3.1, 4.4]
            >> st.isi()
                [2.1, 1., 1.3]
        
        See also
            cv_isi
        """
        return numpy.diff(self.spike_times)

    def mean_rate(self, t_start=None, t_stop=None):
        """ 
        Returns the mean firing rate between t_start and t_stop, in Hz
        
        Inputs:
            t_start - in ms. If not defined, the one of the SpikeTrain object is used
            t_stop  - in ms. If not defined, the one of the SpikeTrain object is used
        
        Examples:
            >> spk.mean_rate()
                34.2
        """
        if (t_start == None) & (t_stop == None):
            t_start = self.t_start
            t_stop  = self.t_stop
            idx     = self.spike_times
        else:
            if t_start == None: 
                t_start = self.t_start
            else:
                t_start = max(self.t_start, t_start)
            if t_stop == None: 
                t_stop=self.t_stop
            else:
                t_stop = min(self.t_stop, t_stop)
            idx = numpy.where((self.spike_times >= t_start) & (self.spike_times <= t_stop))[0]
        return 1000.*len(idx)/(t_stop-t_start)

    def cv_isi(self):
        """
        Return the coefficient of variation of the isis.
        
        cv_isi is the ratio between the standard deviation and the mean of the ISI
          The irregularity of individual spike trains is measured by the squared
        coefficient of variation of the corresponding inter-spike interval (ISI)
        distribution normalized by the square of its mean.
          In point processes, low values reflect more regular spiking, a
        clock-like pattern yields CV2= 0. On the other hand, CV2 = 1 indicates
        Poisson-type behavior. As a measure for irregularity in the network one
        can use the average irregularity across all neurons.
        
        http://en.wikipedia.org/wiki/Coefficient_of_variation
        
        See also
            isi, cv_kl
            
        """
        isi = self.isi()
        if len(isi) > 0:
            return numpy.std(isi)/numpy.mean(isi)
        else:
            logging.debug("Warning, a CV can't be computed because there are not enough spikes")
            return numpy.nan

    def cv_kl(self, bins=100):
        """
        Provides a measure for the coefficient of variation to describe the
        regularity in spiking networks. It is based on the Kullback-Leibler
        divergence and decribes the difference between a given
        interspike-interval-distribution and an exponential one (representing
        poissonian spike trains) with equal mean.
        It yields 1 for poissonian spike trains and 0 for regular ones.
        
        Reference:
            http://incm.cnrs-mrs.fr/LaurentPerrinet/Publications/Voges08fens
        
        Inputs:
            bins - the number of bins used to gather the ISI
        
        Examples:
            >> spklist.cv_kl(100)
                0.98
        
        See also:
            cv_isi
            
        """
        isi = self.isi() / 1000.
        if len(isi) < 2:
            logging.debug("Warning, a CV can't be computed because there are not enough spikes")
            return numpy.nan
        else:
            if newnum:
                proba_isi, xaxis = numpy.histogram(isi, bins=bins, normed=True, new=True)
                xaxis = xaxis[:-1]
            else:
                proba_isi, xaxis = numpy.histogram(isi, bins=bins, normed=True)
            proba_isi /= numpy.sum(proba_isi)
            bin_size = xaxis[1]-xaxis[0]
            # differential entropy: http://en.wikipedia.org/wiki/Differential_entropy
            KL = - numpy.sum(proba_isi * numpy.log(proba_isi+1e-16)) + numpy.log(bin_size)
            KL -= -numpy.log(self.mean_rate()) + 1.
            CVkl=numpy.exp(-KL)
            return CVkl
    

    def fano_factor_isi(self):
        """ 
        Return the fano factor of this spike trains ISI.
        
        The Fano Factor is defined as the variance of the isi divided by the mean of the isi
        
        http://en.wikipedia.org/wiki/Fano_factor
        
        See also
            isi, cv_isi
        """
        isi = self.isi()
        if len(isi) > 0:
            fano = numpy.var(isi)/numpy.mean(isi)
            return fano
        else: 
            raise Exception("No spikes in the SpikeTrain !")

    def time_axis(self, time_bin=10):
        """
        Return a time axis between t_start and t_stop according to a time_bin
        
        Inputs:
            time_bin - the bin width
            
        Examples:
            >> st = SpikeTrain(range(100),0.1,0,100)
            >> st.time_axis(10)
                [ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90n 100]
        
        See also
            time_histogram
        """
        if newnum:
            axis = numpy.arange(self.t_start, self.t_stop+time_bin, time_bin)
        else:
            axis = numpy.arange(self.t_start, self.t_stop, time_bin)
        return axis

    def raster_plot(self, t_start=None, t_stop=None, display=True, kwargs={}):
        """
        Generate a raster plot with the SpikeTrain in a subwindow of interest,
        defined by t_start and t_stop. 
        
        Inputs:
            t_start - in ms. If not defined, the one of the SpikeTrain object is used
            t_stop  - in ms. If not defined, the one of the SpikeTrain object is used
            display - if True, a new figure is created. Could also be a subplot
            kwargs  - dictionary contening extra parameters that will be sent to the plot 
                      function
        
        Examples:
            >> z = subplot(221)
            >> st.raster_plot(display=z, kwargs={'color':'r'})
        
        See also
            SpikeList.raster_plot
        """
        if t_start is None: t_start = self.t_start
        if t_stop is None:  t_stop = self.t_stop
        spikes = numpy.extract((self.spike_times >= t_start) & (self.spike_times <= t_stop), self.spike_times)
        subplot = get_display(display)
        if not subplot or not HAVE_PYLAB:
            print PYLAB_ERROR
            return
        else:
            if len(spikes) > 0:
                xlabel = "Time (ms)"
                ylabel = "Neurons #"
                set_labels(subplot, xlabel, ylabel)
                subplot.scatter(spikes,numpy.ones(len(spikes)),**kwargs)
                pylab.draw()

    def time_offset(self, offset):
        """
        Add an offset to the SpikeTrain object. t_start and t_stop are
        shifted from offset, so does all the spike times.
         
        Inputs:
            offset - the time offset, in ms
        
        Examples:
            >> spktrain = SpikeTrain(arange(0,100,10))
            >> spktrain.time_offset(50)
            >> spklist.spike_times
                [  50.,   60.,   70.,   80.,   90.,  100.,  110.,  
                120.,  130.,  140.]
        """
        self.t_start += offset
        self.t_stop  += offset
        self.spike_times += offset

    def time_slice(self, t_start, t_stop):
        """ 
        Return a new SpikeTrain obtained by slicing between t_start and t_stop. The new 
        t_start and t_stop values of the returned SpikeTrain are the one given as arguments
        
        Inputs:
            t_start - begining of the new SpikeTrain, in ms.
            t_stop  - end of the new SpikeTrain, in ms.

        Examples:
            >> spk = spktrain.time_slice(0,100)
            >> spk.t_start
                0
            >> spk.t_stop
                100
        """
        spikes = numpy.extract((self.spike_times >= t_start) & (self.spike_times <= t_stop), self.spike_times)
        return SpikeTrain(spikes, t_start, t_stop)


    def time_histogram(self, time_bin=10, normalized=True):
        """
        Bin the spikes with the specified bin width. The first and last bins
        are calculated from `self.t_start` and `self.t_stop`.
        
        Inputs:
            time_bin   - the bin width for gathering spikes_times
            normalized - if True, the bin values are scaled to represent firing rates
                         in spikes/second, otherwise otherwise it's the number of spikes 
                         per bin.
        
        Examples:
            >> st=SpikeTrain(range(0,100,5),0.1,0,100)
            >> st.time_histogram(10)
                [200, 200, 200, 200, 200, 200, 200, 200, 200, 200]
            >> st.time_histogram(10, normalized=False)
                [2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
        
        See also
            time_axis
        """
        bins = self.time_axis(time_bin)
        if newnum:
            hist, edges = numpy.histogram(self.spike_times, bins, new=newnum)
        else:
            hist, edges = numpy.histogram(self.spike_times, bins)
        if normalized and isinstance(time_bin, int): # what about normalization if time_bin is a sequence?
            hist *= 1000.0/time_bin
        return hist

    def relative_times(self):
        """
        Rescale the spike times to make them relative to t_start.
        
        Note that the SpikeTrain object itself is modified, t_start 
        is substracted to spike_times, t_start and t_stop
        """
        if self.t_start != 0:
            self.spike_times -= self.t_start
            self.t_stop      -= self.t_start
            self.t_start      = 0.0

    def distance_victorpurpura(self, spktrain, cost=0.5):
        """
        Function to calculate the Victor-Purpura distance between two spike trains.
        See J. D. Victor and K. P. Purpura,
            Nature and precision of temporal coding in visual cortex: a metric-space
            analysis.,
            J Neurophysiol,76(2):1310-1326, 1996
        
        Inputs:
            spktrain - the other SpikeTrain
            cost     - The cost parameter. See the paper for more information
        """
        nspk_1      = len(self)
        nspk_2      = len(spktrain)
        if cost == 0: 
            return abs(nspk_1-nspk_2)
        elif cost > 1e9 :
            return nspk_1+nspk_2
        scr = numpy.zeros((nspk_1+1,nspk_2+1))
        scr[:,0] = numpy.arange(0,nspk_1+1)
        scr[0,:] = numpy.arange(0,nspk_2+1)
            
        if nspk_1 > 0 and nspk_2 > 0:
            for i in xrange(1,nspk_1+1):
                for j in xrange(1,nspk_2+1):
                    scr[i,j] = min(scr[i-1,j]+1,scr[i,j-1]+1)
                    scr[i,j] = min(scr[i,j],scr[i-1,j-1]+cost*abs(self.spike_times[i-1]-spktrain.spike_times[j-1]))
        return scr[nspk_1,nspk_2]


    def distance_kreuz(self, spktrain, dt=0.1):
        """
        Function to calculate the Kreuz/Politi distance between two spike trains
        See  Kreuz, T.; Haas, J.S.; Morelli, A.; Abarbanel, H.D.I. & Politi, A. 
            Measuring spike train synchrony. 
            J Neurosci Methods, 165:151-161, 2007

        Inputs:
            spktrain - the other SpikeTrain
            dt       - the bin width used to discretize the spike times
        
        Examples:
            >> spktrain.KreuzDistance(spktrain2)
        
        See also
            VictorPurpuraDistance
        """
        N              = (self.t_stop-self.t_start)/dt
        vec_1          = numpy.zeros(N, float)
        vec_2          = numpy.zeros(N, float)
        result         = numpy.zeros(N, float)
        idx_spikes     = numpy.array(self.spike_times/dt,int)
        previous_spike = 0
        if len(idx_spikes) > 0:
            for spike in idx_spikes[1:]:
                vec_1[previous_spike:spike] = (spike-previous_spike)*dt
                previous_spike = spike
        idx_spikes     = numpy.array(spktrain.spike_times/dt,int)
        previous_spike = 0
        if len(idx_spikes) > 0:
            for spike in idx_spikes[1:]:
                vec_2[previous_spike:spike] = (spike-previous_spike)*dt
                previous_spike = spike
        idx = numpy.where(vec_1 < vec_2)[0]
        result[idx] = vec_1[idx]/vec_2[idx] - 1
        idx = numpy.where(vec_1 > vec_2)[0]
        result[idx] = -vec_2[idx]/vec_1[idx] + 1
        return numpy.sum(numpy.abs(result))/len(result)
    
    

    ####################################################################
    ### TOO SPECIFIC METHOD ?
    ### Better documentation
    ####################################################################
    def tuning_curve(self, var_array, normalized=False, method='sum'):
        """
        Calculate a firing-rate tuning curve with respect to some variable.
        Assumes that some variable, such as stimulus orientation, varies
        throughout the recording. The values taken by this variable should be
        supplied in a numpy array `var_array`. The spike train is binned
        according to the number of values in `var_array`, e.g., if there are
        N values, the spikes are binned with a bin width
            (`self.t_stop`-`self.t_start`)/N
        so that each bin is associated with a particular value of the variable
        in `var_array`.

        The return value is a dictionary whose keys are the distinct values in
        `val_array`. The values in the dict depend on the arguments `method` and
        `normalized`.

        If `normalized` is False, the responses (bin values) are spike counts,
        if True, they are firing rates.
        If `method` == "max", each value is the maximum response for a given
        value of the variable.
        If `method` == "sum", each value is the summed response...
        If `method` == "mean", ... you get the idea.

        (If someone can rewrite this more clearly, please do so!)
        """
        binwidth = (self.t_stop - self.t_start)/len(var_array)
        time_histogram = self.time_histogram(binwidth, normalized=normalized)
        assert len(time_histogram) == len(var_array)
        tuning_curve = {}
        counts = {}
        for k, x in zip(var_array, time_histogram):
            if not tuning_curve.has_key(k):
                tuning_curve[k] = 0
                counts[k] = 0
            if method in ('sum', 'mean'):
                tuning_curve[k] += x
                counts[k] += 1
            elif method == 'max':
                tuning_curve[k] = max(x, tuning_curve[k])
            else:
                raise Exception()
        if method == 'mean':
            for k in tuning_curve.keys():
                tuning_curve[k] = tuning_curve[k]/float(counts[k])
        return tuning_curve


    ####################################################################
    ### TOO SPECIFIC METHOD ?
    ### Better documentation
    ####################################################################
    def frequency_spectrum(self, time_bin):
        """
        Returns the frequency spectrum of the time histogram together with the
        frequency axis.
        """
        hist       = self.time_histogram(time_bin, normalized=False)
        freq_spect = analysis.simple_frequency_spectrum(hist)
        freq_bin   = 1000.0/self.duration() # Hz
        freq_axis  = numpy.arange(len(freq_spect)) * freq_bin    
        return freq_spect, freq_axis    


    ####################################################################
    ### TOO SPECIFIC METHOD ?
    ### Better documentation
    ####################################################################
    def f1f0_ratio(self, time_bin, f_stim):
        """
        Returns the F1/F0 amplitude ratio where the input stimulus frequency is
        f_stim.
        """
        freq_spect = self.frequency_spectrum(time_bin)[0]
        F0 = freq_spect[0]
        freq_bin = 1000.0/self.duration()
        try:
            F1 = freq_spect[int(round(f_stim/freq_bin))]
        except IndexError:
            errmsg = "time_bin (%f) too large to estimate F1/F0 ratio for an input frequency of %f" % (time_bin, f_stim)
            errmsg += "\nFrequency_spectrum: %s" % freq_spect
            raise Exception(errmsg)
        return F1/F0
        


class SpikeList(object):
    """
    SpikeList(spikes, id_list, t_start=None, t_stop=None, dims=None)
    
    Return a SpikeList object which will be a list of SpikeTrain objects.

    Inputs:
        spikes  - a list of (id,time) tuples (id being in id_list)
        id_list - the list of the ids of all recorded cells (needed for silent cells)
        t_start - begining of the SpikeList, in ms. If None, will be infered from the data
        t_stop  - end of the SpikeList, in ms. If None, will be infered from the data
        dims    - dimensions of the recorded population, if not 1D population
    
    t_start and t_stop are shared for all SpikeTrains object within the SpikeList
    
    Examples:
        >> sl = SpikeList([(0, 0.1), (1, 0.1), (0, 0.2)], range(2))
        >> type( sl[0] )
            <type SpikeTrain>
    
    See also
        load_spikelist
    """
    #######################################################################
    ## Constructor and key methods to manipulate the SpikeList objects   ##
    #######################################################################
    def __init__(self, spikes, id_list, t_start=None, t_stop=None, dims=None):
        """
        Constructor of the SpikeList object

        See also
            SpikeList, load_spikelist
        """
        self.t_start     = t_start
        self.t_stop      = t_stop
        self.dimensions  = dims
        self.spiketrains = {}
        
        ##### First method :-)
        ## Lists are still the best. Nevertheless, now that we have very good
        ## loading times, we may start to think about using a generator instead
        ## of having to allocate all the memory for the list...
        
        import bisect
        spikes.sort()
        if len(spikes) > 0:
            ind, time = zip(*spikes)
        else:
            ind, time = [], []
        for id in id_list:
            beg = bisect.bisect_left(ind,id)
            end = bisect.bisect_right(ind,id)
            self.spiketrains[id] = SpikeTrain(time[beg:end], self.t_start, self.t_stop)
        
        if len(self) > 0 and (self.t_start is None or self.t_stop is None):
            self.__calc_startstop()

    def id_list(self):
        """ 
        Return the list of all the cells ids contained in the
        SpikeList object
        
        Examples
            >> spklist.id_list()
                [0,1,2,3,....,9999]
        """
        return numpy.array(self.spiketrains.keys())

    def copy(self):
        """
        Return a copy of the SpikeList object
        """
        spklist = SpikeList([], [], self.t_start, self.t_stop, self.dimensions)
        for id in self.id_list():
            spklist.append(id, self.spiketrains[id])
        return spklist

    def __calc_startstop(self):
        """
        t_start and t_stop are shared for all neurons, so we take min and max values respectively.
        TO DO : check the t_start and t_stop parameters for a SpikeList. Is it commun to
        all the spikeTrains within the spikelist or each spikelistes do need its own.
        """
        if len(self) > 0:
            if self.t_start is None:
                start_times = numpy.array([self.spiketrains[idx].t_start for idx in self.id_list()])
                self.t_start = numpy.min(start_times)
                logging.debug("Warning, t_start is infered from the data : %f" %self.t_start)
                for id in self.spiketrains.keys():
                    self.spiketrains[id].t_start = self.t_start
            if self.t_stop is None:
                stop_times = numpy.array([self.spiketrains[idx].t_stop for idx in self.id_list()])
  &nb