Part 2

import numpy as np
import scipy.linalg as slin
import fractions
import math
from math import factorial, e, ceil

# Probability of Waiting
def prob_waiting(traffic_intesity, num_agents):
    x = ((traffic_intesity ** num_agents) / factorial(round(num_agents))) * num_agents / (num_agents - traffic_intesity)
    y = 1

    for i in range(round(num_agents)):
        y += (traffic_intesity ** i) / factorial(i)

    return x / (y + x)

# Erlang C formula as rho lambda function
erlang_c = lambda rho, s: (
    (s * (rho ** s) / math.factorial(s)) /
    ((s * (rho ** s) / math.factorial(s)) + sum((rho ** n) / math.factorial(n) for n in range(s)))
    if s > 0 else 0
)


def custom_function(mu, P0, lam, s, t, rho):
    a = rho / s
    exp_term1 = math.exp(-mu * t)
    factorial_term = (lam / mu) ** s / (math.factorial(s) * (1 - a))
    exp_term2 = 1 - math.exp(-mu * (s - 1 - lam / mu) * t)
    denominator = s - 1 - lam / mu
    return exp_term1 * (1 + P0 * factorial_term * exp_term2 / denominator)

np.set_printoptions(formatter={'all':lambda x: str(fractions.Fraction(x).limit_denominator())})

s = 3
mu = 1/8
u = lambda i: i*mu if i <3 else 3*mu
lambda_ = 1/3
rho = lambda_/mu
a = rho/s

pi0 = np.exp(-lambda_/mu)

A = np.zeros((8,8))

if __name__ == "__main__":
    for i in range(8):
        if(i == 0):
            A[i][0] = -lambda_
            A[i][i+1] = lambda_

        if(i == 7):
            A[i][0] = 0

        else:
            A[i][i-1] = u(i)
            A[i][i] = -(u(i) + lambda_)
            A[i][i+1] = lambda_

    P0 = 0
    for i in range(s):
        P0 += rho**i/math.factorial(i)

    P0 += rho**s/math.factorial(s)*1/(1-a)
    P0 = P0**(-1)

    P = slin.expm(A*24)
    print(f"pi0: {P0}")
    print(f"Question 7 probability matrix: \n{P}\n - with the final answer being 140737/993654")
    L0 = (P0/(1+2+3)) * (lambda_/mu)**3 * ((lambda_/u(3))/ (1- lambda_/u(3))**2)
    L = L0 + lambda_/mu
    print(f"Average number of patients: {L}")

    S8 = A[:-1, :-1]
    U = np.linalg.inv(-S8)
    alpha = [1,0,0,0,0,0,0]
    EX8 = alpha@U@np.identity(7)
    print(f"Total expected time until more than 6 guests in question 8: {np.sum(EX8)}")
    U2 = U@U
    print(U2.shape)
    alpha2 = np.array(alpha)*2
    EX2 = alpha2@U2@np.identity(7)
    print(f"We have that E[X^2]: {np.sum(EX2)} and E[X]^2: {np.sum(EX8)**2}, such that V[X] = {np.sum(EX2)-np.sum(EX8)**2}")




    #mu, P0, lam, s, t, rho):
    print(custom_function(mu,P0,lambda_,s,1,rho))

    GT = lambda t: np.exp(-mu*t)*(1
                                  +
                                  (rho**s*P0*(1-np.exp(-mu*(s - 1 - rho)*t)))
                                  /
                                  (math.factorial(s)*(1-a)*(s - 1 - rho))
                                  )



    FW = lambda t: 1 - erlang_c(rho,s) * np.exp(-mu*(s-rho)*t)
    #TODO: I stedet for !waiting så prøv at sæt pi0,pi1,pi2 ind.
    FT = lambda t: (1-erlang_c(rho,s))*(1-np.exp(-mu*t)) + FW(t)

    print(GT(24))

    #M/G/1 queue values - Question 11
    E_S = 1/10
    V_S = 1/4800
    lambda_mg = 9
    rho_mg = lambda_mg*E_S
    E_S2 = V_S + E_S ** 2
    W = (lambda_mg * E_S2) / (2 * (1 - rho_mg))
    T = W + E_S
    k = T/3
    E_PT = lambda_mg * W + rho_mg
    print(f"k is equal to: {k}\n")
    print(f"Average number of morbidly obese patients: {k*E_PT}")

    #Question 12 MGF
    alpha_gamma = 4/3
    chi_gamma = 4.5
    E_J3 = (alpha_gamma*(alpha_gamma + 1)*(alpha_gamma +2))/(chi_gamma ** 3)
    E_S3 = (5/60)**3 + E_J3/1000

    #We split up the terms to make it more readable
    VPT_1 = (lambda_mg**3 * E_S3)/(3*(1-rho_mg))
    VPT_2 = (((lambda_mg ** 2 * E_S2))/(2*(1-rho_mg)))**2
    VPT_3 = ( lambda_mg**2 * (3 - 2 * rho_mg) * E_S2 )/(2*(1 - rho_mg))  + rho_mg*(1-rho_mg)
    VPT = VPT_1 + VPT_2 + VPT_3

    print(f"Answers to question 12\n E[J^3]: {E_J3} \n E[S^3]: {E_S3} \n V[P_T]: {VPT} \n V[M(PT)]: {k**2 *VPT}")