% Master file to produce accuracy plots
clear all
rng(1337)
curdict=cd;
addpath([curdict '/glmnet_matlab'])
addpath([curdict '/Graphical Lasso'])
addpath([curdict '/Functions'])
if isempty(gcp('nocreate')) == 1
    parpool([4 28],'IdleTimeout',240);
end
dbstop if error
pctRunOnAll warning('off','all')

%%
%Model Specifications
FinalTableReg='';
Type = 'Network';%

target = 'FEVD';

Rule = 0.2;

Regularization.A = 'elastic-net';
Regularization.Sigma = 'Sample-Cov';

Regularization2.A = 'Areg';
Regularization2.Sigma = 'GLASSO';

Regularization3.A = 'Areg';
Regularization3.Sigma = 'Threshold';

deltas = [logspace(-4,0,5) 1.25:0.25:4];  %logspace(-1,0,5)
F = 100;                 %Parameters for elastic Net
Na = [100];              %Dimension. Put array in here if output table is wished
T = [400:-10:0]+Na;     %Timepoints. Put array in here if output table is wished
R = 500;
H = 10;
%% Main simulation
clear relOracleReg OracleNormReg norms D
t = 0;
for Ta = T
    clear D Areg
    t = t+1;
    %create Data
    dimension = [Na Ta+1 R];
    density = 1/sqrt(Na);
    
    [Y, real] = DGP(dimension,Type,density);
    
    %% estimate VAR parameters
    tic
    [regCurve(t,:), olsCurve(t,:), tpfpReg(:,:,t), tpfpOLS(:,:,t),Areg] = estimationBias(Y,F,Regularization,deltas,real,target,H);
    disp([Regularization.A '_'  Regularization.Sigma ' with T=' num2str(Ta) ' took ' num2str(toc/60) ' minutes'])
    
    %% Sigma
    tic
    if Ta <= Na*1.1
        selection = randi(R,min([50 R]),1);
        Ya = Y(:,:,selection);
    else
        selection = 1:R;
        Ya = Y;
    end
    [regCurve2(t,:), ~,tpfpReg2(:,:,t), ~, ~] = estimationBias(Ya,Areg(:,:,selection),Regularization2,deltas,real,target,H);
    disp([Regularization2.A '_' Regularization2.Sigma ' with T=' num2str(Ta) ' took ' num2str(toc/60) ' minutes'])
    
    %% Both
    tic
    
    [regCurve3(t,:), ~,tpfpReg3(:,:,t), ~, ~] = estimationBias(Y,Areg,Regularization3,deltas,real,target,H);
    disp([Regularization3.A '_' Regularization3.Sigma ' with T=' num2str(Ta) ' took ' num2str(toc/60) ' minutes'])
    
end
%% Plotting section
% Rename Legend
if strcmp(Regularization2.A,'Areg')
    Regularization2.A = 'elastic-net';
end
if strcmp(Regularization3.A,'Areg')
    Regularization3.A = 'elastic-net';
end
% smooth curves if nessecary
smoothing=0;
if smoothing==1
    regTemp=regCurve2;
        timePoints=1:size(regCurve2,1);
    regTemp(timePoints(end)+1,:)=regTemp(timePoints(end),:);
    for t=timePoints(T<200)
    regCurve2(t,1)=0.5*regTemp(t,1)+0.5*(regTemp(t-1,1)+regTemp(t+1,1));
    end
end
%% Bias Variance Plot
close all

fig1=figure;
fig1.PaperPositionMode = 'manual';
fig1.PaperOrientation = 'landscape';
fig1.PaperSize = [ 10  5 ];
fig1.PaperPosition = [ 0 0 1 1];

breakplot = max(abs(regCurve3(:,1))./abs(regCurve2(:,1)));

Tile = tiledlayout(1,2);
Tile.TileSpacing = 'none';

Titles={'Bias','Frob-Norm','Frobenius Norm to $1_N$','','Variance'};
yTitles={'Average difference to true FEVD (log-scale)','','','','Average Variance of an entry (log-scale)'};

