Changeset 215

Timestamp:

10/16/08 13:45:00 (3 years ago)

Author:

apdavison

Message:

Merged cleanup branch changes r198:214 into the trunk.

Location:

trunk

Files:

• INSTALL (modified) (1 diff)
• README (modified) (1 diff)
• examples/single_neuron/simple_single_neuron.py (modified) (1 diff)
• src/__init__.py (modified) (1 diff)
• src/analysis.py (modified) (1 diff)
• src/io.py (copied) (copied from branches/cleanup/src/io.py)
• src/plotting.py (modified) (3 diffs)
• src/signals.py (modified) (31 diffs)
• src/spikes.py (deleted)

trunk/INSTALL

@@ -1 @@
-NeuroTools Insallation:
+NeuroTools Insallation
+======================
 
 First download the package files:
-svn co https://neuralensemble.kip.uni-heidelberg.de/svn/NeuroTools/trunk NeuroTools
-cd NeuroTools
 
-As root
+        svn co https://neuralensemble.org/svn/NeuroTools/trunk NeuroTools
+        cd NeuroTools
 
-$ python setup.py install
+Install as root (if you want a global install)
 
-or for those without root access, to a local location, something like:
+        $ python setup.py install
 
-# python setup.py install --prefix=$HOME/opt/mystuff
+or for those without root access, install to a writable location, something like:
+
+        # python setup.py install --prefix=$HOME/opt/mystuff
 
 Then you need to add the location: 
 
-$HOME/opt/mystuff/lib/python2.4/site-packages/ 
+$HOME/opt/mystuff/lib/python2.5/site-packages/ 
 
 to your PYTHON_PAH or (in python) your sys.path
 
-Note, here lib is replaced by lib64 on 64bit systems, and python2.4 is
+Note, here lib is replaced by lib64 on 64bit systems, and python2.5 is
 (obviously) replaced by your python version.
 
