% Updated by 07/01/2022
% Tasks: compare LR and GET
%% calculate LR
clear
S = 10000;             % number of samples
klist = [2 4 8];       % dimension of data
Nlist = [100, 400, 1600];    % number of observations

for nk = 1:size(klist,2)
    k = klist(nk);
    
    for nN = 1:size(Nlist,2)
        
        N = Nlist(nN);
        
        LR = [];
        filename=strcat('LR_N',num2str(N),'_k',num2str(k));
        mu = zeros(1,k); Sigma = eye(k);
        gs = GlobalSearch('Display','off');
        t = 3;      
        parfor s = 1:S  
            
            rng(s,'twister')

            % generate data from H0
            data = [];
            while size(data,1)<N
                X = mvnrnd(mu,Sigma,t*N);
                u = rand(t*N,1);  

                delete = u>1/2;
                X(delete,:) = [];
                data = [data;X];
            end
            data = data(1:N,:);

            muhat = mean(data);
            L00(s) = concentrated_llhSN(data,muhat,zeros(1,k));
            X0 = [zeros(1,k);eye(k)/N^(1/3)];

            nLR = @(a) -concentrated_llhSN(data, a(1:k), a(k+1:2*k));
            nLg = []; para = [];
            for j = 1: size(X0,1)
                problem = createOptimProblem('fmincon','x0',[muhat X0(j,:)],'objective',nLR,...
                    'lb',-10*ones(1,2*k),'ub',10*ones(1,2*k));
                [para(j,:),nLg(j)] = run(gs,problem);
            end

            [nLu(s),j] = min(nLg);
            para_u(s,:) = para(j,:);

        end

        LR = 2*(-nLu'-L00');

        save(filename)

    end
end

%% Figure 1
clear; clc; close all;

temp = load('Results/size_power_N100_k2_m137.mat','GET');
GET_N100k2 = temp.GET(:,1);
temp = load('Results/size_power_N400_k2_m137.mat','GET');
GET_N400k2 = temp.GET(:,1);
temp = load('Results/size_power_N1600_k2_m137.mat','GET');
GET_N1600k2 = temp.GET(:,1);
temp = load('Results/size_power_N100_k4_m137.mat','GET');
GET_N100k4 = temp.GET(:,1);
temp = load('Results/size_power_N400_k4_m137.mat','GET');
GET_N400k4 = temp.GET(:,1);
temp = load('Results/size_power_N1600_k4_m137.mat','GET');
GET_N1600k4 = temp.GET(:,1);

temp = load('Results/LR_N100_k2.mat'); LR_N100k2 = temp.LR;
temp = load('Results/LR_N400_k2.mat'); LR_N400k2 = temp.LR;
temp = load('Results/LR_N1600_k2.mat'); LR_N1600k2 = temp.LR;
temp = load('Results/LR_N100_k4.mat'); LR_N100k4 = temp.LR;
temp = load('Results/LR_N400_k4.mat'); LR_N400k4 = temp.LR;
temp = load('Results/LR_N1600_k4.mat'); LR_N1600k4 = temp.LR;

fuente = 12;
figure(1) %% K =2
scatter(LR_N100k2,GET_N100k2,9,'MarkerEdgeColor',[0.2,0.2,.6],...
              'MarkerFaceColor',[1,1,1],...
              'LineWidth',1);

hold on
scatter(LR_N400k2,GET_N400k2,8,'Marker','p','MarkerEdgeColor',[0.6,0.2,.2],...
              'MarkerFaceColor',[0.6,0.2,.2],...
              'LineWidth',1);
set(gco,'fontname','times','fontsize',fuente);
set(gca,'linewidth',1,'fontname','times','fontsize',fuente);
set(gco,'fontname','times','fontsize',fuente);

hold on
scatter(LR_N1600k2,GET_N1600k2,7,'Marker','+','MarkerEdgeColor',[0.2,0.4,.2],...
              'MarkerFaceColor',[0.2,0.4,.2],...
              'LineWidth',1.5);%{'color'},{[0.2,0.2,.6]})
set(gco,'fontname','times','fontsize',fuente);
set(gca,'linewidth',1,'fontname','times','fontsize',fuente);
set(gco,'fontname','times','fontsize',fuente);