for k = [1 5]
    if breakplot <= 0.01 && k == 1
        ax = nexttile(Tile);
        semilogy(T,olsCurve(1:end,k),T,abs(regCurve(:,k)),':',T,abs(regCurve2(:,k)),'-.','LineWidth',1.5)
        ax.FontSize = 10;
        xticks([100:100:max(T)]);
        title(Titles{k})
        ylabel(yTitles{k})
        ax.YGrid='on';
        
        legend('OLS + Sample Cov',[Regularization.A ' + ' Regularization.Sigma], [Regularization2.A ' + ' Regularization2.Sigma],...
            [Regularization3.A ' + ' Regularization3.Sigma] , ['indicates that the' newline 'true line is negative' ],'FontSize',8,'location','northeast')
        legend boxoff
        hold on
        shifter = 10^(-floor(log10(breakplot))-1);
        semilogy(T, abs(regCurve3(:,k))*shifter, '--','LineWidth',1.5)
%                 semilogy(T, abs(regCurve3(:,k)), '--','LineWidth',1.5)

        py1 = [min(olsCurve(:,k))]*ones(1,2);
        py2 = 1.5*[max(regCurve3(:,k))*10^(-floor(log10(breakplot))-1)]*ones(1,2);
        
        Ylabels = ax.YTickLabel;
        kk = 0; multLines = 4;
        ax.Clipping = 'off'; 
        for ii = linspace(py2(1),py1(1),multLines*log10(shifter))
            kk = kk + 1;
            if rem(kk,multLines) == 1
                yline(ii,'-','LineWidth',0.5,'Color',ones(1,3)*0.15,'Alpha',0.15)
                ax.YTick = sort([ax.YTick ii]);
            else
                yline(ii,':','LineWidth',0.5,'Color',ones(1,3)*0.15,'Alpha',0.15)
            end
        end
        
        px = [min(T)-10 max(T)];
        fill([px flip(px)],[py1 flip(py2)],'w','EdgeColor','none');
        px = [100-5 100+10];
        
        plot(px,py1,'k','LineWidth',0.75);
        plot(px,py2,'k','LineWidth',0.75);

        for ii = 1:length(ax.YTick)
            if ax.YTick(ii) < py2(1)
                newTick = -log10(shifter) + log10(ax.YTick(ii));
                ax.YTickLabel{ii} = ['10^{' num2str(newTick) '}'];
                counter = 0;
            elseif ax.YTick(ii) > py1(1)
                distance2last = length(ax.YTick) - ii;
                ax.YTickLabel{ii} = Ylabels{length(Ylabels)-distance2last};
            else
                counter = counter +1;
                ax.YTickLabel{ii} = ['10^{' num2str(newTick + counter) '}'];
            end
        end
        hold off
    else
        ax = nexttile(Tile);
        
        semilogy(T,olsCurve(1:end,k),T,abs(regCurve(:,k)),':',T,abs(regCurve2(:,k)),'-.',T,abs(regCurve3(:,k)),'--','LineWidth',1.5)
        ax.FontSize = 10;
        xticks([100:100:max(T)]);
        title(Titles{k})
        ylabel(yTitles{k})
        ax.YGrid='on';
        
    end
    if k==1
        legend('OLS + Sample Cov',[Regularization.A ' + ' Regularization.Sigma], [Regularization2.A ' + ' Regularization2.Sigma],...
            [Regularization3.A ' + ' Regularization3.Sigma] ,'FontSize',8,'location','northeast')
        legend boxoff
    end
    xlabel('T')
    box off
    xlim([min(T) max(T)])
    
end

print('-dpdf', './Graphs/BiasOLS.pdf','-fillpage')
% print('./Graphs/BiasOLS.png','-dpng')

%% TP FN ROC Plot
fig2=figure;
fig2.PaperPositionMode = 'manual';
fig2.PaperOrientation = 'landscape';
fig2.PaperSize = [ 45  15 ];
fig2.PaperPosition = [ 0 0 1 1];
Ylims=[ 0 1; 0 1]; YTicks=[0 50 100; 0 50 100; 0 0.5 1]; CutValue = round((1-density)*length(tpfpOLS),0);
for k=1:2
    subplot(1,3,k)
    plot(T,squeeze(tpfpOLS(k,CutValue,:)),T,squeeze(tpfpReg(k,CutValue,:)),':',...
        T,squeeze(tpfpReg2(k,CutValue,:)),'k-.',T,squeeze(tpfpReg3(k,CutValue,:)),'--','LineWidth',2) %
    ax = gca;
    ax.FontSize = 18;
    xticks([100:100:max(T)]);
    yticks([0.1:0.1:1])
    axis([min(T)-3 max(T)+3 Ylims(k,:)]);
    box off
end
legend('OLS + Sample Cov',[Regularization.A ' + ' Regularization.Sigma], [Regularization2.A ' + ' Regularization2.Sigma],...
    [Regularization3.A ' + ' Regularization3.Sigma] ,'FontSize',8,'location','southeast')
