Changeset 462

Timestamp:

09/02/10 15:42:45 (17 months ago)

Author:

pierre

Message:

Add a gamma_generator function into stgen, and simplify all the positions index in SpikeTrains? objects with dimensions

Location:

trunk/src

Files:

• signals/spikes.py (modified) (7 diffs)
• stgen.py (modified) (2 diffs)

trunk/src/signals/spikes.py

@@ -915 +915 @@
     #def __setslice__(self, i, j):
 
-
-
     def __setitem__(self, id, spktrain):
         assert isinstance(spktrain, SpikeTrain), "A SpikeList object can only contain SpikeTrain objects"
@@ -1072 +1070 @@
         
         Examples:
-            >> spklist.id_list()
+            >> spklist.id_list
                 [0,2,5]
             >> spklist.complete(arange(5))
-            >> spklist.id_list()
+            >> spklist.id_list
                 [0,1,2,3,4]
         """
@@ -1712 +1710 @@
             return id-offset
         if len(self.dimensions) == 2:
-            x = (id-offset) % self.dimensions[1]
-            y = self.dimensions[0] - 1 - numpy.floor((id-offset)/self.dimensions[1]).astype(int)
+            x = (id-offset) % self.dimensions[0]
+            y = ((id-offset)/self.dimensions[0]).astype(int)
             return (x,y)
 
@@ -1735 +1733 @@
         assert len(position) == len(tuple(self.dimensions)), "position does not have the correct shape !"
         if len(self.dimensions) == 1:
-            return position+offset
+            return position + offset
         if len(self.dimensions) == 2:
-            return (self.dimensions[0] - 1 - position[1])*self.dimensions[1] + position[0] + offset
+            return position[1]*self.dimensions[1] + position[0] + offset
 
 
@@ -1786 +1784 @@
             activity_map = numpy.zeros(self.dimensions, float)
             rates        = spklist.mean_rates()
-            id_offset = min(self.id_list())
+            #id_offset    = min(self.id_list())
+            #x,y          = spklist.id2position(spklist.id_list(), id_offset)
+            x,y          = spklist.id2position(spklist.id_list())            
+            #j,i = x, self.dimensions[0] - 1 - y
             for count, id in enumerate(spklist.id_list()):
-                x,y = spklist.id2position(id, id_offset)
-                j,i = x, self.dimensions[0] - 1 -y
-                activity_map[i,j] = rates[count]
+                #activity_map[i[count],j[count]] = rates[count]                            
+                activity_map[x[count],y[count]] = rates[count]
             if not subplot or not HAVE_PYLAB or not HAVE_MATPLOTLIB:
                 return activity_map
@@ -2247 +2247 @@
             time      = time[sort_idx]
             pos       = pos[sort_idx]
+            x,y       = spk.id2position(pos)
             max_idx   = len(time)-1
             logging.info('Making movie %s - this make take a while' % output)
@@ -2254 +2255 @@
                 while (t_start < t_stop):
                     activity_map = numpy.zeros(spk.dimensions)
-                    while ((time[idx] < t_start + time_bin) and (idx < max_idx)):                        
-                        x,y = spk.id2position(pos[idx])
-                        j,i = x, self.dimensions[0] - 1 -y
-                        activity_map[i,j] += 1
+                    while ((time[idx] < t_start + time_bin) and (idx < max_idx)):                                                
+                        #j,i = x, self.dimensions[0] - 1 -y
+                        activity_map[x[idx],y[idx]] += 1
                         idx += 1
                     im.set_array(activity_map)

trunk/src/stgen.py

@@ -255 +255 @@
             return spikes
 
-
+    def gamma_generator(self, a, b, t_start=0.0, t_stop=1000.0, array=False,debug=False):
+        """
+        Returns a SpikeTrain whose spikes are a realization of a Poisson process
+        with the given rate (Hz) and stopping time t_stop (milliseconds).
+
+        Note: t_start is always 0.0, thus all realizations are as if 
+        they spiked at t=0.0, though this spike is not included in the SpikeList.
+
+        Inputs:
+            a,b     - the parameters of the gamma process
+            t_start - the beginning of the SpikeTrain (in ms)
+            t_stop  - the end of the SpikeTrain (in ms)
+            array   - if True, a numpy array of sorted spikes is returned,
+                      rather than a SpikeTrain object.
+
+        Examples:
+            >> gen.poisson_generator(50, 0, 1000)
+            >> gen.poisson_generator(20, 5000, 10000, array=True)
+         
+        See also:
+            inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
+        """
+
+        #number = int((t_stop-t_start)/1000.0*2.0*rate)
+        
+        # less wasteful than double length method above
+        n = (t_stop-t_start)/1000.0*(a*b)
+        number = numpy.ceil(n+3*numpy.sqrt(n))
+        if number<100:
+            number = min(5+numpy.ceil(2*n),100)
+        
+        if number > 0:
+            isi = self.rng.gamma(a, b, number)*1000.0
+            if number > 1:
+                spikes = numpy.add.accumulate(isi)
+            else:
+                spikes = isi
+        else:
+            spikes = numpy.array([])
+
+        spikes+=t_start
+        i = numpy.searchsorted(spikes, t_stop)
+
+        extra_spikes = []
+        if i==len(spikes):
+            # ISI buf overrun
+            
+            t_last = spikes[-1] + self.rng.gamma(a, b, 1)[0]*1000.0
+
+            while (t_last<t_stop):
+                extra_spikes.append(t_last)
+                t_last += self.rng.gamma(a, b, 1)[0]*1000.0
+            
+            spikes = numpy.concatenate((spikes,extra_spikes))
+
+            if debug:
+                print "ISI buf overrun handled. len(spikes)=%d, len(extra_spikes)=%d" % (len(spikes),len(extra_spikes))
+
+
+        else:
+            spikes = numpy.resize(spikes,(i,))
+
+        if not array:
+            spikes = SpikeTrain(spikes, t_start=t_start,t_stop=t_stop)
+
+
+        if debug:
+            return spikes, extra_spikes
+        else:
+            return spikes
+            
     def inh_poisson_generator(self, rate, t, t_stop, array=False):
         """
@@ -395 +465 @@
             # find index in "t" time
             t_i = numpy.searchsorted(t[t_i:],ps[i],'right')-1+t_i
-            if rn[i]<gamma_hazard((ps[i]-t_last)/1000.0,a[t_i],b[t_i])/rmax:
+            if rn[i]<gamma_hazard_scipy((ps[i]-t_last)/1000.0,a[t_i],b[t_i])/rmax:
                 # keep spike
                 t_last = ps[i]