Changeset 307 for trunk/src/stgen.py

Timestamp:

11/08/08 10:58:36 (4 years ago)

Author:

emuller

Message:

Fixed inh_gamma_generator hazard function numerics
Fixed poisson_process ISI buffer overrun
Added stgen module docstring
Added more tests
Added inh_gamma_psth.py example

Files:

• trunk/src/stgen.py (modified) (13 diffs)

trunk/src/stgen.py

@@ +1 @@
+"""
+NeuroTools.stgen
+================
+
+A collection of tools for stochastic process generation.
+
+
+Classes
+-------
+
+StGen - Object to generate stochastic processes of various kinds
+        and return them as SpikeTrain or AnalogSignal objects.
+
+
+Functions
+---------
+
+shotnoise_fromspikes - Convolves the provided spike train with shot decaying exponential.
+
+gamma_hazard - Compute the hazard function for a gamma process with parameters a,b.
+"""
+
+
 # TODO : needs refactoring
 # TODO: needs python gsl wrapper / convert to numpy
@@ -9 +32 @@
 
 
-def gamma_hazard(x, a, b, dt=1e-4):
+def gamma_hazard_scipy(x, a, b, dt=1e-4):
     """
     Compute the hazard function for a gamma process with parameters a,b
@@ -21 +44 @@
     """
 
-    # TODO: this algorithm has numerical problems
+    # This algorithm is presently not used by
+    # inh_gamma_generator as it has numerical problems
     # Try: 
     # plot(stgen.gamma_hazard(arange(0,1000.0,0.1),10.0,1.0/50.0))
@@ -39 +63 @@
     else:
         return val
+
+
+def gamma_hazard(x, a, b, dt=1e-4):
+    """
+    Compute the hazard function for a gamma process with parameters a,b
+    where a and b are the parameters of the gamma PDF:
+    y(t) = x^(a-1) \exp(-x/b) / (\Gamma(a)*b^a)
+
+    Inputs:
+        x   - in units of seconds
+        a   - dimensionless
+        b   - in units of seconds
+    """
+
+    # Used by inh_gamma_generator
+
+    # Ideally, I would like to see an implementation which does not depend on RPy
+    # but the gamma_hazard_scipy above using scipy exhibits numerical problems, as it does not
+    # support directly returning the log.
+
+    if not check_dependency('rpy'):
+        raise ImportError("gamma_hazard requires RPy (http://rpy.sourceforge.net/)")
+
+    from rpy import r
+
+    # scipy.special.gammaincc has numerical problems
+    #Hpre = -log(scipy.special.gammaincc(a,(x-dt)/b))
+    #Hpost = -log(scipy.special.gammaincc(a,(x+dt)/b))
+
+    # reverting to the good old r.pgamma
+    Hpre = -r.pgamma(x-dt,shape=a,scale=b,lower=False,log=True)
+    Hpost = -r.pgamma(x+dt,shape=a,scale=b,lower=False,log=True)
+    val =  0.5*(Hpost-Hpre)/dt
+
+    return val
+
+
     
 
@@ -45 +106 @@
     def __init__(self, rng=None, seed=None):
         """ 
-        Spike Train Generator
-        Object to generate spiking random processes with various statistics 
-        (inhomogeneous poisson, inhomogeneous gamma, etc.) as SpikeTrain objects.
+        Stochastic Process Generator
+        ============================
+
+        Object to generate stochastic processes of various kinds
+        and return them as SpikeTrain or AnalogSignal objects.
+      
 
         Inputs:
-            rng - Seed for the random number generator. Can be None, or 
+            rng - The random number generator state object (optional). Can be None, or 
                   a numpy.random.RandomState object, or an object with the same 
                   interface.
+
+            seed - A seed for the rng (optional).
 
         If rng is not None, the provided rng will be used to generate random numbers, 
@@ -60 +126 @@
         Examples:
             >> x = StGen()
+
+
+
+        StGen Methods:
+
+        Spiking point processes:
+        ------------------------
+ 
+        poisson_generator - homogeneous Poisson process
+        inh_poisson_generator - inhomogeneous Poisson process (time varying rate)
+        inh_gamma_generator - inhomogeneous Gamma process (time varying a,b)
+
+        Continuous time processes:
+        --------------------------
+
+        OU_generator - Ohrnstein-Uhlenbeck process
+        
+
+        See also:
+          shotnoise_fromspikes
+
         """
 
@@ -75 +162 @@
 
 
-    def poisson_generator(self, rate, t_start=0.0, t_stop=1000.0, array=False):
+    def poisson_generator(self, rate, t_start=0.0, t_stop=1000.0, array=False,debug=False):
         """
         Returns a SpikeList whose spikes are a realization of a Poisson process
@@ -98 +185 @@
         """
 
-        number = int((t_stop-t_start)/1000.0*2.0*rate)
+        #number = int((t_stop-t_start)/1000.0*2.0*rate)
+        
+        # less wasteful than double length method above
+        n = (t_stop-t_start)/1000.0*rate
+        number = numpy.ceil(n+3*numpy.sqrt(n))
+        if number<100:
+            number = min(5+numpy.ceil(2*n),100)
+        
         if number > 0:
             isi = self.rng.exponential(1.0/rate, number)*1000.0
@@ -111 +205 @@
         i = numpy.searchsorted(spikes, t_stop)
 
+        extra_spikes = []
         if i==len(spikes):
