<jupyter_start><jupyter_text>Weighted Erlang Case Throughout this document, the following packages are required:<jupyter_code>import numpy as np
import scipy
import math
from scipy.stats import binom, erlang, poisson
from scipy.optimize import minimize
import csv

with open('phase/SCV_1.00_minima_n5.csv','r') as csvfile:
    
    reader = csv.reader(csvfile)
    old_minima = list(reader)
#     for row in reader:
#         print(','.join(row))
type(eval(old_minima[0][0])[99])<jupyter_output><empty_output><jupyter_text>Weighted Erlang distributionIn this section, we assume that the service time $B$ equals w.p. $p\in[0,1]$ an Erlang-distributed random variable with $K$ exponentially distributed phases, each of them having mean $\mu^{-1}$, and with probability $1-p$ an Erlang-distributed random variable with $K+1$ exponentially distributed phases, again with mean $\mu^{-1}$:\begin{align*}B \stackrel{\text{d}}{=} \sum_{i=1}^{K}X_i + X_{K+1}\mathbb{1}_{\{U > p\}},\end{align*}where $X_i \stackrel{iid}{\sim} \text{Exp}(\mu)$ and $U\sim\text{Unif}[0,1]$.<jupyter_code>def SCV_to_params(SCV):
    
    # weighted Erlang case
    if SCV <= 1:
        K = math.floor(1/SCV)
        p = ((K + 1) * SCV - math.sqrt((K + 1) * (1 - K * SCV))) / (SCV + 1)
        mu = K + (1 - p)
    
        return K, p, mu
    
    # hyperexponential case
    else:
        p = 0.5 * (1 + np.sqrt((SCV - 1) / (SCV + 1)))
        mu = 1 # 1 / mean
        mu1 = 2 * p * mu
        mu2 = 2 * (1 - p) * mu
        
        return p, mu1, mu2<jupyter_output><empty_output><jupyter_text>We keep the following parameters fixed:<jupyter_code>n = 5
omega = 0.5 ###
Delta = 0.01
K, p, mu = SCV_to_params(0.99)
t_MAX = int(5/Delta)
M = 500

epsilon = 0.005<jupyter_output><empty_output><jupyter_text>The recursion of the dynamic program is given as follows. For $i=1,\dots,n-1$, $k=1,\dots,i$, and $m\in\mathbb{N}_0$,\begin{align*}\xi_i(k,m) &= \inf_{t\in \mathbb{N}_0}\Big(\omega \bar{f}^{\circ}_{k,m\Delta}(t\Delta) + (1-\omega)\bar{h}^{\circ}_{k,m\Delta} +\sum_{\ell=2}^{k}\sum_{j=0}^{t}\bar{q}_{k\ell,mj}(t)\xi_{i+1}(\ell,j) +P^{\downarrow}_{k,m\Delta}(t\Delta)\xi_{i+1}(1,0) +P^{\uparrow}_{k,m\Delta}(t\Delta)\xi_{i+1}(k+1,m+t)\Big),\end{align*}whereas, for $k=1,\dots,n$ and $m\in \mathbb{N}_0$,\begin{align*}\xi_n(k,m) = (1-\omega)\bar{h}^{\circ}_{k,m\Delta}.\end{align*}<jupyter_code>from functools import lru_cache
# # compute the poisson pmfs and cdfs
# poisson_pmf = np.zeros((K+1,t_MAX**2+1), dtype=np.float64)
# poisson_cdf = np.zeros((K+1,t_MAX**2+1), dtype=np.float64)

# for z in range(K+1):
#     for t_index in range(t_MAX+1):
#         t = t_index*Delta
#         poisson_pmf[z][t_index] = poisson.pmf(z, mu*t)
#         poisson_cdf[z][t_index] = poisson.cdf(z, mu*t)

