Changeset 486

Timestamp:

07/14/11 18:32:11 (10 months ago)

Author:

mpereira

bzr:base-revision:

svn-v4:400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk:485

bzr:committer:

Michael Pereira <michael.fsp@gmail.com>

bzr:file-ids:

debian/rules 351@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fdebian%2Frules
doc/signals.txt 273@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fdoc%2Fsignals.txt
examples/retina/retina.py 77@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fexamples%2Fretina%2Fretina.py
examples/single_neuron/simple_single_neuron.py 78@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fexamples%2Fsingle_neuron%2Fsimple_single_neuron.py
src/analysis.py 106@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fanalysis.py
src/datastore/__init__.py 213@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fdatastore%2F__init__.py
src/datastore/django_orm/__init__.py 213@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fdatastore%2Fdjango_orm%2F__init__.py
src/datastore/django_orm/models.py 213@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fdatastore%2Fdjango_orm%2Fmodels.py
src/datastore/interface.py 213@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fdatastore%2Finterface.py
src/datastore/shelve_ds.py 213@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fdatastore%2Fshelve_ds.py
src/facets/mixedutils.py 284@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Ffacets%2Fmixedutils.py
src/plotting.py 106@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fplotting.py
src/signals/intervals.py 390@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fsignals%2Fintervals.py
src/signals/spikes.py 375@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fsignals%2Fspikes.py
src/spike2/__init__.py 214@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fspike2%2F__init__.py
src/spike2/sonpy/README 214@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fspike2%2Fsonpy%2FREADME
src/spike2/sonpy/__init__.py 214@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fspike2%2Fsonpy%2F__init__.py
src/spike2/sonpy/_marker.py 214@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fspike2%2Fsonpy%2F_marker.py
src/spike2/sonpy/_waveform.py 214@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fspike2%2Fsonpy%2F_waveform.py
src/spike2/sonpy/son.py 214@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Fspike2%2Fsonpy%2Fson.py
src/utilities/__init__.py 225@400bb3d0-8d2e-0410-a2c0-f692ac833539:trunk%2Fsrc%2Futilities%2F__init__.py

bzr:mapping-version:

v4

bzr:merge:

michael.fsp@gmail.com-20110714162651-g03c2jgqof2dangy

bzr:repository-uuid:

400bb3d0-8d2e-0410-a2c0-f692ac833539

bzr:revision-id:

michael.fsp@gmail.com-20110714163140-v3qk0357xh8xn2zv

bzr:revno:

430

bzr:revprop:branch-nick:

trunk

bzr:root:

trunk

bzr:text-revisions:

debian/rules michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
doc/signals.txt michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
examples/retina/retina.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
examples/single_neuron/simple_single_neuron.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/analysis.py michael.fsp@gmail.com-20110714162651-g03c2jgqof2dangy
src/datastore/__init__.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/datastore/django_orm/__init__.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/datastore/django_orm/models.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/datastore/interface.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/datastore/shelve_ds.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/facets/mixedutils.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/plotting.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/signals/intervals.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/signals/spikes.py michael.fsp@gmail.com-20110714162651-g03c2jgqof2dangy
src/spike2/__init__.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/spike2/sonpy/README michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/spike2/sonpy/__init__.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/spike2/sonpy/_marker.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/spike2/sonpy/_waveform.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/spike2/sonpy/son.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah
src/utilities/__init__.py michael.fsp@gmail.com-20110526110242-ofx6nnkfmr7delah

bzr:timestamp:

2011-07-14 18:31:40.723999977 +0200

bzr:user-agent:

bzr2.3.1+bzr-svn1.0.5dev0

svn:original-date:

2011-07-14T16:31:40.724000Z

Message:

major improvements of SpikeTrain.instantaneous_rate()

Location:

trunk

Files:

• debian/rules (modified) (1 prop)
• doc/signals.txt (modified) (1 prop)
• examples/retina/retina.py (modified) (1 prop)
• examples/single_neuron/simple_single_neuron.py (modified) (1 prop)
• src/analysis.py (modified) (4 diffs)
• src/datastore/__init__.py (modified) (1 prop)
• src/datastore/django_orm/__init__.py (modified) (1 prop)
• src/datastore/django_orm/models.py (modified) (1 prop)
• src/datastore/interface.py (modified) (1 prop)
• src/datastore/shelve_ds.py (modified) (1 prop)
• src/facets/mixedutils.py (modified) (1 prop)
• src/plotting.py (modified) (1 prop)
• src/signals/intervals.py (modified) (3 diffs)
• src/signals/spikes.py (modified) (5 diffs)
• src/spike2/__init__.py (modified) (1 prop)
• src/spike2/sonpy/README (modified) (1 prop)
• src/spike2/sonpy/__init__.py (modified) (1 prop)
• src/spike2/sonpy/_marker.py (modified) (1 prop)
• src/spike2/sonpy/_waveform.py (modified) (1 prop)
• src/spike2/sonpy/son.py (modified) (1 prop)
• src/utilities/__init__.py (modified) (1 prop)

