Hello list:

On the previous topic I asked about using simecol for discrete event, process-based simulation.  Included bellow is the code for the simulation of a simple M:M:1 - GD:oo:oo queue using three approaches (activity based, event based and process based).  The first two scripts are R code, the last is Python using the SimPy library.

Please note the following:
1) The activity-based version is the slowest of all.  Furthermore, as the ocurrence of arrivals and services (both Poisson process) is determined by Bernoulli trials, this approach is very inexact.  No doubt this simulation model could be implemented by the individual based simecol objects (more time-efficiently?)

2) The event based version is more efficient, as the clock jumps from event to event.  However, the programmer has to take care of all the book-keeping details (event scheduling).  Consider for example how one would have to re-write the code completely if more servers are used, there is queue reneging behavior, etc.

3) the process based approach is the simplest for the programmer.  It would be nice if something like SimPy could be implemented in R.  It occurs to me that the fork R-package, which handles processes, could be used for this end.

My questions are:

Is this process based approach way too far out of the simecol scope (simecol objects being essentially based on time iterators, solvers and such)?
Is the process-based approach useful for the implementation of other ecological models? 
Or considering the event-based approach, could such programs be vectorized? Is there a way to implement solvers in simecol that advance the clock to the nearest future event?

Thanks and Regards,

jlrp 

##=====================================================================
## Simulation of a simple M:M:1 - GD:oo:oo queue
## Activity based approach
## JosÃ© Romero - sep 4, 2010
##=====================================================================

#---------------------    
# initial parameters
#---------------------

tmax <- 10000  v        #maximum time of simulation
alpha <- 6            #arrival rate
lambda <- 8            #service rate
queue <- NULL        #client queue, initially empty
service_times <- NULL        #record of service times
times <- numeric(0)        #record of total time for each state
step <- 0.001            #time-steps
p_alpha <- alpha*step        #probability of arrival in one time step
p_lambda <- lambda*step        #probability of departure in one time step
clock <- 0


#---------------------    
# function definitions
#---------------------

not_empty <- function(queue) return(!(length(queue)==0))

empty <- function(queue) return (length(queue)==0)

Update_clock <- function() {
    state <- length(queue) + 1
    if (is.na(times[state])) times[state] <<- 0
    (sample(size=1,c(0,1),prob=c(1-p_alpha,p_alpha))==1)
    times[state] <<- times[state] + step
    clock <<- clock + step    
    }

#----------------------------------------------------------------------
# Main simulation loop:
#----------------------------------------------------------------------

for (i in 1:(tmax/step)) {
    arrival <- (sample(size=1,c(0,1),prob=c(1-p_alpha,p_alpha))==1)
    if (arrival) queue <- c(queue,clock)
    dispatch <- (sample(size=1,c(0,1),prob=c(1-p_lambda,p_lambda))==1) & 
            (not_empty(queue))
    if (dispatch) {
        arrival_time <- queue[1]
        queue <- queue[-1]
        service_times <- c(service_times,clock - arrival_time)
    }
    Update_clock()
}
    

#----------------------------------------------------------------------
# Simulation results:
#----------------------------------------------------------------------

#first convert all NA's in times to 0
times[is.na(times)] <- 0

#then divide times by clock to obtain probability distribution
times <- times / clock

#last compare to theoretical state distribution
#and plot results
    max_x <- length(times)-1    
    theoretical <- dgeom(x=0:max_x,prob=1- alpha/lambda)
    plot(x=rep(0:max_x,2),y=c(times,theoretical),col=c(rep("blue",max_x+1),
        rep("red",max_x+1)))
    mean(service_times)


##=====================================================================
## Simulation of a simple M:M:1 - GD:oo:oo queue
## Event based approach
## JosÃ© Romero - sep 4, 2010
##=====================================================================

#---------------------    
# initial parameters
#---------------------

tmax <- 10000            #maximum time of simulation
alpha <- 6            #arrival rate
lambda <- 8            #service rate
queue <- NULL            #client queue, initially empty
service_times <- NULL        #record of service times
times <- numeric(0)        #record of total time for each state
clock <- 0


#---------------------    
# function definitions
#---------------------

not_empty <- function(queue) return(!(length(queue)==0))

empty <- function(queue) return (length(queue)==0)

Update_clock <- function(next_event_time) {
    time_dif <- next_event_time - clock
    state <- length(queue) + 1
    if (is.na(times[state])) times[state] <<- 0
    times[state] <<- times[state] + time_dif
    clock <<- next_event_time    
    }

Arrival <- function(queue) {
    queue <<- c(queue, clock)
    arrival_time <- clock + rexp(n=1,rate=alpha)
    return(arrival_time)
    }

Dispatch <- function(queue) {
    arrival_time <- queue[1]
    queue <<- queue[-1]
    return(arrival_time)
    }

#----------------------------------------------------------------------
# Main simulation loop:
#----------------------------------------------------------------------

next_arrival <- rexp(1,rate=alpha)

while (clock<tmax) {
    
    if (empty(queue)) next_dispatch <- next_arrival+rexp(n=1,rate=lambda)
    
    while (next_arrival<=next_dispatch) {
        Update_clock(next_arrival)        
        next_arrival <- Arrival(queue)     
    }
    while (not_empty(queue) & (next_dispatch<next_arrival)) {
        Update_clock(next_dispatch)
        service_time <- clock - Dispatch(queue)
        service_times <- c(service_times,service_time)    
        if (not_empty(queue)) next_dispatch <- clock + rexp(n=1,rate=lambda)
    }
}

#----------------------------------------------------------------------
# Simulation results:
#----------------------------------------------------------------------

#first convert all NA's in times to 0
times[is.na(times)] <- 0

#then divide times by clock to obtain probability distribution
times <- times / clock

#last compare to theoretical state distribution
#and plot results
    max_x <- length(times)-1    
    theoretical <- dgeom(x=0:max_x,prob=1- alpha/lambda)
    plot(x=rep(0:max_x,2),y=c(times,theoretical),col=c(rep("blue",max_x+1),
        rep("red",max_x+1)))
    mean(service_times)


##=====================================================================
## Simulation of a simple M:M:1 - GD:oo:oo queue
## Process based approach
## Jose Romero - sep 4, 2010
##=====================================================================

#---------------------    
# initial parameters
#---------------------

from SimPy.Simulation import *
from random import *


class Client(Process):
    served_clients = []
    def go(self,arrival_time,myServer):
        yield request,self,myServer
        yield hold,self,expovariate(lambd=8)
        yield release,self,myServer
                
class Client_source(Process):
    def go(self,finish,myServer):
        while(now()<finish):
            c=Client()
            activate(c,c.go(arrival_time=now(),myServer=myServer))
            yield hold,self,expovariate(lambd=6)
        
initialize()
endtime=1000
server=Resource(capacity=1,monitored=True,monitorType=Monitor)
client_source=Client_source()
activate(client_source,client_source.go(finish=endtime,myServer=server))
simulate(until=endtime)
print server.waitMon


 

	[[alternative HTML version deleted]]

