Changeset 207

Timestamp:

09/22/08 17:03:25 (4 months ago)

Author:

pierre

Message:

Implement some distances for my personal work, the Victor Purpura one and the Kreuz one, for the SpikeLists objects

Files:

• branches/cleanup/src/io.py (modified) (6 diffs)
• branches/cleanup/src/signals.py (modified) (16 diffs)

branches/cleanup/src/io.py

@@ -14 +14 @@
-            raise Exception("The file %s does not exists !" %filename)
-    
+    
+    ### Function to load/save data with the pickle package. Much more faster than a simple
+    ### plain text file. Useful if you don't care about reading the spike files
-    def pickle_load(self):
-        return cPickle.load(file(self.filename,"r"))
@@ -20 +23 @@
-        return cPickle.dump(object, file(filename, "w"))
-    
+    
+    ### Function to load/save data in a text format.
-    def txt_load(self, sepchar = "\t", skipchar = "#"):
-        myfile = open(self.filename, "r", DEFAULT_BUFFER_SIZE)
@@ -79 +84 @@
-    
-    def get_spikelist(self, id_list=None, dt=None, t_start=None, t_stop=None, dims=None, label=None):
-        
-        # First, we try to see if the data are stored in a pickled way (MUCH FASTER)
-        try :
@@ -90 +94 @@
-                logging.debug("dt is infered from the file header")
-                dt = self.header['dt']
-
-
-
-
-
-
+
+
+
+
+
+
+
-            if dims is None and 'dimensions' in self.header:
-                dims = self.header['dimensions']
@@ -121 +126 @@
-
-
-TxtAnalogSignalLoader
+AnalogSignalIO
-    """
-    Convenient object to load txt file and process them in order
@@ -144 +149 @@
-        return 0
-    
-SpikeList
+AnalogSignal
-        header = self.parse_header()
-        

branches/cleanup/src/signals.py

@@ -166 +166 @@
-        
-        Inputs:
-    
+
-        
-        Examples:
@@ -394 +394 @@
-            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))
-
-    ####################################################################
@@ -746 +807 @@
-        
-    
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
@@ -876 +945 @@
-    def mean_rate_std(self, t_start=None, t_stop=None):
-        """
-Mean firing rate averaged
+firing rates
-        between t_start and t_stop
-
@@ -897 +966 @@
-    def mean_rates(self, t_start=None, t_stop=None):
-        """ 
-
+firing 
-
-        Inputs:
@@ -1122 +1191 @@
-
-
-
+spklist=None, 
-        """
-        Function to generate an array of cross correlations computed
@@ -1129 +1198 @@
-        Inputs:
-            pairs    - Could be an int, and then N random pairs of cells will be selected, 
-
+it 
-            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
@@ -1141 +1212 @@
-        if type(pairs) == int:
-            spk1 = self.id_slice(pairs)
-elf
+pklist
-            N    = pairs
-        else:
@@ -1147 +1218 @@
-            N      = len(pairs[:,0])
-            spk1   = self.id_slice(pairs[:,0])
-elf
+pklist
-        hist_1 = spk1.spike_histogram(time_bin)
-        hist_2 = spk2.spike_histogram(time_bin)
@@ -1180 +1251 @@
-            pylab.draw()
-
+    
+    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 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 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):
+        """
+        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(numpy.sum(firing_rate, axis=0)/self.N)
-
-    def mean_rate_covariance(self, spikelist, time_bin):
-        """
-
-
-
+
+
+
+
+
+
+
+
+
-        """
-        if not isinstance(spikelist, SpikeList):
@@ -1193 +1345 @@
-        frate_1 = self.firing_rate(time_bin)
-        frate_1 = numpy.sum(frate_1, axis=0)/self.N
-SpkL
+spikel
-        frate_2 = numpy.sum(frate_2, axis=0)/spikelist.N
-      
-    
+
+
-        return cov
-
@@ -1433 +1585 @@
-
-
+
+
-def load_spikelist(filename, id_list=None, dt = None, t_start=None, t_stop=None, dims=None, label=None):
-    """
@@ -1451 +1605 @@
-
-    If dims, dt, t_start, t_stop or id_list are None, they will be infered from either the data or from the header.
-
-
-
-
-
-
+
-    """
-
-
+
-    return spike_loader.get_spikelist(id_list, dt, t_start, t_stop, dims, label)
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
@@ -1498 +1659 @@
-            raise("Incompatible time interval for the creation of the AnalogSignal")
-
+    def __getslice__(self, i, j):
+        """
+        Return a sublist of the spike_times vector of the SpikeTrain
+        """
+        return self.signal[i:j]
+
+    def duration(self):
+        """
+        Return the duration of the SpikeTrain
+        """
+        return self.t_stop - self.t_start
+
-    def __str__(self):
-        return str(self.signal)
@@ -1509 +1682 @@
-
-class AnalogSignalList(object):
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    def __init__(self, signals, id_list, dt=None, t_start=None, t_stop=None):
-        self.id_list = id_list
-
+      
-        self.t_start = t_start
-        self.t_stop  = t_stop
-
+     
-        self.analog_signals = {}
-        for id in id_list:
@@ -1542 +1739 @@
-        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):