'''

Problem 4.2




'''
import numpy as np
import MarkovLib as ml


def PjMM1(j, lmbda , mu ):
    rho = lmbda / mu
    return rho ** j * (1 - rho)

def PjMM1N(j, lmbda , mu , N ):
    rho = lmbda / mu
    if rho == 1:
        return 1 / (1 + N)
    else:
        return rho ** j * (1 - rho) / (1 - rho ** (N + 1))

def L_MM1(lmbda , mu ):
    rho = lmbda / mu
    return rho / (1 - rho)


def W_MM1(lmbda , mu ):
    rho = lmbda / mu
    return lmbda * rho / (1 - rho)


def PjMMNN(j , lmbda , mu , N ):
    if j > N: 
        return 0
    else:
        rho = lmbda / mu
        t = 1 
        h = t
        for i in range(1, N+1):
            t *=   rho / i
            h +=  t
        h = 1 / h
        for i in range(1,j+1):
            h *=   rho / i
        return h

def PjMMCN(j , lmbda , mu , C , N ):
    A = np.zeros([N+1,N+1])
    for i in range(N):  
        A[i, i + 1] = lmbda
        if i+1 > C:
            A[i + 1, i] = C * mu
        else:
            A[i + 1, i] = (i+1) * mu
    P =  ml.P_asymptotic(A)
    
    if j > N: 
        return 0
    else:
        return P[j]
 
def L_MMCN(lmbda , mu , C , N ):
    e = 0
    for j in range(N+1):
        e += j*PjMMCN(j, lmbda, mu, C, N)
    return e

def W_MMCN(lmbda , mu , C , N ):
    return lmbda * L_MMCN(lmbda, mu, C, N)

lmbda = 4  
mu = 60/25
C = 3 
N = 3 
print("\nN = 3 = maximum number of customers")
print("j   P_j(Direct formula)  P_j(Markov) ")
for j in range(N+1):
    print(j ,"     ",round(PjMMNN(j , lmbda , mu ,  N ),4), "            ",\
          round(PjMMCN(j , lmbda , mu , C , N ),4))

print("E(number in the system) = ", L_MMCN(lmbda,mu,C,N))
print("E(time in the system) = " , W_MMCN(lmbda,mu,C,N))
print("Rate of lost customers = ",PjMMCN(N,lmbda,mu,C,N)*lmbda )

N = 5
print("\nN = 5 = maximum number of customers")
print("j P_j(Markov) ")
for j in range(N+1):
    print(j , round(PjMMCN(j , lmbda , mu , C , N ),4))

print("E(number in the system) = ", L_MMCN(lmbda,mu,C,N))
print("E(time in the system) = " , W_MMCN(lmbda,mu,C,N))
print("Rate of lost customers = ",PjMMCN(N,lmbda,mu,C,N)*lmbda )


print("\nReduction in rate of lost customers = " , \
      PjMMCN(3,lmbda,mu,C,3)*lmbda -PjMMCN(5,lmbda,mu,C,5)*lmbda )

A = np.zeros([N+1,N+1])
for i in range(N):  
    A[i, i + 1] = lmbda
    if i+1 > C:
        A[i + 1, i] = C * mu
    else:
        A[i + 1, i] = (i+1) * mu
        
initState = 0
P_11 = ml.P_timeDependent(A,11-10,0.1,initState) 
P_15 = ml.P_timeDependent(A,15-10,0.1,initState)        
P_18 = ml.P_timeDependent(A,18-10,0.1,initState)
 
print("\nTime dependent solution\nj P_j(11:00) P_j(15:00) P_j(18:00)")
print("- ---------- ---------- ----------")
for j in range(N+1):
    print(j ," ", f"{P_11[j]:.4f}","   ",f"{P_15[j]:.4f}","   ",f"{P_18[j]:.4f}")
print("\nSteady state reach long before 15:00")