%
% These examples are based on the package documentation:
% http://www.ctan.org/tex-archive/macros/latex/contrib/minted
%
\documentclass{article}

\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage{minted}
\usepackage{courier}

\begin{document}

\title{Source Code of Numerical Experimentation}
\author{}
\maketitle

\section{Date Generating Process}
For better reproducibility, we use the following code to generate the Poisson arrivals to simulate the performances of different admission policies. The scaling factor $N$ ranges from $\{100, 200, \ldots, 1500\}$, and we set the time horizon to be quadratic in $N$ for large $N$ and linear in $N$ with a large multiplication factor when $N$ is small.

\begin{minted}{python}
import numpy as np
import os

# Parameters
Lambda = np.array([3, 2, 5])
mu_inv = np.array([2, 1, 3])
n = 3   # number of customer types
d = 2   # number of customer types

for N in Nvec:
    K = 20 * N;
    epochN = max(N, 1000);

    ArrivalTimes = np.empty((n, K))
    ServiceTimes = np.empty((n, K))
    customers = range(n)
    Time_vec = np.zeros(epochN)

    # Create folder if not exist
    folder = "Data/N=" + str(N)
    if not os.path.exists(folder):
        os.makedirs(folder)

    # Generate data
    for k in range(epochN):
        for i in customers:
            InterArrivalTimes = np.random.exponential(1/N/Lambda[i], K)
            ArrivalTimes[i] = np.cumsum(InterArrivalTimes)
            ServiceTimes[i] = np.random.exponential(mu_inv[i], K)
        T = np.amin(ArrivalTimes[:, -1])
        if k == 0:
            Time_vec[k] = T
        else:
            Time_vec[k] = Time_vec[k-1] + T
        ArrivalTimes = np.where(ArrivalTimes<T, ArrivalTimes, 0)
        # Combine arrival times and service times
        TimeChain = np.array([])
        for i in customers:
            Chain = ArrivalTimes[i, :]
            Chain = Chain[Chain != 0]
            # first row arrival time, second row departure time, 
            # third row type
            Chain = np.vstack((Chain, Chain+ServiceTimes[i, 0:Chain.size], 
                    (i+1)*np.ones((1, Chain.size))))
            if TimeChain.size:
                TimeChain = np.hstack((TimeChain, Chain))
            else:
                TimeChain = Chain
        # sort by arrival time
        TimeChain = TimeChain[:, TimeChain[0, :].argsort()]   

        # Save data of this epoch
        file_name = "Data/N=" + str(N) + "/Epoch-" + str(k) + '.npz'
        np.savez(file_name, TimeChain=TimeChain,T=T)
\end{minted}

For Network 2 and Network 3, we respectively adjust the parameters in the code to:
\begin{minted}{python}
# Parameters
Lambda = np.array([2, 3, 5, 1, 6, 2, 3])
mu_inv = np.array([1, 3, 2, 3, 5, 4, 2])
n = 7   # number of customer types
d = 4   # number of customer types
\end{minted}
and
\begin{minted}{python}
# Parameters
Lambda = np.array([2, 2, 2, 2, 3, 4])
mu_inv = np.array([1, 1, 1, 2, 3, 4])
n = 6   # number of customer types
d = 3   # number of resource types
\end{minted}



\section{Simulation of Policies}
After generating the arrivals (saved as ".npz" files) using the above code, we use the following code (and variant) to simulate the average reward of sequential policies with respect to the scaling factor $N \in \{100, 200, \ldots, 1500\}$. The following code piece simulates the total reward collected under the corrected headcount threshold policy (CHT) with scaling factor $N=100$ in Network 3 of Section 7 of the paper.
We run the same simulation code with $N\in\{100, 200, 1500\}$ and save the simulation results to plot the scaling curve.

\begin{minted}{python}
import numpy as np

#########################
# Scaling factor
N = 100
#########################
# Parameters
LP = 7.5
n = 6
d = 3
Lambda = np.array([2, 2, 2, 2, 3, 4])
mu_inv = np.array([1, 1, 1, 2, 3, 4])
q = np.array([2, 2, 2])
c = np.array([20, 20, 20])
r = np.array([1, 1, 1, 4, 6, 8])

epochN = max(N, 1000)
Time_vec = np.zeros(epochN)
AverageReward_CHT = np.zeros(epochN)
postState_CHT = np.zeros(n)

