Changeset 451

Show
Ignore:
Timestamp:
08/23/10 13:47:57 (18 months ago)
Author:
mpereira
Message:

Added make_kernel() to analysis.py

Location:
trunk
Files:
14 added
3 modified

Legend:

Unmodified
Added
Removed

  • trunk/setup.py

    r400 r451  
    1515                'NeuroTools.datastore', 
    1616                'NeuroTools.datastore.django_orm', 
     17                'NeuroTools.optimize', 
    1718               ], 
    1819    package_data={'NeuroTools': ['doc/*.txt', 'README']}, 

  • trunk/src/analysis.py

    r321 r451  
    9292            return k 
    9393         
     94def make_kernel(form, sigma, time_stamp_resolution, direction=1, kernel=None, 
     95                m_idx=None): 
     96    """ 
     97    Constructs a numeric linear convolution kernel of basic shape to be used 
     98    for data smoothing (linear low pass filtering) and firing rate estimation 
     99    from single trial or trial-averaged spike trains. 
     100     
     101    Exponential and alpha kernels may also be used to represent postynaptic 
     102    currents / potentials in a linear (current-based) model. 
     103     
     104    Adapted from original script written by Martin P. Nawrot for the FIND MATLAB 
     105    toolbox. 
     106    See FIND - a unified framework for neural data analysis, 
     107        Meier R, Egert U, Aertsen A, Nawrot MP; Neural Netw. 2008 Oct; 
     108        21(8):1085-93.  
     109         
     110        Nawrot M, Aertsen A, Rotter S (1999) 
     111        Single-trial estimation of neuronal firing rates - From single neuron 
     112        spike trains to population activity. 
     113        J Neurosci Meth 94: 81-92 
     114         
     115    Inputs: 
     116        form                  - kernel form (string) Currently implemented forms 
     117                                are BOX - boxcar, TRI - triangle, GAU - gaussian 
     118                                , EPA - epanechnikov, EXP - exponential, ALP - 
     119                                alpha function. 
     120                                EXP and ALP are aymmetric kernel forms and 
     121                                assume optional parameter 'direction' 
     122        sigma                 - standard deviation of the distribution 
     123                                associated with kernel shape. This parameter 
     124                                defines the time resolution of the kernel 
     125                                estimate and makes different kernels comparable 
     126                                (cf. [1] for symetric kernels). This is used 
     127                                here as an alternative definition to the cut-off 
     128                                frequency of the associated linear filter. 
     129        time_stamp_resolution - temporal resolution of input and output in SI 
     130                                units. E.g. value 0.001 means 1ms time 
     131                                resolution 
     132        direction             - Asymmetric kernels have two possible directions. 
     133                                The values are -1 or 1, default is 1. The 
     134                                definition here is that for direction = 1 the 
     135                                kernel represents the impulse response function 
     136                                of the linear filter. Default value is 1. 
     137     
     138    Outputs: 
     139        kernel  - array of kernel. The length of this array is always an odd 
     140                  number to represent symmetric kernels such that the center bin 
     141                  coincides with the median of the numeric array, i.e for a 
     142                  triangle, the maximum will be at the center bin with equal 
     143                  number of bins to the right and to the left. 
     144        norm    - for rate estimates. The kernel vector is normalized such that 
     145                  the sum of all entries equals unity sum(kernel)=1. When 
     146                  estimating rate functions from discrete spike data (0/1) the 
     147                  additional parameter 'norm' allows for the normalization to 
     148                  rate in spikes per second. 
     149                  Use: rate = norm * scipy.signal.lfilter(kernel, 1, spike_data) 
     150        m_idx   - index of the numerically determined median (center of gravity) 
     151                  of the kernel function 
     152    """ 
     153    assert form.upper() in ('BOX','TRI','GAU','EPA','EXP','ALP'), "form must \ 
     154    be one either 'BOX','TRI','GAU','EPA','EXP' or 'ALP'!" 
     155     
     156    assert direction in (1,-1), "direction must be either 1 or -1" 
     157     
     158    norm = 1./time_stamp_resolution 
     159 
     160    if form.upper() == 'BOX': 
     161        w = 2.0 * sigma * numpy.sqrt(3) 
     162        width = 2 * numpy.floor(w / 2.0 / time_stamp_resolution) + 1  # always odd number of bins 
     163        height = 1. / width 
     164        kernel = numpy.ones((1, width)) * height  # area = 1 
     165         
     166    elif form.upper() == 'TRI': 
     167        w = 2 * sigma * numpy.sqrt(6) 
     168        halfwidth = numpy.floor(w / 2.0 / time_stamp_resolution) 
     169        trileft = numpy.arange(1, halfwidth + 2) 
     170        triright = numpy.arange(halfwidth, 0, -1)  # odd number of bins 
     171        triangle = numpy.append(trileft, triright) 
     172        kernel = triangle / triangle.sum()  # area = 1 
     173         
     174    elif form.upper() == 'EPA': 
     175        w = 2.0 * sigma * numpy.sqrt(5) 
     176        halfwidth = numpy.floor(w / 2.0 / time_stamp_resolution) 
     177        base = numpy.arange(-halfwidth, halfwidth + 1) 
     178        parabula = base**2 
     179        epanech = parabula.max() - parabula  # inverse parabula 
     180        kernel = epanech / epanech.sum()  # area = 1 
     181         
     182    elif form.upper() == 'GAU': 
     183        SI_sigma = sigma / 1000.0 
     184        w = 2.0 * sigma * 2.7  # > 99% of distribution weight 
     185        halfwidth = numpy.floor(w / 2.0 / time_stamp_resolution)  # always odd  
     186        base = numpy.arange(-halfwidth, halfwidth + 1) / 1000.0 * ( 
     187            time_stamp_resolution) 
     188        g = numpy.exp(-(base**2) / 2.0 / SI_sigma**2) / SI_sigma / numpy.sqrt( 
     189            2.0 * numpy.pi) 
     190        kernel = g / g.sum() 
     191         
     192    elif form.upper() == 'ALP': 
     193        SI_sigma = sigma / 1000.0 
     194        w = 5.0 * sigma 
     195        alpha = numpy.arange(1, (2.0 * numpy.floor( 
     196            (w / time_stamp_resolution / 2.0)) + 1) + 1) / 1000.0 * \ 
     197        time_stamp_resolution 
     198        alpha = (2.0 / SI_sigma**2) * alpha * numpy.exp(-alpha * numpy.sqrt(2) \ 
     199                                                        / SI_sigma) 
     200        kernel = alpha / alpha.sum()  # normalization 
     201        if direction == -1: 
     202            kernel = numpy.flipud(kernel) 
     203             
     204    elif form.upper() == 'EXP': 
     205        SI_sigma = sigma / 1000.0 
     206        w = 5.0 * sigma 
     207        expo = numpy.arange(1, (2.0 * numpy.floor(w / time_stamp_resolution / ( 
     208            2.0)) + 1) + 1) / 1000.0 * time_stamp_resolution 
     209        expo = numpy.exp(-expo / SI_sigma) 
     210        kernel = expo / expo.sum() 
     211        if direction == -1: 
     212            kernel = numpy.flipud(kernel) 
     213     
     214    kernel = kernel.ravel() 
     215    m_idx = numpy.nonzero(kernel.cumsum() >= 0.5)[0].min() 
     216     
     217    return kernel, norm, m_idx 
     218 
    94219def simple_frequency_spectrum(x): 
    95220    """ 

  • trunk/test/test_analysis.py

    r321 r451  
    22Unit tests for the NeuroTools.analysis module 
    33""" 
     4import numpy 
     5import scipy.io 
     6import unittest 
     7from numpy import pi, sin 
    48 
    59from NeuroTools import analysis 
    6 import numpy, unittest 
    7 from numpy import pi, sin 
     10 
    811 
    912# test simple_frequency_spectrum 
    … …  
    1114 
    1215class AnalysisTest(unittest.TestCase): 
     16 
     17    def setUp(self): 
     18        #The following are used for the make_kernel testcases 
     19        self.box = scipy.io.loadmat('analysis/make_kernel/box.mat') 
     20        self.tri = scipy.io.loadmat('analysis/make_kernel/tri.mat') 
     21        self.epa = scipy.io.loadmat('analysis/make_kernel/epa.mat') 
     22        self.gau = scipy.io.loadmat('analysis/make_kernel/gau.mat') 
     23        self.alp = scipy.io.loadmat('analysis/make_kernel/alp.mat') 
     24        self.exp = scipy.io.loadmat('analysis/make_kernel/exp.mat') 
     25        self.alp_reversed = scipy.io.loadmat( 
     26            'analysis/make_kernel/alp_reversed.mat') 
     27        self.exp_reversed = scipy.io.loadmat( 
     28            'analysis/make_kernel/exp_reversed.mat') 
    1329 
    1430    def testSimpleFrequencySpectrum(self): 
    … …  
    3652        z=analysis.ccf(a,a) 
    3753        assert z[len(z)/2] == 1 
     54     
     55    def testMakeKernelBox(self): 
     56        true_kernel = self.box['kernel'].ravel() 
     57        true_norm = self.box['norm'].ravel()[0] 
     58        true_m_idx = self.box['m_idx'].ravel()[0] - 1 
     59         
     60        kernel, norm, m_idx = analysis.make_kernel('BOX', 1, 0.001) 
     61         
     62        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     63                                                ) 
     64        self.assertEqual(true_norm, norm, 3) 
     65        self.assertEqual(true_m_idx, m_idx, 3) 
     66         
     67    def testMakeKernelTri(self): 
     68        true_kernel = self.tri['kernel'].ravel() 
     69        true_norm = self.tri['norm'].ravel()[0] 
     70        true_m_idx = self.tri['m_idx'].ravel()[0] - 1 
     71         
     72        kernel, norm, m_idx = analysis.make_kernel('TRI', 1, 0.001) 
     73         
     74        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     75                                                ) 
     76        self.assertEqual(true_norm, norm, 3) 
     77        self.assertEqual(true_m_idx, m_idx, 3) 
     78         
     79    def testMakeKernelEpa(self): 
     80        true_kernel = self.epa['kernel'].ravel() 
     81        true_norm = self.epa['norm'].ravel()[0] 
     82        true_m_idx = self.epa['m_idx'].ravel()[0] - 1 
     83         
     84        kernel, norm, m_idx = analysis.make_kernel('EPA', 1, 0.001) 
     85         
     86        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     87                                                ) 
     88        self.assertEqual(true_norm, norm, 3) 
     89        self.assertEqual(true_m_idx, m_idx, 3) 
     90         
     91    def testMakeKernelGau(self): 
     92        true_kernel = self.gau['kernel'].ravel() 
     93        true_norm = self.gau['norm'].ravel()[0] 
     94        true_m_idx = self.gau['m_idx'].ravel()[0] - 1 
     95         
     96        kernel, norm, m_idx = analysis.make_kernel('GAU', 1, 0.001) 
     97         
     98        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     99                                                ) 
     100        self.assertEqual(true_norm, norm, 3) 
     101        self.assertEqual(true_m_idx, m_idx, 3) 
     102         
     103    def testMakeKernelAlp(self): 
     104        true_kernel = self.alp['kernel'].ravel() 
     105        true_norm = self.alp['norm'].ravel()[0] 
     106        true_m_idx = self.alp['m_idx'].ravel()[0] - 1 
     107         
     108        kernel, norm, m_idx = analysis.make_kernel('ALP', 1, 0.001) 
     109         
     110        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     111                                                ) 
     112        self.assertEqual(true_norm, norm, 3) 
     113        self.assertEqual(true_m_idx, m_idx, 3) 
    38114 
     115    def testMakeKernelExp(self): 
     116        true_kernel = self.exp['kernel'].ravel() 
     117        true_norm = self.exp['norm'].ravel()[0] 
     118        true_m_idx = self.exp['m_idx'].ravel()[0] - 1 
     119         
     120        kernel, norm, m_idx = analysis.make_kernel('EXP', 1, 0.001) 
     121         
     122        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     123                                                ) 
     124        self.assertEqual(true_norm, norm, 3) 
     125        self.assertEqual(true_m_idx, m_idx, 3) 
     126     
     127    def testMakeKernelAlpReversed(self): 
     128        """Same as testAlp but with direction = -1 
     129        """ 
     130        true_kernel = self.alp_reversed['kernel'].ravel() 
     131        true_norm = self.alp_reversed['norm'].ravel()[0] 
     132        true_m_idx = self.alp_reversed['m_idx'].ravel()[0] - 1 
     133         
     134        kernel, norm, m_idx = analysis.make_kernel('ALP', 1, 0.001, 
     135                                                  direction = -1) 
     136        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     137                                                ) 
     138        self.assertEqual(true_norm, norm, 3) 
     139        self.assertEqual(true_m_idx, m_idx, 3) 
     140         
     141    def testMakeKernelExpReversed(self): 
     142        """Same as testExp but with direction = -1 
     143        """ 
     144        true_kernel = self.exp_reversed['kernel'].ravel() 
     145        true_norm = self.exp_reversed['norm'].ravel()[0] 
     146        true_m_idx = self.exp_reversed['m_idx'].ravel()[0] - 1 
     147         
     148        kernel, norm, m_idx = analysis.make_kernel('EXP', 1, 0.001, direction = -1) 
     149         
     150        numpy.testing.assert_array_almost_equal(true_kernel, kernel, decimal = 3 
     151                                                ) 
     152        self.assertEqual(true_norm, norm, 3) 
     153        self.assertEqual(true_m_idx, m_idx, 3) 
     154         
    39155if __name__ == "__main__": 
    40156    unittest.main()