% Experiment 1

d_list = [20 50 100];
s_list = [5 10 20];
p_list = [0.1 0.3 0.5 0.7 0.9];

for s = s_list
    for p = p_list
        simulation1(100, s, p, 30);
    end
end

for d = d_list
    for p = p_list
        simulation1(d, 10, p, 30);
    end
end

res = readtable('simulation1.csv');
res = res{:,:};

figure;
tiledlayout(1,2);
nexttile
d = 100;
s = 5; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
plot(p, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
s = 10; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
plot(p, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
s = 20; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
plot(p, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('p', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13); 
ylim([0, 1]); xlim([0.1, 0.9]); xticks([0.1 0.3 0.5 0.7 0.9]);
title('Exact recovery rates when \it d = 100')

nexttile
d = 100;
s = 5; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
bar_lambda = X(:, 5); max2 = X(:, 6); max3 = X(:, 7);
mu1 = X(:, 8); mu2 = X(:, 9); mu3 = X(:, 10);
rescaled = min([mu1*sqrt(log(s)), ...
    (bar_lambda.*mu2)./(max2*max(s^(2.5), s*sqrt(s*(d-s)))), ...
    (bar_lambda.*mu3)./(max3*sqrt(log(d-s)))],[],2);
rescaled = sqrt(p./(1-p)).*rescaled;
plot(rescaled, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
s = 10; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
bar_lambda = X(:, 5); max2 = X(:, 6); max3 = X(:, 7);
mu1 = X(:, 8); mu2 = X(:, 9); mu3 = X(:, 10);
rescaled = min([mu1*sqrt(log(s)), ...
    (bar_lambda.*mu2)./(max2*max(s^(2.5), s*sqrt(s*(d-s)))), ...
    (bar_lambda.*mu3)./(max3*sqrt(log(d-s)))],[],2);
rescaled = sqrt(p./(1-p)).*rescaled;
plot(rescaled, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
s = 20; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
bar_lambda = X(:, 5); max2 = X(:, 6); max3 = X(:, 7);
mu1 = X(:, 8); mu2 = X(:, 9); mu3 = X(:, 10);
rescaled = min([mu1*sqrt(log(s)), ...
    (bar_lambda.*mu2)./(max2*max(s^(2.5), s*sqrt(s*(d-s)))), ...
    (bar_lambda.*mu3)./(max3*sqrt(log(d-s)))],[],2);
rescaled = sqrt(p./(1-p)).*rescaled;
plot(rescaled, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('Rescaled', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13);
title('Exact recovery rates when \it d = 100')

lgd = legend({'s=5', 's=10', 's=20'},...
    'location', 'southwest', 'FontSize', 13);
lgd.Layout.Tile = 'south';

figure;
tiledlayout(1,2);
nexttile
s = 10;
d = 20; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
plot(p, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
d = 50; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
plot(p, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
d = 100; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
plot(p, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('p', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13); 
ylim([0, 1]); xlim([0.1, 0.9]); xticks([0.1 0.3 0.5 0.7 0.9]);
title('Exact recovery rates when \it s = 10')

nexttile
s = 10;
d = 20; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
bar_lambda = X(:, 5); max2 = X(:, 6); max3 = X(:, 7);
mu1 = X(:, 8); mu2 = X(:, 9); mu3 = X(:, 10);
rescaled = min([mu1*sqrt(log(s)), ...
    (bar_lambda.*mu2)./(max2*max(s^(2.5), s*sqrt(s*(d-s)))), ...
    (bar_lambda.*mu3)./(max3*sqrt(log(d-s)))],[],2);
rescaled = sqrt(p./(1-p)).*rescaled;
plot(rescaled, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
d = 50; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
bar_lambda = X(:, 5); max2 = X(:, 6); max3 = X(:, 7);
mu1 = X(:, 8); mu2 = X(:, 9); mu3 = X(:, 10);
rescaled = min([mu1*sqrt(log(s)), ...
    (bar_lambda.*mu2)./(max2*max(s^(2.5), s*sqrt(s*(d-s)))), ...
    (bar_lambda.*mu3)./(max3*sqrt(log(d-s)))],[],2);
rescaled = sqrt(p./(1-p)).*rescaled;
plot(rescaled, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
d = 100; X = res((res(:,1)==d)&(res(:,2)==s), :);
p = X(:, 3); er_rate = X(:, 4);
bar_lambda = X(:, 5); max2 = X(:, 6); max3 = X(:, 7);
mu1 = X(:, 8); mu2 = X(:, 9); mu3 = X(:, 10);
rescaled = min([mu1*sqrt(log(s)), ...
    (bar_lambda.*mu2)./(max2*max(s^(2.5), s*sqrt(s*(d-s)))), ...
    (bar_lambda.*mu3)./(max3*sqrt(log(d-s)))],[],2);
rescaled = sqrt(p./(1-p)).*rescaled;
plot(rescaled, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('Rescaled', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13);
title('Exact recovery rates when \it s = 10')

lgd = legend({'d=20', 'd=50', 'd=100'},...
    'location', 'southwest', 'FontSize', 13);
lgd.Layout.Tile = 'south';

% Experiment 2

spec_gap_list = [10 30 50];
normal_sd_list = [0.1 0.3 0.5];
rho_list = [0.1 0.01];
p_list = [0.1 0.3 0.5 0.7 0.9];

for spec_gap = spec_gap_list
    for normal_sd = normal_sd_list
        for rho = rho_list
            for p = p_list
                simulation2(spec_gap, normal_sd, rho, p, 30);
            end
        end
    end
end

res = readtable('simulation2.csv');
res = res{:,:};

figure;
tiledlayout(1,3);
nexttile
normal_sd = 0.1;
spec_gap = 50; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
spec_gap = 30; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
spec_gap = 10; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('p', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13); 
ylim([0, 1]); xlim([0.1, 0.9]); xticks([0.1 0.3 0.5 0.7 0.9]);
title('Exact recovery rates when \it\sigma_{normal} = 0.1')

nexttile
normal_sd = 0.3;
spec_gap = 50; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
spec_gap = 30; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
spec_gap = 10; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('p', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13); ylim([0, 1]); xlim([0.1, 0.9]);
xticks([0.1 0.3 0.5 0.7 0.9]);
title('Exact recovery rates when \it\sigma_{normal} = 0.3')

nexttile
normal_sd = 0.5;
spec_gap = 50; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-o', 'LineWidth', 1, 'Color', [1 0.3 0])
hold on
spec_gap = 30; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-^', 'LineWidth', 1, 'Color', [0.2 0.7 0.2])
hold on
spec_gap = 10; 
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.1), :);
er_rate1 = X(:, 5);
X = res((res(:,1)==spec_gap)&(res(:,2)==normal_sd)&(res(:,3)==0.01), :);
er_rate2 = X(:, 5);
p = X(:, 4); er_rate = max(er_rate1, er_rate2);
plot(p, er_rate, '-x', 'LineWidth', 1, 'Color', [0 0.3 1])
hold off
xlabel('p', 'FontSize', 13); ylabel('Exact recovery rate', 'FontSize', 13); ylim([0, 1]); xlim([0.1, 0.9]);
xticks([0.1 0.3 0.5 0.7 0.9]);
title('Exact recovery rates when \it\sigma_{normal} = 0.5')

lgd = legend({'\lambda_1 - \lambda_2 =50', '\lambda_1 - \lambda_2 =30', '\lambda_1 - \lambda_2 =10'},...
    'location', 'southwest', 'FontSize', 13);
lgd.Layout.Tile = 'south';