Changeset 431

Timestamp:

08/26/09 15:10:23 (3 years ago)

Author:

pierre

Message:

Add the interval notion to the AnalogSignal? object, now debugged. Start to fix the AnalogSignalList? objects...

Location:

branches/interval/src/signals

Files:

• analogs.py (modified) (13 diffs)
• intervals.py (modified) (6 diffs)
• spikes.py (modified) (2 diffs)

branches/interval/src/signals/analogs.py

@@ -31 +31 @@
 """
 
-import os, re, numpy
+import os, re, numpy, copy
 from NeuroTools import check_dependency, check_numpy_version, analysis
 from NeuroTools.io import *
@@ -78 +78 @@
     def __init__(self, signal, dt, t_start=0, t_stop=None, interval=None):
         #logging.debug("Creating AnalogSignal. len(signal)=%d, dt=%g, t_start=%g, t_stop=%s" % (len(signal), dt, t_start, t_stop))
-        self.signal  = signal
-        self.dt      = float(dt)
+        self.signal   = signal
+        self.dt       = float(dt)
         self.interval = self.create_intervals(t_start, t_stop, interval)
         self.init_times()
@@ -88 +88 @@
         if self.t_start >= self.t_stop:
             raise Exception("Incompatible time interval for the creation of the AnalogSignal. t_start=%s, t_stop=%s" % (self.t_start, self.t_stop))
+
+        self.times  = self.interval.time_axis()
+        idx         = self.interval.idx_slice_times(self.times)
+        self.signal = self.signal[idx]
 
     def create_intervals(self, t_start=None, t_stop=None, interval=None):
@@ -95 +99 @@
             # interval is fully defined by the user
             if type(interval) == Interval:
-                interval_out = interval.copy()
+                interval_out = interval.copy()
             else :
-                interval_out = Interval(interval) 
+                interval_out = Interval(interval)
         else:
             if t_start is None and t_stop is None :
@@ -116 +120 @@
                     else :
                         interval_out = Interval(t_start)
+        
         return interval_out
 
@@ -135 +140 @@
         Return the duration of the SpikeTrain
         """
-        return self.t_stop - self.t_start
+        return self.interval.total_duration()
 
     def __str__(self):
@@ -156 +161 @@
         Return a copy of the AnalogSignal object
         """
-        return AnalogSignal(self.signal, self.dt, self.t_start, self.t_stop)
+        return copy.deepcopy(self)
 
     def time_axis(self, normalized=False):
@@ -169 +174 @@
             norm = 0.
 
-        time_pavement = numpy.arange(self.t_start-norm, self.t_stop-norm, self.dt)
-        return SpikeTrain(time_pavement, interval=self.interval).spike_times
+        time_pavement = numpy.arange(self.t_start, self.t_stop, self.dt) - norm
+        return self.interval.slice_times(time_pavement)
     
     def time_offset(self, offset):
@@ -239 +244 @@
         assert t_stop <= self.t_stop
         assert t_stop > t_start
-        
-        t = self.time_axis()
-        i_start = int((t_start-self.t_start)/self.dt)
-        i_stop = int((t_stop-self.t_start)/self.dt)
-        signal = self.signal[i_start:i_stop]
-
-        return AnalogSignal(signal, self.dt, t_start, t_stop)
+        interval = Interval((t_start, t_stop))
+        return self.interval_slice(interval)
 
     def interval_slice(self, interval):
@@ -261 +261 @@
             time_slice
         """
