clc
format compact
disp('Model 2 simulation results analysis ...')


flgReport = 0
K = {'100', '200', '500', '1000', '2000'};

fname_seg1  = 'm2_monte_carlo_simulation_results_n'; %fbaccala2001a_ex3_cdefg_n100_alpha1p_m1000
fname_seg2  = '_alpha1p_m1000';
fname_total = 'm2_monte_carlo_simulation_results_total.mat';

flgFTL = 0;
if exist('flgMatTotalLoaded') == 1,
    flgFTL = strcmp(flgMatTotalLoaded,fname_total);
end;

disp('Be patient, it will take a while ...')
disp('If you want to see further detail, set variable flgReport = 1') 
disp('')

%% ======================
if (~flgFTL) & (exist(fname_total) == 2),
    Ktemp = K;
    disp(['loading ' fname_total ' ...'])
    eval(['load ' fname_total])
    K = Ktemp; clear Ktemp
elseif ~flgFTL,
    
    save tmp.mat K fname_seg1 fname_seg2 fname_total flgFTL
    clear all
    load tmp.mat
    delete tmp.mat
    
    nPlotPoints = 32;
    %=====================
    %True structure;
    trueStruct = ...
        [NaN 0  0  0  0;
        1  NaN 0  0  0;
        1   0 NaN 0  0;
        1   0  0 NaN 1;
        0   0  0  1 NaN];
    
    [mm,mm]=size(trueStruct);
    
    %True Reachability map;
    trueReachab = ...
        [NaN 0  0  0  0;
        1  NaN 0  0  0;
        1   0 NaN 0  0;
        1   0  0 NaN 1;
        1   0  0  1 NaN];
    
    for kFile = 1:length(K),
        skF = int2str(kFile);
        disp(['loading ' fname_seg1 K{kFile} fname_seg2 '.mat ...'])
        
        eval(['load ' fname_seg1 K{kFile} fname_seg2]);
        
        eval(['iPdc' skF '=zeros(1000,mm,mm,32);'])
        eval(['pValuesiPdc' skF '=zeros(1000,mm,mm,32);'])
        eval(['pValuesiPdcT' skF '=zeros(1000,mm,mm,32);'])
        
        eval(['iDtf' skF '=zeros(1000,mm,mm,32);'])
        eval(['SS' skF '=zeros(1000,mm,mm,32);'])
        
        eval(['pValuesGct' skF '=zeros(1000,mm,mm);'])
        eval(['pValuesGctT' skF '=zeros(1000,mm,mm);'])
        eval(['pValuesiGct' skF '=zeros(1000,mm,mm);'])
        eval(['pValmvgc' skF ' = zeros(1000,mm,mm);'])
        eval(['pValmvgcT' skF ' = zeros(1000,mm,mm);'])
        eval(['pValmvgcun' skF ' = zeros(1000,mm,mm);'])
        eval(['sigmvgcT' skF ' = zeros(1000,mm,mm);'])
        eval(['sigmvgcun' skF ' = zeros(1000,mm,mm);'])
        eval(['Fmvgc' skF ' = zeros(1000,mm,mm);'])
        eval(['Fmvgcun' skF ' = zeros(1000,mm,mm);'])
        
        for k = 1:1000,
            c = C{k}; d = D{k}; e = E{k}; f = F{k}; g = G{k};
            for r = 1:mm
                for s = 1:mm
                    eval(['iPdc' skF '(k,r,s,:) = getCij(c.pdc,r,s);'])
                    eval(['pValuesiPdc' skF '(k,r,s,:) = getCij(c.pvalues,r,s);'])
                    eval(['pValuesiPdcT' skF '(k,r,s,:) = getCij(g.pvalues,r,s);'])
                    eval(['iDtf' skF '(k,r,s,:) = getCij(d.dtf,r,s);'])
                    eval(['pValuesiDtf' skF '(k,r,s,:) = getCij(d.pvalues,r,s);'])
                    eval(['SS' skF '(k,r,s,:) = getCij(c.SS,r,s);'])
                end;
            end;
            eval(['pValuesGct' skF '(k,:,:) = c.pValue_gct;'])
            eval(['pValuesGctT' skF '(k,:,:) = g.pValue_gct;'])
            eval(['pValuesiGct' skF '(k,:,:) = c.pValue_igct;'])
            eval(['pValmvgc' skF '(k,:,:) = e.pval;'])
            eval(['pValmvgcT' skF '(k,:,:) = f.pval;'])
            eval(['pValmvgcun' skF '(k,:,:) = e.pvalun;'])
            eval(['sigmvgc' skF '(k,:,:) = e.sig;'])
            eval(['sigmvgcT' skF '(k,:,:) = f.sig;'])
            eval(['sigmvgcun' skF '(k,:,:) = e.sigun;'])
            eval(['Fmvgc' skF '(k,:,:) = e.F;'])
            eval(['Fmvgcun' skF '(k,:,:) = e.Fun;'])
        end;
        eval(['pValuesGct' skF '(pValuesGct' skF '==0)=eps;'])
        eval(['pValmvgc' skF '(pValmvgc' skF ' == 0) = eps;'])
        eval(['pValuesGctT' skF '(pValuesGctT' skF '==0)=eps;'])
        eval(['pValmvgcT' skF '(pValmvgcT' skF ' == 0) = eps;'])
        eval(['pValmvgcun' skF '(pValmvgcun' skF ' == 0) = eps;'])
    end;
    clear C D E F G c d e f g
    flgMatTotalLoaded = fname_total;
    eval(['save ' fname_total])
