Changeset 215

Timestamp:

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

Author:

apdavison

Message:

Merged cleanup branch changes r198:214 into the trunk.

Files:

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

trunk/INSTALL

@@ -1 @@
-
+
+
-
-First download the package files:
-svn co https://neuralensemble.kip.uni-heidelberg.de/svn/NeuroTools/trunk NeuroTools
-cd NeuroTools
-
-
+
+
-
-
+
-
-
+
-
-
+
+
+
-
-Then you need to add the location: 
-
-4
+5
-
-to your PYTHON_PAH or (in python) your sys.path
-
-4
+5
-(obviously) replaced by your python version.
-
-
-
-
+
+
+
+
+
-and voila!
-

trunk/README

@@ -1 @@
-
+

trunk/examples/single_neuron/simple_single_neuron.py

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

trunk/src/__init__.py

@@ -1 @@
-pikes', 'signals', 'stgen', 'utilities
+ignals', '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 @@
-
-
-
+
-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 @@
-"""
-
-
-
+
+
+
+
-try :
-    import pylab
-
-
-
-
+
+
+
+
+
+
+
+
+
-
-class SpikeTrain(object):
-    """
-
+
+
-
-    Event times are given in a list (sparse representation) in milliseconds.
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    """
-    # 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):
-        """
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-
-        """
-
-        if dt is not None:
-<=
+ <= 
-                raise ValueError("dt must be greater than zero")
-            #self.spike_times *= dt
-
-
-
-
-
+
+
+
-        self.spike_times = numpy.array(spike_times, 'float')
-
@@ -69 +71 @@
-        # the spikes with t >= t_start
-        if self.t_start is not None:
-
-
-
-
+
-
-        # 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]
-
-    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...
+
-        # 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:
-
+ 
-                self.t_start = 0
-            if self.t_stop is None:
@@ -111 +110 @@
-
-        if self.t_start >= self.t_stop :
-for the creation of the SpikeTrain. t_start=%s, t_stop=
+: t_start = %s, 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):
-        """
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
-        """
-        spike_times = self.spike_times.copy()
-
-
-
+
+
-            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):
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
-    def mean_rate(self, t_start=None, t_stop=None):
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            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.
-
-
+
+
+
+
-        """
-        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):
-
+
+
+
+
+
+
+
+
-        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 !")
-
-
+=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)
-
-color='b'
+display=True, kwargs={}
-        """
-        Generate a raster plot with the SpikeTrain in a subwindow of interest,
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        return SpikeTrain(spikes, self.dt, t_start, t_stop)
-
-
+
+
-        """
-        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 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):
-        """
-
-
+
+
+
+
+
-        if self.t_start != 0:
-            self.spike_times -= self.t_start
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    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.
-        """
-
+      
-        freq_spect = analysis.simple_frequency_spectrum(hist)
-
-
+  
+ 
-        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
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-
-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 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   ##
-    #######################################################################
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-