-        signal_list = []
-
-        if type(interval) is not Interval :
-            interval = Interval(interval)
-
-        for itv in interval.interval_data :
-            signal_list.append(self.signal[(itv[0]-self.t_start)/self.dt:(itv[1]-self.t_start)/self.dt])
-
-        return AnalogSignal(signal_list, self.dt, interval)
-
-    def threshold_detection(self, threshold=None, format=None,sign='above'):
+        interval = interval.intersect(self.interval)
+        idx      = interval.idx_slice_times(self.times)
+        return AnalogSignal(self.signal[idx], self.dt, interval=interval)
+
+    def threshold_detection(self, threshold=None, format=None, sign='above'):
         """
         Returns the times when the analog signal crosses a threshold.
@@ -353 +347 @@
         # recalculate everything into timesteps, is more stable against rounding errors
         # and subsequent cutouts with different sizes
-        events = numpy.floor(numpy.array(events)/self.dt)
-        t_min_l = numpy.floor(t_min/self.dt)
-        t_max_l = numpy.floor(t_max/self.dt)
-        t_start = numpy.floor(self.t_start/self.dt)
-        t_stop = numpy.floor(self.t_stop/self.dt)
-        
-        for spike in events:
-            if ((spike-t_min_l) >= t_start) and ((spike+t_max_l) < t_stop):
-                spike = spike - t_start
-                if average:
-                    result += self.signal[(spike-t_min_l):(spike+t_max_l)]
-                else:
-                    result.append(self.signal[(spike-t_min_l):(spike+t_max_l)])
-                Nspikes += 1
+        t_min_l  = numpy.floor(t_min/self.dt)
+        t_max_l  = numpy.floor(t_max/self.dt)
+        result   = numpy.zeros((t_min_l+t_max_l), numpy.float32)
+        events_interval = copy.copy(self.interval)
+        events_interval.offset_start(t_min, colapse=True)
+        events_interval.offset_stop(-t_max, colapse=True)
+        events = events_interval.slice_times(events)
+        assert len(events) > 0, "the PSTH windows [event-t_min, event+t_max] should be included within the spike train interval"
+        for ev in events:
+            ev = numpy.floor((ev - self.interval.t_start())/self.dt)
+            if average:
+                result += self.signal[(spike-t_min_l):(spike+t_max_l)]
+            else:
+                result.append(self.signal[(spike-t_min_l):(spike+t_max_l)])
+            Nspikes += 1
         if average:
             result = result/Nspikes
@@ -415 +410 @@
         assert (t_min >= 0) and (t_max >= 0), "t_min and t_max should be greater than 0"
         assert len(events) > 0, "events should not be empty and should contained at least one element"
-        
+        events_interval = copy.copy(self.interval)
+        events_interval.offset_start(t_min, colapse=True)
+        events_interval.offset_stop(-t_max, colapse=True)
+        events = events_interval.slice_times(events)
+        assert len(events) > 0, "the PSTH windows [event-t_min, event+t_max] should be included within the spike train interval"
         result = {}
-        for index, spike in enumerate(events):
-            if ((spike-t_min) >= self.t_start) and ((spike+t_max) < self.t_stop):
-                spike = spike - self.t_start
-                t_start_new = (spike-t_min)
-                t_stop_new = (spike+t_max)
-                result[index] = self.time_slice(t_start_new, t_stop_new)
+        for index, ev in enumerate(events):
+            t_start_new = (ev-t_min)
+            t_stop_new  = (ev+t_max)
+            result[index] = self.time_slice(t_start_new, t_stop_new)
         return result
-        
-            
+
+
+
 class AnalogSignalList(object):
     """
@@ -537 +535 @@
         if len(self) > 0:
             errmsgs = []
-            print numpy.shape(val.signal)
             val = val.interval_slice(self.interval)
-            print numpy.shape(val.signal)
             if not val.interval.is_equal(self.interval) :
                 raise Exception("the provided signal doesn't cover the AnalogSignalList interval")

branches/interval/src/signals/intervals.py

@@ -5 +5 @@
     from interval import * 
 
-import numpy
+import numpy, copy
 
 class Interval(object):
@@ -50 +50 @@
     def __iter__(self):
         return iter(self.sub_intervals)
+    
+    def __contains__(self, other):
+        return numpy.all(numpy.any(x.inf <= y.inf and y.sup <= x.sup for x in self.interval_data) for y in other)
     
     def __getslice__(self, i, j):
@@ -102 +105 @@
         Return a copy of the SpikeTrain object
         """
-        return Interval(self.sub_intervals)
+        return copy.deepcopy(self)
 
     def offset_start(self, shift, from_stop=False, colapse=False) :
@@ -163 +166 @@
 
     def slice_times(self, times) :
-        if times.__class__.__name__ == 'SpikeTrain' :
-            times = numpy.array(times.spike_times)
-        else :
-            times = numpy.array(times)
         spikes_selector = numpy.zeros(len(times), dtype=numpy.bool)
+        times = numpy.array(times)
         if HAVE_INTERVAL:
             for itv in self.interval_data :
@@ -174 +174 @@
             spikes_selector = (times >= self.t_start()) & (times <= self.t_stop())
         return numpy.extract(spikes_selector, times)
+    
+    def time_axis(self, dt=0.1):
+        result = numpy.array([], float)
+        if HAVE_INTERVAL:
+            for itv in self.interval_data :
+                result = numpy.concatenate((result, numpy.arange(itv[0],itv[1],dt)))
+        else:
+            return numpy.arange(self.t_start(),self.t_stop(),dt)
+        return result
 
     def idx_slice_times(self, times):
         spikes_selector = numpy.zeros(len(times), dtype=numpy.bool)
         if HAVE_INTERVAL:
-            for itv in self.interval_data :
-                spikes_selector = spikes_selector + (times > itv[0])*(times <= itv[1])
+            for itv in self.interval_data:
+                spikes_selector = spikes_selector + (times >= itv[0])*(times < itv[1])
         else:
             spikes_selector = (times >= self.t_start()) & (times <= self.t_stop())
@@ -186 +195 @@
     def is_equal(self, itv) :
         return self.sub_intervals == itv.sub_intervals
-
 
     def intersect(self, itv) :

branches/interval/src/signals/spikes.py

@@ -197 +197 @@
         Return a copy of the SpikeTrain object
         """
-        return SpikeTrain(self.spike_times, interval=self.interval)
+        return copy.deepcopy(self)
 
 
@@ -869 +869 @@
         """
 
-        spklist = SpikeList([], [], None, None, self.interval, self.dimensions)
-        for id in self.id_list():
-            spklist.append(float(id), self.spiketrains[float(id)])
-        return spklist
+        copy.deepcopy(self)
 
     def create_intervals(self, t_start=None, t_stop=None, interval=None):