else
    disp(['Using previously loaded data from ' flgMatTotalLoaded ' ...'])
end;

%% =================
    %True structure;
    trueStruct = ...
        [NaN 0  0  0  0;
        1  NaN 0  0  0;
        1   0 NaN 0  0;
        1   0  0 NaN 1;
        0   0  0  1 NaN];
    
    [mm,mm]=size(trueStruct);
    
    %True Reachability map;
    trueReachab = ...
        [NaN 0  0  0  0;
        1  NaN 0  0  0;
        1   0 NaN 0  0;
        1   0  0 NaN 1;
        1   0  0  1 NaN];

M = length(K);
N = 1000; %Number of Monte Carlo simulations
[mm,mm]=size(trueStruct);
%if ~((exist('trueStruct') == 1) & (exist('FP') == 1)),
nFreq = 32;

alpha = 0.01;
FP = NaN*zeros(mm,mm,M); FPT = FP; FN = NaN*zeros(mm,mm,M); FNT = FP;
FPmv = FP; FNmv = FP; FPmvT = FP; FNmvT = FP; FPunmv = FP; FNunmv = FP;
FNipdc = FP; FPipdc = FP; FNidtf = FP; FPidtf = FP;
FNmvadj = FP; FPmvadj = FP; FNmvTadj = FP; FPmvTadj = FP;
FPipdcf = NaN*zeros(mm,mm,M,nFreq); FNipdcf = FPipdcf;

for kFile = 1:M,
    skF = int2str(kFile);
    for r = 1:mm,
        for s = 1:mm,
            for i=1:N,
                eval(['ptemp(i)= max(-log10(pValuesiPdc' skF ...
                    '(i,r,s,:)));'])
                eval(['ptempDtf(i)= max(-log10(pValuesiDtf' skF ...
                    '(i,r,s,:)));'])
                eval(['ptempT(i)= max(-log10(pValuesiPdcT' skF ...
                    '(i,r,s,:)));'])
                eval(['pValuesiPdcmin' skF '(i,r,s) = min(pValuesiPdc' ...
                                       skF '(i,r,s,:));'])
                eval(['pValuesiPdcTmin' skF '(i,r,s)= min(pValuesiPdcT' ...
                                        skF '(i,r,s,:));'])
                eval(['pValuesiDtfmin' skF '(i,r,s) = min(pValuesiDtf' ...
                                       skF '(i,r,s,:));'])
            end;
            if trueStruct(r,s) == 0,
                eval(['FP(r,s,kFile)=sum(-log10(pValuesGct' skF ...
                                                         '(:,r,s)) > 2);'])
                
                eval(['FPmv(r,s,kFile)=sum(-log10(pValmvgc' skF ...
                                                         '(:,r,s)) > 2);'])
                eval(['FPT(r,s,kFile)=sum(-log10(pValuesGctT' skF ...
                                                         '(:,r,s)) > 2);'])
                eval(['FPmvT(r,s,kFile)=sum(-log10(pValmvgcT' skF ...
                                                         '(:,r,s)) > 2);'])
                
                % ====== FP with p-values adjustment for cMVGC =========
                eval(['FPmvadj(r,s,kFile)=sum(sigmvgc' skF '(:,r,s));'])
                eval(['FPmvTadj(r,s,kFile)=sum(sigmvgcT' skF '(:,r,s));'])
                
                eval(['FPunmv(r,s,kFile)=sum(-log10(pValmvgcun' skF ...
                                                         '(:,r,s)) > 2);'])
                FPipdc(r,s,kFile)=sum(ptemp > 2);
                FPidtf(r,s,kFile)=sum(ptempDtf > 2);
                 for freq=1:nFreq,
                    eval(['FPipdcf(r,s,kFile,freq) = sum(-log10(pValuesiPdc' ...
                                                skF '(:,r,s,freq)) > 2);'])
                end;               
                
            elseif trueStruct(r,s) == 1,
                eval(['FN(r,s,kFile)=sum(-log10(pValuesGct' skF ...
                                                        '(:,r,s)) <= 2);'])
                eval(['FNmv(r,s,kFile)=sum(-log10(pValmvgc' skF ...
                                                        '(:,r,s)) <= 2);'])
                eval(['FNT(r,s,kFile)=sum(-log10(pValuesGctT' skF ...
                                                        '(:,r,s)) <= 2);'])
                eval(['FNmvT(r,s,kFile)=sum(-log10(pValmvgcT' skF ...
                                                        '(:,r,s)) <= 2);'])
                
                % ===== FN with p-values adjustment for cMVGC ==========
                eval(['FNmvadj(r,s,kFile)=sum(sigmvgc' skF ...
                                                         '(:,r,s) == 0);'])
                eval(['FNmvTadj(r,s,kFile)=sum(sigmvgcT' skF ...
                                                         '(:,r,s) == 0);'])
                
                eval(['FNunmv(r,s,kFile)=sum(-log10(pValmvgcun' skF ...
                                                        '(:,r,s)) <= 2);'])
                FNipdc(r,s,kFile)=sum(ptemp <= 2);
                FNipdcT(r,s,kFile)=sum(ptempT <= 2);
                for freq=1:nFreq,
                    eval(['FNipdcf(r,s,kFile,freq) = sum(-log10(pValuesiPdc' ...
                                               skF '(:,r,s,freq)) <= 2);'])
                end;
                
            end
            if trueReachab(r,s) == 0,
                FPidtf(r,s,kFile)=sum(ptempDtf > 2);
            elseif trueReachab(r,s) == 1,
                FNidtf(r,s,kFile)=sum(ptempDtf <= 2);
            end
        end;
    end;
    