trunk/src/analysis.py

@@ -143 +143 @@
                sua2. Thus positive values of int relate to events of sua2 that
                lead events of sua1. Units are the same as the input arrays.
-        int_ - predictior: accumulated differences based on the prediction
+        int_ - predictor: accumulated differences based on the prediction
                method. The length of int_ is n_pred * length(int).  Units are
                the same as the input arrays.
@@ -216 +216 @@
     
     # Plot the results if display=True   
-    subplot = get_display(display)
-    if not subplot or not HAVE_PYLAB:
+    if display:
+        subplot = get_display(display)
+        if not subplot or not HAVE_PYLAB:
+            return int, int_, norm
+        else:
+            # Plot the cross-correlation
+            try:
+                counts, bin_edges = numpy.histogram(int, **kwargs)
+                edge_distances = numpy.diff(bin_edges)
+                bin_centers = bin_edges[1:] - edge_distances/2
+                counts = counts / norm
+                xlabel = "Time"
+                ylabel = "Cross-correlation coefficient"
+                #NOTE: the x axis corresponds to the upper edge of each bin
+                subplot.plot(bin_centers, counts, label='cross-correlation', color='b')
+                if predictor is None:
+                    set_labels(subplot, xlabel, ylabel)
+                    pylab.draw()
+                elif predictor is 'shuffle':            
+                    # Plot the predictor
+                    norm_ = norm * n_pred
+                    counts_, bin_edges_ = numpy.histogram(int_, **kwargs)
+                    counts_ = counts_ / norm_
+                    subplot.plot(bin_edges_[1:], counts_, label='predictor')
+                    subplot.legend()
+                    pylab.draw()
+            except ValueError:
+                print "There are no correlated events within the selected lag"\
+                " window of %s" % lag
+    else:
         return int, int_, norm
-    else:
-        # Plot the cross-correlation
-        try:
-            counts, bin_edges = numpy.histogram(int, **kwargs)
-            counts = counts / norm
-            xlabel = "Time"
-            ylabel = "Cross-correlation coefficient"
-            #NOTE: the x axis corresponds to the upper edge of each bin
-            subplot.plot(bin_edges[1:], counts, label='cross-correlation')
-            if predictor is None:
-                set_labels(subplot, xlabel, ylabel)
-                pylab.draw()
-            elif predictor is 'shuffle':            
-                # Plot the predictor
-                norm_ = norm * n_pred
-                counts_, bin_edges_ = numpy.histogram(int_, **kwargs)
-                counts_ = counts_ / norm_
-                subplot.plot(bin_edges_[1:], counts_, label='predictor')
-                subplot.legend()
-                pylab.draw()
-        except ValueError:
-            print "There are no correlated events within the selected lag"\
-            " window of %s" % lag
+    
 def _dict_max(D):
     """
@@ -254 +260 @@
             return k
         
-def make_kernel(form, sigma, time_stamp_resolution, direction=1, kernel=None,
-                m_idx=None):
+def make_kernel(form, sigma, time_stamp_resolution, direction=1):
     """
     Constructs a numeric linear convolution kernel of basic shape to be used
@@ -312 +317 @@
         m_idx   - index of the numerically determined median (center of gravity)
                   of the kernel function
+                  
+    Further comments:
+    
+    Assume matrix X of n spike trains represented as binary vector (0/1).
+    
+    To obtain single trial rate function of trial one should use
+        r = norm * numpy.convolve(sua, kernel)
+    To obtain trial-averaged spike train one should use
+        r_avg = norm * numpy.convolve(sua, np.mean(X,1))
+        
+        where X is an array of shape (l,n), where n is the number of trials and
+        l the length of each trial
+
+    Note that filter assumes a causal filtering. To represent the resulting
+    firing rate function at center of gravity (acausal) use
+
+        r=r(length(kernel):end);    -> discard first length(kernel)-1
+                                       elements which considered zero entries
+                                       (border effect)
+        t=((1:length(r))+m_idx)*TimeStampResolution -> in seconds
+        plot(t,r)   -> in second time scale
     """
     assert form.upper() in ('BOX','TRI','GAU','EPA','EXP','ALP'), "form must \
-    be one either 'BOX','TRI','GAU','EPA','EXP' or 'ALP'!"
+    be one of either 'BOX','TRI','GAU','EPA','EXP' or 'ALP'!"
     
     assert direction in (1,-1), "direction must be either 1 or -1"

trunk/src/signals/intervals.py

@@ -4 +4 @@
 
 if HAVE_INTERVAL:
-    from interval import * 
+    from interval import *
 
 import numpy
@@ -31 +31 @@
             self.end_times   = end_times
             # write the intervals to an interval object (pyinterval)
-            test = isinstance(start_times, int) or isinstance(start_times, float)
+            scalar_types = (int, float, numpy.float, numpy.float32, numpy.float64, numpy.int, numpy.int8,
+                            numpy.int16, numpy.int32, numpy.int64)
+            test = isinstance(start_times, scalar_types)
             if test:
-                start_times = [start_times]
-            test = isinstance(end_times, int) or isinstance(end_times, float)
+                self.start_times = [self.start_times]
+            test = isinstance(end_times, scalar_types)
             if test:
-                end_times = [end_times]
-            if len(start_times) != len(end_times) :
+                self.end_times = [self.end_times]
+            if len(self.start_times) != len(self.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])))
@@ -45 +47 @@
             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
+            self.start_times = [start_times]
+            self.end_times   = [end_times]
 
     def intersect(self, itv) :

trunk/src/signals/spikes.py

@@ -27 +27 @@
 
 import os, re, numpy
+import scipy.signal
 import logging
 from NeuroTools import check_dependency, check_numpy_version
@@ -85 +86 @@
     #######################################################################
     def __init__(self, spike_times, t_start=None, t_stop=None):
+        #TODO: add information about sampling rate at time of creation
+        
         """
         Constructor of the SpikeTrain object
