#lang scribble/doc
@(require scribble/manual
          scribblings/icons
          (for-label scheme/base
                     plot/plot
                     "../simulation-with-graphics.ss"))

@title[#:tag "simpy"]{SimPy and Simulation Collection Example Comparisons}

@local-table-of-contents[]

This chapter had some example from the SimPy web site and equivalent examples in the simulation collection.

@section{Cellphone}

@subsection{Simpy Example---Cellphone}

@verbatim{
#!/usr/bin/env python
""" Simulate the operation of a BCC cellphone system using SimPy

"""
from __future__ import generators   # not needed for Python 2.3+
from SimPy.Simulation import *
from random import Random,expovariate

class Generator(Process):
 """ generates a sequence of calls """
 def __init__(self, maxN, lam):
     global busyEndTime
     Process.__init__(self)
     self.name = "generator"
     self.maxN = maxN
     self.lam = lam
     self.iatime = 1.0/lam
     self.rv = Random(gSeed)
     busyEndTime = now() # simulation start time

 def execute(self):
     for i in range(self.maxN):
         j = Job(i)
         activate(j,j.execute())
         yield hold,self,self.rv.expovariate(lam)
     self.trace("WARNING generator finished -- ran out of jobs")

 def trace(self,message):
     if GTRACING: print "%7.4f \t%s"%(self.time(), message)


class Job(Process):
    """ instances of the Job class represent calls arriving at
    random at the chnnels.
    """
    def __init__(self,i):
        Process.__init__(self)
        self.i = i
        self.name = "Job"+`i`

    def execute(self):
        global busyStartTime,totalBusyVisits,totalBusyTime
        global Nfree,busyEndTime,Jrv
        self.trace("arrived ")
        if Nfree == 0: self.trace("blocked and left")
        else:
             self.trace("got a channel")
             Nfree -=  1
             if Nfree == 0:
                 self.trace("start busy period======")
                 busyStartTime = now()
                 totalBusyVisits += 1
                 interIdleTime = now() - busyEndTime
             yield hold,self,Jrv.expovariate(mu)
             self.trace("finished")
             if Nfree == 0:
                 self.trace("end   busy period++++++")
                 busyEndTime = now()
                 busy = now() - busyStartTime
                 self.trace("         busy  = %9.4f"%(busy,))
                 totalBusyTime +=busy
             Nfree += 1

    def trace(self,message):
         if TRACING: print "%7.4f \t%s %s "%(now(), message , self.name)


class Statistician(Process):
     """ observes the system at intervals """
     def __init__(self,Nhours,interv,gap):
         Process.__init__(self)
         self.Nhours = Nhours
         self.interv = interv
         self.gap = gap
         self.name="Statistician"

     def execute(self):
         global busyStartTime,totalBusyTime,totalBusyVisits
         global hourBusyTime,hourBusyVisits
         for i in range(self.Nhours):
             yield hold,self,self.gap
             totalBusyTime = 0.0
             totalBusyVisits = 0
             if Nfree == 0: busyStartTime = now()
             yield hold,self,self.interv
             if Nfree == 0: totalBusyTime += now()-busyStartTime
             if STRACING:  print "%7.3f %5d"%(totalBusyTime,totalBusyVisits)
             m.tally(totalBusyTime)
             bn.tally(totalBusyVisits)
         print("Busy Time:   mean = %10.5f var= %10.5f"%(m.mean(),m.var()))
         print("Busy Number: mean = %10.5f var= %10.5f"%(bn.mean(),bn.var()))


     def trace(self,message):
         if STRACING: print "%7.4f \t%s"%(self.time(), message)

totalBusyVisits = 0
totalBusyTime   = 0.0
NChannels =  4  # number of channels in the cell
Nfree   = NChannels
maxN    = 1000
gSeed   = 11111111
JrvSeed = 3333333
lam = 1.0
mu = 0.6667
meanLifeTime = 1.0/mu
TRACING  = 0
GTRACING = 0
STRACING = 1
Nhours  =  10
interv = 60.0    # monitor
gap    = 15.0    # monitor
print "lambda    mu      s  Nhours interv  gap"
print "%7.4f %6.4f %4d %4d    %6.2f %6.2f"%(lam,mu,NChannels,Nhours,interv,gap)

m = Monitor()
bn=Monitor()
Jrv = Random(JrvSeed)
s = Statistician(Nhours,interv,gap)
initialize()
g = Generator(maxN, lam)
activate(g,g.execute())
activate(s,s.execute())
simulate(until=10000.0)}

The output from the program follows.

@verbatim{
lambda    mu      s  Nhours interv  gap
 1.0000 0.6667    4   10     60.00  15.00
  3.274     8
  0.444     2
  3.649     8
  0.652     3
  2.277     6
  1.372     9
  0.411     4
  0.067     2
  2.910     8
  5.073     9
Busy Time:   mean =    2.01297 var=    2.55814
Busy Number: mean =    5.90000 var=    7.49000}

@subsection{Simulation Collection Example---Cellphone}

@schememod[
scheme/base
(code:comment " Cellphone simulation from PySim")

(require (planet williams/simulation/simulation))
(require (planet williams/science/random-distributions))

(code:comment " Parameters")

(define Nchannels 4)
(define maxN 1000)
(define lam 1.0)
(define mu 0.6667)
(define meanLifeTime (/ 1.0 mu))
(define Nhours 10)
(define interv 60.0)
(define gap 15.0)

(code:comment " Globals")

(define Nfree Nchannels)
(define totalBusyVisits 0)
(define totalBusyTime 0.0)
(define busyStartTime 0.0)
(define busyEndTime 0.0)


(code:comment " Statistics")

(define m #f)
(define bn #f)

(define-process (generator n)
  (do ((i 0 (+ i 1)))
      ((= i n) (void))
    (schedule now (job i))
    (wait (random-exponential (/ 1.0 lam)))))

(define-process (job i)
  (when (> Nfree 0)
    (set! Nfree (- Nfree 1))
    (when (= Nfree 0)
      (code:comment " Start busy period")
      (set! busyStartTime (current-simulation-time))
      (set! totalBusyVisits (+ totalBusyVisits 1)))
    (work (random-exponential (/ 1.0 mu)))
    (when (= Nfree 0)
      (code:comment " End busy period")
      (set! busyEndTime (current-simulation-time))
      (set! totalBusyTime (+ totalBusyTime
                             (- (current-simulation-time)
                                busyStartTime))))
    (set! Nfree (+ Nfree 1))))

(define-process (statistician)
  (do ((i 0 (+ i 1)))
      ((= i Nhours) (void))
    (code:comment " wait the specified gap time")
    (wait gap)
    (code:comment " initialize")
    (set! totalBusyTime 0.0)
    (set! totalBusyVisits 0)
    (when (= Nfree 0)
      (set! busyStartTime (current-simulation-time)))
    (code:comment " wait the specified interval time")
    (wait interv)
    (code:comment " trace busy time and busy visits")
    (when (= Nfree 0)
      (set! totalBusyTime (+ totalBusyTime
                             (- (current-simulation-time)
                                busyStartTime))))
    (printf "~a: busy time = ~a; busy visits = ~a~n"
            (current-simulation-time) totalBusyTime totalBusyVisits)
    (code:comment " Tally statistics")
    (set-variable-value! m totalBusyTime)
    (set-variable-value! bn totalBusyVisits))
  (code:comment " Print final statistics")
  (printf "Busy time:   mean = ~a, variance = ~a~n"
          (variable-mean m) (variable-variance m))
  (printf "Busy number: mean = ~a, variance = ~a~n"
          (variable-mean bn) (variable-variance bn)))

(define (run-simulation)
  (with-new-simulation-environment
   (code:comment " Print parameters")
   (printf "lambda = ~a~n" lam)
   (printf "mu     = ~a~n" mu)
   (printf "s      = ~a~n" Nchannels)
   (printf "Nhours = ~a~n" Nhours)
   (printf "interv = ~a~n" interv)
   (printf "gap    = ~a~n" gap)
   (code:comment " Initialize global variables")
   (set! Nfree Nchannels)
   (set! totalBusyVisits 0)
   (set! totalBusyTime 0.0)
   (set! busyEndTime 0.0)
   (code:comment " Create statistics")
   (set! m (make-variable))
   (tally (variable-statistics m))
   (set! bn (make-variable))
   (tally (variable-statistics bn))
   (code:comment " Create initial processes and start the execution")
   (schedule (at 0.0) (generator maxN))
   (schedule (at 0.0) (statistician))
   (schedule (at 10000.0) (stop-simulation))
   (start-simulation)))

(run-simulation)]

The output from the program follows.

@verbatim{
lambda = 1.0
mu     = 0.6667
s      = 4
Nhours = 10
interv = 60.0
gap    = 15.0
75.0: busy time = 3.19230597325965; busy visits = 11
150.0: busy time = 2.4780498393357675; busy visits = 7
225.0: busy time = 3.5391222097609045; busy visits = 12
300.0: busy time = 2.1031429717027947; busy visits = 6
375.0: busy time = 1.3915494625956057; busy visits = 3
450.0: busy time = 2.9838869629721785; busy visits = 6
525.0: busy time = 3.8242483424334637; busy visits = 8
600.0: busy time = 2.7636093265058435; busy visits = 9
675.0: busy time = 1.73398034270042; busy visits = 8
750.0: busy time = 1.957177682358406; busy visits = 9
Busy time:   mean = 2.596707311362503, variance = 0.5790891717346138
Busy number: mean = 7.9, variance = 6.089999999999996}

@section{Jackson}

@subsection{Simpy Example---Jackson}

@verbatim{
#!/usr/bin/env python
"""
  Simulation of a jackson queue network

"""
from __future__ import generators  # not needed in Python 2.3+
from SimPy.Simulation import *
from random import Random,expovariate,uniform

def sum(P):
    """ Sum of the real vector P """
    sumP = 0.0
    for i in range(len(P)): sumP += P[i]
    return (sumP)



def choose2dA(i,P,g):
    """  return a random choice from a set j = 0..n-1
       with probs held in list of lists P[j] (n by n)
       using row i
       g = random variable
       call:  next = choose2d(i,P,g)
    """
    U = g.random()
    sumP = 0.0
    for j in range(len(P[i])):  # j = 0..n-1
        sumP +=  P[i][j]
        if U < sumP: break
    return(j)

## -----------------------------------------------
class Msg(Process):
    """a message"""
    def __init__(self,i,node):
        Process.__init__(self)
        self.i = i
        self.node = node

    def execute(self):
        """ executing a message """
        global noInSystem
        startTime = now()
        noInSystem += 1
        ##print "DEBUG noInSystm = ",noInSystem
        NoInSystem.accum(noInSystem)
        self.trace("Arrived node  %d"%(self.node,))
        while self.node <> 3:
            yield request,self,processor[self.node]
            self.trace("Got processor %d"%(self.node,))
            st = Mrv.expovariate(1.0/mean[self.node])
            Stime.observe(st)
            yield hold,self,st
            yield release,self,processor[self.node]
            self.trace("Finished with %d"%(self.node,))
            self.node = choose2dA(self.node,P,Mrv)
            self.trace("Transfer to   %d"%(self.node,))
        TimeInSystem.observe(now()-startTime)
        self.trace(    "leaving       %d"%(self.node,),noInSystem)
        noInSystem -= 1
        NoInSystem.accum(noInSystem)

    def trace(self,message,nn=0):
        if MTRACING: print "%7.4f %3d %10s %3d"%(now(),self.i, message,nn)


## -----------------------------------------------
class Generator(Process):
 """ generates a sequence of msgs """
 def __init__(self, rate,maxT,maxN):
     Process.__init__(self)
     self.name = "Generator"
     self.rate = rate
     self.maxN = maxN
     self.maxT = maxT
     self.g = Random(11335577)
     self.i = 0

 def execute(self):
     while (now() < self.maxT)  & (self.i < self.maxN):
         self.i+=1
         p = Msg(self.i,startNode)
         activate(p,p.execute())
         ## self.trace("starting "+p.name)
         yield hold,self,self.g.expovariate(self.rate)
     self.trace("generator finished with "+`self.i`+" ========")
     self.stopSim()

 def trace(self,message):
     if GTRACING: print "%7.4f \t%s"%(now(), message)

## -----------------------------------------------
# generator:
GTRACING = 0
# messages:
rate = 0.5
noInSystem = 0
MTRACING = 0
startNode = 0
maxNumber = 10000
maxTime = 100000.0
Mrv = Random(77777)
TimeInSystem=Monitor()
NoInSystem= Monitor()
Stime = Monitor()

processor = [Resource(1),Resource(1),Resource(1)]
mean=[1.0, 2.0, 1.0]
P = [[0, 0.5, 0.5, 0],[0,0,0.8,0.2], [0.2,0,0,0.8]]

initialize()
g = Generator(rate,maxTime,maxNumber)
activate(g,g.execute())
simulate(until=100.0)

print "Mean numberm= %8.4f"%(NoInSystem.timeAverage(),)
print "Mean delay  = %8.4f"%(TimeInSystem.mean(),)
print "Mean stime  = %8.4f"%(Stime.mean(),)
print "Total jobs  = %8d"%(TimeInSystem.count(),)
print "Total time  = %8.4f"%(now(),)
print "Mean rate   = %8.4f"%(TimeInSystem.count()/now(),)}

The following is the output from the program.

@verbatim{
Mean numberm=   1.5412
Mean delay  =   3.9369
Mean stime  =   1.0332
Total jobs  =       38
Total time  = 100.0000
Mean rate   =   0.3800}

@subsection{Simulation Collection Example---Jackson}

@schememod[
scheme/base
;;; Jackson Network from PySim

(require (planet williams/simulation/simulation-with-graphics))
(require (planet williams/science/random-distributions))

;;; Parameters

(define n 3)
(define m #(1.0 2.0 1.0))
(define p (vector
           (make-discrete #(0.0 0.5 0.5 0.0))
           (make-discrete #(0.0 0.0 0.8 0.2))
           (make-discrete #(0.2 0.0 0.0 0.8))))
(define np '(1 1 1))
(define mu (/ 1.0 1.5))
(define max-messages 10000)

;;; Global Variables (set in main)

(define n-arrivals 0)
(define n-in-system #f)
(define time-in-system #f)
(define processor #f)

;;; Data Collection

(define-process (message i)
  (let ((arrive-time (current-simulation-time))
        (work-duration 0.0))
    (set-variable-value!
     n-in-system (+ (variable-value n-in-system) 1))
    (let loop ((node 0))               ; start node is always 0
      (when (< node n)
        (with-resource ((vector-ref processor node))
          (work (random-exponential (vector-ref m node))))
        (loop (random-discrete (vector-ref p node)))))
    (set-variable-value!
     n-in-system (- (variable-value n-in-system) 1))
    (set-variable-value!
     time-in-system (- (current-simulation-time)
                       arrive-time))))

(define (generator)
  (do ((i 0 (+ i 1)))
      ((= i max-messages) (void))
    (set! n-arrivals (+ n-arrivals 1))
    (schedule now (message i))
    (wait (random-exponential (/ mu)))))

(define (main)
  (with-new-simulation-environment
    ;; Initialize global variables
    (set! n-arrivals 0)
    (set! n-in-system (make-variable 0))
    (accumulate (variable-history n-in-system))
    (set! time-in-system (make-variable))
    (tally (variable-statistics time-in-system))
    (set! processor (list->vector
                     (map make-resource np)))
    ;;
    (schedule (at 0.0) (generator))
    (schedule (at 5000.0) (stop-simulation))
    (start-simulation)
    (printf "Mean number in system = ~a~n" (variable-mean n-in-system))
    (printf "Mean delay in system  = ~a~n" (variable-mean time-in-system))
    (printf "Total time run        = ~a~n" (current-simulation-time))
    (printf "Total jobs arrived    = ~a~n" n-arrivals)
    (printf "Total jobs completed  = ~a~n" (variable-n time-in-system))
    (printf "Average arrival rate  = ~a~n" (/ n-arrivals
                                              (current-simulation-time)))
    (write-special (history-plot (variable-history n-in-system)))
    (newline)))

(main)]

The following is the output from the program.

@verbatim{
Mean number in system = 14.726063392545544
Mean delay in system  = 21.56635336305308
Total time run        = 5000.0
Total jobs arrived    = 3419
Total jobs completed  = 3413
Average arrival rate  = 0.6838}

@image["scribblings/images/jackson.png"]