end;


trueStruct

trueReachab

%% ======== Boxplots 1 GCT p-values =======================================
str_var = 'GCT';
str_varname = 'pValuesGct';
str_bartitle = 'Model 2';
str_suptitle = 'GCT (Wald)';
xboxplot_var;

%% ======== Boxplots 2 iPDC p-values ======================================
str_var = 'iPDC';
str_varname = 'pValuesiPdc';
str_bartitle = 'Model 2';
str_suptitle = 'iPDC';
xboxplot_var;

%% ======== Boxplots 3 cMVGC p-values =====================================
str_var = 'cMVGC';
str_varname = 'pValmvgc';
str_model = 'Model 2';
str_suptitle = 'cMVGC';
xboxplot_var;

%% ======== Boxplots 4 iDTF p-values ======================================
str_var = 'iDTF';
str_varname = 'pValuesiDtf';
str_bartitle = 'Model 2';
str_suptitle = 'iDTF';
xboxplot_var;

%% ======== Scatterplots GCT x iPDC =======================================
str_varX = 'pValuesGct';
str_FNX = ''; % GCT: ''; iPDC: 'ipdc'; SethGC: 'mv'
str_measureX = 'GCT';
str_varY = 'pValuesiPdc';
str_FNY = 'ipdc'; % see above
str_measureY = 'iPDC';
str_model = 'Model 2';
flgAdj = 0;

str_suptitle = 'GCT x iPDC';

xscatterplot_vars;

%% ======== Scatterplots GCT x iDTF =======================================
str_varX = 'pValuesGct';
str_FNX = ''; % GCT: ''; iPDC: 'ipdc'; SethGC: 'mv'
str_measureX = 'GCT';
str_varY = 'pValuesiDtf';
str_FNY = 'idtf'; % see above
str_measureY = 'iDTF';
str_model = 'Model 2';
flgAdj = 0;

str_suptitle = 'GCT x iDTF';

xscatterplot_vars;

%% ======== Scatterplots iPDC x iDTF ======================================
str_varX = 'pValuesiPdc';
str_FNX = 'ipdc'; % GCT: ''; iPDC: 'ipdc'; SethGC: 'mv'
str_measureX = 'iPDC';
str_varY = 'pValuesiDtf';
str_FNY = 'idtf'; % see above
str_measureY = 'iDTF';
str_model = 'Model 2';
flgAdj = 0;

str_suptitle = 'iPDC x iDTF';

xscatterplot_vars;

%% ======== Scatterplots GCT x cMVGT ======================================
str_varX = 'pValuesGct';
str_FNX = '';
str_measureX = 'GCT';
str_varY = 'pValmvgc';
str_FNY = 'mv';
str_measureY = 'MVGC';
str_model = 'Model 2';
flgAdj = 1;

xscatterplot_vars

%% ======== Scatterplots GCT x cMVGT True Model Order =====================
str_varX = 'pValuesGctT';
str_FNX = 'T';
str_measureX = 'GCT';
str_varY = 'pValmvgcT';
str_FNY = 'mvT';
str_measureY = 'MVGC';
str_model = 'Model 2 with True order';
flgAdj = 1;

xscatterplot_vars

if ~exist('m2_analysis_results.mat'),
    save m2_analysis_results
end;