# # # methode 1
# # B_sfs = np.zeros(t_MAX**2+1)
# # for t_index in range(t_MAX**2+1):
# # #     B_sfs[t_index] = B_sf(t_index*Delta)
# #     B_sfs[t_index] = poisson_cdf[K-1][t_index] + (1 - p) * poisson_pmf[K][t_index]
    
# # methode 2
# # B_sfs = np.zeros(t_MAX**2+1)
# poisson_pmfs = [np.zeros(t_MAX**2+1) for z in range(K+1)]
# poisson_cdfs = [np.zeros(t_MAX**2+1) for z in range(K+1)]

# # for t_index in range(t_MAX**2+1):
# #     t = t_index*Delta
# #     cdf, pmf = poisson.cdf(K-1, mu*t), poisson.pmf(K, mu*t)
    
# #     poisson_pmfs[K][t_index] = pmf
# #     B_sfs[t_index] = cdf + (1 - p) * pmf

# for z in range(n):
#     for t_index in range(t_MAX**2+1):
#         t = t_index*Delta
        
#         poisson_pmfs[z][t_index] = poisson.pmf(z, mu*t)
#         poisson_cdfs[z][t_index] = poisson.cdf(z, mu*t)

# B_sfs = poisson_cdfs[K-1] + (1 - p) * poisson_pmfs[K]
# B_sfs[-240000]
# B_sfs_alt[-240000]
# poisson_pmfs[0][-240000]
# print(np.mean(B_sfs),np.std(B_sfs))
# print(np.mean(B_sf),np.std(B_sf))

# poisson_pmf
### TRANSITION PROBABILITIES ###

# 1. No client has been served before time t.
@lru_cache(maxsize=128)
def B_sf(t):
    """The survival function P(B > t)."""
#     try:
#         t_index = int(t/Delta)
#         return B_sf[t_index]
#     except:
    return poisson.cdf(K-1, mu*t) + (1 - p) * poisson.pmf(K, mu*t)
        

@lru_cache(maxsize=128)
def P_up(k, u, t):
    """Computes P(N_t- = k | N_0 = k, B_0 = u)."""
#     try:
#         u_index, t_index = int(u/Delta), int(t/Delta)
#         return B_sf[u_index+t_index] / B_sf[u_index]
#     except:
    return B_sf(u+t) / B_sf(u)


# 2. All clients have been served before time t.
@lru_cache(maxsize=128)
def gamma(z, t):
    """Computes P(Z_t = z | B > t)."""
#     try:
#         t_index = int(t/Delta)
#         gamma_circ = poisson_pmfs[z-1][t_index]

#         if z == K + 1:
#             gamma_circ *= (1 - p)

#         return gamma_circ / B_sfs[t_index]
# #     except:
    gamma_circ = poisson.pmf(z-1, mu*t)

    if z == K + 1:
        gamma_circ *= (1 - p)

    return gamma_circ / B_sf(t)

@lru_cache(maxsize=128)
def P_k0(k, z, t):
    """Computes P(N_t- = 0 | N_0 = k, X_0 = z)."""
    if z <= K:
        return sum([binom.pmf(m, k, 1-p) * erlang.cdf(t, k*K-z+1+m, scale=1/mu) for m in range(k+1)])
    elif z == K + 1:
        return sum([binom.pmf(m, k-1, 1-p) * erlang.cdf(t, (k-1)*K+1+m, scale=1/mu) for m in range(k)])

# @lru_cache(maxsize=128)
def P_down(k, u, t):
    """Computes P(N_t- = 0 | N_0 = k, B_0 = u)."""
    return sum([gamma(z, u) * P_k0(k, z, t) for z in range(1, K+2)])


# 3. Some (but not all) clients have been served before time t.
@lru_cache(maxsize=128)
def psi(v, t, k, l):
    """
    Computes P(t-v < Erl(k,mu) < t, Erl(k,mu) + Erl(l-k,mu) > t),
    where Erl(k,mu) and Erl(l-k,mu) are independent.
    """
