Changeset 453

Timestamp:

08/23/10 14:25:44 (18 months ago)

Author:

mpereira

Message:

Added crosscorrelate() to analysis.py

Location:

trunk

Files:

• src/analysis.py (modified) (3 diffs)
• src/signals/spikes.py (modified) (1 diff)
• test/analysis/crosscorrelate/spike_data (added)
• test/test_analysis.py (modified) (3 diffs)

trunk/src/analysis.py

@@ -22 +22 @@
 import os, numpy
 from NeuroTools import check_dependency
-
 
 def arrays_almost_equal(a, b, threshold):
@@ -45 +44 @@
 
     """
-    assert(x.ndim == y.ndim, "Inconsistent shape !")
+    assert x.ndim == y.ndim, "Inconsistent shape !"
 #    assert(x.shape == y.shape, "Inconsistent shape !")
     if axis is None:
@@ -79 +78 @@
     return iFxy/varxy
 
+from NeuroTools.plotting import get_display, set_labels
+
+HAVE_PYLAB = check_dependency('pylab')
+
+
+def crosscorrelate(sua1, sua2, lag=None, n_pred=1, predictor=None,
+                   display=False, kwargs={}):
+    """
+    Calculates the cross-correlation between two vectors containing event times.
+    Returns (int, int_, norm). See below for details.
+    
+    Adapted from original script written by Martin P. Nawrot for the FIND MATLAB
+    toolbox.
+    See FIND - a unified framework for neural data analysis,
+        Meier R, Egert U, Aertsen A, Nawrot MP; Neural Netw. 2008 Oct;
+        21(8):1085-93. 
+     
+    Inputs:
+        sua1      - array of event times. Can be either a column/row vector.
+        sua2      - array of event times. Can be either a column/row vector.
+                    If sua2 == sua1 the result is the
+                    autocorrelogram.
+        lag       - the max. lag for which relative event timing is considered
+                    with a max. difference of +/- lag. A default lag is computed
+                    from the inter-event interval of the longer of the two sua
+                    arrays 
+        n_pred    - number of surrogate compilations for the predictor. This
+                    influences the total length of the predictor output array
+                    int_
+        predictor - string array determines the type of bootstrap predictor to
+                    be used:
+                        shuffle - shuffles inter-event intervals of the longer
+                                  input array and calculates relative
+                                  differences with the shorter input array.
+                                  n_pred determines the number of repeated
+                                  shufflings, resulting differences are pooled
+                                  from all repeated shufflings
+        display   - if True the corresponding plots will be displayed. If False,
+                    int, int_ and norm will be returned.
+                    when display = True and n_pred > 1, the averaged predictor
+                    will be plotted.
+        kwargs    - arguments to be passed to numpy.histogram.
+
+    Outputs:
+        int  - accumulated differences of events in sua1 minus the events in
+               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
+               method. The length of int_ is n_pred * length(int).  Units are
+               the same as the input arrays.
+        norm - normalization factor used to scale the bin heights in int and
+               int_. int/norm and int_/norm correspond to the linear
+               correlation coefficient.
+    
+    Examples:
+        >> crosscorrelate(numpy_array1, numpy_array2)
+        >> crosscorrelate(spike_train1.spike_times, spike_train2.spike_times)
+        >> crosscorrelate(spike_train1.spike_times, spike_train2.spike_times,
+                          lag = 150.0)
+        >> crosscorrelate(spike_train1.spike_times, spike_train2.spike_times,
+                          display=True, kwargs={'bins':100})
+            
+    See also:
+        ccf
+    """    
+    assert predictor is 'shuffle' or predictor is None, "predictor must be \
+    either None or 'shuffle'. Other predictors are not yet implemented."
+    
+    #Check whether sua1 and sua2 are SpikeTrains or arrays
+    sua = []
+    for x in (sua1, sua2):
+        if x.ndim == 1:
+            sua.append(x)
+        elif x.ndim == 2 and (x.shape[0] == 1 or x.shape[1] == 1):
+            sua.append(x.ravel())
+        else:
+            raise TypeError("sua1 and sua2 must be either instances of the \
+                            SpikeTrain class or column/row vectors")
+    sua1 = sua[0]
+    sua2 = sua[1]
+    
+    if sua1.size < sua2.size:
+        if lag is None:
+            lag = numpy.ceil(10*numpy.mean(numpy.diff(sua1)))
+        reverse = False
+    else:
+        if lag is None:
+            lag = numpy.ceil(20*numpy.mean(numpy.diff(sua2)))
+        sua1, sua2 = sua2, sua1
+        reverse = True
+            
+    #construct predictor
+    if predictor is 'shuffle':
+        isi = numpy.diff(sua2)
+        sua2_ = numpy.array([])
+        for ni in xrange(1,n_pred+1):
+            idx = numpy.random.permutation(isi.size-1)
+            
+            sua2_ = numpy.append(sua2_, numpy.add(numpy.insert(
+                (numpy.cumsum(isi[idx])), 0, 0), sua2.min() + (
+                numpy.random.exponential(isi.mean()))))
+            
+    #calculate cross differences in spike times
+    int = numpy.array([])
+    int_ = numpy.array([])
+    for k in xrange(0, sua1.size):
+        int = numpy.append(int, sua1[k] - sua2[numpy.nonzero(
+            (sua2 > sua1[k] - lag) & (sua2 < sua1[k] + lag))])
+    if predictor == 'shuffle':
+        for k in xrange(0, sua1.size):
+            int_ = numpy.append(int_, sua1[k] - sua2_[numpy.nonzero(
+                (sua2_ > sua1[k] - lag) & (sua2_ < sua1[k] + lag))])
+    
+    if reverse is True:
+        int = -int
+        int_ = -int_
+        
+    norm = numpy.sqrt(sua1.size * sua2.size)
+    
+    # Plot the results if display=True
+    subplot = get_display(display)
+    if not subplot or not HAVE_PYLAB:
+        return int, int_, norm
+    else:
+        # Plot the cross-correlation
+        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()
 
 def _dict_max(D):

trunk/src/signals/spikes.py

@@ -28 +28 @@
 import os, re, numpy
 import logging
-from NeuroTools import check_dependency, check_numpy_version, analysis
+from NeuroTools import check_dependency, check_numpy_version
+try:
+    from NeuroTools import analysis
+except ImportError:
+    pass
+
 from NeuroTools.io import *
 from NeuroTools.plotting import get_display, set_axis_limits, set_labels, SimpleMultiplot, progress_bar

trunk/test/test_analysis.py

@@ -7 +7 @@
 from numpy import pi, sin
 
-from NeuroTools import analysis
-
+from NeuroTools import analysis, signals
 
 # test simple_frequency_spectrum
@@ -27 +26 @@
         self.exp_reversed = scipy.io.loadmat(
             'analysis/make_kernel/exp_reversed.mat')
+        
+        #Used for the testcases of the crosscorrelate function
+        spk = signals.load_spikelist('analysis/crosscorrelate/spike_data')
+        self.spk0 = spk[0].spike_times
+        self.spk1 = spk[1].spike_times
 
     def testSimpleFrequencySpectrum(self):
@@ -153 +157 @@
         self.assertEqual(true_m_idx, m_idx, 3)
         
+    def testCrosscorrelateNoLag(self):
+            int, int_, norm = analysis.crosscorrelate(self.spk0, self.spk1)
+            #The following are output was generated with the FIND MATLAB toolbox
+            matlab_int = numpy.loadtxt('analysis/crosscorrelate/out_matlab_int')
+            numpy.testing.assert_array_almost_equal(int, matlab_int,
+                                                    decimal=3)
+            #The int_ output has a random component and for this reason the test
+            # cases are not as trivial
+            
+    def testCrosscorrelateLag100(self):
+        """Test case with lag within the length of the input array
+        """
+        int, int_, norm, = analysis.crosscorrelate(self.spk0, self.spk1,
+                                                   lag=100.0)
+        matlab_int = numpy.loadtxt('analysis/crosscorrelate/out_matlab_int_lag_100')
+        numpy.testing.assert_array_almost_equal(int, matlab_int, decimal = 3)
+            
+    def testCrosscorrelateLag500(self):
+        """Test case with lag is higher than the trial length
+        """
+        int, int_, norm = analysis.crosscorrelate(self.spk0, self.spk1,
+                                                  lag=500.0)
+        matlab_int = numpy.loadtxt('analysis/crosscorrelate/out_matlab_int_lag_500')
+        numpy.testing.assert_array_almost_equal(int, matlab_int, decimal = 3)
+            
 if __name__ == "__main__":
     unittest.main()