import random, math, pylab, heapq 
                # heapq provides the Heap Queue structure which allows to manipulate efficiently sorted list.
                # Here we use a list of copies of the system, sorted according to their next time of evolution.

c=.25           # c = creation rate 0->1 ; 1-c = annihilation rate 1->0
N=100000         # population size
populat=[(0.,0,0) for n in range(N)]
                # initial state of the population : 
                # each member of (or "copy" in) the population is described by a pair (time,state) 
                #   time  = next time at which it will evolve
                #   state = 0 or 1 = empty or occupied 
                #   K = observed activity starting from time 0
t,tmax=0.,40.    # initial and maximal time

def escaperate(n):  
    if n==0:
        rate=c
    else:
        rate=(1.-c)
    return rate

def P1(t):  #theoretical probability to be in state 1 at time t
    return c-c*math.exp(-t)

heapq.heapify(populat)  # orders the population into a Heap Queue ; useful if the initial state of the population is not ordered in time
while t<tmax:
    (t,state,K)=heapq.heappop(populat) # pops and yields the first element of populat, which is always the next to evolve
    Deltat=random.expovariate(escaperate(state)) # interval until next evolution
    newstate=(t+Deltat,1-state,K+1)      
    heapq.heappush(populat,newstate)

print 'Theoretical value of the probability to be in state 1 at time t=',tmax, ': P1(',tmax,')=',P1(tmax)
states=[state for (time,state,K) in populat]
print 'Proportion of copies of the system being   in state 1 at time t=',tmax, ': f1(',tmax,')=', 1-sum(states)*1./N

pylab.title('histogram for the state')
data=[1-state for (time,state,K) in populat]
pylab.hist(data,bins=2,normed=True)
pylab.show()

pylab.title('histogram for the dynamical activity')
data=[K for (time,state,K) in populat]
pylab.hist(data,bins=20,normed=True)
pylab.show()