Changeset 315

Timestamp:

11/10/08 19:17:08 (2 months ago)

Author:

LaurentPerrinet

Message:

Unit Test of signals.py was not working / fixed the computation of CV_KL / some clean-up on the example with a simple neuron

Files:

• trunk/examples/single_neuron/CRF_neuron_vs_signal.py (modified) (6 diffs)
• trunk/examples/single_neuron/fiber.param (modified) (1 diff)
• trunk/examples/single_neuron/simple_single_neuron.py (modified) (7 diffs)
• trunk/src/signals.py (modified) (7 diffs)
• trunk/test/test_signals.py (modified) (4 diffs)

trunk/examples/single_neuron/CRF_neuron_vs_signal.py

@@ -23 +23 @@
-
-# this is not mandatory but just a "easy_install progressbar" away
-
+3 
-import progressbar # see http://projects.scipy.org/pipermail/scipy-dev/2008-January/008200.html
-
@@ -30 +30 @@
-N_exp_noise = 10
-N_exp_snr = 25
-
-
+
+
+
+
+
-
-name = sys.argv[0].split('.')[0] # name of the current script withpout the '.py' part
@@ -40 +43 @@
-
-    # calculates the dimension of the parameter space
-
+s
-
-    # creates results array with size of parameter space dimension
@@ -49 +52 @@
-    pbar=progressbar.ProgressBar(widgets=[name, " ", progressbar.Percentage(), ' ',
-            progressbar.Bar(), ' ', progressbar.ETA()], maxval=numpy.prod(results_dim))
-
+s
-        params = myFibers.params
-        params.update(experiment) # updates what changed in the dictionary
@@ -55 +58 @@
-        data = myFibers.run(params,verbose=False)
-        # calculating the index in the parameter space
-
+s
-        # put the data at the right position in the results array
-        CRF[index] = data.mean_rate()#
@@ -65 +68 @@
-results.close()
-
-#print CRFsr
-#numpy.array(p.noise_std._values),numpy.array(p.snr._values),
-pylab.plot(CRF.transpose()) #color = (sin(2*pi*noise_list)**2,cos(2*pi*noise_list)**2,1))

trunk/examples/single_neuron/fiber.param

@@ -7 +7 @@
-    "connectseed": 12345789,
-  },
-1
-8
+2
+2
-  "weight": 1.0,
-5
+10
-}

trunk/examples/single_neuron/simple_single_neuron.py

