Changeset 408

Timestamp:

07/16/09 11:50:11 (3 years ago)

Author:

pierre

Message:

Start to change all the NeuroTools object to use a more general concept of the interval instead of t_start, t_stop. Will help a lot the analyze of the data when objects have several trials or stimulus presented during a single run.

Location:

branches/interval/src/signals

Files:

• intervals.py (modified) (4 diffs)
• spikes.py (modified) (45 diffs)

branches/interval/src/signals/intervals.py

@@ -22 +22 @@
     """
     
-    def __init__(self, start_times, end_times) :
+    def __init__(self, sub_intervals):
         """
         Constructor of the Interval object.
 
         """
-        if HAVE_INTERVAL:
-            self.start_times = start_times
-            self.end_times   = end_times
-            # write the intervals to an interval object (pyinterval)
-            test = isinstance(start_times, int) or isinstance(start_times, float)
-            if test:
-                start_times = [start_times]
-            test = isinstance(end_times, int) or isinstance(end_times, float)
-            if test:
-                end_times = [end_times]
-            if len(start_times) != len(end_times) :
-                raise Exception("There sould be an equal number of starts and stops")
-            self.interval_data = interval(*numpy.transpose(numpy.array([start_times,end_times])))
-        else:
-            test = isinstance(start_times, int) or isinstance(start_times, float)
-            assert test, "Interval package not present, start_times should be a number !"
-            test = isinstance(end_times, int) or isinstance(end_times, float)
-            assert test, "Interval package not present, end_times should be a number !"
-            self.start_times = start_times
-            self.end_times   = end_times
+        if not type(sub_intervals) is list:
+            sub_intervals = list((sub_intervals,))
+        
+        self._check_interval(sub_intervals)
+        self.sub_intervals = sub_intervals
+        self.start_times   = numpy.array(self.sub_intervals)[:,0]
+        self.end_times     = numpy.array(self.sub_intervals)[:,1]
+
+        if HAVE_INTERVAL:
+            self.interval_data = interval(*numpy.transpose(numpy.array([self.start_times, self.end_times])))
+
+    def _check_interval(self, sub_intervals):
+        if not HAVE_INTERVAL:
+            assert len(sub_intervals) == 1, "Interval package not present, interval must be (t_start, t_stop) !"
+        for item in sub_intervals:
+            assert (len(item) == 2), "Intervals must be a list of tuple (t_start, t_stop) !"
+            assert item[1] > item[0], "Intervals must have tuple with t_start < t_stop !"
 
     def intersect(self, itv) :
-        self.interval_data = self.interval_data & itv.interval_data
+        result = Interval(self.sub_intervals)
+        result.interval_data = result.interval_data & itv.interval_data
+        return result
 
     def union(self, itv) :
-        self.interval_data = self.interval_data | itv.interval_data
-
+        if HAVE_INTERVAL:
+            result = Interval(self.sub_intervals)
+            result.interval_data = result.interval_data & itv.interval_data
+            result.interval_data = result.interval_data | itv.interval_data
+        else:
+            result = Interval(self.sub_intervals)
+            result.start_times   = min(result.start_times, spiketrain.start_times)
+            result.end_times     = max(result.end_times, spiketrain.end_times)
+        return result
+    
     def __str__(self):
         return str(self.interval_data)
@@ -60 +67 @@
         return shape(self.interval_data)[0]
     
+    def __iter__(self):
+        return iter(self.sub_intervals)
+    
     def __getslice__(self, i, j):
         """