#     try:
#     t_index = int(t/Delta)
#     return sum([poisson_pmfs[j][t_index] * binom.sf(j-k, j, v/t) for j in range(k, l)])
#     except:
    return sum([poisson.pmf(j, mu*t) * binom.sf(j-k, j, v/t) for j in range(k, l)])


@lru_cache(maxsize=128)
def q(k, l, z, v, t):
    """Computes P(N_t = l, B_t < v | N_0 = k, Z_0 = z)."""

    q = 0

    if z <= K:
        for m in range(k-l+2):
            I_klmz = (k - l + 1) * K - z + m + 1
            E = p * psi(v, t, I_klmz, I_klmz+K) + (1 - p) * psi(v, t, I_klmz, I_klmz+K+1)
            q += binom.pmf(m, k-l+1, 1-p) * E

    elif z == K + 1:
        for m in range(k-l+1):
            I_klm = (k - l) * K + m + 1
            E = p * psi(v, t, I_klm, I_klm+K) + (1 - p) * psi(v, t, I_klm, I_klm+K+1)
            q += binom.pmf(m, k-l, 1-p) * E

    return q

def trunc(x, y):
    """Truncates x to the interval [0,y]."""
    return min(max(0, x), y)

@lru_cache(maxsize=128)
def q_bar(k, l, m, j, t):
    """Approximates P(N_{t*Delta} = l, B_{t*Delta} \in d(j*Delta) | N_0 = k, B_0 = m * Delta)."""

    lower, upper = trunc((j-0.5)*Delta, t*Delta), trunc((j+0.5)*Delta, t*Delta)
    q_bar = sum([gamma(z, m*Delta) * (q(k, l, z, upper, t*Delta) - q(k, l, z, lower, t*Delta)) for z in range(1, K+2)])

    return q_bar
### OBJECTIVE/RISK FUNCTION ###

@lru_cache(maxsize=128)
def f(k, t):
#     try:
#     t_index = int(t/Delta)
#     print(t_index,k-1)
#     return (1 - poisson_cdfs[k-1][t_index]) * t - (1 - poisson_cdfs[k][t_index]) * k / mu
#     except:
    return poisson.sf(k-1, mu*t) * t - poisson.sf(k, mu*t) * k / mu

@lru_cache(maxsize=128)
def f_bar(k, z, t):  # TODO: deze kan sneller door wat op te slaan
    if z <= K:
        return sum([binom.pmf(m, k, 1 - p) * f(k*K+1-z+m, t) for m in range(k+1)])
    elif z == K + 1:
        return sum([binom.pmf(m, k-1, 1 - p) * f((k-1)*K+1+m, t) for m in range(k)])

@lru_cache(maxsize=128)
def f_circ(k, u, t):
    return sum([gamma(z, u) * f_bar(k, z, t) for z in range(1, K+2)])

@lru_cache(maxsize=128)
def h_bar(k, z):

    if k == 1:
        return 0
    elif z <= K: ################### TODO
        return ((k - 1) * (K + 1 - p) + 1 - z) / mu
    elif z == K + 1:
        return ((k - 2) * (K + 1 - p) + 1) / mu

@lru_cache(maxsize=128)
def h_circ(k, u):
    return sum([gamma(z, u) * h_bar(k, z) for z in range(1, K+2)])
### COST ###
def cost_we(t, i, k, m):
    """Computes (approximately) the cost when t is the next interarrival time."""
    
    m = int(round(m))
    
    cost = omega * f_circ(k, m*Delta, t*Delta) + (1 - omega) * h_circ(k, m*Delta)
    cost += P_down(k, m*Delta, t*Delta) * xi_we(i+1, 1, 0)
    cost += P_up(k, m*Delta, t*Delta) * xi_we(i+1, k+1, m+t)

    for l in range(2, k+1):
        for j in range(t+1):
            cost += q_bar(k, l, m, j, t) * xi_we(i+1, l, j)
    
    return cost