legend boxoff

Tvalue = length(T);
subplot(1,3,3)
plot(tpfpOLS(2,:,Tvalue),tpfpOLS(1,:,Tvalue),tpfpReg(2,:,Tvalue),tpfpReg(1,:,Tvalue),':',...
    tpfpReg2(2,:,Tvalue),tpfpReg2(1,:,Tvalue),'k-.',tpfpReg3(2,:,Tvalue),tpfpReg3(1,:,Tvalue),'--','LineWidth',2)
hold on
plot([0 1],[0 1],'k')
hold off
ax = gca;
ax.FontSize = 18;
xticks([0 0.5 1]);
yticks([0 0.5 1])
box off
print('-dpdf', './Graphs/tpfpROC.pdf','-fillpage')
%% Acc, ROC, AUROC Plot
fig = figure;
fig.PaperPositionMode = 'manual';
fig.PaperOrientation = 'landscape';
fig.PaperSize = [ 12  4 ];
fig.PaperPosition = [ 0 0 1 1];
subplot(1,3,1)
 CutValue = round(0.9*length(tpfpOLS),0);

 
Accuracy(1,:)=(squeeze(tpfpOLS(1,CutValue,:))*0.1*Na^2+squeeze(tpfpOLS(1,CutValue,:))*0.9*Na^2)/Na^2;
Accuracy(2,:)=(squeeze(tpfpReg(1,CutValue,:))*0.1*Na^2+squeeze(tpfpReg(1,CutValue,:))*0.9*Na^2)/Na^2;
Accuracy(3,:)=(squeeze(tpfpReg2(1,CutValue,:))*0.1*Na^2+squeeze(tpfpReg2(1,CutValue,:))*0.9*Na^2)/Na^2;
Accuracy(4,:)=(squeeze(tpfpReg3(1,CutValue,:))*0.1*Na^2+squeeze(tpfpReg3(1,CutValue,:))*0.9*Na^2)/Na^2;
plot(T,Accuracy(1,:),T,Accuracy(2,:),':',...
    T,Accuracy(3,:),'-.',T,Accuracy(4,:),'--','LineWidth',1.5)
box off
xlabel('T')
ylabel('Accurarcy')
title('Accuracy')

subplot(1,3,2)
[~,Tvalue]=min(abs(150-T));
plot(tpfpOLS(2,:,Tvalue),tpfpOLS(1,:,Tvalue),tpfpReg(2,:,Tvalue),tpfpReg(1,:,Tvalue),':',...
    tpfpReg2(2,:,Tvalue),tpfpReg2(1,:,Tvalue),'-.',tpfpReg3(2,:,Tvalue),tpfpReg3(1,:,Tvalue),'--','LineWidth',1.5)
hold on
plot([0 1],[0 1],'k')
hold off
title(['ROC (T=' num2str(T(Tvalue)) ')'])
xlabel('False positive rate (FPR)')
ylabel('True positive rate (TPR)')

box off

subplot(1,3,3)
for t=1:size(tpfpOLS,3)
    AUROC(1,t) = -trapz(tpfpOLS(2,:,t),fliplr(tpfpOLS(1,:,t)));
    AUROC(2,t) = -trapz(tpfpReg(2,:,t),fliplr(tpfpReg(1,:,t)));
    AUROC(3,t) = -trapz(tpfpReg2(2,:,t),fliplr(tpfpReg2(1,:,t)));
    AUROC(4,t) = -trapz(tpfpReg3(2,:,t),fliplr(tpfpReg3(1,:,t)));
end
plot(T,AUROC(1,:),T,AUROC(2,:),':',...
    T,AUROC(3,:),'-.',T,AUROC(4,:),'--','LineWidth',1.5)
legend('OLS + Sample Cov',[Regularization.A ' + ' Regularization.Sigma], [Regularization2.A ' + ' Regularization2.Sigma],...
    [Regularization3.A ' + ' Regularization3.Sigma] ,'FontSize',8,'location','southeast')
legend boxoff
title('AUROC')
ylabel('Area under ROC')

xlabel('T')
box off

print('-dpdf', './Graphs/accAUROC.pdf','-fillpage')

%% Save results in mat file
% only save if run for 200 reps or more
if R > 199
    clear A Y Ya real fig fig1 poolname Tile ax
    close all
    save(['./Output/NormResults_MCr_' num2str(R) '_' date])
end