Changeset 323

Timestamp:

11/12/08 17:00:08 (2 months ago)

Author:

emuller

Message:

Added adapting markov process (1d and 2d) to stgen

Files:

• trunk/src/stgen.py (modified) (2 diffs)
• trunk/test/test_stgen.py (modified) (1 diff)

trunk/src/stgen.py

@@ -21 +21 @@
-"""
-
-
-# TODO : needs refactoring
-# TODO: needs python gsl wrapper / convert to numpy
-# TODO make it generate spiketrains?
-
-from NeuroTools import check_dependency
@@ -392 +388 @@
-    else:
-        inh_gamma_generator = _inh_gamma_generator_unimplemented
+
+
+
+    def _inh_adaptingmarkov_generator_python(self, a, bq, tau, t, t_stop, array=False):
+
+        """
+        Returns a SpikeList whose spikes are an inhomogeneous
+        realization (dynamic rate) of the so-called adapting markov
+        process (see references). The implementation uses the thinning
+        method, as presented in the references.
+
+        This is the 1d implementation, with no relative refractoriness.
+        For the 2d implementation with relative refractoriness, 
+        see the inh_2dadaptingmarkov_generator.
+
+        Inputs:
+            a,bq    - arrays of the parameters of the hazard function where a[i] and bq[i] 
+                     will be active on interval [t[i],t[i+1]]
+            tau    - the time constant of adaptation (in milliseconds).
+            t      - an array specifying the time bins (in milliseconds) at which to 
+                     specify the rate
+            t_stop - length of time to simulate process (in ms)
+            array  - if True, a numpy array of sorted spikes is returned,
+                     rather than a SpikeList object.
+
+        Note: 
+            - t_start=t[0]
+
+            - a is in units of Hz.  Typical values are available 
+              in Fig. 1 of Muller et al 2007, a~5-80Hz (low to high stimulus)
+
+            - bq here is taken to be the quantity b*q_s in Muller et al 2007, is thus
+              dimensionless, and has typical values bq~3.0-1.0 (low to high stimulus)
+
+            - tau_s has typical values on the order of 100 ms
+
+
+        References:
+
+        Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier 
+        Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories
+        Neural Comput. 2007 19: 2958-3010.
+        
+        Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
+
+        Examples:
+        
+        See also:
+            inh_poisson_generator, inh_gamma_generator, inh_2dadaptingmarkov_generator
+
+        """
+
+        from numpy import shape
+
+        if shape(t)!=shape(a) or shape(a)!=shape(bq):
+            raise ValueError('shape mismatch: t,a,b must be of the same shape')
+
+        # get max rate and generate poisson process to be thinned
+        rmax = numpy.max(a)
+        ps = self.poisson_generator(rmax, t_start=t[0], t_stop=t_stop, array=True)
+
+        isi = numpy.zeros_like(ps)
+        isi[1:] = ps[1:]-ps[:-1]
+        isi[0] = ps[0] #-0.0 # assume spike at 0.0
+
+        # return empty if no spikes
+        if len(ps) == 0:
+            return SpikeTrain(numpy.array([]), t_start=t[0],t_stop=t_stop)
+
+
+        # gen uniform rand on 0,1 for each spike
+        rn = numpy.array(self.rng.uniform(0, 1, len(ps)))
+
+        # instantaneous a,bq for each spike
+
+        idx=numpy.searchsorted(t,ps)-1
+        spike_a = a[idx]
+        spike_bq = bq[idx]
+
+        keep = numpy.zeros(shape(ps),bool)
+
+        # thin spikes
+
+        i = 0
+        t_last = 0.0
+        t_i = 0
+        # initial adaptation state is unadapted, i.e. large t_s
+        t_s = 1000*tau
+
+        while(i<len(ps)):
+            # find index in "t" time, without searching whole array each time
+            t_i = numpy.searchsorted(t[t_i:],ps[i],'right')-1+t_i
+
+            # evolve adaptation state
+            t_s+=isi[i]
+
+            if rn[i]<a[t_i]*numpy.exp(-bq[t_i]*numpy.exp(-t_s/tau))/rmax:
+                # keep spike
+                keep[i] = True
+                # remap t_s state
+                t_s = -tau*numpy.log(numpy.exp(-t_s/tau)+1)
+            i+=1
+
+
+        spike_train = ps[keep]
+
+        if array:
+            return spike_train
+
+        return SpikeTrain(spike_train, t_start=t[0],t_stop=t_stop)
+
+
+    # use slow python implementation for the time being
+    # TODO: provide optimized C/weave implementation if possible
+
+        
+    inh_adaptingmarkov_generator = _inh_adaptingmarkov_generator_python
+
+
+    def _inh_2Dadaptingmarkov_generator_python(self, a, bq, tau_s, tau_r, qrqs, t, t_stop, array=False):
+
+        """
+        Returns a SpikeList whose spikes are an inhomogeneous
+        realization (dynamic rate) of the so-called 2D adapting markov
+        process (see references).  2D implies the process has two
+        states, an adaptation state, and a refractory state, both of
+        which affect its probability to spike.  The implementation
+        uses the thinning method, as presented in the references.
+
+        For the 1d implementation, with no relative refractoriness,
+        see the inh_adaptingmarkov_generator.
+
+        Inputs:
+            a,bq    - arrays of the parameters of the hazard function where a[i] and bq[i] 
+                     will be active on interval [t[i],t[i+1]]
+            tau_s    - the time constant of adaptation (in milliseconds).
+            tau_r    - the time constant of refractoriness (in milliseconds).
+            qrqs     - the ratio of refractoriness conductance to adaptation conductance.
+                       typically on the order of 200.
+            t      - an array specifying the time bins (in milliseconds) at which to 
+                     specify the rate
+            t_stop - length of time to simulate process (in ms)
+            array  - if True, a numpy array of sorted spikes is returned,
+                     rather than a SpikeList object.
+
+        Note: 
+            - t_start=t[0]
+
+            - a is in units of Hz.  Typical values are available 
+              in Fig. 1 of Muller et al 2007, a~5-80Hz (low to high stimulus)
+
+            - bq here is taken to be the quantity b*q_s in Muller et al 2007, is thus
+              dimensionless, and has typical values bq~3.0-1.0 (low to high stimulus)
+
+            - qrqs is the quantity q_r/q_s in Muller et al 2007, 
+              where a value of qrqs = 3124.0nS/14.48nS = 221.96 was used.
+
+            - tau_s has typical values on the order of 100 ms
+            - tau_r has typical values on the order of 2 ms
+
+
+        References:
+
+        Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier 
+        Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories
+        Neural Comput. 2007 19: 2958-3010.
+        
+        Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
+
+        Examples:
+        
+        See also:
+            inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
+
+        """
+
+        from numpy import shape
+
+        if shape(t)!=shape(a) or shape(a)!=shape(bq):
+            raise ValueError('shape mismatch: t,a,b must be of the same shape')
+
+        # get max rate and generate poisson process to be thinned
+        rmax = numpy.max(a)
+        ps = self.poisson_generator(rmax, t_start=t[0], t_stop=t_stop, array=True)
+
+        isi = numpy.zeros_like(ps)
+        isi[1:] = ps[1:]-ps[:-1]
+        isi[0] = ps[0] #-0.0 # assume spike at 0.0
+
+        # return empty if no spikes
+        if len(ps) == 0:
+            return SpikeTrain(numpy.array([]), t_start=t[0],t_stop=t_stop)
+
+
+        # gen uniform rand on 0,1 for each spike
+        rn = numpy.array(self.rng.uniform(0, 1, len(ps)))
+
+        # instantaneous a,bq for each spike
+
+        idx=numpy.searchsorted(t,ps)-1
+        spike_a = a[idx]
+        spike_bq = bq[idx]
+
+        keep = numpy.zeros(shape(ps),bool)
+
+        # thin spikes
+
+        i = 0
+        t_last = 0.0
+        t_i = 0
+        # initial adaptation state is unadapted, i.e. large t_s
+        t_s = 1000*tau_s
+        t_r = 1000*tau_s
+
+        while(i<len(ps)):
+            # find index in "t" time, without searching whole array each time
+            t_i = numpy.searchsorted(t[t_i:],ps[i],'right')-1+t_i
+
+            # evolve adaptation state
+            t_s+=isi[i]
+            t_r+=isi[i]
+
+            if rn[i]<a[t_i]*numpy.exp(-bq[t_i]*(numpy.exp(-t_s/tau_s)+qrqs*numpy.exp(-t_r/tau_r)))/rmax:
+                # keep spike
+                keep[i] = True
+                # remap t_s state
+                t_s = -tau_s*numpy.log(numpy.exp(-t_s/tau_s)+1)
+                t_r = -tau_r*numpy.log(numpy.exp(-t_r/tau_r)+1)
+            i+=1
+
+
+        spike_train = ps[keep]
+
+        if array:
+            return spike_train
+
+        return SpikeTrain(spike_train, t_start=t[0],t_stop=t_stop)
+
+
+    # use slow python implementation for the time being
+    # TODO: provide optimized C/weave implementation if possible
+
+        
+    inh_2Dadaptingmarkov_generator = _inh_2Dadaptingmarkov_generator_python
+
+
+
+
+
-
-

trunk/test/test_stgen.py

@@ -227 +227 @@
-
-
+    def testInhAdaptingMarkov(self):
+
+
+        stg = stgen.StGen()
+
+        from numpy import array
+
+        t = array([0.0,5000.0])
+        # using parameters from paper
+        a = array([23.18,47.24])
+        bq = array([0.10912,0.09794])*14.48
+
+        # expected mean firing rates for 2D case
+        alpha = array([7.60, 10.66])
+
+        tau = 110.0
+        t_stop = 10000.0
+
+        st = stg.inh_adaptingmarkov_generator(a,bq,tau,t,t_stop)
+        assert isinstance(st,signals.SpikeTrain)
+
+        st = stg.inh_adaptingmarkov_generator(a,bq,tau,t,t_stop,array=True)
+        assert isinstance(st, numpy.ndarray)
+
+        assert st[-1] < t_stop
+        assert st[0] > t[0]
+
+
+    def testInh2DAdaptingMarkov(self):
+
+
+        stg = stgen.StGen()
+
+        from numpy import array
+
+        t = array([0.0,10000.0])
+        # using parameters from paper
+        a = array([23.18,47.24])
+        bq = array([0.10912,0.09794])*14.48
+
+        # expected mean firing rates for 2D case
+        alpha = array([7.60, 10.66])
+
+        tau_s = 110.0
+        tau_r = 1.97
+        qrqs = 221.96
+        t_stop = 20000.0
+
+        st = stg.inh_2Dadaptingmarkov_generator(a,bq,tau_s,tau_r,qrqs,t,t_stop)
+        assert isinstance(st,signals.SpikeTrain)
+
+        st = stg.inh_2Dadaptingmarkov_generator(a,bq,tau_s,tau_r,qrqs,t,t_stop,array=True)
+        assert isinstance(st, numpy.ndarray)
+
+        assert st[-1] < t_stop
+        assert st[0] > t[0]
+
+        spikes1 = len(st[st<10000])
+        spikes2 = len(st[st>10000])
+        
+        # should be approximately alpha[0]*10.0 spikes
+        assert numpy.clip(spikes1,60,100)==spikes1
+
+        # should be approximately alpha[1]*10.0 spikes
+        assert numpy.clip(spikes2,80,140)==spikes2
+
+
-
-