def xi_we(i, k, m):
    """Implements the Weighted Erlang Case."""
    
#     print("begin",i,k,m)
    
    # truncate time in service m
    if m > M:
        m = M
    
    if xi_matrix[i-1][k-1][m-1]:  # retrieve stored value
        pass
    elif i == n:
        xi_matrix[i-1][k-1][m-1] = (1 - omega) * h_circ(k, m*Delta)        
    else:
        
        if m >= 3 and xi_matrix[i-1][k-1][m-2] and xi_matrix[i-1][k-1][m-3]:
#             print(i,k,m)
#             print(xi_matrix[i-1][k-1])
            if abs(xi_matrix[i-1][k-1][m-2] - xi_matrix[i-1][k-1][m-3]) < epsilon:
                xi_matrix[i-1][k-1][m-1:] = [xi_matrix[i-1][k-1][m-2]] * (M - (m - 1))
                minima[i-1][k-1][m-1:] = [minima[i-1][k-1][m-2]] * (M - (m - 1))
            
                print(i,k,m,"break")
                return xi_matrix[i-1][k-1][m-1]
            
                ### TODO: fill all coming values with current cost & minimum
                    
        t_guess = eval(old_minima[i-1][k-1])[m-1]
        cost_guess = cost_we(t_guess, i, k, m)

        t_new = t_guess

        while True:
#             print(i, k, m, t_guess, cost_guess)

            t_new -= 1
            cost_new = cost_we(t_new, i, k, m)
            if cost_new < cost_guess:
                t_guess = t_new
                cost_guess = cost_new
            else:
                break

        xi_matrix[i-1][k-1][m-1] = cost_guess
        minima[i-1][k-1][m-1] = t_guess

        print("eind", i, k, m, t_guess, cost_guess)

    return xi_matrix[i-1][k-1][m-1]
xi_matrix = [[[None for m in range(t_MAX)] for k in range(i+1)] for i in range(n)]
minima = [[[None for m in range(t_MAX)] for k in range(i+1)] for i in range(n)]

# %%cython inline
# print("Cost:", xi_we(1, 1, 0))
for i in np.arange(n,0,-1):
    for k in range(1,i+1):
        print("i =",i,"k =",k)
        for m in range(t_MAX):
            xi_we(i,k,m)
eind 14 1 0 69 0.3465760690660555
eind 14 2 86 168 1.0258562318250255
eind 14 2 85 168 1.0258562318250237
eind 13 1 0 86 0.7756027183751549
eind 14 3 279 267 1.659409653583836
eind 14 2 0 168 1.0258562318250244
eind 14 2 1 168 1.0258562318250255
eind 14 2 2 168 1.0258562318250244
14 2 3 break
eind 14 3 278 267 1.659409653583836
eind 13 2 88 191 1.477601341587333
import csv

with open('phase/SCV_0.90_minima.csv','w', newline='') as myfile:
    out = csv.writer(myfile)
    out.writerows(minima)

    with open('phase/SCV_0.90_minima.csv','w', newline='') as myfile:
    out = csv.writer(myfile)
    out.writerows(minima)

xi_matrix[14][1]
## t op 0.01 afronden
## SCV op 0.05 afronden
## m op 0.01 afronden
# u_max = 10**2 ### IDEE: Gamma's opslaan
# gamma_matrix = [[None for z in range(K+1)] for u in range(u_max)]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19166220317800392, 'fraction_numerical': 0.05284770576682735, 'mean_word_length': 2.636157091853294, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '443303', 'license': 'MIT', 'n_tokens_mistral': 5353, 'n_tokens_neox': 4987, 'n_words': 1142, 'path': '.ipynb_checkpoints/Phase-Type Case-checkpoint.ipynb', 'repo_name': 'Roshanmahes/Appointment-Scheduling', 'size': '11202'}