hold on
h2 = line([00,40],[0,40]);
set(h2,'color','k');
plot([0:0.1:40],[0:0.1:40],'k','LineWidth',1);
legend('{\itn} = 100','{\itn} = 400','{\itn} = 1,600','45º','Location','SouthEast')

xlabel('LR')
ylabel('GET')
xlim([0 40])
ylim([0 40])

figure(2) %% K =2
scatter(LR_N100k4,GET_N100k4,9,'MarkerEdgeColor',[0.2,0.2,.6],...
              'MarkerFaceColor',[1,1,1],...
              'LineWidth',1);

hold on
scatter(LR_N400k4,GET_N400k4,8,'Marker','p','MarkerEdgeColor',[0.6,0.2,.2],...
              'MarkerFaceColor',[0.6,0.2,.2],...
              'LineWidth',1);
set(gco,'fontname','times','fontsize',fuente);
set(gca,'linewidth',1,'fontname','times','fontsize',fuente);
set(gco,'fontname','times','fontsize',fuente);

hold on
scatter(LR_N1600k4,GET_N1600k4,7,'Marker','+','MarkerEdgeColor',[0.2,0.4,.2],...
              'MarkerFaceColor',[0.2,0.4,.2],...
              'LineWidth',1.5);%{'color'},{[0.2,0.2,.6]})
set(gco,'fontname','times','fontsize',fuente);
set(gca,'linewidth',1,'fontname','times','fontsize',fuente);
set(gco,'fontname','times','fontsize',fuente);

hold on
h2 = line([00,40],[0,40]);
set(h2,'color','k');
plot([0:0.1:40],[0:0.1:40],'k','LineWidth',1);

xlabel('LR')
ylabel('GET')
xlim([0 40])
ylim([0 40])

S = 10000;
% Gaussian ranks of GET
[~,IN100k2] = sort(GET_N100k2); [~,IN100k2] = sort(IN100k2); GGN100k2 = norminv(IN100k2/(S+1));
[~,IN400k2] = sort(GET_N400k2); [~,IN400k2] = sort(IN400k2); GGN400k2 = norminv(IN400k2/(S+1));
[~,IN1600k2] = sort(GET_N1600k2); [~,IN1600k2] = sort(IN1600k2); GGN1600k2 = norminv(IN1600k2/(S+1));
[~,IN100k4] = sort(GET_N100k4); [~,IN100k4] = sort(IN100k4); GGN100k4 = norminv(IN100k4/(S+1));
[~,IN400k4] = sort(GET_N400k4); [~,IN400k4] = sort(IN400k4); GGN400k4 = norminv(IN400k4/(S+1));
[~,IN1600k4] = sort(GET_N1600k4); [~,IN1600k4] = sort(IN1600k4); GGN1600k4 = norminv(IN1600k4/(S+1));
% Gaussian ranks of LR
[~,IN100k2] = sort(LR_N100k2); [~,IN100k2] = sort(IN100k2); GLN100k2 = norminv(IN100k2/(S+1));
[~,IN400k2] = sort(LR_N400k2); [~,IN400k2] = sort(IN400k2); GLN400k2 = norminv(IN400k2/(S+1));
[~,IN1600k2] = sort(LR_N1600k2); [~,IN1600k2] = sort(IN1600k2); GLN1600k2 = norminv(IN1600k2/(S+1));
[~,IN100k4] = sort(LR_N100k4); [~,IN100k4] = sort(IN100k4); GLN100k4 = norminv(IN100k4/(S+1));
[~,IN400k4] = sort(LR_N400k4); [~,IN400k4] = sort(IN400k4); GLN400k4 = norminv(IN400k4/(S+1));
[~,IN1600k4] = sort(LR_N1600k4); [~,IN1600k4] = sort(IN1600k4); GLN1600k4 = norminv(IN1600k4/(S+1));
% 
[corr(GGN100k2,GLN100k2),corr(GGN400k2,GLN400k2),corr(GGN1600k2,GLN1600k2)]
[corr(GGN100k4,GLN100k4),corr(GGN400k4,GLN400k4),corr(GGN1600k4,GLN1600k4)]