% The results should be compared to the PRCC results section in
% Supplementary Material D and Table D.1 for different N (specified by
% "runs" in the script below
clear all;
close all;

rng(1);
% Sample size N
runs=100;

% LHS MATRIX 
Parameter_settings_LHS;

% LHS_Call(xmin,xmean,xmax,xsd,nsample,distrib,threshold)
r_LHS = LHS_Call(1, r, 30, 0 ,runs,'unif'); 
rcw_LHS = LHS_Call(1, Nc, 5, 0 ,runs,'unif'); 
rdw_LHS = LHS_Call(1, Nd, 20, 0 ,runs,'unif'); 
m_LHS = LHS_Call(5, m, 20, 0 ,runs,'unif'); 
T1_LHS = LHS_Call(0.01, T1, 0.5 , 0, runs,'unif'); 
T2_LHS = LHS_Call(0.01, T2, 0.5, 0 ,runs,'unif'); 
bc1_LHS = LHS_Call(7, bc1, 25, 0 ,runs,'unif'); 
bc2_LHS = LHS_Call(7, bc2, 25, 0 ,runs,'unif'); 
bw1_LHS = LHS_Call(2, bw1, 25, 0, runs,'unif'); 
bw2_LHS = LHS_Call(2, bw2, 25, 0 ,runs,'unif'); 

% LHS MATRIX and PARAMETER LABELS
LHSmatrix=[r_LHS rcw_LHS rdw_LHS m_LHS T1_LHS T2_LHS bc1_LHS bc2_LHS bw1_LHS bw2_LHS];
D_lhs = zeros(1,runs); W_lhs = zeros(1,runs); 
for x=1:runs %Run solution x times choosing different values
    f=@ODE_LHS;
    x
    %LHSmatrix(x,:);
    options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
    [t,y]=ode45(@(t,y)f(t,y,LHSmatrix,x,runs),tspan,y0,options); 
    
    % Save only the outputs at the end points (equilibrium):
    % MORE EFFICIENT
    D_lhs(1,x)=y(end,1)/LHSmatrix(x,3); % diner fraction
    W_lhs(1,x)=y(end,2); % waiter fraction

end
% Save the workspace
save Model_LHS.mat;
%CALCULATE PRCC
alpha=0.01;
[prccD, signD, sign_labelD]=PRCC(LHSmatrix,D_lhs,1:length(time_points),PRCC_var,alpha);
[prccW, signW, sign_labelW]=PRCC(LHSmatrix,W_lhs,1:length(time_points),PRCC_var,alpha);

% figure; hist(D_lhs)
% figure; hist(W_lhs)
figure; scatter(D_lhs,W_lhs)
%PRCC_PLOT(LHSmatrix,V_T,1:length(time_points),PRCC_var,'y')