beta = 1; p = 10; n = 100; 

seed = 1; randn('state',seed); 

sigma = 1; 

iter = 50000; 

alpha_arr = [0.25:0.25:10] / 100; 
n_alpha = length(alpha_arr); 
s_alpha_arr = zeros(n_alpha,1); 
s_tag_alpha = zeros(n_alpha,1); 


%% find s_alpha by inverting the TW distribution
if beta==1  load TW_beta1.mat
else load TW_beta2.mat
end

[mu_np sigma_np] = KN_mu_sigma(n,p,beta);   mu_np = mu_np / n; sigma_np = sigma_np / n; 

for j=1:n_alpha
    alpha = alpha_arr(j); 
    idx = min ( find(TW_s > 1-alpha) ); 
    s_alpha_arr(j) = x(idx);
end

t_alpha = mu_np + s_alpha_arr * sigma_np;            


for j=1:n_alpha
    alpha = alpha_arr(j); 
    t_tilde_alpha(j) = TW_trace_ratio_threshold(p,n,beta,alpha); 
end


P_L1 = zeros(n_alpha,1); 
P_U   = zeros(n_alpha,1); 
P_U_tilde=zeros(n_alpha,1); 

for counter = 1:iter
    
          %% generate data
            X = zeros(n,p); 
            if beta==1     % add noise
                X = X + sigma * randn(n,p); % real noise
            else
                X = X + 1/sqrt(2) * sigma * complex(randn(n,p),randn(n,p)); 
            end
            
                S = 1/n * X' * X;   
                S = (S + S') / 2 ; 
                lambda = eig(S); %eigs(S,'lm',options);            
                lambda = sort(lambda,'descend'); 

                detect_results = lambda(1) > t_alpha ; 
    
                P_L1 = P_L1 + detect_results;
                
                detect_results = lambda(1)/(1/p * sum(lambda)) > t_alpha; 
                
                P_U = P_U + detect_results; 
                
                detect_results = lambda(1)/(1/p * sum(lambda)) > t_tilde_alpha; 
                
                P_U_tilde = P_U_tilde + detect_results; 

    if mod(counter,500)==0
%%        
        fsize = 20; 
        figure(5); clf; set(gca,'fontsize',fsize); 

        plot(alpha_arr*100,alpha_arr*100,'r-'); hold on; 
        plot(alpha_arr*100,P_L1*100/counter,'b.-');
        hold on ; 
        plot(alpha_arr*100,P_U*100/counter,'mo-'); 
        plot(alpha_arr*100,P_U_tilde*100/counter,'gs-'); 
        title(['p = ' num2str(p) ', n = ' num2str(n) ', \beta = ' num2str(beta)]); 
        grid on; 
        axis([0 max(alpha_arr*100) 0 max(alpha_arr*100)]); 
        h = legend('\mbox{$\alpha$}','\mbox{$\Pr[\ell_1 > t(\alpha)$]}','\mbox{$\Pr[U>t(\alpha)]$}', ...
            '\mbox{$\Pr[U>\tilde{t}(\alpha)]$}',2); 
        set(h,'interpreter','latex'); 
        set(h,'fontsize',fsize+2); 
        
        xlabel('\alpha (in percents)'); 
        
        drawnow; 
        
    end
end

