% Updated by 07/01/2022
% Tasks: calculate the size and power 
%% calculate test statistics
clear
S = 10000;             % number of samples
klist = [2 4 8];       % dimension of data
Nlist = [100, 400, 1600];    % number of observations
skew = 0.6;            % skewness coefficient in H_a1 and H_a2

[m, dif] = fmincon(@(m) abs(m_skew(m)-skew),0);

for nk = 1:size(klist,2)
    nk
    k = klist(nk);
    
    for nN = 1:size(Nlist,2)
        nN
        clearvars -except S m k Nlist klist nN nk skew

        N = Nlist(nN); % number of observations

        mu= zeros(1,k); Omega = eye(k); omega = diag(diag(Omega));

        A(1,:) = zeros(1,k);                                            % H0
        A(2,:) = sqrt(m/2/k)*ones(1,k);                                 % H_a1
        A(3,:) = [sqrt(m/5/k)*ones(1,k/2)  2*sqrt(m/5/k)*ones(1,k/2)];  % H_a2

        filename=strcat('Results/size_power_N',num2str(N),'_k',num2str(k),'_m',num2str(round(m*10)));
        gs = GlobalSearch('Display','off','NumTrialPoints',3000);

        X0 = eye(k); % get k initial values for GET

        t=3;
        for na = 1:size(A,1)
            a = A(na,:);    

            parfor s = 1:5
                
                rng(s,'twister')

                % generate data from H0/ H_a1/ H_a2
                data=[];
                while size(data,1)<N
                    X = mvnrnd(mu,Omega,t*N);
                    u = rand(t*N,1);  
                    delete = u>(normcdf((a*omega^-1*(X-ones(t*N,1)*mu)')'));
                    X(delete,:) = [];
                    data = [data;X];
                end
                data = data(1:N,:);    

                % calculate the test statistic of GET
                GETk = [];
                for i = 1:k
                    problem = createOptimProblem('fmincon','x0',X0(:,i),'objective',@(a) -GETstat(data,a),...
                        'nonlcon',@normalization);
                    [~,GETk(i)] = run(gs,problem);
                end
                GET(s,na) = max(-GETk);
                GETall(s,:,na) = -GETk;

                % calculate the test statistic of Margins
                margin(s,na) = Marginstat(data);

                % calculate the test statistic of AA
                lm_AA(s,na) = LMAAstat(data);

                % calculate the test statistic of GMM
                GMM(s,na) = GMMstat(data);

            end           
        end
%         save(filename)
    end
end

% %% Table 1
% clear
% 
% Nlist = [100 400];
% klist = [2 4 8];
% Table_1 = [];
% 
% for nk = 1:size(klist,2)
%     
%     k = klist(nk);
%     sizepower_k = [];
%     for nN = 1:size(Nlist,2)
%         N = Nlist(nN);  
% 
%         filename=strcat('Results/size_power_N',num2str(N),'_k',num2str(k),'_m137');
%         load(filename,'S','GMM','GET','lm_AA','margin')
%         cvGMM = GMM(:,1);
%         cvGET = GET(:,1);
%         for i = 1:3
%             sizepower(1,(i-1)*3+1:i*3) = mean(GET(1:S,i)>quantile(cvGET,[0.99 0.95 0.9]));
%             sizepower(2,(i-1)*3+1:i*3) = mean(lm_AA(1:S,i)>chi2inv([0.99 0.95 0.9],k));
%             sizepower(3,(i-1)*3+1:i*3) = mean(GMM(1:S,i)>quantile(cvGMM,[0.99 0.95 0.9]));
%             sizepower(4,(i-1)*3+1:i*3) = mean(margin(1:S,i)>chi2inv([0.99 0.95 0.9],k));
%         end
%         sizepower_k = [sizepower_k sizepower*100];
%         clearvars -except Nlist klist nk nN k N sizepower_k Table_1
%     end
%     Table_1 = [Table_1; sizepower_k];
% end

