function SIM = TRI_CY_run_the_sim(lambda_scope, lambda_count, T_scope, ...
    T_count, slots, d_mu_scope, d_mu_count, d_var_scope, ...
    d_var_count, prob_list)

% TRI_CY_run_the_sim - a simulator getting sim result about
% priority method
% On input:
%     lambda_scope (1 * 2 vector): arrival rate of ED urgent cases
%     lambda_count(int): how many lambda the user want to separate
%     T_scope (1 * 2 vector): max time for simulation period
%     T_count(int): how many T the user want to separate
%     number_bins(int): how many bins that user want (just support 2 for
%     now)
%     slots (vector): time bin for priority
%     d_mu_scope (1 * 2 vector): mean duration of surgery
%     d_mu_count(int): how many d_mu the user want to separate
%     d_var_scope (1 * 2 vector): variance of duration of surgery
%     d_var_count(int): how many d_var the user want to separate
%     priority_count(int): how many priority the user want to separate
% On output:
%     SIM (nx6 array): Simulation about priority method
%       col 1: lambda
%       col 2: T
%       col 3: d_mu
%       col 4: d_var
%       col 5: priority method delay patients number
%       col 6: 1C1S method delay patients number
% Call:
%      lambda_scope = [0.1, 3];
%      T_scope = [24, 192];
%      slots = [0, 3; 3, 12; 12, 24];
%      d_mu_scope = [0.5, 3];
%      d_var_scope = [0, 1];
%      prob_list = [0.05; 0.6; 0.35];
%      SIM = TRI_CY_run_the_sim(lambda_scope, 30, T_scope, 8, slots,
%      d_mu_scope, 6, d_var_scope, 11, prob_list);
% Author:
%     Chengyu Yang
%     UU
%     Fall 2023
%

lambda_min = lambda_scope(1, 1);
lambda_max = lambda_scope(1, 2);
T_min = T_scope(1, 1);
T_max = T_scope(1, 2);
d_mu_min = d_mu_scope(1, 1);
d_mu_max = d_mu_scope(1, 2);
d_var_min = d_var_scope(1, 1);
d_var_max = d_var_scope(1, 2);

lambda_list = linspace(lambda_min, lambda_max, lambda_count);
T_list = linspace(T_min, T_max, T_count);
d_mu_list = linspace(d_mu_min, d_mu_max, d_mu_count);
d_var_list = linspace(d_var_min, d_var_max, d_var_count);

priority_prob_list = prob_list;
SIM = [];

for lambda = lambda_list
    for T = T_list
        for d_mu = d_mu_list
            for d_var = d_var_list
                for priority_prob = priority_prob_list

                    U_SORT_delay_number_min = 0;
                    U_SORT_delay_number_fair = 0;
                    % U_SORT_delay_number_combine = 0;
                    FCFS_delay_number = 0;
                    Total_Patient_number = 0;

                    rng('default');
                    TB = slots;

                    for times = 1:10

                        patients = TRI_CY_arrivals(lambda,T,priority_prob,d_mu,d_var,0.5);

                        p_queue = TRI_CY_p_queue_method(patients,TB, 0, 10);
                        U_SORT_EMR = TRI_CY_p_queue_to_EMR(patients, p_queue, TB);
                        U_SORT_delay_number_min = U_SORT_delay_number_min + numel(find(U_SORT_EMR(:, 6) > 0));

                        p_queue = TRI_CY_p_queue_method(patients,TB, 1, 10);
                        U_SORT_EMR = TRI_CY_p_queue_to_EMR(patients, p_queue, TB);
                        U_SORT_delay_number_fair = U_SORT_delay_number_fair + numel(find(U_SORT_EMR(:, 6) > 0));


                        FCFS_EMR = TRI_CY_sim_1C1S(patients, slots, 0, 24, 100);
                        FCFS_delay_number = FCFS_delay_number + numel(find(FCFS_EMR(:, 6) > 0));
                        
                        [row_count, ~] = size(patients);
                        Total_Patient_number = Total_Patient_number +  row_count;

                    end

                    delay_percentage = [U_SORT_delay_number_min,U_SORT_delay_number_fair, FCFS_delay_number] / Total_Patient_number;
                    
                    SIM = [SIM; lambda, T, d_mu, d_var, delay_percentage, Total_Patient_number];
                    
                end
            end
        end
    end
end