@@ -21 +21 @@
-import sys, os, tempfile
-# choosing the simulator
-brian
+nest2
-# the link to read SpikeList files with NeuroTools
-from NeuroTools.signals import load_spikelist
@@ -27 +27 @@
-from NeuroTools.parameters import ParameterSet
-
-def define(N=100):
-    """
-    defines the parameter set the first time
-
-    """
-    #these are fixed
-    simulation_params = ParameterSet({'dt'        : 0.1,# discretization step in simulations (ms)
-                'simtime'   : 40000*0.1,      # float; (ms)
-                'syn_delay' : 1.0,         # float; (ms)
-                'kernelseed' : 4321097,       # array with one element per thread
-                'connectseed' : 12345789  # seed for random generator(s) used during simulation
-                })
-    # these may change
-    P = ParameterSet({'simulation': simulation_params,
-                    'N'         : N,
-                    'noise_std' : 8.0,       # (nA??) standard deviation of the internal noise
-                    'snr'       : 1.0,      # (nA??) size of the input signal
-                    'weight'    : 1.0 },
-                    label="fiber_params")
-
-    P.save('file://' + os.getcwd() + '/fiber.param')
-
-class FiberChannel(object):
@@ -55 +34 @@
-    """
-    def __init__(self,N=100):
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-        # keep all parameters just change the size
-        self.params.N = N
-        print self.params.pretty()
-
@@ -79 +62 @@
-        N = params.N
-        #dc_generator
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-        # target population
-        output = sim.Population(N , sim.IF_cond_exp)
@@ -101 +80 @@
-        sim.Projection(noise, output, conn)
-
-#
-#
+
+
-
-        output.record()
@@ -118 +97 @@
-        Timer.reset()  # start timer on construction
-
-=
-
+ = 
+#
-        output.printSpikes(output_filename)#
-        output_DATA = load_spikelist(output_filename,N,
@@ -138 +117 @@
-        return output_DATA
-
-
-
-if __name__ == '__main__':
-
-
+
-    spikes = myFibers.run(myFibers.params)
-    spikes.raster_plot()

trunk/src/signals.py

@@ -262 +262 @@
-    def mean_rate(self, t_start=None, t_stop=None):
-        """ 
-
+s
-        
-        Inputs:
@@ -302 +302 @@
-        
-        See also
-
+
+
-        """
-        isi = self.isi()
@@ -310 +311 @@
-            logging.debug("Warning, a CV can't be computed because there are not enough spikes")
-            return numpy.nan
+
+
+
+    def cv_kl(self, bins=100):
+        """
+        Provides a measure for the coefficient of variation to describe the
+        regularity in spiking networks. It is based on the Kullback-Leibler
+        divergence and decribes the difference between a given
+        interspike-interval-distribution and an exponential one (representing
+        poissonian spike trains) with equal mean.
+        It yields 1 for poissonian spike trains and 0 for regular ones.
+        
+        Reference:
+            http://incm.cnrs-mrs.fr/LaurentPerrinet/Publications/Voges08fens
+        
+        
+        Examples:
+            >> spklist.cv_kl()
+                0.98
+        
+        See also:
+            cv_isi
+            
+        """
+        isi = self.isi() / 1000.
+        if newnum:
+            proba_isi, xaxis = numpy.histogram(isi, bins=bins, normed=True, new=True)
+            #xaxis = xaxis[:len(xaxis)-1]
+        else:
+            proba_isi, xaxis = numpy.histogram(isi, bins=bins, normed=True)
+
+        proba_isi /= numpy.sum(proba_isi)
+        bin_size = xaxis[1]-xaxis[0]
+
+        # differential entropy: http://en.wikipedia.org/wiki/Differential_entropy
+        KL = - numpy.sum(proba_isi * numpy.log(proba_isi+1e-16)) + numpy.log(bin_size)
+        KL -= -numpy.log(self.mean_rate()) + 1.
+        CVkl=numpy.exp(-KL)
+        
+        return CVkl
+    
-
-    def fano_factor_isi(self):
@@ -1124 +1166 @@
-            pylab.draw()
-
-
+    
-    def cv_isi(self, float_only=False):
-        """
-cv coefficients for all the SpikeTrains objects
+CV coefficients for each SpikeTrains object
-        within the SpikeList. Return NaN when not enough spikes are present
-        
@@ -1141 +1183 @@
-
-        See also:
-
-
-
+
+
+
+
-        for id in self.id_list():
-.append(self.spiketrains[id].cv_isi()
-
+[id] = self.spiketrains[id].cv_isi(
+
-        if float_only:
-            cvs_isi = numpy.extract(numpy.logical_not(numpy.isnan(cvs_isi)),cvs_isi)
@@ -1152 +1195 @@
-
-
+    def cv_kl(self, bins = 50):
+        """
+        Return the list of all the CV coefficients for each SpikeTrains object
+        within the SpikeList.
+        
+        See also:
+            cv_isi_hist, cv_local, cv_isi, SpikeTrain.cv_kl
+            
+        """
+        cvs_kl = numpy.empty(len(self.id_list()))
+        for id in self.id_list():
+            cvs_kl[id] = self.spiketrains[id].cv_kl(bins = bins)
+        return cvs_kl
-
-    def cv_isi_hist(self, bins=50, display=False, kwargs={}):
@@ -1188 +1244 @@
-
-
-
-    def cv_kl(self):
-        """
-        Provides another (new?) measure for the coefficient of variation to describe the
-        regularity in spiking networks. It is based on the Kullback-Leibler
-        divergence and decribes the difference between a given
-        interspike-interval-distribution and an exponential one (representing
-        poissonian spike trains) with equal mean.
-        It yields approx. 1 for poissonian spike trains and 0 for regular ones.
-        
-        Attention: Contrary to the normal cv which is calculated for each
-        single neuron and then averaged, this function needs to synchronously
-        consider the spike times of all neurons.
-        
-        Examples:
-            >> spklist.cv_kl()
-                0.98
-        
-        See also:
-            cv_isi, cv_isi_hist, cv_local
-        """
-        
-        isis=numpy.concatenate(self.isi())
-        mISI=numpy.mean(isis)
-        norm=len(isis)
-        min=numpy.min(isis)
-        max=numpy.max(isis)
-        bins = numpy.arange(min,max,1)
-        #print 'CV-KL   norm=',norm,' #bins=',len(bins),' ev=',mISI
-        hist = numpy.histogram(isis,bins)
-        dist = hist[0].astype(float)/norm
-        testsum = numpy.sum(dist)
-        if testsum>1.1 or testsum<0.9 :
-            print 'CV-KL   testsum=',testsum
-        H=0.
-        for j in range(len(dist)):
-            if dist[j] > 0. :
-                H = H - (dist[j]*numpy.log(dist[j]))
-        KL=0.
-        KL=-H+numpy.log(mISI)+1.
-        CVkl=numpy.exp(-KL)
-        
-        return CVkl
-
-

trunk/test/test_signals.py

@@ -133 +133 @@
-        assert 0.9 < spk1.cv_isi() < 1.1
-        self.assertEqual(spk2.cv_isi(), 0)
+
+
+    def testCvKL(self):
+        poisson_param = 1./10 # 1 / firing_frequency
+        isi           = numpy.random.exponential(poisson_param, 10000)
+        poisson_times = numpy.cumsum(isi)*1000. # To convert the spikes_time in ms
+        spk1 = signals.SpikeTrain(poisson_times)
+        assert 0.9 < spk1.cv_kl(bins = 1000) < 1.1
+        # does not depend on bin size
+        assert 0.9 < spk1.cv_kl(bins = 100) < 1.1
+        # does not depend on time
+        poisson_param = 1./4
+        isi           = numpy.random.exponential(poisson_param, 10000)
+        poisson_times = numpy.cumsum(isi)*1000. # To convert the spikes_time in ms
+        spk1 = signals.SpikeTrain(poisson_times)
+        assert 0.9 < spk1.cv_kl() < 1.1
+        spk2 = signals.SpikeTrain(range(10), t_stop=10)
+        self.assertEqual(spk2.cv_isi(), 0)
-    
-    def testHistogram(self):
-        poisson_param = 1./40
-
+0
-        poisson_times = numpy.cumsum(isi)*1000. # To convert the spikes_time in ms
-        spk = signals.SpikeTrain(poisson_times)
@@ -191 +209 @@
-        for idx in xrange(nb_cells):
-            param   = 1./frequencies[idx]
-
+0
-            pspikes = numpy.cumsum(isi)*1000. # To convert the spikes_time in ms
-            for spike in pspikes: 
@@ -279 +297 @@
-    
-    def testCVKL(self):
-self.spk.cv_kl(
+numpy.mean(self.spk.cv_kl()
-    
-    def testCVLocal(self):
@@ -370 +388 @@
-        cc1 = self.spk.pairwise_cc_zero([(i,i) for i in xrange(5)], time_bin=0.1)
-        cc2 = self.spk.pairwise_cc_zero(5, time_bin=0.1)
+        #print cc1, cc2
+
-        assert (0 <= cc1 <= 1) and (cc1 > cc2)
-