@@ -90 +100 @@
         """
         return Interval(self.start_times, self.end_times, self.t_start, self.t_stop)
+
+
+    def offset_start(self, shift, from_stop=False) :
+        """
+        Shift all stop times of the interval by shift (ms). If from_start is
+        true, use start times as the reference to compute the new stop times  
+        """
+        interval_array = numpy.array(self.interval_data)
+
+        if from_stop == False :
+            if (interval_array[:,0] + shift > interval_array[:,1]).any() :
+                raise Exception("the shift is too big to preserve start/stop temporal order")
+            interval_array[:,0] += shift
+        else :
+            if shift > 0 :
+                raise Exception("the shift should be negative to preserve start/stop temporal order")
+            interval_array[:,0] = interval_array[:,1] + shift
+
+        self.interval_data = interval(*interval_array)
+
+
+    def offset_stop(self, shift, from_start=False) :
+        """
+        Shift all stop times of the interval by shift (ms). If from_start is
+        true, use start times as the reference to compute the new stop times  
+        """
+        interval_array = numpy.array(self.interval_data)
+
+        if from_start == False :
+            if (interval_array[:,1] + shift < interval_array[:,0]).any() :
+                raise Exception("the shift is too big to preserve start/stop temporal order")
+            interval_array[:,1] += shift
+        else :
+            if shift < 0 :
+                raise Exception("the shift should be positive to preserve start/stop temporal order")
+            interval_array[:,1] = interval_array[:,0] + shift
+
+        self.interval_data = interval(*interval_array)
+
+    def offset_full(self, shift) :
+        """
+        Shift the whole interval by 'shift' ms
+        """
+        self.interval_data += shift
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
     def offset(self, start=None, end=None) :
@@ -144 +213 @@
         return numpy.extract(spikes_selector, times)
 
+    def is_equal(self, itv) :
+        return self.sub_intervals == itv.sub_intervals
 
 

branches/interval/src/signals/spikes.py

@@ -49 +49 @@
 
 newnum = check_numpy_version()
-
+INTERVAL_WARNING = "WARNING: t_start, t_stop are deprecated ! Use a syntax like interval=(t_start, t_stop) instead"
 
 class SpikeTrain(object):
@@ -76 +76 @@
     ## Constructor and key methods to manipulate the SpikeTrain objects  ##
     #######################################################################
-    def __init__(self, spike_times, t_start=None, t_stop=None):
+    def __init__(self, spike_times, interval=None):
         """
         Constructor of the SpikeTrain object
@@ -83 +83 @@
             SpikeTrain
         """
-
-        self.t_start     = t_start
-        self.t_stop      = t_stop
-        self.spike_times = numpy.array(spike_times, numpy.float32)
-
-        # If t_start is not None, we resize the spike_train keeping only
-        # the spikes with t >= t_start
-        if self.t_start is not None:
-            self.spike_times = numpy.extract((self.spike_times >= self.t_start), self.spike_times)
-
-        # If t_stop is not None, we resize the spike_train keeping only
-        # the spikes with t <= t_stop
-        if self.t_stop is not None:
-            self.spike_times = numpy.extract((self.spike_times <= self.t_stop), self.spike_times)
+        if interval is not None:
+            self.interval = interval
+        else:
+            try:
+                t_start = min(spike_times)
+            except Exception:
+                print "Error in guessing t_start (first spike), spikes may be empty !"
+                t_start = 0
+            try:
+                t_stop  = max(spike_times)
+            except Exception:
+                print "Error in guessing t_stop (last spike), spikes may be empty !"
+                t_stop = 0.1
+            self.interval = Interval((t_start, t_stop))
+
+        self.spike_times = self.interval.slice_times(numpy.array(spike_times, numpy.float32))
 
         # We sort the spike_times. May be slower, but is necessary by the way for quite a 
@@ -107 +109 @@
         # 1 element  : t_start = time, t_stop = time + 0.1
         # several    : t_start = min(time), t_stop = max(time)
-        
-        size = len(self.spike_times)
-        if size == 0:
-            if self.t_start is None: 
-                self.t_start = 0
-            if self.t_stop is None:
-                self.t_stop  = 0.1
-        elif size == 1: # spike list may be empty
-            if self.t_start is None:
-                self.t_start = self.spike_times[0]
-            if self.t_stop is None:
-                self.t_stop = self.spike_times[0] + 0.1
-        elif size > 1:
-            if self.t_start is None:
-                self.t_start = numpy.min(self.spike_times)
-            if numpy.any(self.spike_times < self.t_start):
-                raise ValueError("Spike times must not be less than t_start")
-            if self.t_stop is None:
-                self.t_stop = numpy.max(self.spike_times)
-            if numpy.any(self.spike_times > self.t_stop):
-                raise ValueError("Spike times must not be greater than t_stop")
-
-        if self.t_start >= self.t_stop :
-            raise Exception("Incompatible time interval : t_start = %s, t_stop = %s" % (self.t_start, self.t_stop))
-        if self.t_start < 0:
-            raise ValueError("t_start must not be negative")
-        if numpy.any(self.spike_times < 0):
-            raise ValueError("Spike times must not be negative")
+
+
 
     def __str__(self):
@@ -152 +128 @@
         Return the time parameters of the SpikeTrain (t_start, t_stop)
         """
-        return (self.t_start, self.t_stop)
+        return (self.interval.t_start(), self.interval.t_stop())
     
     def is_equal(self, spktrain):