# Simulation epochs
for k in range(epochN):
    # Load data
    file_name = "Data/N=" + str(N) + "/Epoch-" + str(k) + '.npz'
    data = np.load(file_name)
    TimeChain = data['TimeChain']
    T = data['T']

    if k == 0:
        Time_vec[k] = T
    else:
        Time_vec[k] = Time_vec[k-1] + T

    arrivalN = TimeChain.shape[1]       # total number of arrivals

    # Initialization
    N_hat = N * q - np.log(N) * c
    if k == 0:
        # Corrected Headcount Threshold Policy (CHT)
        X_CHT = np.zeros(n)
        preN_CHT = 0
        preChain_CHT = np.array([])
    else:
        # Corrected Headcount Threshold Policy (CHT)
        X_CHT = postState_CHT
        preN_CHT = int(np.sum(X_CHT))
        preChain_CHT = postChain_CHT

    # Simulations
    TotalReward_CHT = 0
    Decisions_CHT = np.zeros(arrivalN)

    for i in range(arrivalN):
        arriving_time = TimeChain[0, i]
        arriving_type = int(TimeChain[2, i])

        #########################
        # Corrected Headcount Threshold Policy (CHT)
        # Update headcounts
        if preChain_CHT.size:
            Truncation = np.hstack((preChain_CHT, TimeChain[:, 0:i]))
        else:
            Truncation = TimeChain[:, 0:i]
        if Truncation.size > 0:
            departure_vec = (Truncation[1, :] > arriving_time)
            accept_vec = np.hstack([np.ones(preN_CHT), Decisions_CHT[0:i]])
            index = departure_vec * accept_vec
            type_vec = Truncation[2, index > 0]
            for j in range(n):
                X_CHT[j] = np.count_nonzero(type_vec == j + 1)
        #########################
        # Admission control
        # Y are headcounts
        Y_a = X_CHT[0] + X_CHT[2] + X_CHT[3]
        Y_b = X_CHT[0] + X_CHT[1] + X_CHT[4]
        Y_c = X_CHT[1] + X_CHT[2] + X_CHT[5]
        # Z are headcount deviations
        Z_a = N_hat[0] - Y_a
        Z_b = N_hat[1] - Y_b
        Z_c = N_hat[2] - Y_c

        if arriving_type == 1:
            if (Z_a > 0) and (Y_b < N * q[1]):
                decision = 1
            else:
                decision = 0
        elif arriving_type == 2:
            if (Z_b > 0) and (Y_c < N * q[2]):
                decision = 1
            else:
                decision = 0
        elif arriving_type == 3:
            if (Z_c > 0) and (Y_a < N * q[0]):
                decision = 1
            else:
                decision = 0
        elif arriving_type == 4:
            if Y_a < N*q[0]:
                decision = 1
            else:
                decision = 0
        elif arriving_type == 5:
            if Y_b < N*q[1]:
                decision = 1
            else:
                decision = 0
        elif arriving_type == 6:
            if Y_c < N*q[2]:
                decision = 1
            else:
                decision = 0

        Decisions_CHT[i] = decision
        TotalReward_CHT += r[arriving_type - 1] * decision

    # Update reward vector
    if k == 0:
        AverageReward_CHT[k] = TotalReward_CHT / T
    else:
        AverageReward_CHT[k] = (AverageReward_CHT[k - 1] * Time_vec[k - 1] 
                                + TotalReward_CHT) / Time_vec[k]

    # Update remaining chain for next epoch
    # Corrected Headcount Threshold Policy (CHT)
    index = (TimeChain[1, :] > T) * (Decisions_CHT == 1)
    if k == 0:
        postChain_CHT = TimeChain[:, index > 0]
    else:
        postChain_CHT = np.hstack((preChain_CHT[:, preChain_CHT[1, :] > T], 
                        TimeChain[:, index > 0]))
    postChain_CHT[0:2, :] = postChain_CHT[0:2, :] - T
    for j in range(n):
            postState_CHT[j] = np.count_nonzero(postChain_CHT[2, :] == j+1)

# Save results
file_name = "N=" + str(N) + '.npz'
np.savez(file_name, Time_vec=Time_vec,AverageReward_CHT=AverageReward_CHT, 
        LP=LP, N=N, c=c, epochN=epochN)
\end{minted}

The simulation of the other two networks in Section 7 follows the same fashion with minimum code change (parameters and policy actions).
The above code will save the cumulative average rewards collected at different time points along the entire time horizon, which are then used to create plots in the paper. 

\end{document}