from heapq import heappush, heappop
import random
from scipy.special import factorial

class Stats:
    def __init__(self, sim):
        self.sim = sim
        self.waiting_time = 0.0
        self.servicing_time = 0.0
        
    def avg_wait_time(self):
        return self.waiting_time / self.sim.t
    
    def avg_service_time(self):
        return self.servicing_time / self.sim.t
    
    def avg_system_time(self):
        return (self.waiting_time + self.servicing_time) / self.sim.t
    
    def prb_empty(self):
        c = self.sim.servers
        rho = self.sim.arrival_rate / self.sim.service_rate
        if self.sim.servers == 1:
            return 1.0 - rho
        return 1 / (sum(rho**n/factorial(n) for n in range(0, c)) +
                    rho**c/factorial(c) * 1/(1-rho/c))
    
    def exp_wait_time(self):
        c = self.sim.servers
        if c == 1:
            return self.exp_system_time() - self.exp_service_time()
        rho = self.sim.arrival_rate / self.sim.service_rate
        return self.prb_empty()*rho**(c+1) / (factorial(c-1)*(c-rho)**2)
    
    def exp_service_time(self):
        c = self.sim.servers
        rho = self.sim.arrival_rate / self.sim.service_rate
        if c == 1:
            return rho
        return self.exp_system_time() - self.exp_wait_time()
    
    def exp_system_time(self):
        c = self.sim.servers
        rho = self.sim.arrival_rate / self.sim.service_rate
        if c == 1:
            return rho / (1 - rho)
        return self.exp_wait_time() + rho
    
    def show(self):
        summaries = [
            'avg_wait_time',
            'exp_wait_time',
            'avg_service_time',
            'exp_service_time',
            'avg_system_time',
            'exp_system_time',
        ]
        for s in summaries:
            print('{} {:.2f}'.format(s, getattr(self, s)()))

class Sim:
    def __init__(self, arrival_rate, service_rate, servers, debug=False):
        self.arrival_rate = arrival_rate
        self.service_rate = service_rate
        self.servers = servers
        self.debug = debug
        
        self.t = 0.0
        self.events = []
        self.waiting = 0
        self.servicing = 0
        
        self.stats = Stats(self)
        
    def schedule_arrival(self):
        delay = random.expovariate(self.arrival_rate)
        heappush(self.events, (self.t + delay, 'arrival'))
        
    def service(self):
        assert self.waiting > 0
        self.waiting -= 1
        self.servicing += 1
        delay = random.expovariate(self.service_rate)
        heappush(self.events, (self.t + delay, 'serviced'))
        
    def process_event(self):
        t, event = heappop(self.events)
        if self.debug:
            print(t, event)
        self.stats.waiting_time += (t - self.t) * self.waiting
        self.stats.servicing_time += (t - self.t) * self.servicing
        self.t = t
        if event == 'arrival':
            self.waiting += 1
            self.schedule_arrival()
            if self.servicing < self.servers:
                self.service()
        elif event == 'serviced':
            assert self.servicing > 0
            self.servicing -= 1
            if self.waiting > 0:
                self.service()
                
    def simulate(self, duration):
        self.schedule_arrival()
        while self.t < duration:
            self.process_event()
        return self.stats

Sim(arrival_rate=0.6, service_rate=1.0, servers=1).simulate(1e6).show()

avg_wait_time 0.90
exp_wait_time 0.90
avg_service_time 0.60
exp_service_time 0.60
avg_system_time 1.50
exp_system_time 1.50

Sim(arrival_rate=0.6, service_rate=0.5, servers=3).simulate(1e6).show()

avg_wait_time 0.09
exp_wait_time 0.09
avg_service_time 1.20
exp_service_time 1.20
avg_system_time 1.29
exp_system_time 1.29