@@ -165 +141 @@
             time_parameters()
         """
-        test = (self.time_parameters() == spktrain.time_parameters())
+        test = (self.interval.is_equal(spktrain.interval))
         return numpy.all(self.spike_times == spktrain.spike_times) and test
     
@@ -172 +148 @@
         Return a copy of the SpikeTrain object
         """
-        return SpikeTrain(self.spike_times, self.t_start, self.t_stop)
+        return SpikeTrain(self.spike_times, interval=self.interval)
 
 
@@ -179 +155 @@
         Return the duration of the SpikeTrain
         """
-        return self.t_stop - self.t_start
+        return self.interval.total_duration()
     
     
@@ -202 +178 @@
         self.spike_times = numpy.insert(self.spike_times, self.spike_times.searchsorted(spiketrain.spike_times), \
                                         spiketrain.spike_times)
-        self.t_start     = min(self.t_start, spiketrain.t_start)
-        self.t_stop      = max(self.t_stop, spiketrain.t_stop)
+        self.interval = self.interval.union(spiketrain.interval)
 
     def format(self, relative=False, quantized=False):
@@ -244 +219 @@
         """
         
-        return SpikeTrain(self.spike_times+jitter*(numpy.random.normal(loc=0.0,scale=1.0,size=self.spike_times.shape[0])),t_start=self.t_start,t_stop=self.t_stop)
+        spk = self.spike_times+jitter*(numpy.random.normal(loc=0.0,scale=1.0,size=self.spike_times.shape[0]))
+        return SpikeTrain(spk, interval=self.interval)
         
 
@@ -266 +242 @@
         return numpy.diff(self.spike_times)
 
-    def mean_rate(self, t_start=None, t_stop=None):
+    def mean_rate(self, interval=None):
         """ 
         Returns the mean firing rate between t_start and t_stop, in Hz
         
         Inputs:
-            t_start - in ms. If not defined, the one of the SpikeTrain object is used
-            t_stop  - in ms. If not defined, the one of the SpikeTrain object is used
+            interval - An interval object
         
         Examples:
@@ -278 +253 @@
                 34.2
         """
-        if (t_start == None) & (t_stop == None):
-            t_start = self.t_start
-            t_stop  = self.t_stop
-            idx     = self.spike_times
-        else:
-            if t_start == None: 
-                t_start = self.t_start
-            else:
-                t_start = max(self.t_start, t_start)
-            if t_stop == None: 
-                t_stop=self.t_stop
-            else:
-                t_stop = min(self.t_stop, t_stop)
-            idx = numpy.where((self.spike_times >= t_start) & (self.spike_times <= t_stop))[0]
-        return 1000.*len(idx)/(t_stop-t_start)
+        if interval:
+            new_interval = interval.intersect(self.interval)
+            idx          = new_interval.slice_times(self.spike_times)
+        else:
+            idx          = self.spike_times
+            new_interval = self.interval
+        return 1000.*len(idx)/new_interval.total_duration()
 
     def cv_isi(self):
@@ -396 +363 @@
         """
         if newnum:
-            axis = numpy.arange(self.t_start, self.t_stop+time_bin, time_bin)
-        else:
-            axis = numpy.arange(self.t_start, self.t_stop, time_bin)
+            axis = numpy.arange(self.interval.t_start(), self.interval.t_stop()+time_bin, time_bin)
+        else:
+            axis = numpy.arange(self.interval.t_start(), self.interval.t_stop(), time_bin)
         return axis
 
-    def raster_plot(self, t_start=None, t_stop=None, interval=None, display=True, kwargs={}):
+    def raster_plot(self, interval=None, display=True, kwargs={}):
         """
         Generate a raster plot with the SpikeTrain in a subwindow of interest,
@@ -407 +374 @@
         
         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
+            interval - An interval object
             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 
@@ -420 +386 @@
             SpikeList.raster_plot
         """
-        if t_start is None: t_start = self.t_start
-        if t_stop is None:  t_stop = self.t_stop
         
         if interval is None:
-            interval = Interval(t_start, t_stop)
+            interval = self.interval
         
         spikes  = interval.slice_times(self.spike_times)
@@ -454 +418 @@
                 120.,  130.,  140.]
         """
-        self.t_start += offset
-        self.t_stop  += offset
+        self.interval.offset_full(offset)
         self.spike_times += offset
 
@@ -474 +437 @@
                 100
         """
-        spikes = numpy.extract((self.spike_times >= t_start) & (self.spike_times <= t_stop), self.spike_times)
-        return SpikeTrain(spikes, t_start, t_stop)
+        interval = Interval((t_start, t_stop)) 
+        return SpikeTrain.interval_slice(interval)
 
     def interval_slice(self, interval):
@@ -492 +455 @@
                 100
         """