-Developpers of NeuroTools may be intersted in having their lastly updated version from SVN used by python (and therefore on python's path). A solution is to symbolically link the src folder to a folder included in the path:
-cd my_local_site-packages_folder
-ln -s where_I_check-out_neuralensemble/NeuroTools/trunk/src NeuroTools
+Developpers of NeuroTools may be interested in having their lastly updated version from SVN used by python (and therefore on python's path). A solution is to symbolically link the src folder to a folder included in the path:
+
+        cd my_local_site-packages_folder
+        ln -s where_I_check-out_neuralensemble/NeuroTools/trunk/src NeuroTools
+
 and voila!
 

trunk/README

@@ -1 @@
-NeuroTools SVN : see https://facets-vm3.kip.uni-heidelberg.de/trac/NeuroTools/wiki
+For more info see the NeuroTools trac http://neuralensemble.org/trac/NeuroTools

trunk/examples/single_neuron/simple_single_neuron.py

@@ -23 +23 @@
 import pyNN.nest2 as sim
 # the link to read SpikeList files with NeuroTools
-from NeuroTools.spikes import loadSpikeList
+from NeuroTools.signals import loadSpikeList
 # using parameters utility
 from NeuroTools.parameters import ParameterSet

trunk/src/__init__.py

@@ -1 @@
-__all__ = ['analysis', 'benchmark', 'parameters', 'plotting', 'sandbox', 'spikes', 'signals', 'stgen', 'utilities']
+__all__ = ['analysis', 'benchmark', 'parameters', 'plotting', 'sandbox', 'signals', 'stgen', 'utilities', 'io']
 
 """from tables import __version__

trunk/src/analysis.py

@@ -9 +9 @@
 import os, numpy
 
+def ccf(x, y, axis=None):
+    """Computes the cross-correlation function of two series `x` and `y`.
+        Note that the computations are performed on anomalies (deviations from
+        average).
+        Returns the values of the cross-correlation at different lags.
+        Lags are given as [0,1,2,...,n,n-1,n-2,...,-2,-1] (not any more)
+        
+        :Parameters:
+            `x` : 1D MaskedArray
+                Time series.
+            `y` : 1D MaskedArray
+                Time series.
+            `axis` : integer *[None]*
+                Axis along which to compute (0 for rows, 1 for cols).
+                If `None`, the array is flattened first.
+    """
+    assert(x.ndim == y.ndim, "Inconsistent shape !")
+#    assert(x.shape == y.shape, "Inconsistent shape !")
+    if axis is None:
+        if x.ndim > 1:
+            x = x.ravel()
+            y = y.ravel()
+        npad = x.size + y.size
+        xanom = (x - x.mean(axis=None))
+        yanom = (y - y.mean(axis=None))
+        Fx = numpy.fft.fft(xanom, npad, )
+        Fy = numpy.fft.fft(yanom, npad, )
+        iFxy = numpy.fft.ifft(Fx.conj()*Fy).real
+        varxy = numpy.sqrt(numpy.inner(xanom,xanom) * numpy.inner(yanom,yanom))
+    else:
+        npad = x.shape[axis] + y.shape[axis]
+        if axis == 1:
+            if x.shape[0] != y.shape[0]:
+                raise ValueError, "Arrays should have the same length!"
+            xanom = (x - x.mean(axis=1)[:,None])
+            yanom = (y - y.mean(axis=1)[:,None])
+            varxy = numpy.sqrt((xanom*xanom).sum(1) * (yanom*yanom).sum(1))[:,None]
+        else:
+            if x.shape[1] != y.shape[1]:
+                raise ValueError, "Arrays should have the same width!"
+            xanom = (x - x.mean(axis=0))
+            yanom = (y - y.mean(axis=0))
+            varxy = numpy.sqrt((xanom*xanom).sum(0) * (yanom*yanom).sum(0))
+        Fx = numpy.fft.fft(xanom, npad, axis=axis)
+        Fy = numpy.fft.fft(yanom, npad, axis=axis)
+        iFxy = numpy.fft.ifft(Fx.conj()*Fy,n=npad,axis=axis).real
+    # We juste turn the lags into correct positions:
+    iFxy = numpy.concatenate((iFxy[len(iFxy)/2:len(iFxy)],iFxy[0:len(iFxy)/2]))
+    return iFxy/varxy
 
-def record(output, cfilename = 'SpikeTrainPlay.wav', fs=44100, enc = 'pcm26'):
-    """ record the 'sound' produced by a neuron. Takes a spike train as the
-    output.
-
-    >>> record(my_spike_train)
-
-    """
-
-
-    # from the spike list
-    simtime_seconds = (output.t_stop - output.t_start)/1000.
-    #time = numpy.linspace(0, simtime_seconds , fs*simtime_seconds)
-    (trace,time) = numpy.histogram(output.spike_times*1000., fs*simtime_seconds)
-
-
-    # TODO convolve with proper spike...
-    spike = numpy.ones((fs/1000.,)) # one ms
-
-    trace = numpy.convolve(trace, spike, mode='same')#/2.0
-    trace /= numpy.abs(trace).max() * 1.1
-
-    from scikits.audiolab import wavwrite
-    wavwrite(trace, cfilename, fs = fs, enc = enc)
-
-def play(output):
-    """
-    plays a spike list to the audio output
-
-    play(spike_list) where spike_list is a spike_list object
-
-    see playing_with_simple_single_neuron.py for a sample use
-
-    >>> play(my_spike_train)
-
-    TODO: make it possible to play multiple spike trains in stereo
-    """
-
-
-    from tempfile import mkstemp
-    fd, cfilename   = mkstemp('SpikeListPlay.wav')
-    try:
-        record(output, cfilename)
-        import pyaudio
-        import wave
-
-        chunk = 1024
-        wf = wave.open(cfilename, 'rb')
-        p = pyaudio.PyAudio()
-
-        # open stream
-        stream = p.open(format =
-                        p.get_format_from_width(wf.getsampwidth()),
-                        channels = wf.getnchannels(),
-                        rate = wf.getframerate(),
-                        output = True)
-
-        # read data
-        data = wf.readframes(chunk)
-
-        # play stream
-        while data != '':
-            stream.write(data)
-            data = wf.readframes(chunk)
-
-        stream.close()
-        p.terminate()
-    except:
-        print "Error playing the SpikeTrain "
-        # finally
-        os.remove(cfilename)
-
-    # Python 2.4 compatibility
-    # finally:
-    os.remove(cfilename)
 
 def _dict_max(D):

trunk/src/plotting.py

@@ -1 +1 @@
 
-import numpy #, pylab
-import sys
+import numpy, pylab, sys
 from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
 from matplotlib.figure import Figure
@@ -38 +37 @@
 
 
-def raster_plot(spike_list,output=None):# limits of the plot
-    import pylab
-    DATA=spike_list.as_list_id_list_time()
-    pylab.plot(DATA[1],DATA[0],'.')
-    pylab.ylabel('neuron ID')
-    pylab.xlabel('time (s)')
-    pylab.axis([spike_list.t_start, spike_list.t_stop, 0, spike_list.N])
-    if not(output==None):
-        pylab.savefig(output)
-    #else:
-        #pylab.show()
-
-
-
-
 def set_frame(ax,boollist,linewidth=2):
     assert len(boollist) == 4
@@ -65 +49 @@
             if draw:
                 ax.add_line(side)
+
+
 
 class SimpleMultiplot(object):

trunk/src/signals.py

@@ -5 +5 @@
 """
 
-import numpy
-import os
+import numpy, os, logging, re
+from NeuroTools import analysis
+from NeuroTools import io
+
 try :
     import pylab
-except Exception:
-    print "Warning: pylab not present"
-from NeuroTools import analysis
-import logging
+    ENABLE_PLOTS     = True
+except ImportError:
+    ENABLE_PLOTS     = False
+    MATPLOTLIB_ERROR = """
+    Matplotlib not detected so plots are disabled. 
+    To turn on plots, please install the Matplotlib package
+    """
+    print MATPLOTLIB_ERROR
+
 
 class SpikeTrain(object):
     """
-    This class defines a spike train as a list of the events times.
+    SpikeTrain(spikes_times, dt=None, 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.
 
-    >>> s1 = SpikeTrain([0.0, 0.1, 0.2, 0.5])
-    >>> s2 = SpikeTrain(numpy.array([0, 1, 2, 5]), dt=0.1)
-    >>> assert all(s1.spike_times == s2.spike_times)
-    >>> print s1.format(relative=True)
-    [ 0.   0.1  0.1  0.3]
-    >>> s1.isi()
-    array([ 0.1,  0.1,  0.3])
-    >>> s1.mean_rate()
-    8.0
-    >>> s1.cv_isi()
-    0.565685424949
+    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)
+        dt          - time precision of the spike events
+
+    Examples:
+        >>> s1 = SpikeTrain([0.0, 0.1, 0.2, 0.5])
+        >>> s2 = SpikeTrain([0, 1, 2, 5], dt=0.1)
+        >>> assert all(s1.spike_times == s2.spike_times)
+        >>> s1.isi()
+        array([ 0.1,  0.1,  0.3])
+        >>> s1.mean_rate()
+        8.0
+        >>> s1.cv_isi()
+        0.565685424949
     """
-    # TODO in the definition, should a spike train be ordered?
-
-    # TODO : should we handle different population shapes differently?
 
     #######################################################################
@@ -41 +51 @@
     def __init__(self, spike_times, dt=None, t_start=None, t_stop=None):
         """
-        `spike_times` is a list/numpy array of spike times (in milliseconds)
-
-        TODO: We proposed initially that :
-        If `dt` is specified, the time values should be ints,
-        and will be multiplied by `dt`, otherwise time values should be floats.
-        Markus suggested that the internal representation should be integer and
-        that analysis methods also.
-
-        If `t_start` and `t_stop` are not specified, they are inferred from the data.
-
-        the format adopted for the representation is relative=False, quantized=False
+        Constructor of the SpikeTrain object
+        
+        See also
+            SpikeTrain
 
         """
 
         if dt is not None:
-            if dt<=0:
+            if dt <= 0:
                 raise ValueError("dt must be greater than zero")
             #self.spike_times *= dt
 
-        self.dt = dt
-        self.t_start = t_start
-        self.t_stop  = t_stop
-
+        self.dt          = dt
+        self.t_start     = t_start
+        self.t_stop      = t_stop
         self.spike_times = numpy.array(spike_times, 'float')
 
@@ -69 +71 @@
         # the spikes with t >= t_start
         if self.t_start is not None:
-            idx = numpy.where(self.spike_times >= self.t_start)[0]
-            self.spike_times = self.spike_times[idx]
-            assert numpy.all(self.spike_times >= self.t_start), \
-              "Spike times out of range (t_start=%s, min(spike_times)=%s" % (self.t_start,self.spike_times.min())
+            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:
-            idx = numpy.where(self.spike_times <= self.t_stop)[0]
-            self.spike_times = self.spike_times[idx]
-            assert numpy.all(self.spike_times <= self.t_stop), \
-                "Spike times out of range (t_stop=%s, max(spike_times)=%s" % (self.t_stop,self.spike_times.max())
-
+            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)
         # 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
@@ -91 +90 @@
         size = len(spike_times)
         if size == 0:
-            if self.t_start is None:
+            if self.t_start is None: 
                 self.t_start = 0
             if self.t_stop is None:
@@ -111 +110 @@
 
         if self.t_start >= self.t_stop :
-            raise Exception("Incompatible time interval for the creation of the SpikeTrain. t_start=%s, t_stop=%s" % (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")
@@ -122 +121 @@
     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 copy(self):
+        """
+        Return a copy of the SpikeTrain object
+        """
+        return SpikeTrain(self.spike_times, self.dt, self.t_start, self.t_stop)
+    
+    def __getdisplay__(self,display):
+        """
+        Return a pylab object with a plot() function to draw the plots.
+        
+        Inputs:
+            display - if True, a new figure is created. Otherwise, if display is a
+                      subplot object, this object is returned.
+        """
+        if display is False:
+            return None
+        elif display is True:
+            pylab.figure()
+            return pylab
+        else:
+            return display
+    
+    def __labels__(self, subplot, xlabel, ylabel):
+        """
+        Function to put some labels on a plot
+        
+        Inputs:
+            subplot - the targeted plot
+            xlabel  - a string for the x label
+            ylabel  - a string for the y label
+        """
+        if hasattr(subplot, 'xlabel'):
+            subplot.xlabel(xlabel)
+            subplot.ylabel(ylabel)
+        else:
+            subplot.set_xlabel(xlabel)
+            subplot.set_ylabel(ylabel)
 
     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.concatenate((self.spike_times, spiketrain.spike_times))
+        self.spike_times = numpy.sort(self.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):
         """
-        a function to format a spike train from one format to another.
-        outputs a numpy array
-
+        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()
-        spike_times.sort()
-
-        if relative and len(spike_times)>0:
+
+        if relative and len(spike_times) > 0:
             spike_times[1:] = spike_times[1:] - spike_times[:-1]
 
@@ -145 +223 @@
         return spike_times
 
+
+
     #######################################################################
     ## Analysis methods that can be applied to a SpikeTrain object       ##
     #######################################################################
+    
     def isi(self):
-        # TODO this needs some thinking to know how to handle the border, in particular the
-        # first spike and t_start
-        return self.format(relative=True, quantized=False)[1:]
-
-    # Returns the mean firing rate of the SpikeTrain
+        """
+        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):
-        """ Mean firing rate between t_start and t_stop in Hz
-
-        NOTE: avoided calling where when defaults settings t_start=self.t_start, t_stop=self.t_stop
-        """
-        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
-            if t_stop==None: t_stop=self.t_stop
+        """ 
+        Return 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
+        """
+        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
+            if t_stop == None: t_stop=self.t_stop
             idx = numpy.where((self.spike_times >= t_start) & (self.spike_times <= t_stop))[0]
         return 1000.*len(idx)/(t_stop-t_start)
 
-    # Returns the cv of the isi of the SpikeTrain
     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)
@@ -181 +273 @@
         Poisson-type behavior. As a measure for irregularity in the network one
         can use the average irregularity across all neurons.
-
-        TODO: is it useful to get the std of CV?
+        
+        See also
+            isi
+
         """
         isi = self.isi()
-        #TODO return an exception if mean.isi= 0
         if len(isi) > 0:
             return numpy.std(isi)/numpy.mean(isi)
@@ -192 +285 @@
 
     def fano_factor_isi(self):
-        """ returns the fano factor of this spike trains ISI (see SpikeList.fano_factor)"""
+        """ 
+        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
+        
+        See also
+            isi, cv_isi
+        """
         isi = self.isi()
         if len(isi) > 0:
-            #firing_rate = self.time_histogram(time_bin,False)      
             fano = numpy.var(isi)/numpy.mean(isi)
             return fano
@@ -201 +300 @@
             raise Exception("No spikes in the SpikeTrain !")
 
-    def time_axis(self, time_bin):
+    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, 90]
+        
+        See also
+            time_histogram
+
         """
         return numpy.arange(self.t_start, self.t_stop, time_bin)
 
-    def raster_plot(self, t_start=None, t_stop=None, color='b'):
+    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. Those are not the global t_start and t_stop
-        of the SpikeTrain objects. If not defined, we use the one of the SpikeTrain
-        object
-        """
-        if t_start is None:
-            t_start = self.t_start
-        if t_stop is None:
-            t_stop = self.t_stop
-        idx = numpy.where((self.spike_times >= t_start) & (self.spike_times <= t_stop))[0]
-        spikes = self.spike_times[idx]
-        if len(spikes) > 0:
-            pylab.figure()
-            pylab.scatter(spikes,numpy.ones(len(spikes)), c=color)
-            pylab.xlabel("Time (ms)", size="x-large")
-            pylab.ylabel("Neuron #", size="x-large")
-
-    # Method to sort a SpikeTrain
-    def sort_by_time(self):
-        self.spike_time = sort(self.spike_time)
-
-    def subSpikeTrain(self, t_start, t_stop):
-        """ Returns a spike train sliced between t_start and t_stop
-        t_start and t_stop may either be single values or sequences of start
-        and stop times.
-        """
-        if hasattr(t_start, '__len__'):
-            assert len(t_start) == len(t_stop)
-            mask = False
-            for t0,t1 in zip(t_start, t_stop):
-                mask = mask | (self.spike_times >= t0) & (self.spike_times <= t1)
-            t_start = t_start[0]
-            t_stop = t_stop[-1]
-        else:
-            mask = (self.spike_times >= t_start) & (self.spike_times <= t_stop)
-        idx = numpy.where(mask)[0]
-        spikes = self.spike_times[idx]
+        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 = self.__getdisplay__(display)
+        if not subplot or not ENABLE_PLOTS:
+            print MATPLOTLIB_ERROR
+            return
+        else:
+            if len(spikes) > 0:
+                xlabel = "Time (ms)"
+                ylabel = "Neurons #"
+                self.__labels__(subplot, xlabel, ylabel)
+                subplot.scatter(spikes,numpy.ones(len(spikes)),**kwargs)
+                pylab.draw()
+
+    def add_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.add_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
+        
+        Inputs:
+            t_start - begining of the new SpikeTrain, in ms.
+            t_stop  - end of the new SpikeTrain, in ms.
+
+        """
+        spikes = numpy.extract((self.spike_times >= t_start) & (self.spike_times <= t_stop), self.spike_times)
         return SpikeTrain(spikes, self.dt, t_start, t_stop)
 
-    def time_histogram(self, time_bin , normalized=True):
+
+    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`.
-        If `normalized` is True, the bin values are scaled so as to represent
-        firing rates in spikes/second, otherwise it's the number of spikes per bin.
-        """
-        if hasattr(time_bin, '__len__'):
-            bins = time_bin
-        else:
-            bins = self.time_axis(time_bin)
-        
-        (hist, lower_edges ) = numpy.histogram(self.spike_times, bins)
+        
+        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)
+        hist, lower_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
@@ -267 +413 @@
     def relative_times(self):
         """
-        Rescale spike times to make them relative to t_start.
-        t_start becomes 0."""
+        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
-
+            self.t_stop      -= self.t_start
+            self.t_start      = 0.0
+
+    def VictorPurpuraDistance(self, spktrain, cost):
+        """
+        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 informations
+        """
+        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 KreuzDistance(self, spktrain):
+        """
+        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, 2007, 165, 151-161
+
+        Inputs:
+            spktrain - the other SpikeTrain
+        """
+        N              = (self.t_stop-self.t_start)/self.dt
+        vec_1          = numpy.zeros(N)
+        vec_2          = numpy.zeros(N)
+        result         = numpy.zeros(N)
+        idx_spikes     = self.spike_times/self.dt
+        previous_spike = 0
+        for spike in idx_spikes:
+            vec_1[previous_spike:spike] = (spike-previous_spike)*self.dt
+            previous_spike = spike
+        vec_1[previous_spike-1:N] = vec_1[previous_spike-2]
+        idx_spikes     = spktrain.spike_times/self.dt
+        previous_spike = 0
+        for spike in idx_spikes:
+            vec_2[previous_spike:spike] = (spike-previous_spike)*self.dt
+            previous_spike = spike
+        vec_2[previous_spike-1:N] = vec_2[previous_spike-2]
+        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))
+    
+    
+    def record(self, filename = 'SpikeTrain.wav', fs=44100):
+        """ 
+        Fancy tool, not really useful but still...
+        Record the 'sound' produced by a SpikeTrain. Need the scikits.audiolab
+        package
+        
+        Inputs:
+            filename - the name of the output file (should have .wav extension)
+            fs       - the sampling rate
+        
+        Examples:
+            >>> spktrain.record("my_sound.wav")
+        """
+        # from the spike list
+        simtime_seconds = (self.t_stop - self.t_start)/1000.
+        (trace,time) = numpy.histogram(self.spike_times*1000., fs*simtime_seconds)
+        # TODO convolve with proper spike...
+        spike = numpy.ones((fs/1000.,)) # one ms
+        trace = numpy.convolve(trace, spike, mode='same')#/2.0
+        trace /= numpy.abs(trace).max() * 1.1
+        try:
+            from scikits.audiolab import wavwrite
+        except ImportError:
+            print "You need the scikits.audiolab package to produce sounds !"
+        wavwrite(trace, filename, fs = fs, enc = 'pcm26')
+
+    ####################################################################
+    ### TOO SPECIFIC METHOD ?
+    ### Better documentation
+    ####################################################################
     def tuning_curve(self, var_array, normalized=False, method='sum'):
         """
@@ -320 +562 @@
         return tuning_curve
 
+
+    ####################################################################
+    ### TOO SPECIFIC METHOD ?
+    ### Better documentation
+    ####################################################################
     def frequency_spectrum(self, time_bin):
         """
@@ -325 +572 @@
         frequency axis.
         """
-        hist = self.time_histogram(time_bin, normalized=False)
+        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    
+        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):
         """
@@ -346 +598 @@
             raise Exception(errmsg)
         return F1/F0
-    
-    def merge(self, spiketrain, relative_times=False):
-        """
-        Add the spike times from `spiketrain` to this `SpikeTrain`s spike list.
-        """
-        if relative_times:
-            self.relative_times()
-            spiketrain.relative_times()
-        self.spike_times = numpy.concatenate((self.spike_times, spiketrain.spike_times))
-        self.spike_times.sort()
-        self.t_start = min(self.t_start, spiketrain.t_start)
-        self.t_stop = max(self.t_stop, spiketrain.t_stop)
-        
+        
+
 
 class SpikeList(object):
     """
-    A SpikeList object is a list of SpikeTrain objects.
-
-    >>> sl = SpikeList(3, [(0, 0.1), (1, 0.1), (0, 0.2)])
-    >>> type( sl.spiketrains[0] )
-    <type SpikeTrain>
-    >>> sl.spiketrains[0].spike_times
-    array([ 0.1,  0.2])
-    >>> sl.as_ids_times()
-    (array([0,0,1]), array([0.1,0.2,0.1]))
-    >>> sl.as_ids_times(relative=True)
-    (array([0,1,0]), array([0.1,0.0,0.1]))
-    >>> sl.as_ids_times(quantized=0.01)
-    (array([0,0,1]), array([10,20,10]))
-    >>> sl.as_list_id_time()
-    [(0,0.1), (0,0.2), (1,0.1)]
-    >>> sl.as_id_list_times()
-    [(0, array([0.1, 0.2])), (1, array([0.1]))]
-    >>> sl.as_time_list_ids()
-    [(0.1, [0,1]), (0.2, [0])]
-    >>> sl.as_2byN_array()
-    >>> array([[0,0,1],[0.1,0.2,0.1]])
+    SpikeList(spikes, id_list, dt=None, 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)
+        is_list - the list of the ids of all recorded cells (needed for silent cells)
+        dt      - if dt is specified, time values should be floats
+        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
+    
+    dt, t_start and t_stop are shared for all SpikeTrains object within the SpikeList
+    
+    Examples:
+        >>> sl = SpikeList(3, [(0, 0.1), (1, 0.1), (0, 0.2)])
+        >>> type( sl[0] )
+        >>> <type SpikeTrain>
+    
+    See also
+        loadSpikeList
 
     """
@@ -388 +629 @@
     ## Constructor and key methods to manipulate the SpikeList objects   ##
     #######################################################################
-    def __init__(self, spikes, id_list, dt=None, t_start=None, t_stop=None, label=''):
-        """
-        `spikes` is a list/tuple of (id,t) tuples (id in id_list)
-        `id_list` is the list of ids of all recorded cells (needed for silent cells)
-        If `dt`is specified, the time values should be ints,
-        and will be multiplied by `dt`, otherwise time values should be floats.
-        If `t_start` and `t_stop` are not specified, they are inferred from the data.
-
-        dt, t_start and t_stop are shared for all SpikeTrains
-
-        """
-        self.id_list = id_list
-        self.t_start = t_start
-        self.t_stop = t_stop
-        self.dt = dt
-        self.label = label
-        # transform spikes in a spike array
+    def __init__(self, spikes, id_list, dt=None, t_start=None, t_stop=None, dims=None):
+        """
+        Constructor of the SpikeList object
+
+        See also
+            SpikeList, loadSpikeList
+        """
+        self.t_start     = t_start
+        self.t_stop      = t_stop
+        self.dt          = dt
+        self.dimensions  = dims
         self.spiketrains = {}
+
         for id in id_list:
             self.spiketrains[id] = []
         for id,time in spikes:
-            if id in id_list: #id_list can be a subset of the list of recorded neurons
+            if id in self.id_list(): #id_list can be a subset of the list of recorded neurons
                 self.spiketrains[id].append(time)
 
         # writing as a list of SpikeTrains
         for id,spikes in self.spiketrains.items(): #
-            self.spiketrains[id] = SpikeTrain(spike_times=spikes, dt=self.dt, t_start=self.t_start, t_stop=self.t_stop)
+            self.spiketrains[id] = SpikeTrain(spikes, self.dt, 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 N(self):
-        return len(self.id_list)
+    def id_list(self):
+        """ 
+        Return the list of all the cells ids contained in the
+        SpikeList object
+        """
+        return numpy.array(self.spiketrains.keys())
+
+    def copy(self):
+        """
+        Return a copy of the SpikeList object
+        """
+        # Maybe not optimal, should be optimized
+        spklist = SpikeList([], [], self.dt, 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):
@@ -429 +680 @@
         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 = start_times.min()
+                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])
-                self.t_stop  = stop_times.max()
+                stop_times = numpy.array([self.spiketrains[idx].t_stop for idx in self.id_list()])
+                self.t_stop  = numpy.max(stop_times)
+                logging.debug("Warning, t_stop  is infered from the data : %f" %self.t_stop)
                 for id in self.spiketrains.keys():
                     self.spiketrains[id].t_stop = self.t_stop
@@ -441 +694 @@
             raise Exception("No SpikeTrains")
 
-    def __getitem__(self, i):
-        return self.spiketrains[i]
-
-    def __setitem__(self, i, val):
-        assert isinstance(val, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects"
-        self.spiketrains[i] = val
-        self.id_list.append(i)
+    def __getitem__(self, id):
+        return self.spiketrains[id]
+    
+    def __getslice__(self, i, j):
+        """
+        Return a new SpikeList object with all the ids between i and j
+        """
+        ids = numpy.where((self.id_list() > i) and (self.id_list() < j))
+        return self.id_slice(ids)
+    
+    #def __setslice__(self, i, j):
+
+    def __setitem__(self, id, spktrain):
+        assert isinstance(spktrain, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects"
+        self.spiketrains[id] = spktrain
         self.__calc_startstop()
 
@@ -456 +717 @@
         return len(self.spiketrains)
 
-    def append(self, id, spiketrain):
+    def __sub_id_list(self, sub_list=None):
+        if sub_list == None:
+            return self.id_list()
+        if type(sub_list) == int:
+            return numpy.random.permutation(self.id_list())[0:sub_list]
+        if type(sub_list) == list:
+            return sub_list
+    
+    def __getdisplay__(self,display):
+        """
+        Return a pylab object with a plot() function to draw the plots.
+        
+        Inputs:
+            display - if True, a new figure is created. Otherwise, if display is a
+                      subplot object, this object is returned.
+        """
+        if display is False:
+            return None
+        elif display is True:
+            pylab.figure()
+            return pylab
+        else:
+            return display
+    
+    def __labels__(self, subplot, xlabel, ylabel):
+        """
+        Function to put some labels on a plot
+        
+        Inputs:
+            subplot - the targeted plot
+            xlabel  - a string for the x label
+            ylabel  - a string for the y label
+        """
+        if hasattr(subplot, 'xlabel'):
+            subplot.xlabel(xlabel)
+            subplot.ylabel(ylabel)
+        else:
+            subplot.set_xlabel(xlabel)
+            subplot.set_ylabel(ylabel)
+
+    def append(self, id, spktrain):
         """
         Add a SpikeTrain object to the SpikeList
-        """
-        assert isinstance(spiketrain, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects"
-        if id in self.id_list:
-            raise Exception("Id already present in SpikeList.Use setitem instead()")
-        else:
-            self.spiketrains[id] = spiketrain
-            self.id_list.append(id)
-        self.t_start = min(self.t_start, spiketrain.t_start) or spiketrain.t_start # in case self.t_start is None
-        self.t_stop = max(self.t_stop, spiketrain.t_stop)
-
+        
+        Inputs:
+            id       - the id of the new cell
+            spktrain - the SpikeTrain object representing the new cell
+        
+        The SpikeTrain object is sliced according to the t_start and t_stop times
+        of the SpikeLlist object
+        
+        Examples
+            >>> st=SpikeTrain(range(0,100,5),0.1,0,100)
+            >>> spklist.append(999, st)
+            >>> spklist[999]
+        
+        See also
+            concatenate, __setitem__
+        """
+        assert isinstance(spktrain, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects"
+        if id in self.id_list():
+            raise Exception("id %d already present in SpikeList. Use __setitem__ (spk[id]=...) instead()" %id)
+        else:
+            self.spiketrains[id] = spktrain.time_slice(self.t_start, self.t_stop)
+            
     def get_time_parameters(self):
         """
-        Returns the time parameters of the SpikeList (t_start, t_stop, dt)
+        Return the time parameters of the SpikeList (t_start, t_stop, dt)
         """
         return (self.t_start, self.t_stop, self.dt)
 
-    # Same as for the SpikeTrain object
     def time_axis(self, time_bin):
+        """
+        Return a time axis between t_start and t_stop according to a time_bin
+        
+        Inputs:
+            time_bin - the bin width
+        
+        See also
+            spike_histogram
+        """
         return numpy.arange(self.t_start, self.t_stop, time_bin)
 
-    def concatenate(self, SpikeList_list):
-        """
-        Concatenation of a list of SpikeLists to the current SpikeList
-        """
+    def concatenate(self, spklists):
+        """
+        Concatenation of SpikeLists to the current SpikeList.
+        
+        Inputs:
+            spklists - could be a single SpikeList or a list of SpikeLists
+        
+        The concatenated SpikeLists must have similar (t_start, t_stop, dt), and
+        they can't shared similar cells. All their ids have to be different.
+        
+        See also
+            append, merge, __setitem__
+        """
+        if isinstance(spklists, SpikeList):
+            spklists = [spklists]
         # We check that Spike Lists have similar time_axis
-        sl_= SpikeList_list[0]
-        for sl in SpikeList_list:
+        for sl in spklists:
             if not sl.get_time_parameters() == self.get_time_parameters():
                 raise Exception("Spike Lists should have similar time_axis")
-        for sl in SpikeList_list:
-            for id in sl.id_list:
+        for sl in spklists:
+            for id in sl.id_list():
                 self.append(id, sl.spiketrains[id])
 
@@ -500 +831 @@
         """
         for id, spiketrain in spikelist.spiketrains.items():
-            if id in self.id_list:
+            if id in self.id_list():
                 self.spiketrains[id].merge(spiketrain, relative_times)
             else:
@@ -508 +839 @@
                 
     
-    def idsubSpikeList(self, id_list):
-        """
-        Generate a new SpikeList truncated from a particular sublist of cells
-        """
-        # We check what are the elements that are in self.id_list and not in
-        # id_list. We remove such elements from the SpikeList
-        new_SpkList = SpikeList([], [], self.dt, self.t_start, self.t_stop)
-        if isinstance(id_list,int):
-            id_list = numpy.random.permutation(self.id_list)[0:id_list]
+    def id_slice(self, id_list):
+        """
+        Return a new SpikeList obtained by selecting particular ids
+        
+        Inputs:
+            id_list - Can be an integer (and then N random cells will be selected)
+                      or a sublist of the current ids
+        
+        The new SpikeList inherits the time parameters (t_start, t_stop, dt)
+        
+        Examples:
+            >>> spklist.id_list
+            >>> [830, 1959, 1005, 416, 1011, 1240, 729, 59, 1138, 259]
+            >>> new_spklist = spklist.id_slice(5)
+            >>> new_spklist.id_list
+            >>> [1011, 729, 1138, 416, 59]
+
+        See also
+            time_slice
+        """
+        new_SpkList = SpikeList([], [], self.dt, self.t_start, self.t_stop, self.dimensions)
+        id_list = self.__sub_id_list(id_list)
         for id in id_list:
             try:
                 new_SpkList.append(id, self.spiketrains[id])
             except Exception:
-                print "Item %s is not in the source SpikeList" %id
+                print "id %d is not in the source SpikeList" %id
         return new_SpkList
 
-    def timesubSpikeList(self, t_start, t_stop):
-        """
-        Generate a new SpikeList truncated from particular time boundaries
-        returns a new SpikeList
-
-        """
-        new_SpkList = SpikeList([], [], self.dt, t_start, t_stop)
-        for id in self.id_list:
-            new_SpkList.append(id, self.spiketrains[id].subSpikeTrain(t_start, t_stop))
+    def time_slice(self, t_start, t_stop):
+        """
+        Return a new SpikeList obtained by slicing between t_start and t_stop
+        
+        Inputs:
+            t_start - begining of the new SpikeTrain, in ms.
+            t_stop  - end of the new SpikeTrain, in ms.
+        
+        See also
+            id_slice
+        """
+        new_SpkList = SpikeList([], [], self.dt, t_start, t_stop, self.dimensions)
+        for id in self.id_list():
+            new_SpkList.append(id, self.spiketrains[id].time_slice(t_start, t_stop))
         new_SpkList.__calc_startstop()
         return new_SpkList
+    
+    def add_offset(self, offset):
+        """
+        Add an offset to the whole SpikeList object. t_start and t_stop are
+        shifted from offset, so does all the SpikeTrain.
+         
+        Inputs:
+            offset - the time offset, in ms
+        
+        Examples:
+            >>>spklist.t_start
+            >>> 1000
+            >>>spklist.add_offset(50)
+            >>>spklist.t_start
+            >>> 1050
+        """
+        self.t_start += offset
+        self.t_stop  += offset
+        for id in self.id_list():
+            self.spiketrains[id].add_offset(offset)
+    
+    def first_spike_time(self):
+        """
+        Get the time of the first real spike in the SpikeList
+        """
+        first_spike = self.t_stop
+        is_empty    = True
+        for id in self.id_list():
+            if len(self.spiketrains[id]) > 0:
+                is_empty = False
+                if self.spiketrains[id].spike_times[0] < first_spike:
+                    first_spike = self.spiketrains[id].spike_times[0]
+        if is_empty:
+            raise Exception("No spikes can be found in the SpikeList object !")
+        else:
+            return first_spike
+    
+    def last_spike_time(self):
+        """
+        Get the time of the last real spike in the SpikeList
+        """
+        last_spike = self.t_start
+        is_empty    = True
+        for id in self.id_list():
+            if len(self.spiketrains[id]) > 0:
+                is_empty = False
+                if self.spiketrains[id].spike_times[-1] > last_spike:
+                    last_spike = self.spiketrains[id].spike_times[-1]
+        if is_empty:
+            raise Exception("No spikes can be found in the SpikeList object !")
+        else:
+            return last_spike
+    
+    
+    def select_ids(self, criteria=None):
+        """
+        Return the list of all the cells in the SpikeList that will match the criteria
+        expressed with the following syntax. 
+        
+        Inputs : 
+            criteria - a string that can be evaluated on a SpikeTrain object, where the 
+                       SpikeTrain should be named ``cell''.
+        
+        Exemples:
+            >>> spklist.select_ids("cell.mean_rate() > 0") (all the active cells)
+            >>> spklist.select_ids("cell.mean_rate() == 0") (all the silent cells)
+            >>> spklist.select_ids("len(cell.spiketimes) > 10")
+            >>> spklist.select_ids("mean(cell.isi()) < 1")
+        """
+        selected_ids = []
+        for id in self.id_list():
+            cell = self.spiketrains[id]
+            if eval(criteria):
+                selected_ids.append(id)
+        return selected_ids
+
+    
+    def save(self, filename, method="text"):
+        """
+        Save the SpikeList in a text or binary file
+        
+        Inputs:
+            filename - name of the output file
+            method   - "text"   -> a classical human readible text file
+                       "pickle" -> binary backup, much more faster to load/save, 
+                                   but no longer human readible
+                       "hdf5"   -> a compress an optimized structure, not
+                                   implemented yet
+        """
+        spike_loader = io.SpikeListIO(filename)
+        if method == "pickle":
+            spike_loader.save_spikelist_pickle(self, filename)
+        if method == "text":
+            spike_loader.save_spikelist_txt(self, filename)
 
 
@@ -542 +985 @@
 
     
-    def isi(self, nbins=100, display=False):
+    def isi(self):
         """
         Return the list of all the isi vectors for all the SpikeTrains objects
-        within the SpikeList. If display is True, then it plots the distribution
-        of the ISI
+        within the SpikeList.
+        
+        See also:
+            isi_hist
         """
         isis = []
-        for idx,id in enumerate(self.id_list):
+        for id in self.id_list():
             isis.append(self.spiketrains[id].isi())
-        if not display:
-            return isis
-        else:
-            ISI = numpy.array([])
-            for idx in xrange(self.N()):
-                ISI = numpy.concatenate((ISI,isis[idx]))
-            values, xaxis = numpy.histogram(ISI, nbins, normed=1)
-            pylab.figure()
-            pylab.plot(xaxis, values)
-            pylab.xlabel("Inter Spike Interval (ms)", size="x-large")
-            pylab.ylabel("% of Neurons", size="x-large")
-
-    def isi_hist(self, bins):
-        """Return the histogram of the ISI.
-        
-        bins may either be an integer, giving the number of bins (between the
-        min and max of the data) or a list/array containing the lower edges of
-        the bins.
-        
-        Returns a tuple, (histogram values, bin edges)
-        """
-        isis = numpy.concatenate(self.isi())
-        return numpy.histogram(isis, bins=bins)
-
-    
-    def cv_isi(self, nbins=100, display=False):
+        return isis
+
+
+    def isi_hist(self, bins=50, display=False, kwargs={}):
+        """
+        Return the histogram of the ISI.
+        
+        Inputs:
+            bins    - the number of bins (between the min and max of the data) 
+                      or a list/array containing the lower edges of the bins.
+            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)
+            >>> spklist.isi_hist(10, display=z, kwargs={'color':'r'})
+
+        See also:
+            isi
+        """
+        isis          = numpy.concatenate(self.isi())
+        values, xaxis = numpy.histogram(isis, bins=bins, normed=True)
+        subplot       = self.__getdisplay__(display)
+        if not subplot or not ENABLE_PLOTS:
+            return values, xaxis
+        else:
+            xlabel = "Inter Spike Interval (ms)"
+            ylabel = "% of Neurons"
+            self.__labels__(subplot, xlabel, ylabel)
+            subplot.plot(xaxis, values, **kwargs)
+            pylab.draw()
+
+    
+    def cv_isi(self):
         """
         Return the list of all the cv coefficients for all the SpikeTrains objects
-        within the SpikeList. If display is True, then it plots the distribution
-        of the CVs
+        within the SpikeList.
+        
+        See also:
+            cv_isi_hist
         """
         cvs_isi = []
-        for idx,id in enumerate(self.id_list):
+        for id in self.id_list():
             isi = self.spiketrains[id].isi()
             if len(isi) > 1:
                 cvs_isi.append(numpy.std(isi)/numpy.mean(isi))
-        if not display:
+        if len(cvs_isi) > 0:
             return cvs_isi
         else:
-            CV = numpy.array([])
-            for idx in xrange(len(cvs_isi)):
-                CV = numpy.concatenate((CV,[cvs_isi[idx]]))
-            values, xaxis = numpy.histogram(CV, nbins, normed=1)
-            pylab.figure()
-            pylab.plot(xaxis, values)
-            pylab.xlabel("Inter Spike Interval CV", size="x-large")
-            pylab.ylabel("% of Neurons", size="x-large")
-
-    def cv_isi_hist(self, bins):
-        cvs = numpy.array(self.cv_isi())
-        return numpy.histogram(cvs, bins=bins)
-
-    def time_axis(self, time_bin):
-        return numpy.arange(self.t_start, self.t_stop, time_bin)
+            raise Exception("No isi can be computed in the SpikeList !")
+
+
+    def cv_isi_hist(self, bins=50, display=False, kwargs={}):
+        """
+        Return the histogram of the cv coefficients.
+        
+        Inputs:
+            bins    - the number of bins (between the min and max of the data) 
+                      or a list/array containing the lower edges of the bins.
+            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)
+            >>> spklist.cv_isi_hist(10, display=z, kwargs={'color':'r'})
+        
+        See also:
+            cv_isi
+        """
+        cvs           = numpy.array(self.cv_isi())
+        values, xaxis = numpy.histogram(cvs, bins=bins, normed=True)
+        subplot       = self.__getdisplay__(display)
+        if not subplot or not ENABLE_PLOTS:
+            return values, xaxis
+        else:
+            xlabel = "Inter Spike Interval (ms)"
+            ylabel = "% of Neurons"
+            self.__labels__(subplot, xlabel, ylabel)
+            subplot.plot(xaxis, values, **kwargs)
+            pylab.draw()
+
 
     def mean_rate(self, t_start=None, t_stop=None):
         """
-        Return the mean firing rate averaged accross all SpikeTrains
-
-        see mean_rates
+        Return the mean firing rate averaged accross all SpikeTrains between t_start and t_stop.
+        
+        Inputs:
+            t_start - begining of the selected area to compute mean_rate, in ms
+            t_stop  - end of the selected area to compute mean_rate, in ms
+        
+        If t_start or t_stop are not defined, those of the SpikeList are used
+        
+        Examples:
+            >>> spklist.mean_rate()
+            >>> 12.63
+        
+        See also
+            mean_rates, mean_rate_std
         """
         return numpy.mean(self.mean_rates(t_start, t_stop))
-    
+
+
     def mean_rate_std(self, t_start=None, t_stop=None):
         """
-        Std deviation of the Mean firing rate averaged accross all SpikeTrains
-
-        see mean_rate
+        Standard deviation of the firing rates accross all SpikeTrains 
+        between t_start and t_stop
+
+        Inputs:
+            t_start - begining of the selected area to compute std(mean_rate), in ms
+            t_stop  - end of the selected area to compute std(mean_rate), in ms
+        
+        If t_start or t_stop are not defined, those of the SpikeList are used
+        
+        Examples:
+            >>> spklist.mean_rate_std()
+            >>> 13.25
+
+        See also
+            mean_rate, mean_rates
         """
         return numpy.std(self.mean_rates(t_start, t_stop))
-    
+
+
     def mean_rates(self, t_start=None, t_stop=None):
-        """ returns a vector of the size of id_list giving the mean rate for each neuron
-
-        see SpikeTrain.mean_rate
+        """ 
+        Returns a vector of the size of id_list giving the mean firing rate for each neuron
+
+        Inputs:
+            t_start - begining of the selected area to compute std(mean_rate), in ms
+            t_stop  - end of the selected area to compute std(mean_rate), in ms
+        
+        If t_start or t_stop are not defined, those of the SpikeList are used
+        
+        See also
+            mean_rate, mean_rate_std
         """
         rates = []
-        for id in self.id_list:
+        for id in self.id_list():
             rates.append(self.spiketrains[id].mean_rate(t_start, t_stop))
-
         return rates
     
-    def rate_distribution(self, nbins=25, normalize=True, display=False):
+    
+    def rate_distribution(self, nbins=25, normalize=True, display=False, kwargs={}):
         """
         Return a vector with all the mean firing rates for all SpikeTrains.
-        If display is True, then it plots the distribution of the rates
-        """
-        #rates = numpy.zeros(self.N(), float)
-        #for idx,id in enumerate(self.id_list):
-        #    rates[idx] = self.spiketrains[id].mean_firing_rate()
-        rates = self.mean_rates()
-        if not display:
+        
+        Inputs:
+            bins    - the number of bins (between the min and max of the data) 
+                      or a list/array containing the lower edges of the bins.
+            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
+        
+        See also
+            mean_rate, mean_rates
+        """
+        rates   = self.mean_rates()
+        subplot = self.__getdisplay__(display)
+        if not subplot or not ENABLE_PLOTS:
             return rates
         else:
-            values, xaxis = numpy.histogram(rates, nbins, normed=1)
-            pylab.plot(xaxis,values)
-            pylab.ylabel("% of Neurons", size="x-large")
-            pylab.xlabel("Average Firing Rate (Hz)", size="x-large")
-
-    def spike_histogram(self, time_bin, normalized=False, display=False):
+            values, xaxis = numpy.histogram(rates, nbins, normed=True)
+            xlabel = "Average Firing Rate (Hz)"
+            ylabel = "% of Neurons"
+            self.__labels__(subplot, xlabel, ylabel)
+            subplot.plot(xaxis, values, **kwargs)
+            pylab.draw()
+            
+
+
+    def spike_histogram(self, time_bin, normalized=False, display=False, kwargs={}):
         """
         Generate an array with all the spike_histograms of all the SpikeTrains
-        objects within the SpikeList. If display is True, then it plots the
-        mean firing rate of the all population along time
-        """
-        if hasattr(time_bin, '__len__'):
-            nbins = len(time_bin)
-        else:
-            nbins = numpy.ceil((self.t_stop-self.t_start)/float(time_bin))
-        spike_hist = numpy.zeros((self.N(), nbins), float)
-        logging.debug("nbins = %d" % nbins)
-        for idx,id in enumerate(self.id_list):
-            try:
-                spike_hist[idx,:] = self.spiketrains[id].time_histogram(time_bin,normalized)
-            except ValueError:
-                print spike_hist[idx,:].shape, self.spiketrains[id].time_histogram(time_bin,normalized).shape
-                print nbins, self.t_start, self.t_stop
-                print self.spiketrains[id].t_start, self.spiketrains[id].t_stop
-                raise
-        if display:
-            pylab.figure()
-            pylab.plot(self.time_axis(time_bin),sum(spike_hist)/self.N())
-            pylab.ylabel("Mean Number of Spikes per bin", size="x-large")
-            pylab.xlabel("Time (ms)", size="x-large")
-        return spike_hist
-
-    def firing_rate(self, time_bin, display=False):
-        """
-        Calculate firing rate traces (in Hz) from arrays of spike times.
-
-        Spike times and binwidth should are in milliseconds.
-        Returns a tuple (list_of_firing_rate_traces,bins)
-
-        >>> pylab.plot(output[0].time_axis(dt),sum(output.firing_rate(dt)))
-        """
-        return self.spike_histogram(time_bin, normalized=True, display=display)
-
-    # Compute the Fano Factor of the population. Need to be checked
+        objects within the SpikeList.
+        
+        Inputs:
+            time_bin   - the time bin used to gather the data
+            normalized - if True, the histogram are in Hz (spikes/second), otherwise they are
+                         in spikes/bin
+            display    - if True, a new figure is created. Could also be a subplot. The averaged
+                         spike_histogram over the whole population is then plotted
+            kwargs     - dictionary contening extra parameters that will be sent to the plot 
+                         function
+        
+        See also
+            firing_rate, time_axis
+        """
+        nbins      = self.time_axis(time_bin)
+        N          = len(self)
+        spike_hist = numpy.zeros((N, len(nbins)), float)
+        subplot    = self.__getdisplay__(display)
+        for idx,id in enumerate(self.id_list()):
+            spike_hist[idx,:] = self.spiketrains[id].time_histogram(time_bin, normalized)
+        if not subplot or not ENABLE_PLOTS:
+            return spike_hist
+        else:
+            if normalized:
+                ylabel = "Firing rate (Hz)"
+            else:
+                ylabel = "Spikes per bin"
+            xlabel = "Time (ms)"
+            self.__labels__(subplot, xlabel, ylabel)
+            subplot.plot(self.time_axis(time_bin),numpy.sum(spike_hist, axis=0)/N,**kwargs)
+            pylab.draw()
+            
+
+    def firing_rate(self, time_bin, display=False, kwargs={}):
+        """
+        Generate an array with all the instantaneous firing rates along time (in Hz) 
+        of all the SpikeTrains objects within the SpikeList.
+        
+        Inputs:
+            time_bin   - the number of bins (between the min and max of the data) 
+                         or a list/array containing the lower edges of the bins.
+            display    - if True, a new figure is created. Could also be a subplot. The averaged
+                         spike_histogram over the whole population is then plotted
+            kwargs     - dictionary contening extra parameters that will be sent to the plot 
+                         function
+        
+        See also
+            spike_histogram, time_axis
+        """
+        return self.spike_histogram(time_bin, normalized=True, display=display, kwargs=kwargs)
+
+
     def fano_factor(self, time_bin):
+        """
+        Compute the Fano Factor of the population activity.
+        
+        Inputs:
+            time_bin   - the number of bins (between the min and max of the data) 
+                         or a list/array containing the lower edges of the bins.
+        
+        The Fano Factor is computed as the variance of the averaged activity divided by its
+        mean
+        
+        See also
+            spike_histogram, firing_rate
+        """
         firing_rate = self.spike_histogram(time_bin)
-        firing_rate = sum(firing_rate)
-        fano = numpy.var(firing_rate)/numpy.mean(firing_rate)
+        firing_rate = numpy.sum(firing_rate,axis=0)
+        fano        = numpy.var(firing_rate)/numpy.mean(firing_rate)
         return fano
-    
+
+
     def fano_factors_isi(self):
-        """ returns a list containing the fano factors for each neuron"""
+        """ 
+        Return a list containing the fano factors for each neuron
+        
+        See also
+            isi, isi_cv
+        """
         fano_factors = []
-        for id in self.id_list:
+        for id in self.id_list():
             try:
                 fano_factors.append(self.spiketrains[id].fano_factor_isi())
@@ -707 +1260 @@
 
     
-    def id2position(self, id, dims):
-        """
-        Allow to translate a gid (ranging from 0 to N*M) into
-        a position on a N*M grid. dims should be given as a tuple
-        (N,M)
-        """
-        # Then we translate it in 2D
-        if len(dims) == 1:
+    def id2position(self, id):
+        """
+        Return a position (x,y) from an id if the cells are aranged on a
+        grid of size dims, as defined in the dims attribute of the SpikeList object
+        
+        Inputs:
+            id - the id of the cell
+        
+        The 'dimensions' attribute of the SpikeList must be defined
+        
+        See also
+            activity_map, activity_movie
+        """
+        if self.dimensions is None:
+            raise Exception("Dimensions of the population are not defined ! Set spikelist.dims")
+        if len(self.dimensions) == 1:
             return id
-        if len(dims) == 2:
-            x = numpy.floor(id/dims[0])
-            y = id % dims[1]
+        if len(self.dimensions) == 2:
+            x = numpy.floor(id/self.dimensions[0])
+            y = id % self.dimensions[1]
             return (x,y)
 
     
-    def activity_map(self, dims, bounds=None, display=False):
-        """
-        Generate a map of the activity during t_start and t_stop.
-        If dims is a tuple, then cells are placed on a grid of size
-        (N,M), else if dims is an array of size (2,nb_cells) with the
-        x (first line) and y (second line) flotting positions of the cells,
-        we generate a scatter plot. bounds is a parameters allowing to specify
-        the range of the colorbar
-        """
-        if isinstance(dims, tuple) or isinstance(dims, list):
-            activity_map = numpy.zeros(dims,float)
-            rates = self.mean_rates()
-            for id in self.id_list:
-                position = self.id2position(id, dims)
+    def activity_map(self, t_start=None, t_stop=None, float_positions=None, display=False, kwargs={}):
+        """
+        Generate a 2D map of the activity averaged between t_start and t_stop.
+        If t_start and t_stop are not defined, we used those of the SpikeList object
+        
+        Inputs:
+            t_start         - if not defined, the one of the SpikeList is used
+            t_stop          - if not defined, the one of the SpikeList is used
+            float_positions - None by default, meaning that the dimensions attribute 
+                              of the SpikeList is used to arange the ids on a 2D grid. 
+                              Otherwise, if the cells have flotting positions, 
+                              float_positions should be an array of size
+                              (2, nb_cells) with the x (first line) and y (second line) 
+                              coordinates of the cells
+            display         - if True, a new figure is created. Could also be a subplot. 
+                              The averaged spike_histogram over the whole population is 
+                              then plotted
+            kwargs          - dictionary contening extra parameters that will be sent 
+                              to the plot function
+
+        The 'dimensions' attribute of the SpikeList is used to turn ids into 2d positions. It should
+        therefore be not empty.
+        
+        Examples:
+            >> spklist.activity_map(0,1000,display=True)
+        
+        See also
+            activity_movie
+        """
+        subplot = self.__getdisplay__(display)
+        
+        if t_start == None: 
+            t_start = self.t_start
+        if t_stop  == None: 
+            t_stop  = self.t_stop
+        if t_start != self.t_start or t_stop != self.t_stop:
+            spklist = self.time_slice(t_start, t_stop)
+        else:
+            spklist = self
+        
+        if float_positions is None:
+            if self.dimensions is None:
+                raise Exception("Dimensions of the population are not defined ! Set spikelist.dims")
+            activity_map = numpy.zeros(self.dimensions, float)
+            rates        = spklist.mean_rates()
+            for id in spklist.id_list():
+                position = spklist.id2position(id)
                 activity_map[position] = rates[id]
-            if not display:
+            if not subplot or not ENABLE_PLOTS:
                 return activity_map
             else:
-                pylab.figure()
-                pylab.imshow(activity_map,interpolation='bicubic')
-                pylab.colorbar()
-                pylab.show()
-                if bounds is not None:
-                    pylab.clim(bounds)
-        elif isinstance(dims, numpy.ndarray):
-            if not len(self.id_list) == len(dims[0]):
-                raise Exception("Error, the number of positions does not match the number of cells in the SpikeList")
-            rates = self.mean_rates()
-            #for id in self.id_list:
-            #    rates.append(self.spiketrains[id].mean_firing_rate())
-            if not display:
+                subplot.imshow(activity_map, **kwargs)
+                subplot.colorbar()
+                pylab.draw()
+        elif isinstance(float_positions, numpy.ndarray):
+            if not len(spklist.id_list) == len(float_positions[0]):
+                raise Exception("Error, the number of flotting positions does not match the number of cells in the SpikeList")
+            rates = spklist.mean_rates()
+            if not subplot or not ENABLE_PLOTS:
                 return rates
             else:
-                x = dims[0,:]
-                y = dims[1,:]
-                pylab.scatter(x,y,c=rates)
-                pylab.colorbar()
-                if bounds is not None:
-                    pylab.clim(bounds)
-
-    def pairwise_correlations(self, pairs, time_bin=1., display=False):
-        """
-        Function to generate an array of cross correlation between computed
-        between pairs of cells within the SpikeTrains. pairs should be therefore
-        a list of (cell_id_1, cell_id_2). If display is True, then if plots the
-        averaged cross correlation over all those pairs
-        """
-        if len(pairs[0]) != len(pairs[1]):
-            raise Exception("Pairs should have the same number of elements")
-        nb_pairs = len(pairs[0])
-        spk1 = self.idsubSpikeList(pairs[0])
-        spk2 = self.idsubSpikeList(pairs[1])
+                x = float_positions[0,:]
+                y = float_positions[1,:]
+                subplot.scatter(x,y,c=rates, **kwargs)
+                subplot.colorbar()
+                pylab.draw()
+
+
+    def pairwise_cc(self, pairs, spklist=None, time_bin=1., averaged=False, display=False, kwargs={}):
+        """
+        Function to generate an array of cross correlations computed
+        between pairs of cells within the SpikeTrains.
+        
+        Inputs:
+            pairs    - Could be an int, and then N random pairs of cells will be selected, 
+                       or it could be a list of tuple (id cell_1, id cell_2)
+            time_bin - The time bin used to gather the spikes
+            spklist  - The other SpikeList object where cells should be taken. By default, 
+                       if None, this is the one calling the function
+            averaged - If true, only the averaged CC among all the pairs is returned (less memory needed)
+            display  - if True, a new figure is created. Could also be a subplot. The averaged
+                       spike_histogram over the whole population is then plotted
+            kwargs   - dictionary contening extra parameters that will be sent to the plot 
+                       function
+        """
+        subplot = self.__getdisplay__(display)
+        
+        if spklist == None:
+            spklist = self
+        
+        if type(pairs) == int:
+            spk1 = self.id_slice(pairs)
+            spk2 = spklist.id_slice(pairs)
+            N    = pairs
+        else:
+            pairs  = numpy.array(pairs)
+            N      = len(pairs[:,0])
+            spk1   = self.id_slice(pairs[:,0])
+            spk2   = spklist.id_slice(pairs[:,1])
         hist_1 = spk1.spike_histogram(time_bin)
         hist_2 = spk2.spike_histogram(time_bin)
-        print spk1, spk2
         length = 2*len(hist_1[0])
-        results = numpy.zeros((nb_pairs,length), float)
-        for idx in xrange(nb_pairs):
+        if not averaged:
+            results = numpy.zeros((nb_pairs,length), float)
+        else:
+            results = numpy.zeros(length, float)
+        for idx in xrange(N):
             # We need to avoid empty spike histogram, otherwise the ccf function
             # will give a nan vector
-            if sum(hist_1[idx]) > 0 and sum(hist_2[idx] > 0):
-                results[idx,:] = ccf(hist_1[idx],hist_2[idx])
-        if (display):
-            results = sum(results)/nb_pairs
-            pylab.figure()
-            xaxis  = time_bin*numpy.arange(-len(results)/2, len(results)/2)
-            pylab.plot(xaxis, results)
-            pylab.xlabel("Time (ms)", size="x-large")
-            pylab.ylabel("Cross Correlation", size="x-large")
-        else:
-            return results
+            if numpy.sum(hist_1[idx]) > 0 and numpy.sum(hist_2[idx] > 0):
+                if not averaged:
+                    results[idx,:] = analysis.ccf(hist_1[idx],hist_2[idx])
+                else:
+                    results += analysis.ccf(hist_1[idx],hist_2[idx])
+        if not subplot or not ENABLE_PLOTS:
+            if not averaged:
+                return results
+            else:
+                return results/N
+        else:
+            if averaged:
+                results = results/N
+            else:
+                results = numpy.sum(results, axis=0)/N
+            xaxis   = time_bin*numpy.arange(-len(results)/2, len(results)/2)
+            xlabel  = "Time (ms)"
+            ylabel  = "Cross Correlation"
+            subplot.plot(xaxis, results)
+            self.__labels__(subplot, xlabel, ylabel)
+            pylab.draw()
+
+
+    def CrossCorrZero(self, pairs, spklist=None, time_bin=1., mode="auto"):
+        #Just the Aertsen Value at 0 delay
+        t_start = float(time_parameters[0])
+        t_stop  = float(time_parameters[1])
+        
+        spk1 = CutSpikeList(f1, time_parameters)
+        spk2 = CutSpikeList(f2, time_parameters)
+    
+        if ((size(spk1)==0) or (size(spk2)==0)):
+            return 0
+        
+        if mode == "auto":
+            if (selected_set==None):
+                selected_set = GetCells(spk1)
+            if isinstance(selected_set,int):
+                all_positions = GetCells(spk1)
+                selected_set  = numpy.random.permutation(all_positions)[0:selected_set]
+            spk1 = CutSpikeList2(spk1,selected_set)
+            spk2 = CutSpikeList2(spk2,selected_set)
+            cells_id = sort(selected_set)
+            for idx in xrange(len(spk1[0])):
+                spk1[1][idx] = numpy.where(cells_id == spk1[1][idx])[0]
+            for idx in xrange(len(spk2[0])):
+                spk2[1][idx] = numpy.where(cells_id == spk2[1][idx])[0]
+            nb_neur   = len(selected_set)
+        
+        if mode == "cross":
+            # If nothing is entered, we perform the cross correlation between the two
+            # whole spike trains
+            if (selected_set==None):
+                cell_1 = GetCells(spk1)
+                cell_2 = GetCells(spk2)
+                N      = min(len(cell_1),len(cell_2))
+                cell_1 = cell_1[0:N]
+                cell_2 = cell_2[0:N]
+            
+            # If we want to select only X pairs of cells to performs the auto
+            # correlation, we select those random pairs.
+            if isinstance(selected_set,int):
+                N = selected_set
+                all_positions = GetCells(spk1)
+                shuffle       = numpy.random.permutation(all_positions)
+                cell_1        = shuffle[0:N]
+                # If we perform a cross correlation within the same spike file,
+                # we prevent the same neuron to be in the two signals, to have a real
+                # non overlapping cross correlation
+                if not(any(equal(spk1,spk2) == False)):
+                    shuffle = numpy.random.permutation(all_positions)
+                    cell_2  = shuffle[0:N]
+                else:
+                    # Otherwise we don't care...
+                    all_positions = GetCells(spk2)
+                    shuffle = numpy.random.permutation(all_positions)
+                    cell_2  = shuffle[0:N]
+            spk1 = CutSpikeList2(spk1,cell_1)
+            spk2 = CutSpikeList2(spk2,cell_2)
+            cells_id_1 = sort(cell_1)
+            cells_id_2 = sort(cell_2)
+            for count, cell in enumerate(cells_id_1):
+                spk1[1][numpy.where(cell == spk1[1])] = count
+            for count, cell in enumerate(cells_id_2):
+                spk2[1][numpy.where(cell == spk2[1])] = count
+            nb_neur = len(cells_id_1)
+        
+        
+        num_bins  = int((t_stop-t_start)/time_bin)+1
+        mat_neur1 = numpy.zeros((num_bins,nb_neur))
+        mat_neur2 = numpy.zeros((num_bins,nb_neur))
+        
+        timespk1 = numpy.array(((spk1[0] - t_start)/time_bin),'int')
+        timespk2 = numpy.array(((spk2[0] - t_start)/time_bin),'int')
+        
+        mat_neur1[timespk1,spk1[1]] = 1
+        mat_neur2[timespk2,spk2[1]] = 1
+        Z = float(sum(sum(mat_neur1*mat_neur2)))
+        
+        N = float(num_bins*nb_neur)
+        X = float(len(spk1[0]))
+        Y = float(len(spk2[0]))
+        
+        X=X/N
+        Y=Y/N
+        Z=Z/N
+        
+        return (Z-X*Y)/sqrt((X*(1-X))*(Y*(1-Y)))
+    
+    def VictorPurpuraDistance(self, id_list, spklist=None, cost=0.5):
+        """
+        Function to calculate the Victor-Purpura distance averaged over N pairs in the SpikeList
+        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:
+            id_list - Could be an int, and then N random pairs will be selected, 
+                      or it could be a list cells ids
+            spklist - The other SpikeList object where cells should be taken. By default, 
+                      if None, this is the one calling the function
+            cost    - The cost parameter. See the paper for more informations. BY default, set to 0.5
+
+        """
+        if spklist == None:
+            spklist = self
+        
+        if type(id_list) == int:
+            spk1 = self.id_slice(id_list)
+            spk2 = spklist.id_slice(spk1.id_list)
+            N    = id_list
+        else:
+            N    = len(id_list)
+            spk1 = self.id_slice(id_list)
+            spk2 = spklist.id_slice(id_list)
+        
+        distance   = 0.
+        
+        for idx in xrange(len(spk1)):
+            distance += spk1[spk1.id_list[idx]].VictorPurpuraDistance(spk2[spk2.id_list[idx]], cost)
+        return distance/N
+    
+    
+    def KreuzDistance(self, id_list, spklist=None):
+        """
+        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, 2007, 165, 151-161
+
+        Inputs:
+            id_list - Could be an int, and then N random pairs will be selected, 
+                      or it could be a list cells ids
+            spklist  - The other SpikeList object where cells should be taken. By default, 
+                       if None, this is the one calling the function
+        """
+        if spklist == None:
+            spklist = self
+        
+        if type(id_list) == int:
+            spk1 = self.id_slice(id_list)
+            spk2 = spklist.id_slice(spk1.id_list)
+            N    = id_list
+        else:
+            N    = len(id_list)
+            spk1 = self.id_slice(id_list)
+            spk2 = spklist.id_slice(id_list)
+        
+        distance   = 0.
+        
+        for idx in xrange(len(spk1)):
+            distance += spk1[spk1.id_list[idx]].KreuzDistance(spk2[spk2.id_list[idx]])
+        return distance/N
 
     def mean_rate_variance(self, time_bin):
         """
-        Function to extract the variance of the firing rate along time,
+        Return the standard deviation of the firing rate along time,
         if events are binned with a time bin.
+        
+        Inputs:
+            time_bin - time bin to bin events
+        
+        See also
+            mean_rate, mean_rates, mean_rate_covariance, firing_rate
         """
         firing_rate = self.firing_rate(time_bin)
-        return numpy.var(sum(firing_rate)/self.N())
-
-    def mean_rate_covariance(self, SpkList, time_bin):
-        """
-        Function to extract the covariance of the firing rate along time,
-        if events are binned with a time bin. We need a second SpikeList
-        object with same time parameters to calculate the covariance
-        """
-        if not isinstance(SpkList, SpikeList):
+        return numpy.var(numpy.sum(firing_rate, axis=0)/len(self))
+
+    def mean_rate_covariance(self, spikelist, time_bin):
+        """
+        Return the covariance of the firing rate along time,
+        if events are binned with a time bin.
+        
+        Inputs:
+            spikelist - the other spikelist to compute covariance 
+            time_bin  - time bin to bin events
+        
+        See also
+            mean_rate, mean_rates, mean_rate_variance, firing_rate
+        """
+        if not isinstance(spikelist, SpikeList):
             raise Exception("Error, argument should be a SpikeList object")
-        if not SpkList.get_time_parameters() == self.get_time_parameters():
+        if not spikelist.get_time_parameters() == self.get_time_parameters():
             raise Exception("Error, both SpikeLists should share common t_start, t_stop and dt")
         frate_1 = self.firing_rate(time_bin)
-        frate_1 = sum(frate_1)/self.N()
-        frate_2 = SpkList.firing_rate(time_bin)
-        frate_2 = sum(frate_2)/SpkList.N()
+        frate_1 = numpy.sum(frate_1, axis=0)/len(self)
+        frate_2 = spikelist.firing_rate(time_bin)
+        frate_2 = numpy.sum(frate_2, axis=0)/len(spikelist)
         N = len(frate_1)
-        cov = sum(frate_1*frate_2)/N-sum(frate_1)*sum(frate_2)/(N*N)
+        cov = numpy.sum(frate_1*frate_2)/N-numpy.sum(frate_1)*numpy.sum(frate_2)/(N*N)
         return cov
-    
-    def raster_plot(self, id_list=None, t_start=None, t_stop=None, colors='b', subplot=None, size=1):
-        """
-        Generate a raster plot of a certain number of cells in the
-        SpikeList object. If id_list is an integer, then N ids will be randomly choosen
-        in id_list. If this is a list, those id will be selected.
-        Raster is made between t_start and t_stop (region of interest, not the global ones)
-        and colors can be a list of color (each for one cells) or a single string
-        to apply the same color to all the raster plots.
-        """
+
+
+    def raster_plot(self, id_list=None, t_start=None, t_stop=None, display=True, kwargs={}):
+        """
+        Generate a raster plot for the SpikeList in a subwindow of interest,
+        defined by id_list, t_start and t_stop. 
+        
+        Inputs:
+            id_list - can be a integer (and then N cells are randomly selected) or a list of ids. If None, 
+                      we use all the ids of the SpikeList
+            t_start - in ms. If not defined, the one of the SpikeList object is used
+            t_stop  - in ms. If not defined, the one of the SpikeList 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)
+            >>> spikelist.raster_plot(display=z, kwargs={'color':'r'})
+        
+        See also
+            SpikeTrain.raster_plot
+        """
+        subplot = self.__getdisplay__(display)
         if id_list == None: 
-            id_list = self.id_list
-            spk = self.idsubSpikeList(id_list)
-        elif id_list != self.id_list:
-            spk = self.idsubSpikeList(id_list)
-            id_list = spk.id_list
-        else:
-            spk = self.idsubSpikeList(id_list)
-
+            id_list = self.id_list()
+            spk = self
+        else:
+            spk = self.id_slice(id_list)
+
+        if t_start is None: t_start = spk.t_start
+        if t_stop is None:  t_stop  = spk.t_stop
+        if t_start != spk.t_start or t_stop != spk.t_stop:
+            spk = spk.time_slice(t_start, t_stop)
+
+        if not subplot or not ENABLE_PLOTS:
+            print MATPLOTLIB_ERROR
+        else:
+            ids, spike_times = spk.convert(format="[ids, times]")
+            idx = numpy.where((spike_times >= t_start) & (spike_times <= t_stop))[0]
+            if len(spike_times) > 0:
+                subplot.scatter(spike_times, ids, **kwargs)
+            xlabel = "Time (ms)"
+            ylabel = "Neuron #"
+            self.__labels__(subplot, xlabel, ylabel)
+            min_id = numpy.min(spk.id_list)
+            max_id = numpy.max(spk.id_list)
+            length = t_stop - t_start
+            subplot.axis([t_start-0.05*length, t_stop+0.05*length, min_id-2, max_id+2])
+            pylab.draw()
+
+
+
+    def activity_movie(self, time_bin=10, t_start=None, t_stop=None, float_positions=None, output="animation.mpg", bounds=(0,5), fps=10, display = True, kwargs={}):
+        """
+        Generate a movie of the activity between t_start and t_stop.
+        If t_start and t_stop are not defined, we used those of the SpikeList object
+        
+        Inputs:
+            time_bin        - time step to bin activity during the movie. 
+                              One frame is the mean rate during time_bin
+            t_start         - if not defined, the one of the SpikeList is used, in ms
+            t_stop          - if not defined, the one of the SpikeList is used, in ms
+            float_positions - None by default, meaning that the dimensions attribute of the SpikeList
+                              is used to arange the ids on a 2D grid. Otherwise, if the cells have
+                              flotting positions, float_positions should be an array of size
+                              (2, nb_cells) with the x (first line) and y (second line) coordinates of
+                              the cells
+            output          - The filename to store the movie
+            bounds          - The common color bounds used during all the movies frame. 
+                              This is a tuple
+                              of values (min, max), in spikes per frame.
+            fps             - The number of frame per second in the final movie
+            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
+
+        The 'dimensions' attribute of the SpikeList is used to turn ids into 2d positions. It should
+        therefore be not empty.
+        
+        Examples:
+            >> spklist.activity_movie(10,0,1000,bounds=(0,5),display=subplot(221),output="test.mpg")
+        
+        See also
+            activity_map
+        """
+        subplot = self.__getdisplay__(display)
         if t_start is None: t_start = self.t_start
-        if t_stop is None:  t_stop = self.t_stop
-
-        if subplot is None:
-            pylab.figure()
-            subplot = pylab
-        if isinstance(colors, str): # this is much, much faster than a different colour per row
-            ids, spike_times = self.as_ids_times()
-            idx = numpy.where((spike_times >= t_start) & (spike_times <= t_stop))[0]
-            spike_times = spike_times[idx]
-            ids = ids[idx]
-            if len(spike_times) > 0:
-                print "Plotting %d points for %s" % (len(spike_times), self.label)
-                subplot.scatter(spike_times, ids, s=size, c=colors)
-        elif len(colors) != len(id_list):
-            raise Exception("The numbers of colors does not match the number of cells!")
-        else: # this is very slow in interactive mode
-            for count,id in enumerate(spk.id_list):
-                st = spk.spiketrains[id]
-                idx = numpy.where((st.spike_times >= t_start) & (st.spike_times <= t_stop))[0]
-                spikes = st.spike_times[idx]
-                if len(spikes) > 0:
-                    subplot.scatter(spikes,count*numpy.ones(len(spikes)), s=size, c=colors[count])
-        xlabel = "Time (ms)"
-        ylabel = "Neuron #"
-        if hasattr(subplot, 'xlabel'):
-            subplot.xlabel(xlabel, size="x-large")
-            subplot.ylabel(ylabel, size="x-large")
-        else:
-            subplot.set_xlabel(xlabel, size="x-large")
-            subplot.set_ylabel(ylabel, size="x-large")
-        subplot.axis([t_start-10, t_stop+10, -1, len(self)+1])
-
-    # Clearly not optimal, but at least this is working... The best way to
-    # do that function would be to go through the sorted list of all the spike times
-    def activity_movie(self, dims, time_bin=10, bounds = (0,20), output="animation.mpg", fps=10):
-        time  = self.t_start
-        files = []
-        while (time < self.t_stop):
-            subSpkList = self.timesubSpikeList(time, time + time_bin)
-            activity_map = subSpkList.activity_map(dims, bounds, display=True)
-            caption = time+time_bin/2
-            pylab.title("Time = %g ms" %caption)
-            fname = "_tmp_%g.jpg" %caption
-            print " Saving Frame", fname
-            pylab.savefig(fname)
-            pylab.close()
-            files.append(fname)
-            time += time_bin
-        print 'Making movie %s - this make take a while' %output
-        command = "mencoder 'mf://_tmp_*.jpg' -mf type=jpg:fps=%d -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o %s" %(fps,output)
-        print command
-        os.system(command)
-        ## cleanup
-        print "Clean up...."
-        for fname in files: os.remove(fname)
-
-    def pw_corr_pearson(self,edge,bins,number_of_neuron_pairs): 
-        """
-        TODO: document/test this function
-
-        """
-        #bins = edge[1]-edge[0]
-        cor = numpy.zeros((number_of_neuron_pairs,))
-        [neuron_ids_array, spike_times_array] = self.as_list_id_list_time()
-        
-        # Pairwise correlation
-        neuron_ids_unique = numpy.unique(neuron_ids_array)
-        
-        if len(neuron_ids_unique)==0:
-            return (0,0)
-        
-        for count in range(number_of_neuron_pairs):
-            # draw two neuron radomly out of the neuron_id_unique pool
-            neuron_1_tmp = neuron_ids_unique[numpy.floor(numpy.random.uniform()*len(neuron_ids_unique))]
-            neuron_2_tmp = neuron_ids_unique[numpy.floor(numpy.random.uniform()*len(neuron_ids_unique))]
-            # get the spike_times of the neurons
-            spike_times_1_tmp = self.spiketrains[neuron_1_tmp].spike_times
-            spike_times_2_tmp = self.spiketrains[neuron_2_tmp].spike_times
-
-            # hist of the spiketrains
-            n1_hist = numpy.histogram(spike_times_1_tmp,bins=bins,range=edge)
-            n2_hist = numpy.histogram(spike_times_2_tmp,bins=bins,range=edge)
+        if t_stop is None:  t_stop  = self.t_stop
+        if not subplot or not ENABLE_PLOTS:
+            print MATPLOTLIB_ERROR
+        else:
+            files        = []
+            activity_map = numpy.zeros(self.dimensions)
+            im           = subplot.imshow(activity_map, **kwargs)
+            im.set_clim(bounds[0],bounds[1])
+            im.colorbar()
+            count     = 0
+            idx       = 0
+            manager   = pylab.get_current_fig_manager()
+            if t_start != self.t_start or t_stop != self.t_stop:
+                spk   = spk.time_slice(t_start, t_stop)
+            else:
+                spk   = self
+            time, pos = spk.convert("times, ids")
+            # We sort the spikes to allow faster process later
+            sort_idx  = time.ravel().argsort()
+            time      = time[sort_idx]
+            pos       = pos[sort_idx]
+
+            if float_positions is None:
+                if self.dimensions is None:
+                    raise Exception("Dimensions of the population are not defined ! Set spikelist.dims")
+            while (t_start < t_stop):
+                activity_map = numpy.zeros(spk.dimensions)
+                while (time[idx] < t_start + time_bin):
+                    addr = spk.id2position(pos[idx])
+                    activity_map[addr] += 1
+                    idx += 1
+                im.set_array(activity_map)
+                subplot.title("time = %d ms" %t_start)
+                manager.canvas.draw()
+                im.set_clim(bounds[0],bounds[1])
+                fname = "_tmp_spikes_%05d.png" %count
+                print " Saving Frame", fname
+                pylab.savefig(fname)
+                files.append(fname)
+                t_start += time_bin
+                count += 1
+            print 'Making movie %s - this make take a while' %output
+            command = "mencoder 'mf://_tmp_*.jpg' -mf type=jpg:fps=%d -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o %s" %(fps,output)
+            print command
+            os.system(command)
+            ## cleanup
+            print "Clean up...."
+            for fname in files: os.remove(fname)
+
+
+    def pairwise_pearson_corrcoeff(self, pairs, spklist=None, time_bin=1.): 
+        """
+        Function to return the mean and the variance of the pearson correlation coefficient. 
+        For more details, see Kumar et al, ....
+        
+        Inputs:
+            pairs    - an int to specify the number of pairs used to estimate the correlation coefficient
+            time_bin - The time bin used to gather the spikes
+        """
+        if spklist == None:
+            spklist = self
+        
+        if type(pairs) == int:
+            spk1 = self.id_slice(pairs)
+            spk2 = spklist.id_slice(pairs)
+            N    = pairs
+        else:
+            pairs  = numpy.array(pairs)
+            N      = len(pairs[:,0])
+            spk1   = self.id_slice(pairs[:,0])
+            spk2   = spklist.id_slice(pairs[:,1])
+        
+        cor          = numpy.zeros(N, float)
+        
+        ## We have to extract only the non silent cells, to avoied problems
+        non_silent = spk1.select_ids("cell.mean_rate() > 0")
+        spk1 = spk1.id_slice(non_silent)
+        non_silent = spk2.select_ids("cell.mean_rate() > 0")
+        spk2 = spk2.id_slice(non_silent)
+        
+        for count in range(N):
+            # draw two neuron randomly out of the neuron_id_unique pool
+            neuron_1 = self.id_list()[numpy.floor(numpy.random.uniform()*len(self))]
+            neuron_2 = self.id_list()[numpy.floor(numpy.random.uniform()*len(self))]
+
+            # get the spike_times of the selected neurons
+            spk_1 = self.spiketrains[neuron_1]
+            spk_2 = self.spiketrains[neuron_2]
+
+            # get the histogram of the spiketrains
+            n1_hist = spk_1.time_histogram(time_bin)
+            n2_hist = spk_2.time_histogram(time_bin)
 
             # correlation
             # TODO: normalize the cor, look in 1.6 in Kumar et al.
             # bruederle: the function corrcoeff actually implements the definition in Kumar 1.6
-            cov = numpy.corrcoef(n1_hist[0],n2_hist[0])[1][0]
-            #print n1_hist[0]
-            #print n2_hist[0]
+            cov = numpy.corrcoef(n1_hist,n2_hist)[1][0]
 
             # bruederle: the expression 'cov' is already, per definition, the pearson correlation coefficient 
@@ -937 +1798 @@
 
         cor_coef_mean = cor.mean()
-        cor_coef_std = cor.std()
-        return cor_coef_mean, cor_coef_std
-
+        cor_coef_std  = cor.std()
+        return (cor_coef_mean, cor_coef_std)
+
+
+    ####################################################################
+    ### TOO SPECIFIC METHOD ?
+    ### Better documentation
+    ####################################################################
     def f1f0_ratios(self, time_bin, f_stim):
         """
@@ -950 +1816 @@
         return f1f0_dict
 
-        # methods for different output formats
-    def as_ids_times(self, relative=False, quantized=False):
-        """Returns (array of ids, array of times)."""
-        spikes = numpy.concatenate([st.spike_times for st in self.spiketrains.itervalues()])
-        ids = numpy.concatenate([id*numpy.ones((len(st.spike_times),)) for id,st in self.spiketrains.iteritems()])
-        assert len(ids) == len(spikes)
-        return ids, spikes
-
-    def as_list_id_time(self, relative=False, quantized=False):
-        """
-        Returns a list/tuple of (id,t) tuples (id in range(0,N)).
-
-        """
-        spikes =[]
-        for id in self.id_list:
-            for time in self.spiketrains[id].spike_times :
-                spikes.append((id, time))
-        return spikes
-
-    def as_list_id_list_time(self, relative=False, quantized=False):
-        """
-        Returns (list of ids, list of times).
-
-        """
-        spikes = self.as_list_id_time()
-        ids, times = [], []
-        for spike in spikes:
-            ids.append(spike[0])
-            times.append(spike[1])
-
-        return [ids, times]
-
-    def as_id_list_times(self, relative=False, quantized=False):
-        pass
-
-    def as_time_list_ids(self, relative=False, quantized=False):
-        pass
-
-    def as_2byN_array(self, relative=False, quantized=False):
-        pass
-
-    def as_spikematrix(self):
-        """ Returns a list giving for each neuron how many spikes per
-        time bin, usually zero everywhere, one if a spike.
-        Useful for plots etc.
-
-        """
-        return self.firing_rate(self.dt)
-
-    def as_pyNN_SpikeArray(self):
-        """
-        to use in conjunction with SpikeSourceArray neurons
-
-        >>> pop = Population((10,),SpikeSourceArray, {'spike_times': sl.as_pyNN_SpikeArray() })
-
-        """
-        spike_array = list([ [] for i in range(self.N())])
-        for spike in spikes:
-            spike_array[spike[0]].append(spike[1])
-
-        return spike_array
-
-"""
-SpikeLists functions
-
-"""
-
-def population2spikelist(output, N, dt , t_start, t_stop):
+
+
+    #######################################################################
+    ## Method to convert the SpikeList into several others format        ##
+    #######################################################################
+
+    def convert(self, format="[times, ids]", relative=False, quantized=False):
+        """
+        Return a new representation of the SpikeList object, in a user designed format.
+            format is an expression containing either the keywords times and ids, 
+            time and id.
+       
+        Inputs:
+            relative -  a boolean to say if a relative representation of the spikes 
+                        times compared to t_start is needed
+            quantized - a boolean to round the spikes_times.
+        
+        Examples: 
+            >>>spk.convert("[times, ids]") will return a list of two elements, the 
+                first one being the array of all the spikes, the second the array of all the
+                corresponding ids
+            >>>spk.convert("[(time,id)]") will return a list of tuples (time, id)
+        
+        See also
+            SpikeTrain.format
+        """
+        is_times = re.compile("times")
+        is_ids   = re.compile("ids")
+        times  = numpy.concatenate([st.format(relative, quantized) for st in self.spiketrains.itervalues()])
+        ids    = numpy.concatenate([id*numpy.ones(len(st.spike_times), int) for id,st in self.spiketrains.iteritems()])
+        if is_times.search(format):
+            if is_ids.search(format):
+                return eval(format)
+            else:
+                raise Exception("You must have a format with [times, ids] or [time, id]")
+        is_times = re.compile("time")
+        is_ids   = re.compile("id")
+        if is_times.search(format):
+            if is_ids.search(format):
+                result = []
+                for id, time in zip(ids, times):
+                    result.append(eval(format))
+            else:
+                raise Exception("You must have a format with [times, ids] or [time, id]")
+            return result
+
+
+
+
+class AnalogSignal(object):
     """
-    TODO sl2pynn
+    AnalogSignal(signal, dt, t_start=None, t_stop=None)
+    
+    Return a AnalogSignal object which will be a analog signal trace
+
+    Inputs:
+        signal  - the vector with the data of the AnalogSignal
+        dt      - the time step between two data points of the sampled analog signal
+        t_start - begining of the signal, in ms. If None, will be set to 0
+        t_stop  - end of the SpikeList, in ms. If None, will be infered from the data
+    
+    Examples:
+        >>> as = AnalogSignal(range(100), dt=0.1, t_start=0, t_stop=10)
+
+    See also
+        AnalogSignalList, load_currenttraces, load_membranetrace, load_conductancetraces
+
     """
-    import os, tempfile
-
-    tmpdir = tempfile.mkdtemp()
-    output_filename=os.path.join(tmpdir,'output.gdf')
-    output.printSpikes(output_filename)#
-    output_DATA = loadSpikeList(output_filename, N= N, dt = dt, t_start=t_start, t_stop=t_stop)
-    os.remove(output_filename)
-    os.rmdir(tmpdir)
-    return output_DATA
-
-def readFile(filename, sepchar = "\t", skipchar = '#'):
-    """
-    Function to read data. Should be optimize to deal with errors in the file
-    that may append sometimes, when nest1 does not save the entire last line.
-    It assumes that the data have been produced by pyNN under the format
-    time (ms) gids (for raster)
-    value (unit) gids (for continuous recordings like Vm, current, conductances)
-    """
-    myfile = open(filename, "r", 10000)
-    contents = myfile.readlines()
-    myfile.close()
-    data = []
-    for line in contents:
-        stripped_line = line.lstrip()
-        if (len(stripped_line) != 0):
-            if (stripped_line[0] != skipchar):
-                items = stripped_line.split(sepchar)
-                data.append([int(float(items[1])), float(items[0])])
-    return(data)
-
-def get_header(filename):
-    cmd = ''
-    f = open(filename, 'r')
-    for line in f.readlines():
-        if line[0] == '#':
-            cmd += line[1:].strip() + ';'
-    f.close()
-    header = {}
-    exec cmd in None, header
-    return header
-
-def loadSpikeList(filename, id_list=None, dt = None, t_start=None, t_stop=None):
-    """
-    Returns a SpikeList object from the tmp file saved by PyNN.
-
-    The input is the filename where pyNN stored the spike list and the discretization step dt.
-    it returns the spike list represented as a list of events. Events are tuples (relative time
-    since last event, neuron_id); neuron_id and time are integers
-
-    All times in milliseconds.
-
-    The PyNN format (with compatible_output=True) is:
-    * comment lines containing header information,
-    * then one line per event:  absolute time in ms, GID
-    (see function printSpikes in nest1.py)
-
-    """
-    if id_list is None: # try to obtain id_list from file header
-        header = get_header(filename)
-        if 'first_id' in header:
-            n_cells = int(header['last_id']) - int(header['first_id']) + 1
-            id_list = n_cells
-        else:
-            raise Exception("id_list not provided and cannot be inferred from file header.")
-    if isinstance(id_list,int): # allows to just specify the number of neurons
-        id_list = range(id_list)
-    spikes = readFile(filename)
-    print "Read %d spikes from %s" % (len(spikes), filename)
-    return SpikeList(spikes, id_list, dt, t_start, t_stop, label=filename)
-
-
-def loadConductanceTraceList(filename, id_list, dt = None, t_start=None, t_stop=None):
-    if isinstance(id_list,int): # allows to just specify the number of neurons
-        id_list = range(id_list)
-    analog_signals = readFile(filename)
-    return ConductanceTraceList(analog_signals, id_list, dt, t_start, t_stop)
-
-
-def loadMembraneTraceList(filename, id_list, dt = None, t_start=None, t_stop=None):
-    if isinstance(id_list,int): # allows to just specify the number of neurons
-        id_list = range(id_list)
-    analog_signals = readFile(filename)
-    return MembraneTraceList(analog_signals, id_list, dt, t_start, t_stop)
-
-
-def loadCurrentTraceList(filename, id_list, dt = None, t_start=None, t_stop=None):
-    if isinstance(id_list,int): # allows to just specify the number of neurons
-        id_list = range(id_list)
-    analog_signals = readFile(filename)
-    return CurrentTraceList(analog_signals, id_list, dt, t_start, t_stop)
-
-
-class AnalogSignal(object):
-
     def __init__(self, signal, dt, t_start=None, t_stop=None):
 
@@ -1130 +1899 @@
         # elements with t <= t_stop
         if self.t_stop is not None:
-            self.signal = self.signal[:(self.t_stop-self.t_start/self.dt)]
-
-        if len(signal) > 0: # spike list may be empty
+            self.signal = self.signal[:(self.t_stop-self.t_start)/self.dt]
+
+        if len(self.signal) > 0: # spike list may be empty
             if self.t_start is None:
                 self.t_start = 0.
@@ -1141 +1910 @@
         #TODO return an exception if self.t_stop < self.t_start (when not empty)
         if self.t_start >= self.t_stop :
-            raise("Incompatible time interval for the creation of the AnalogSignal")
+            raise Exception("Incompatible time interval for the creation of the AnalogSignal")
+
+    def __getslice__(self, i, j):
+        """
+        Return a sublist of the signal vector of the AnalogSignal
+        """
+        return self.signal[i:j]
+
+    def duration(self):
+        """
+        Return the duration of the SpikeTrain
+        """
+        return self.t_stop - self.t_start
 
     def __str__(self):
@@ -1149 +1930 @@
         return len(self.signal)
 
+    def copy(self):
+        """
+        Return a copy of the AnalogSignal object
+        """
+        return AnalogSignal(self.signal, self.dt, self.t_start, self.t_stop)
+
     def time_axis(self):
+        """
+        Return the time axis of the AnalogSignal
+        """
         return numpy.arange(self.t_start, self.t_stop, self.dt)
 
+    def time_slice(self, t_start, t_stop):
+        """ 
+        Return a new AnalogSignal obtained by slicing between t_start and t_stop
+        
+        Inputs:
+            t_start - begining of the new SpikeTrain, in ms.
+            t_stop  - end of the new SpikeTrain, in ms.
+
+        """
+        signal = self.signal[t_start/self.dt:t_stop/self.dt]
+        return AnalogSignal(signal, self.dt, t_start, t_stop)
+        
+
 
 class AnalogSignalList(object):
-
-    def __init__(self, signals, id_list, dt=None, t_start=None, t_stop=None):
-        self.id_list = id_list
-        self.N = len(id_list)
+    """
+    AnalogSignalList(signals, id_list, dt=None, t_start=None, t_stop=None, dims=None)
+    
+    Return a AnalogSignalList object which will be a list of AnalogSignal objects.
+
+    Inputs:
+        signals - a list of the AnalogSignals objects
+        id_list - the list of the ids of all recorded cells (needed for silent cells)
+        dt      - if dt is specified, time values should be floats
+        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
+    
+    dt, t_start and t_stop are shared for all SpikeTrains object within the SpikeList
+    
+    See also
+        load_currentlist load_vmlist, load_conductancelist
+
+    """
+    def __init__(self, signals, id_list, dt=None, t_start=None, t_stop=None, dims=None):
         self.t_start = t_start
         self.t_stop  = t_stop
-        self.dt = dt
+        self.dt      = dt
+        self.dimensions     = dims
         self.analog_signals = {}
         for id in id_list:
@@ -1167 +1987 @@
             self.analog_signals[id].append(value)
 
+        if t_start is None:
+            self.t_start = 0
+
         # writing as a list of SpikeTrains
         for id,signal in self.analog_signals.items(): #TODO: handle missing fields
@@ -1172 +1995 @@
         if t_start is None or t_stop is None:
             self.__calc_startstop()
-            #self.spiketrains.append( SpikeTrain(spike_array[i_neuron], dt=dt ))
+
+    def id_list(self):
+        """ 
+        Return the list of all the cells ids contained in the
+        SpikeList object
+        """
+        return numpy.array(self.analog_signals.keys())
+
+    def copy(self):
+        """
+        Return a copy of the AnalogSignalList object
+        """
+        # Maybe not optimal, should be optimized
+        aslist = AnalogSignalList([], [], self.dt, self.t_start, self.t_stop, self.dimensions)
+        for id in self.id_list():
+            aslist.append(id, self.analog_signals[id])
+        return aslist
 
     def __calc_startstop(self):
@@ -1182 +2021 @@
         if len(self) > 0:
             if self.t_start is None:
-                self.t_start = 0.
-            #if self.t_stop is None:
-            #
+                start_times  = numpy.array([self.analog_signals[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.analog_signals.keys():
+                    self.analog_signals[id].t_start = self.t_start
+            if self.t_stop is None:
+                stop_times  = numpy.array([self.analog_signals[idx].t_stop for idx in self.id_list()])
+                self.t_stop = numpy.max(stop_times)
+                logging.debug("Warning, t_stop  is infered from the data : %f" %self.t_stop)
+                for id in self.analog_signals.keys():
+                    self.analog_signals[id].t_stop = self.t_stop
         else:
             raise Exception("No Analog Signals !")
+    
+    def __getdisplay__(self,display):
+        """
+        Return a pylab object with a plot() function to draw the plots.
+        
+        Inputs:
+            display - if True, a new figure is created. Otherwise, if display is a
+                      subplot object, this object is returned.
+        """
+        if display is False:
+            return None
+        elif display is True:
+            pylab.figure()
+            return pylab
+        else:
+            return display
+    
+    def __labels__(self, subplot, xlabel, ylabel):
+        """
+        Function to put some labels on a plot
+        
+        Inputs:
+            subplot - the targeted plot
+            xlabel  - a string for the x label
+            ylabel  - a string for the y label
+        """
+        if hasattr(subplot, 'xlabel'):
+            subplot.xlabel(xlabel)
+            subplot.ylabel(ylabel)
+        else:
+            subplot.set_xlabel(xlabel)
+            subplot.set_ylabel(ylabel)
 
     def __getitem__(self, i):
@@ -1199 +2078 @@
         return len(self.analog_signals)
 
+    def __sub_id_list(self, sub_list=None):
+        if sub_list == None:
+            return self.id_list()
+        if type(sub_list) == int:
+            return numpy.random.permutation(self.id_list())[0:sub_list]
+        if type(sub_list) == list:
+            return sub_list
+
     def append(self, id, signal):
+        """
+        Add an AnalogSignal object to the AnalogSignalList
+        
+        Inputs:
+            id     - the id of the new cell
+            signal - the AnalogSignal object representing the new cell
+        
+        The AnalogSignal object is sliced according to the t_start and t_stop times
+        of the AnalogSignallist object
+        
+        See also
+            __setitem__
+        """
         assert isinstance(signal, AnalogSignal), "An AnalogSignalList object can only contain AnalogSignal objects"
-        if id in self.id_list:
-            raise Exception("Id already present in AnalogSignalList.Use setitem instead()")
+        if id in self.id_list():
+            raise Exception("Id already present in AnalogSignalList. Use setitem instead()")
         else:
             self.analog_signals[id] = signal
-            self.id_list.append(id)
-            self.N += 1
         self.__calc_startstop()
 
     def time_axis(self):
+        """
+        Return the time axis of the AnalogSignalList object
+        """
         return numpy.arange(self.t_start,self.t_stop,self.dt)
 
-
-
-
-class MembraneTraceList(AnalogSignalList):
-
-    def plot(self, id_list, v_thresh=None):
-        pylab.figure()
-        pylab.hold(1)
-        pylab.ylabel("Analog values [unit]", size="x-large")
-        pylab.xlabel("Time (ms)", size="x-large")
+    def id_slice(self, id_list):
+        """
+        Return a new AnalogSignalList obtained by selecting particular ids
+        
+        Inputs:
+            id_list - Can be an integer (and then N random cells will be selected)
+                      or a sublist of the current ids
+        
+        The new AnalogSignalList inherits the time parameters (t_start, t_stop, dt)
+        
+        See also
+            time_slice
+        """
+        new_AnalogSignalList = AnalogSignalList([], [], self.dt, self.t_start, self.t_stop, self.dimensions)
+        id_list = self.__sub_id_list(id_list)
         for id in id_list:
-            to_be_plot = self.analog_signals[id].signal
-            if v_thresh is not None:
-                to_be_plot = pylab.where(to_be_plot>=v_thresh-0.05,0.0,to_be_plot)
-            pylab.plot(to_be_plot)
-        [x,y] = pylab.xticks()
-        pylab.xticks(x,numpy.array(self.t_start+x*self.dt,'string').tolist())
-
-
-class CurrentTraceList(AnalogSignalList):
-
-    def plot(self, id_list):
-        pylab.figure()
-        pylab.hold(1)
-        pylab.ylabel("Current (nA)", size="x-large")
-        pylab.xlabel("Time (ms)", size="x-large")
-        for id in id_list:
-            pylab.plot(self.analog_signals[id].signal)
-        [x,y] = pylab.xticks()
-        pylab.xticks(x,numpy.array(self.t_start+x*self.dt,'string').tolist())
-
-
-class ConductanceTraceList(AnalogSignalList):
-
-    def plot(self, id_list):
-        pylab.figure()
-        pylab.hold(1)
-        pylab.ylabel("Conductance (nS)", size="x-large")
-        pylab.xlabel("Time (ms)", size="x-large")
-        for id in id_list:
-            pylab.plot(self.analog_signals[id].signal)
-        [x,y] = pylab.xticks()
-        pylab.xticks(x,numpy.array(self.t_start+x*self.dt,'string').tolist())
-
-
-
-
-
-def ccf(x, y, axis=None):
-    """Computes the cross-correlation function of two series `x` and `y`.
-Note that the computations are performed on anomalies (deviations from
-average).
-Returns the values of the cross-correlation at different lags.
-Lags are given as [0,1,2,...,n,n-1,n-2,...,-2,-1] (not any more)
-
-:Parameters:
-    `x` : 1D MaskedArray
-        Time series.
-    `y` : 1D MaskedArray
-        Time series.
-    `axis` : integer *[None]*
-        Axis along which to compute (0 for rows, 1 for cols).
-        If `None`, the array is flattened first.
+            try:
+                AnalogSignalList.append(id, self.analog_signals[id])
+            except Exception:
+                print "id %d is not in the source AnalogSignalList" %id
+        return new_AnalogSignalList
+
+    def time_slice(self, t_start, t_stop):
+        """
+        Return a new AnalogSignalList obtained by slicing between t_start and t_stop
+        
+        Inputs:
+            t_start - begining of the new AnalogSignalList, in ms.
+            t_stop  - end of the new AnalogSignalList, in ms.
+        
+        See also
+            id_slice
+        """
+        new_AnalogSignalList = AnalogSignalList([], [], self.dt, t_start, t_stop, self.dimensions)
+        for id in self.id_list():
+            new_AnalogSignalList.append(id, self.analog_signals[id].time_slice(t_start, t_stop))
+        new_AnalogSignalList.__calc_startstop()
+        return new_AnalogSignalList
+
+    def save(self, filename, method="text"):
+        """
+        Save the AnalogSignal in a text or binary file
+        
+        Inputs:
+            filename - name of the output file
+            method   - "text"   -> a classical human readible text file
+                       "pickle" -> binary backup, much more faster to load/save, 
+                                   but no longer human readible
+                       "hdf5"   -> a compress an optimized structure, not
+                                   implemented yet
+        """
+        as_loader = io.AnalogSignalIO(filename)
+        if method == "pickle":
+            as_loader.save_analogsignal_pickle(self, filename)
+        if method == "text":
+            as_loader.save_analogsignal_txt(self, filename)
+        
+    
+
+class VmList(AnalogSignalList):
+
+    def plot(self, id_list=None, v_thresh=None, display=True, kwargs={}):
+        """
+        Plot all cells in the AnalogSignalList defined by id_list
+        
+        Inputs:
+            id_list - can be a integer (and then N cells are randomly selected) or a 
+                      list of ids. If None, we use all the ids of the SpikeList
+            v_thresh- For graphical purpose, plot a spike when Vm > V_thresh. If None, 
+                      just plot the raw Vm
+            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)
+            >>> aslist.plot(5, v_thresh = -50, display=z, kwargs={'color':'r'})
+        """
+        subplot   = self.__getdisplay__(display)
+        id_list   = self._AnalogSignalList__sub_id_list(id_list)
+        time_axis = self.time_axis()  
+        if not subplot or not ENABLE_PLOTS:
+            print MATPLOTLIB_ERROR
+        else:
+            xlabel = "Time (ms)"
+            ylabel = "Membrane Potential (mV)"
+            self.__labels__(subplot, xlabel, ylabel)
+            for id in id_list:
+                to_be_plot = self.analog_signals[id].signal
+                if v_thresh is not None:
+                    to_be_plot = pylab.where(to_be_plot>=v_thresh-0.05,0.0,to_be_plot)
+                subplot.plot(time_axis, to_be_plot, **kwargs)
+                subplot.hold(1)
+            pylab.draw()
+
+
+class CurrentList(AnalogSignalList):
+
+    def plot(self, id_list=None, v_thresh=None, display=True, kwargs={}):
+        """
+        Plot all cells in the AnalogSignalList defined by id_list
+        
+        Inputs:
+            id_list - can be a integer (and then N cells are randomly selected) or a 
+                      list of ids. If None, we use all the ids of the SpikeList
+            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)
+            >>> aslist.plot(5, display=z, kwargs={'color':'r'})
+        """
+        subplot   = self.__getdisplay__(display)
+        id_list   = self._AnalogSignalList__sub_id_list(id_list)
+        time_axis = self.time_axis()  
+        if not subplot or not ENABLE_PLOTS:
+            print MATPLOTLIB_ERROR
+        else:
+            xlabel = "Time (ms)"
+            ylabel = "Current (nA)"
+            self.__labels__(subplot, xlabel, ylabel)
+            for id in id_list:
+                subplot.plot(time_axis, self.analog_signals[id].signal, **kwargs)
+                subplot.hold(1)
+            pylab.draw()
+
+class ConductanceList(AnalogSignalList):
+
+    def plot(self, id_list=None, v_thresh=None, display=True, kwargs={}):
+        """
+        Plot all cells in the AnalogSignalList defined by id_list
+        
+        Inputs:
+            id_list - can be a integer (and then N cells are randomly selected) or a 
+                      list of ids. If None, we use all the ids of the SpikeList
+            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)
+            >>> aslist.plot(5, display=z, kwargs={'color':'r'})
+        """
+        subplot   = self.__getdisplay__(display)
+        id_list   = self._AnalogSignalList__sub_id_list(id_list)
+        time_axis = self.time_axis()  
+        if not subplot or not ENABLE_PLOTS:
+            print MATPLOTLIB_ERROR
+        else:
+            xlabel = "Time (ms)"
+            ylabel = "Conductance (nS)"
+            self.__labels__(subplot, xlabel, ylabel)
+            for id in id_list:
+                subplot.plot(time_axis, self.analog_signals[id].signal, **kwargs)
+                subplot.hold(1)
+            pylab.draw()
+
+
+
+#############################################################
+## Object Loaders. Functions used to create NeuroTools
+## objects from data generated by pyNN (the most simple form
+## supported right now)
+#############################################################
+
+def load_spikelist(filename, id_list=None, dt = None, t_start=None, t_stop=None, dims=None):
     """
-    assert(x.ndim == y.ndim, "Inconsistent shape !")
-#    assert(x.shape == y.shape, "Inconsistent shape !")
-    if axis is None:
-        if x.ndim > 1:
-            x = x.ravel()
-            y = y.ravel()
-        npad = x.size + y.size
-        xanom = (x - x.mean(axis=None))
-        yanom = (y - y.mean(axis=None))
-        Fx = numpy.fft.fft(xanom, npad, )
-        Fy = numpy.fft.fft(yanom, npad, )
-        iFxy = numpy.fft.ifft(Fx.conj()*Fy).real
-        varxy = numpy.sqrt(numpy.inner(xanom,xanom) * numpy.inner(yanom,yanom))
-    else:
-        npad = x.shape[axis] + y.shape[axis]
-        if axis == 1:
-            if x.shape[0] != y.shape[0]:
-                raise ValueError, "Arrays should have the same length!"
-            xanom = (x - x.mean(axis=1)[:,None])
-            yanom = (y - y.mean(axis=1)[:,None])
-            varxy = numpy.sqrt((xanom*xanom).sum(1) * (yanom*yanom).sum(1))[:,None]
-        else:
-            if x.shape[1] != y.shape[1]:
-                raise ValueError, "Arrays should have the same width!"
-            xanom = (x - x.mean(axis=0))
-            yanom = (y - y.mean(axis=0))
-            varxy = numpy.sqrt((xanom*xanom).sum(0) * (yanom*yanom).sum(0))
-        Fx = numpy.fft.fft(xanom, npad, axis=axis)
-        Fy = numpy.fft.fft(yanom, npad, axis=axis)
-        iFxy = numpy.fft.ifft(Fx.conj()*Fy,n=npad,axis=axis).real
-    # We juste turn the lags into correct positions:
-    iFxy = numpy.concatenate((iFxy[len(iFxy)/2:len(iFxy)],iFxy[0:len(iFxy)/2]))
-    return iFxy/varxy
-
-
+    Returns a SpikeList object from a file. If the file has been generated by PyNN, 
+    a header should be found with following parameters:
+     ---> dims, dt, id of the first cell, id of the last cell. 
+    They must be specified otherwise.  Then the classical PyNN format for text file is:
+     ---> one line per event:  absolute time in ms, GID
+    
+    Inputs:
+        filename - the name of the spike file
+        id_list  - the list of the recorded ids. Can be an int (meaning cells in 
+                   the range (0,..,N)), or a list. 
+        dims     - if the cells were aranged on a 2/3D grid, a tuple with the dimensions
+        dt       - the discretization step, in ms
+        t_start  - begining of the simulation, in ms.
+        t_stop   - end of the simulation, in ms
+
+    If dims, dt, t_start, t_stop or id_list are None, they will be infered from either 
+    the data or from the header.
+    All times are in milliseconds. 
+    The format of the file (text, pickle or hdf5) will be inferred automatically
+    """
+    spike_loader = io.SpikeListIO(filename)
+    return spike_loader.get_spikelist(id_list, dt, t_start, t_stop, dims)
+
+
+
+def load_conductancelist(filename, id_list=None, dt = None, t_start=None, t_stop=None, dims=None):
+    """
+    Returns a ConductanceList object from a file. If the file has been generated by PyNN, 
+    a header should be found with following parameters:
+     ---> dims, dt, id of the first cell, id of the last cell. 
+    They must be specified otherwise.  Then the classical PyNN format for text file is:
+     ---> one line per event:  data value, GID
+    
+    Inputs:
+        filename - the name of the spike file
+        id_list  - the list of the recorded ids. Can be an int (meaning cells in 
+                   the range (0,..,N)), or a list. 
+        dims     - if the cells were aranged on a 2/3D grid, a tuple with the dimensions
+        dt       - the discretization step, in ms
+        t_start  - begining of the simulation, in ms.
+        t_stop   - end of the simulation, in ms
+
+    If dims, dt, t_start, t_stop or id_list are None, they will be infered from either 
+    the data or from the header.
+    All times are in milliseconds. 
+    The format of the file (text, pickle or hdf5) will be inferred automatically
+    """
+    analogsignal_loader = io.AnalogSignalIO(filename)
+    return analogsignal_loader.get_conductancelist(id_list, dt, t_start, t_stop, dims)
+
+
+def load_vmlist(filename, id_list=None, dt=None, t_start=None, t_stop=None, dims=None):
+    """
+    Returns a VmList object from a file. If the file has been generated by PyNN, 
+    a header should be found with following parameters:
+     ---> dims, dt, id of the first cell, id of the last cell. 
+    They must be specified otherwise.  Then the classical PyNN format for text file is:
+     ---> one line per event:  data value, GID
+    
+    Inputs:
+        filename - the name of the spike file
+        id_list  - the list of the recorded ids. Can be an int (meaning cells in 
+                   the range (0,..,N)), or a list. 
+        dims     - if the cells were aranged on a 2/3D grid, a tuple with the dimensions
+        dt       - the discretization step, in ms
+        t_start  - begining of the simulation, in ms.
+        t_stop   - end of the simulation, in ms
+
+    If dims, dt, t_start, t_stop or id_list are None, they will be infered from either 
+    the data or from the header.
+    All times are in milliseconds. 
+    The format of the file (text, pickle or hdf5) will be inferred automatically
+    """
+    analogsignal_loader = io.AnalogSignalIO(filename)
+    return analogsignal_loader.get_vmlist(id_list, dt, t_start, t_stop, dims)
+
+
+
+def load_currentlist(filename, id_list=None, dt = None, t_start=None, t_stop=None, dims=None):
+    """
+    Returns a CurrentList object from a file. If the file has been generated by PyNN, 
+    a header should be found with following parameters:
+     ---> dims, dt, id of the first cell, id of the last cell. 
+    They must be specified otherwise.  Then the classical PyNN format for text file is:
+     ---> one line per event:  data value, GID
+    
+    Inputs:
+        filename - the name of the spike file
+        id_list  - the list of the recorded ids. Can be an int (meaning cells in 
+                   the range (0,..,N)), or a list. 
+        dims     - if the cells were aranged on a 2/3D grid, a tuple with the dimensions
+        dt       - the discretization step, in ms
+        t_start  - begining of the simulation, in ms.
+        t_stop   - end of the simulation, in ms
+
+    If dims, dt, t_start, t_stop or id_list are None, they will be infered from either 
+    the data or from the header.
+    All times are in milliseconds. 
+    The format of the file (text, pickle or hdf5) will be inferred automatically
+    """
+    analogsignal_loader = io.AnalogSignalIO(filename)
+    return analogsignal_loader.get_currentlist(id_list, dt, t_start, t_stop, dims)