@@ -449 +452 @@
         if not subplot or not HAVE_PYLAB:
             print PYLAB_ERROR
-            return
         else:
             if len(spikes) > 0:
@@ -570 +572 @@
             hist = hist.astype(bool).astype(int)
         return hist
-
+    
+    def instantaneous_rate(self, resolution, kernel, norm, t_start=None,
+                           t_stop=None, acausal=True, m_idx=None,
+                           trim=False):
+        """
+        Estimate instantaneous firing rate by kernel convolution.
+        
+        Inputs:
+            resolution  - time stamp resolution of the spike times (ms). the
+                          same resolution will be assumed for the kernel
+            kernel      - kernel function used to convolve with
+            norm        - normalization factor associated with kernel function
+                          (see analysis.make_kernel for details)
+            t_start     - start time of the interval used to compute the firing
+                          rate
+            t_stop      - end time of the interval used to compute the firing
+                          rate (included)
+            acausal     - if True, acausal filtering is used, i.e., the gravity
+                          center of the filter function is aligned with the
+                          spike to convolve
+            m_idx       - index of the value in the kernel function vector that
+                          corresponds to its gravity center. this parameter is
+                          not mandatory for symmetrical kernels but it is
+                          required when assymmetrical kernels are to be aligned
+                          at their gravity center with the event times
+            trim        - if True, only the 'valid' region of the convolved
+                          signal are returned, i.e., the points where there
+                          isn't complete overlap between kernel and spike train
+                          are discarded
+                          NOTE: if True and an assymetrical kernel is provided
+                          the output will not be aligned with [t_start, t_stop]
+        """
+        
+        if t_start is None:
+            t_start = self.t_start
+            
+        if t_stop is None:
+            t_stop = self.t_stop
+        
+        if m_idx is None:
+            m_idx = kernel.size / 2
+        
+        time_vector = numpy.zeros((t_stop - t_start)/resolution + 1)
+        
+        spikes_slice = self.spike_times[(self.spike_times >= t_start) & (
+            self.spike_times <= t_stop)]
+        
+        for spike in spikes_slice:
+            index = (spike - t_start) / resolution
+            time_vector[index] = 1
+        
+        r = norm * scipy.signal.fftconvolve(time_vector, kernel, 'full')
+
+        if acausal is True:
+            if trim is False:
+                r = r[m_idx:-(kernel.size - m_idx)]
+                t_axis = numpy.linspace(t_start, t_stop, r.size)
+                return t_axis, r
+            
+            elif trim is True:
+                r = r[2 * m_idx:-2*(kernel.size - m_idx)]
+                t_start = t_start + m_idx * resolution
+                t_stop = t_stop - ((kernel.size) - m_idx) * resolution
+                t_axis = numpy.linspace(t_start, t_stop, r.size)
+                return t_axis, r
+            
+        if acausal is False:
+            if trim is False:
+                r = r[m_idx:-(kernel.size - m_idx)]
+                t_axis = (numpy.linspace(t_start, t_stop, r.size) +
+                          m_idx * resolution)
+                return t_axis, r
+            
+            elif trim is True:
+                r = r[2 * m_idx:-2*(kernel.size - m_idx)]
+                t_start = t_start + m_idx * resolution
+                t_stop = t_stop - ((kernel.size) - m_idx) * resolution
+                t_axis = (numpy.linspace(t_start, t_stop, r.size) +
+                          m_idx * resolution)
+                return t_axis, r
+    
     def relative_times(self):
         """
@@ -845 +927 @@
         ## using a generator to build the SpikeList...
         
-        if not hasattr(spikes, 'size'): # is an array:
+        if not isinstance(spikes, numpy.ndarray): # is an array:
             spikes = numpy.array(spikes, numpy.float32)
         N = len(spikes)