-            raise RuntimeError("Internal ISI buffer overrun.  Please file a bug report ticket at http://neuralensemble.org/NeuroTools.")
-
-        if array:
-            return numpy.resize(spikes,(i,))
-
-        return SpikeTrain(numpy.resize(spikes,(i,)), t_start=t_start,t_stop=t_stop)
+            # ISI buf overrun
+            
+            t_last = spikes[-1] + self.rng.exponential(1.0/rate, 1)[0]*1000.0
+
+            while (t_last<t_stop):
+                extra_spikes.append(t_last)
+                t_last += self.rng.exponential(1.0/rate, 1)[0]*1000.0
+            
+            spikes = numpy.concatenate((spikes,extra_spikes))
+
+            if debug:
+                print "ISI buf overrun handled. len(spikes)=%d, len(extra_spikes)=%d" % (len(spikes),len(extra_spikes))
+
+
+        else:
+            spikes = numpy.resize(spikes,(i,))
+
+        if not array:
+            spikes = SpikeTrain(spikes, t_start=t_start,t_stop=t_stop)
+
+
+        if debug:
+            return spikes, extra_spikes
+        else:
+            return spikes
 
 
@@ -270 +384 @@
     # TODO: provide optimized C/weave implementation if possible
 
-    inh_gamma_generator = _inh_gamma_generator_python
+
+    def _inh_gamma_generator_unimplemented(self, a, b, t, t_stop, array=False):
+        raise Exception("inh_gamma_generator is disabled as dependency RPy was not found.")
+
+    if check_dependency('rpy'):
+        inh_gamma_generator = _inh_gamma_generator_python
+    else:
+        inh_gamma_generator = _inh_gamma_generator_unimplemented
 
 
@@ -432 +553 @@
     OU_generator = _OU_generator_python2
 
-    # TODO: inhomogeneous OU generator
+    # TODO: optimized inhomogeneous OU generator
 
 
@@ -439 +560 @@
 # TODO fix shotnoise stuff below  ... and write tests
 
-    def shotnoise_generator(self,rate,tau,q,num,t):
-        """ 
-        Generate shotnoise
-
-        Inputs:
-            rate - the rate (in Hz) of the shotnoise
-            tau  - the exponential decay of the synapse
-            g    - the quantal increase for a spike
-
-        quantal-increase-"q"-exponential-decay-"tau" synapse
-        and a poisson spike train of "rate" and "num" """
-
-        g = numpy.zeros(numpy.shape(t),float)
-        for i in xrange(num):
-
-            spikes = poisson_generator(rate,t[-1])
-            dg=exp_conv(spikes,t,q,tau)
-            tmp = add(g,dg,g)
-
-        return g
+def shotnoise_fromspikes(self,rate,tau,q,num,t):
+    """ 
+    Generate shotnoise
+
+    Inputs:
+        rate - the rate (in Hz) of the shotnoise
+        tau  - the exponential decay of the synapse
+        g    - the quantal increase for a spike
+
+    quantal-increase-"q"-exponential-decay-"tau" synapse
+    and a poisson spike train of "rate" and "num" """
+
+    g = numpy.zeros(numpy.shape(t),float)
+    for i in xrange(num):
+
+        spikes = poisson_generator(rate,t[-1])
+        dg=exp_conv(spikes,t,q,tau)
+        tmp = add(g,dg,g)
+
+    return g
 
 
@@ -466 +587 @@
 
 
-def exp_conv(poisson_train,t,q,tau,eps = 1.0e-8):
+def shotnoise_fromspikes(spike_train,q,tau,dt,array=False, eps = 1.0e-8):
     """ 
-    Convolve poisson spike trains with shot decaying exponentials
-    t must be equally spaced arrayrange
-    poisson spike times must all in the range of t
-    otherwise unpredicted results.
-    """
-
-    dt = t[1]-t[0]
+    Convolves the provided spike train with shot decaying exponentials
+    yielding so called shot noise if the spike train is Poisson-like.  
+    Returns an AnalogSignal if array=False, otherwise (shotnoise,t) as numpy arrays. 
+
+   Inputs:
+      spike_train - a SpikeTrain object
+      q - the shot jump for each spike
+      tau - the shot decay time constant in milliseconds
+      dt - the resolution of the resulting shotnoise in milliseconds
+      array - if True, returns (shotnoise,t) as numpy arrays, otherwise an AnalogSignal.
+      eps - a numerical parameter indicating at what value of 
+      the shot kernal the tail is cut.  The default is usually fine.
+
+   Examples:
+      >> stg = stgen.StGen()
+      >> st = stg.poisson_generator(10.0,0.0,1000.0)
+      >> g_e = shotnoise_fromspikes(st,2.0,10.0)
+
+
+   See also:
+      poisson_generator, inh_gamma_generator, OU_generator ...
+   """
+
+    st = spike_train
+
+    t = numpy.arange(st.t_start,st.t_stop,dt)
 
     # time of vanishing significance
-    vs_t = -tau*log(eps/q)
-
-    kern = q*exp(-arrayrange(0.0,vs_t,dt)/tau)
-
-    idx = clip(searchsorted(t,poisson_train),0,len(t)-1)
-
-    a = zeros(shape(t),Float)
-
-    put(a,idx,1.0)
-
-    return convolve(a,kern)[0:len(t)]
+    vs_t = -tau*numpy.log(eps/q)
+
+    kern = q*numpy.exp(-numpy.arange(0.0,vs_t,dt)/tau)
+
+    idx = numpy.clip(numpy.searchsorted(t,poisson_train,'right')-1,0,len(t)-1)
+
+    a = numpy.zeros(shape(t),float)
+
+    a[idx] = 1.0
+
+    y = convolve(a,kern)[0:len(t)]
+
+    result = AnalogSignal(y,dt,t_start=0,t_stop=st.t_stop-st.t_start)
+    result.time_offset(st.t_start)
+    return result
+