%% clear all the enviroment variables
close all
clear all

%% Parameter setting to generate pre-change and post-change time series

d = 20; % the dimension of the parameter space
% in the normal mean detection case, it is the dimension of mean vector
% d can also be seen as the number of sensors in the context of
% multi-sensor detection
Pre_mean = zeros(d,1); % the mean vector before the change-point
w = 100; % w is the window-size used in the procedures

%% Calculate the threshold b in order to make ARL be Gamma
Gamma = 10000; % Gamma is the allowable lower bound for ARL

T = 500; % the length of pre-change sequence
%%

% b_CUSUM = Threshold_normal_mean('CUSUM', T, d, Pre_mean, Gamma, w);
% % b_GLR = Threshold_normal_mean('GLR', T, d, Pre_mean, Gamma, w);
% % b_SRRS = Threshold_normal_mean('SRRS', T, d, Pre_mean, Gamma, w);
% % b_ACM = Threshold_normal_mean('ACM', T, d, Pre_mean, Gamma, w); 
% % b_WLCUSUM = Threshold_normal_mean('WLCUSUM', T, d, Pre_mean, Gamma, w);
% b_WLGLR = Threshold_normal_mean('WLGLR', T, d, Pre_mean, Gamma, w);
% b_WLSRRS = Threshold_normal_mean('WLSRRS', T, d, Pre_mean, Gamma, w);
% b_WLACM = Threshold_normal_mean('WLACM', T, d, Pre_mean, Gamma, w); 
% b_WLACM_L1 = Threshold_normal_mean('WLACM_L1', T, d, Pre_mean, Gamma, w); 
% b_WLSRRS_L1 = Threshold_normal_mean('WLSRRS_L1', T, d, Pre_mean, Gamma, w); 

%% Compare the procedures through the simulated EDDs
% Note that EDD_... is a result vector for 10000 Monte Carlo trials


% This is a library for the threshold already computed. This is used to
% reduce the time for computing it again.
% The following is for d = 20, w = 100:
b_CUSUM = 6.7639; 
b_WLACM = 6.8665;
b_WLGLR = 29.2049;
b_WLSRRS = 6.9289;
b_WLACM_L1 = 8.00;
b_WLSRRS_L1 = 8.03;



% KL divergence between the pre-change and post-change distributions
% KL = 0.5*(Post_mean - Pre_mean)'*(Post_mean - Pre_mean); 

T = 200; % the length of the post-change series


% EDD_GLR = EDD_normal_mean('GLR', Len, d, Post_mean, b_GLR, w);
% EDD_SRRS = EDD_normal_mean('SRRS', Len, d, Post_mean, b_SRRS, w);
% EDD_ACM = EDD_normal_mean('ACM', Len, d, Post_mean, b_ACM, w);
% EDD_WLCUSUM = EDD_normal_mean('WLCUSUM', Len, d, Post_mean, b_WLCUSUM, w);
EDD_CUSUM = EDD_normal_mean('CUSUM', T, d, Pre_mean, Post_mean, b_CUSUM, w);
disp(mean(EDD_CUSUM));
disp(std(EDD_CUSUM));
EDD_WLGLR = EDD_normal_mean('WLGLR', T, d, Pre_mean, Post_mean, b_WLGLR, w);
disp(mean(EDD_WLGLR));
disp(std(EDD_WLGLR));
EDD_WLGLR_V = EDD_normal_mean('WLGLR_V', T, d, Pre_mean, Post_mean, b_WLGLR, w);
disp(mean(EDD_WLGLR_V));
disp(std(EDD_WLGLR_V));
EDD_WLSRRS = EDD_normal_mean('WLSRRS', T, d, Pre_mean, Post_mean, b_WLSRRS, w);
disp(mean(EDD_WLSRRS));
disp(std(EDD_WLSRRS));
EDD_WLACM = EDD_normal_mean('WLACM', T, d, Pre_mean, Post_mean, b_WLACM, w);
disp(mean(EDD_WLACM));
disp(std(EDD_WLACM));
%% Draw plot 1
% This part draws a plot: EDD versus mu;
T = 200;
mu_mesh = [0.3;0.5;0.6;0.7;1;1.5];
l = length(mu_mesh);
EDD_CUSUM_mesh = zeros(l,1); EDDSTD_CUSUM_mesh = zeros(l,1);
EDD_WLGLR_V_mesh = zeros(l,1); EDDSTD_WLGLR_V_mesh = zeros(l,1);
EDD_WLGLR_mesh = zeros(l,1); EDDSTD_WLGLR_mesh = zeros(l,1);
EDD_WLSRRS_mesh = zeros(l,1); EDDSTD_WLSRRS_mesh = zeros(l,1);
EDD_WLACM_mesh = zeros(l,1); EDDSTD_WLACM_mesh = zeros(l,1);
EDD_WLACM_L1_mesh = zeros(l,1); EDDSTD_WLACM_L1_mesh = zeros(l,1);
EDD_WLSRRS_L1_mesh = zeros(l,1); EDDSTD_WLSRRS_L1_mesh = zeros(l,1);


p = 1; % proportion of affected sensors
d1 = round(d*p); 
affected = randsample(d,d1);
r = zeros(d,1); r(affected) = 1;

for i = 1 : l
    Post_mean = mu_mesh(i) * r.* ones(d,1); % the mean vector after the change-point
    [~,EDD_CUSUM_mesh(i),EDDSTD_CUSUM_mesh(i)] = ...
        EDD_normal_mean('CUSUM', T, d, Pre_mean, Post_mean, b_CUSUM, w);
    [~,EDD_WLGLR_mesh(i),EDDSTD_WLGLR_mesh(i)] = ... 
        EDD_normal_mean('WLGLR', T, d, Pre_mean, Post_mean, b_WLGLR, w);
    [~,EDD_WLGLR_V_mesh(i),EDDSTD_WLGLR_V_mesh(i)] = ... 
        EDD_normal_mean('WLGLR_V', T, d, Pre_mean, Post_mean, b_WLGLR, w);
    [~,EDD_WLSRRS_mesh(i), EDDSTD_WLSRRS_mesh(i)] = ... 
        EDD_normal_mean('WLSRRS', T, d, Pre_mean, Post_mean, b_WLSRRS, w);
    [~,EDD_WLSRRS_L1_mesh(i), EDDSTD_WLSRRS_L1_mesh(i)] = ... 
        EDD_normal_mean('WLSRRS', T, d, Pre_mean, Post_mean, b_WLSRRS, w);
    [~,EDD_WLACM_mesh(i), EDDSTD_WLACM_mesh(i)] = ...
        EDD_normal_mean('WLACM', T, d, Pre_mean, Post_mean, b_WLACM, w);
    [~,EDD_WLACM_L1_mesh(i), EDDSTD_WLACM_L1_mesh(i)] = ...
        EDD_normal_mean('WLACM_L1', T, d, Pre_mean, Post_mean, b_WLACM, w);
end
figure; 
co = 1;
plot(mu_mesh(co:end),EDD_CUSUM_mesh(co:end),'--k*');
hold on; plot(mu_mesh(co:end),EDD_WLGLR_mesh(co:end),'--m+')
hold on; plot(mu_mesh(co:end),EDD_WLSRRS_mesh(co:end),'--go')
hold on; plot(mu_mesh(co:end),EDD_WLACM_mesh(co:end),'--rx')
hold on; plot(mu_mesh(co:end),EDD_WLSRRS_L1_mesh(co:end),'-gd')
hold on; plot(mu_mesh(co:end),EDD_WLACM_L1_mesh(co:end),'-rs')

legend('CUSUM','WLGLR, w=100', 'WLSRRS, w=100', 'WLACM, w=100', 'WLSRRS-L1, w=100', 'WLACM-L1, w=100')
xlabel('\mu')
ylabel('EDD')
xlim([mu_mesh(co)-0.01 mu_mesh(end)])
box on;
myboldify;
print -dpdf EDD_mesh_normal_mean_changeMu.pdf

%% Draw plot 2
% This part draws a plot: EDD versus p;
T = 200;
%p_mesh = [0.35;0.5;0.6;0.7;0.9;1];
p_mesh = [0.1;0.2;0.3;0.4;0.5;0.6];
%p_mesh = [0.9];
l = length(p_mesh);
EDD_CUSUM_mesh = zeros(l,1); EDDSTD_CUSUM_mesh = zeros(l,1);
EDD_WLGLR_mesh = zeros(l,1); EDDSTD_WLGLR_mesh = zeros(l,1);
EDD_WLGLR_V_mesh = zeros(l,1); EDDSTD_WLGLR_V_mesh = zeros(l,1);
EDD_WLSRRS_mesh = zeros(l,1); EDDSTD_WLSRRS_mesh = zeros(l,1);
EDD_WLACM_mesh = zeros(l,1); EDDSTD_WLACM_mesh = zeros(l,1);
EDD_WLSRRS_L1_mesh = zeros(l,1); EDDSTD_WLSRRS_L1_mesh = zeros(l,1);
EDD_WLACM_L1_mesh = zeros(l,1); EDDSTD_WLACM_L1_mesh = zeros(l,1);


Time_CUSUM = 0;
Time_WLGLR = 0;
Time_WLGLR_V = 0;
Time_WLACM = 0;
Time_WLSRRS = 0;
Time_WLACM_L1 = 0;
Time_WLSRRS_L1 = 0;
REPS = 1000;
for i = 1 : l
    p = p_mesh(i);
    d1 = round(d*p); 
    affected = randsample(d,d1);
    r = zeros(d,1); r(affected) = 1;
    
    Post_mean = 1 * r.* ones(d,1); % the mean vector after the change-point
    Post_mean = repmat(Post_mean, 1, T);
    % CUSUM
    tic;
    [~,EDD_CUSUM_mesh(i),EDDSTD_CUSUM_mesh(i)] = ...
        EDD_normal_mean('CUSUM', T, d, Pre_mean, Post_mean, b_CUSUM, w);
    Time_CUSUM = Time_CUSUM + toc;
    % WLGLR
    tic;
    [~,EDD_WLGLR_mesh(i),EDDSTD_WLGLR_mesh(i)] = ... 
        EDD_normal_mean('WLGLR', T, d, Pre_mean, Post_mean, b_WLGLR, w);
    Time_WLGLR = Time_WLGLR + toc;
    %WLGLR_V
    tic;
    [~,EDD_WLGLR_V_mesh(i),EDDSTD_WLGLR_V_mesh(i)] = ... 
        EDD_normal_mean('WLGLR_V', T, d, Pre_mean, Post_mean, b_WLGLR, w);
    Time_WLGLR_V = Time_WLGLR_V + toc;
    %WLSRRS
    tic;
    [~,EDD_WLSRRS_mesh(i), EDDSTD_WLSRRS_mesh(i)] = ... 
        EDD_normal_mean('WLSRRS', T, d, Pre_mean, Post_mean, b_WLSRRS, w);
    Time_WLSRRS = Time_WLSRRS + toc;
    %WLACM
    tic;
    [~,EDD_WLACM_mesh(i), EDDSTD_WLACM_mesh(i)] = ...
        EDD_normal_mean('WLACM', T, d, Pre_mean, Post_mean, b_WLACM, w);
    Time_WLACM = Time_WLACM + toc;
    %WLSRRS_L1
    tic;
    [~,EDD_WLSRRS_L1_mesh(i), EDDSTD_WLSRRS_L1_mesh(i)] = ... 
        EDD_normal_mean('WLSRRS', T, d, Pre_mean, Post_mean, b_WLSRRS, w);
    Time_WLSRRS_L1 = Time_WLSRRS_L1 + toc;
    %WLACM_L1
    tic;
    [~,EDD_WLACM_L1_mesh(i), EDDSTD_WLACM_L1_mesh(i)] = ...
        EDD_normal_mean('WLACM_L1', T, d, Pre_mean, Post_mean, b_WLACM, w);
    Time_WLACM_L1 = Time_WLACM_L1 + toc;
end

figure; 
co = 1;
plot(p_mesh(co:end),EDD_CUSUM_mesh(co:end),'--k*');
hold on; plot(p_mesh(co:end),EDD_WLGLR_V_mesh(co:end),'--m+')
hold on; plot(p_mesh(co:end),EDD_WLSRRS_mesh(co:end),'--go')
hold on; plot(p_mesh(co:end),EDD_WLACM_mesh(co:end),'--rx')
hold on; plot(p_mesh(co:end),EDD_WLSRRS_L1_mesh(co:end),'-gd')
hold on; plot(p_mesh(co:end),EDD_WLACM_L1_mesh(co:end),'-rs')
legend('CUSUM','WLGLR, w=100', 'WLSRRS, w=100', 'WLACM, w=100', 'WLSRRS-L1, w=100', 'WLACM-L1, w=100')
xlabel('p')
ylabel('EDD')
xlim([p_mesh(co)-0.01 p_mesh(end)])
box on;
myboldify;
%print -dpdf EDD_mesh_normal_mean_changeP.pdf