-        times           = interval.slice_times(self.spike_times)
-        t_start, t_stop = interval.time_parameters()
-        return SpikeTrain(times, t_start, t_stop)
+        times = interval.slice_times(self.spike_times)
+        return SpikeTrain(times, interval=interval)
         
 
@@ -534 +496 @@
         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
+        padding = -self.interval.t_start()
+        if padding != 0:
+            self.spike_times += padding
+            self.interval.offset_full(padding)
 
     def distance_victorpurpura(self, spktrain, cost=0.5):
@@ -586 +548 @@
             VictorPurpuraDistance
         """
-        N              = (self.t_stop-self.t_start)/dt
+        
+        N              = (self.interval.t_stop()-self.interval.t_start())/dt
         vec_1          = numpy.zeros(N, numpy.float32)
         vec_2          = numpy.zeros(N, numpy.float32)
@@ -730 +693 @@
         hist       = self.time_histogram(time_bin, normalized=False)
         freq_spect = analysis.simple_frequency_spectrum(hist)
-        freq_bin   = 1000.0/self.duration() # Hz
+        freq_bin   = 1000.0/self.interval.total_duration() # Hz
         freq_axis  = numpy.arange(len(freq_spect)) * freq_bin    
         return freq_spect, freq_axis    
@@ -746 +709 @@
         freq_spect = self.frequency_spectrum(time_bin)[0]
         F0 = freq_spect[0]
-        freq_bin = 1000.0/self.duration()
+        freq_bin = 1000.0/self.interval.total_duration()
         try:
             F1 = freq_spect[int(round(f_stim/freq_bin))]
@@ -783 +746 @@
     ## Constructor and key methods to manipulate the SpikeList objects   ##
     #######################################################################
-    def __init__(self, spikes, id_list, t_start=None, t_stop=None, dims=None):
+    def __init__(self, spikes, id_list, interval=None, dims=None):
         """
         Constructor of the SpikeList object
@@ -790 +753 @@
             SpikeList, load_spikelist
         """
-        self.t_start     = t_start
-        self.t_stop      = t_stop
+        #if t_start or t_stop:
+            
+        self.interval    = interval
         self.dimensions  = dims
         self.spiketrains = {}
@@ -808 +772 @@
             spikes = spikes[idx]
             logging.debug("sorted spikes[:10,:] = %s" % str(spikes[:10,:]))
-        
             break_points = numpy.where(numpy.diff(spikes[:, 0]) > 0)[0] + 1
             break_points = numpy.concatenate(([0], break_points))
@@ -815 +778 @@
                 id = spikes[break_points[idx], 0]
                 if id in id_list:
-                    self.spiketrains[id] = SpikeTrain(spikes[break_points[idx]:break_points[idx+1], 1], self.t_start, self.t_stop)
-        
-        if len(self) > 0 and (self.t_start is None or self.t_stop is None):
+                    self.spiketrains[id] = SpikeTrain(spikes[break_points[idx]:break_points[idx+1], 1], interval=self.interval)
+        
+        if len(self) > 0 and (self.interval is None):
             self.__calc_startstop()
 
@@ -847 +810 @@
         """
         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()], numpy.float32)
-                self.t_start = numpy.min(start_times)
-                logging.debug("Warning, t_start is infered from the data : %f" %self.t_start)
+            if self.interval is None:
+                start_times = numpy.array([self.spiketrains[idx].interval.t_start() for idx in self.id_list()], numpy.float32)
+                t_start = numpy.min(start_times)
+                logging.debug("Warning, t_start is infered from the data : %f" %t_start)
+
+                stop_times = numpy.array([self.spiketrains[idx].interval.t_stop() for idx in self.id_list()], numpy.float32)
+                t_stop  = numpy.max(stop_times)
+                logging.debug("Warning, t_stop  is infered from the data : %f" %t_stop)
+
                 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()], numpy.float32)
-                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
+                    self.spiketrains[id].interval = Interval((t_start, t_stop))
         else:
             raise Exception("No SpikeTrains")
@@ -952 +914 @@
             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)
+            self.spiketrains[id] = spktrain.interval_slice(self.interval)
             
     def time_parameters(self):
@@ -958 +920 @@
         Return the time parameters of the SpikeList (t_start, t_stop)
         """
-        return (self.t_start, self.t_stop)
+        return (self.interval.t_start(), self.interval.t_stop())
 
     def time_axis(self, time_bin):
@@ -971 +933 @@
         """
         if newnum:
-            axis = numpy.arange(self.t_start, self.t_stop+time_bin, time_bin)
-        else:
-            axis = numpy.arange(self.t_start, self.t_stop, time_bin)
+            axis = numpy.arange(self.interval.t_start(), self.interval.t_stop() + time_bin, time_bin)
+        else:
+            axis = numpy.arange(self.interval.t_start(), self.interval.t_stop(), time_bin)
         return axis
 
@@ -993 +955 @@
         # We check that Spike Lists have similar time_axis
         for sl in spklists:
-            if not sl.time_parameters() == self.time_parameters():
+            if not sl.interval.is_equal(self.interval):
                 raise Exception("Spike Lists should have similar time_axis")
         for sl in spklists:
@@ -1043 +1005 @@
         missing_ids = id_list.difference(set(self.id_list()))
         for id in missing_ids:
-            self.append(id, SpikeTrain([],self.t_start, self.t_stop))
+            self.append(id, SpikeTrain([],self.interval))
         
     
@@ -1066 +1028 @@
             time_slice, interval_slice
         """
-        new_SpkList = SpikeList([], [], self.t_start, self.t_stop, self.dimensions)
+        new_SpkList = SpikeList([], [], self.interval, self.dimensions)
         id_list = self.__sub_id_list(id_list)
         for id in id_list:
@@ -1086 +1048 @@
             id_slice, interval_slice
         """
-        new_SpkList = SpikeList([], [], 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
+        interval = Interval((t_start, t_stop))
+        return self.interval_slice(interval)
+        
         
     def interval_slice(self, interval):
@@ -1103 +1063 @@
             id_slice, time_slice
         """
-        t_start, t_stop = interval.time_parameters()
-        new_SpkList = SpikeList([], [], t_start, t_stop, self.dimensions)
+        new_SpkList = SpikeList([], [], interval, self.dimensions)
         for id in self.id_list():
             new_SpkList.append(id, self.spiketrains[id].interval_slice(interval))
@@ -1124 +1083 @@
                 1050
         """
-        self.t_start += offset
-        self.t_stop  += offset
+        self.interval.offset_full(offset)
         for id in self.id_list():
             self.spiketrains[id].time_offset(offset)
@@ -1156 +1114 @@
         Get the time of the first real spike in the SpikeList
         """
-        first_spike = self.t_stop
+        first_spike = self.interval.t_stop()
         is_empty    = True
         for id in self.id_list():
@@ -1172 +1130 @@
         Get the time of the last real spike in the SpikeList
         """
-        last_spike = self.t_start
+        last_spike = self.interval.t_start()
         is_empty    = True
         for id in self.id_list():
@@ -1458 +1416 @@
             
 
-    def mean_rate(self, t_start=None, t_stop=None):
+    def mean_rate(self, interval=None):
         """
         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
+            interval - An interval object
         
         If t_start or t_stop are not defined, those of the SpikeList are used
@@ -1478 +1435 @@
 
 
-    def mean_rate_std(self, t_start=None, t_stop=None):
+    def mean_rate_std(self, interval=None):
         """
         Standard deviation of the firing rates accross all SpikeTrains 
@@ -1484 +1441 @@
 
         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
+            interval - An interval object
         
         If t_start or t_stop are not defined, those of the SpikeList are used
@@ -1496 +1452 @@
             mean_rate, mean_rates
         """
-        return numpy.std(self.mean_rates(t_start, t_stop))
-
-
-    def mean_rates(self, t_start=None, t_stop=None):
+        return numpy.std(self.mean_rates(interval))
+
+
+    def mean_rates(self, interval=None):
         """ 
         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
+            interval - An interval object
         
         If t_start or t_stop are not defined, those of the SpikeList are used
@@ -1514 +1469 @@
         rates = []
         for id in self.id_list():
-            rates.append(self.spiketrains[id].mean_rate(t_start, t_stop))
+            rates.append(self.spiketrains[id].mean_rate(interval))
         return rates
     
@@ -2032 +1987 @@
 
 
-    def raster_plot(self, id_list=None, t_start=None, t_stop=None, display=True, kwargs={}):
+    def raster_plot(self, id_list=None, interval=None, display=True, kwargs={}):
         """
         Generate a raster plot for the SpikeList in a subwindow of interest,
@@ -2060 +2015 @@
             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 interval is None:
+            interval = spk.interval
+
+        if not interval.is_equal(spk.interval):
+            spk = spk.interval_slice(interval)
 
         if not subplot or not HAVE_PYLAB: