function [ b ] = Threshold_normal_mean( Name_Procedure, T, d, Pre_mean, Gamma, w,s )
%THRESHOLD_NORMAL_MEAN Summary of this function goes here
%   This function aims to compute the threshold such that the ARL of the
%   procedure is almost Gamma. Here we use the approximation that the
%   threshold b is almost the (exp(-T/Gamma)) quantile of the maximum of
%   the mc Monte Carlo trials

% Number of Monte Carlo trials of running the procedure on a pre-change
% sequence
mc = 100;


switch Name_Procedure
    % CUSUM
    case 'CUSUM'
        % To construct the CUSUM procedure we need an exactly specified
        % post-change mean vector
        Post_mean_CUSUM = 1 * ones(d,1);
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1); stat0 = 0;
            
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            % initialization: compute stat(1)
            x = Signal(:,1);
            stat(1) = max(stat0 + NLL_normal_mean(x, Pre_mean) - ...
                NLL_normal_mean(x, Post_mean_CUSUM), 0);
            for t = 2 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % compute stat(t) from stat(2)
                stat(t) = max(stat(t-1) + NLL_normal_mean(x, Pre_mean) - ...
                    NLL_normal_mean(x, Post_mean_CUSUM), 0);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    % WLGLR
    case 'WLGLR'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            S = zeros(1,w); % The statistic inside the max (over window-size)
            PartialSum = zeros(d,T); % sum of the recent t observations
            
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            % compute stat(1)
            x = Signal(:,1); PartialSum(:,1) = x;
            S(1,1) = x'*x/2; stat(1) = S(1,1);
            
            for t = 2 : T
                % Acquire a new observation
                x = Signal(:,t);
                PartialSum(:,t) = PartialSum(:,t-1) + x;
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    % Using the idea of circle array
                    index = mod(j-1,w)+1;
                    if j == 1
                        y = PartialSum(:,t); % set PartialSum(:,0)=0;
                    else
                        y = PartialSum(:,t) - PartialSum(:,j-1);
                    end
                    S(1,index) = y'*y / (2*(t-j+1));
                end
                stat(t) = max(S);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
        
    % WLACM
    case 'WLACM'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            S = zeros(1,w); eta = zeros(1,w);
            
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = Pre_mean; S(1,index) = 0; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    S(1,index) = S(1,index) + NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                end
                stat(t) = max(S);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    % WLSRRS
    case 'WLSRRS'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            Lambda = zeros(1,w); 
            eta = zeros(1,w); % step size
            
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = Pre_mean; Lambda(1,index) = 1; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    Lambda(1,index) = Lambda(1,index) * exp(NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                end
                stat(t) = log(sum(Lambda));
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));       
    
    % WLACM with projection onto L1-ball
    case 'WLACM_L1'
        % The radius of L1-ball
        %s = 1000;
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            S = zeros(1,w); eta = zeros(1,w);
            
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = Pre_mean; S(1,index) = 0; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    S(1,index) = S(1,index) + NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                    theta(:,index) = ...
                        perform_l1ball_projection(theta(:,index), s);
                end
                stat(t) = max(S);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));        

    % WLSRRS with projection onto L1-ball
    case 'WLSRRS_L1'
        %s = 1000;
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            Lambda = zeros(1,w); 
            eta = zeros(1,w); % step size
            
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = Pre_mean; Lambda(1,index) = 1; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    Lambda(1,index) = Lambda(1,index) * exp(NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                    theta(:,index) = ...
                        perform_l1ball_projection(theta(:,index), s);                    
                end
                stat(t) = log(sum(Lambda));
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    case 'MEI'   
        stat_max = zeros(mc,1);
        a = 1;
        omega = 0.1;
        for i = 1 : mc
            % reset the statistics
            stat = zeros(T,1);
            S = zeros(T,1); % The statistic inside the max (over window-size)
            % generate pre-change sequence with length T
            Signal = randn(d,T);
            
            for t = 1 : T
                for k = max(t-w,1) : t
                    %MLE = 1/mean(Signal(:,k:t),2);
                    temp = mean(Signal(:,k:t),2);
                    mu_hat = zeros(d,1);
                    for p = 1:d
                        if abs(temp(p)) >= omega
                            mu_hat(p) = a*temp(p);
                        else
                            mu_hat(p) = 0;
                        end
                    end
                    S(k) = 0;
                    for l = k : t
                        x = Signal(:,l);
                        S(k) = S(k) + NLL_normal_mean(x,Pre_mean) - ...
                         NLL_normal_mean(x,mu_hat);
                    end
                end
                stat(t) = max(S(max(t-w,1):t));
                    
            end
            stat_max(i) = max(stat);
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    end
end

