#!/usr/bin/env python
from __future__ import generators
from SimPy.Simulation import *
from random import Random,expovariate,uniform
# sip095  erlang-B loss M/H/c/c, M/Ek/c/c, M/M/c/c limited
#        $Author: kgmuller $
# $Revision: 1.1.1.5 $  $Date: 2006/02/02 13:35:45 $
# converted from pyt095 2002 Sept 10
# 2002 9 20 included bcc function directly.
"""
Jobs (e.g messages) arrive randomly at rate $1.0$ per minute at a
  2-server system.  Mean service time is $0.75$ minutes and the
  service-time distribution is (1) exponential, (2) Erlang-5, or (3)
  hyperexponential with $p=1/8$,$m1=2.0$, and $m2=4/7$. However no
  queue is allowed; a job arriving when all the servers are busy is
  rejected.
  
  Develop and run a simulation program to estimate the probability of
  rejection (which, in steady-state, is the same as $p(c)$) Measure
  and compare the probability for each service time distribution.
  Though you should test the program with a trace, running just a few
  jobs, the final runs should be of 10000 jobs without a trace. Stop
  the simulation when 10000 jobs have been generated. Use the random
  number stream 1 for all runs.
"""
__version__ = '\nModel: BCC at 2-server system $Revision: 1.1.1.5 $ $Date: 2006/02/02 13:35:45 $'

dist = ""
def bcc(lam, mu,s):
    """ bcc - blocked customers cleared model

    - returns p[i], i = 0,1,..s.
    - ps = p[s] = prob of blocking
    - lameff = effective arrival rate = lam*(1-ps)

    See Winston 22.11 for Blocked Customers Cleared Model (Erlang B formula)
    """
    rho = lam/mu
    n = range(s+1)
    p = [0]*(s+1)
    p[0] = 1
    sump = 1.0
    for i in n[1:]:
        p[i] = (rho/i)*p[i-1]
        sump = sump + p[i]
    p0 = 1.0/sump
    for i in n:
        p[i] = p[i]*p0
    p0 = p[0]    
    ps = p[s]
    lameff = lam*(1-ps)
    L = rho*(1-ps)
    return {'lambda':lam,'mu':mu,'s':s,
	'p0':p0,'p[i]':p,'ps':ps, 'L':L}

def ErlangVariate(mean,K,g):
    """ Erlang random variate

    mean = mean
    K = shape parameter
    g = rv to be used
    """
    sum = 0.0 ; mu = K/mean
    for i in range(K):
        sum += g.expovariate(mu)
    return (sum)

def HyperVariate(p,m1,m2,g):
    """ Hyperexponential random variate

    p = prob of branch 1
    m1 = mean of exponential, branch 1
    m2 = mean of exponential, branch 2
    g = rv to be used
    """
    if g.random() < p:
        return g.expovariate(1.0/m1)
    else: return g.expovariate(1.0/m2)

def testHyperVariate():
    """ tests the HyerVariate rv generator"""
    ERR=0
    x = (1.0981,1.45546,5.7470156)
    p = 0.0, 1.0 ,0.5
    g = Random(1113355)
    for i in range(3):
        x1 = HyperVariate(p[i],1.0,10.0,g)
        #print p[i], x1
        assert abs(x1 - x[i]) < 0.001,'HyperVariate error'
    

## use bcc instead!!
def erlangB(rho,c):
    """ Erlang's B formula for probabilities in no-queue

    Returns p[n] list
    see also SPlus and R version in que.q mmcK
    que.py has bcc.
    """
    n = range(c+1) ; pn = range(c+1)
    term = 1
    pn[0] = 1
    sum = 1 ; term = 1.0
    i=1
    while i < (c+1):
        term  *= rho/i
        pn[i] = term
        sum += pn[i]
        i += 1
    for i in n: pn[i] = pn[i]/sum
    return(pn)


class JobGen(Process):
    """ generates a sequence of Jobs
    """    
    def __init__(self, JobRate,MaxJob,mu):
         Process.__init__(self)
         self.JobRate = JobRate
         self.MaxJob = MaxJob
         self.g = Random(111333)
         self.i = 0
         self.mu =mu
     
    def execute(self):
         global NoInService, Busy
         for i in range(self.MaxJob):         
             j = Job(i,self.mu)
             activate(j,j.execute(),delay=0.0)
             t = self.g.expovariate(self.JobRate)
             MT.tally(t)
             yield hold,self,t
         self.trace("Job generator finished")

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

class Job(Process):
    """ Jobs that are either accepted or rejected
    """
    g = Random(779977)
    def __init__(self, i, mu):
         Process.__init__(self)
         self.mu = mu
         self.i = i
     
    def execute(self):
        """ Job execution, only if accepted"""
        global NoInService,Busy,dist,NoRejected
        if NoInService < c:
            self.trace("Job %2d accepted b=%1d"%(self.i,Busy))
            NoInService +=1
            if NoInService == c:
                Busy =1
                try: BM.accum(Busy,now())
                except: "accum error BM=",BM
            #yield   hold,self,Job.g.expovariate(self.mu);           dist= "Exponential"
            yield   hold,self,ErlangVariate(1.0/self.mu,5,Job.g);          dist= "Erlang     "
            #yield   hold,self,HyperVariate(1.0/8,m1=2.0,m2=4.0/7,g=Job.g); dist= "HyperExpon "
            NoInService -=1
            Busy =0
            BM.accum(Busy,now())
            self.trace("Job %2d leaving b=%1d"%(self.i,Busy))
        else:
            self.trace("Job %2d REJECT  b=%1d"%(self.i,Busy))
            NoRejected +=1
 
    def trace(self,message):
        if JobTRACING: print "%8.4f \t%s"%(now(), message)


##--------------------------------------------------------------
c = 2
lam = 1.0
mu = 1.0/0.75
p = 1.0/8 ; m1= 2.0;  m2 = 4.0/7.0
K = 5
rho = lam/mu

NoRejected = 0
NoInService = 0
Busy = 0

JobRate = lam
JobMax = 10000

JobTRACING = 0
JobGenTRACING = 0

BM=Monitor()
MT = Monitor()

def main():
    initialize()
    jbg = JobGen(1.0,JobMax,mu)
    activate(jbg,jbg.execute(),0.0)
    ## print now()
    simulate(until=20000.0)
    print "time at the end =",now()
    ##print erlangB(rho,c)
    print "now=",now(), " startTime ",BM.startTime
    print "No Rejected = %d, ratio= %s"%(NoRejected,(1.0*NoRejected)/JobMax)
    print "Busy proportion (%6s) = %8.6f"%(dist,BM.timeAverage(now()),)
    print "Erlang pc (th)                = %8.6f"%(erlangB(rho,c)[c],)

##--------------------------------------------------------------

if __name__ == '__main__':
    print __version__
    ## testHyperVariate()
    main()