% July 2018, Etemadi and Lu
%The observed RSEs of $\widehat{\Delta}_{PES}$ support our estimated RSEs based on Eq. 9
clear all;clf('reset');
data_name={'Ego-facebook';'Wiki-Vote';'AstroPh';'Enron-email';'Brightkite';'Gowalla';'Libimseti';'Digg';...
    'Amazon';'Citeseer';'Web-Google';'Youtube';'Dblp';'As-skitter'; 'FLicker';'Orkut';'Wiki-German';'USpatent';'LiveJournal';... 
    'DBpedia';'Web-Arabic';'Gsh-2015';'MicrosoftAG';'Friendster'};
size_subplot=size(data_name,1);

for i=1: size_subplot
    str1=data_name(i);
    path=char(strcat('our/',str1,'-raw-data-plot-tri.txt'));
    properties=dlmread(path,',',[5 0 5 9]);
    C=properties(1); N=properties(2); M=properties(3); T=properties(4); K=properties(5); CW=properties(6); W=properties(7);KW=properties(9);KWT=properties(10);
    iteration=dlmread(path,',',[6 0 6 0]);
    %iteration=700;
    data=dlmread(path,',',9);
    point_number=1.0*size(data,1)/(iteration+1);
    n=zeros(point_number,1);
    ObRseT=0;
    EstRseT=0;
    AproxEstRseT=0;
    EstRseT1=0;
    p=0;
    for h=1:point_number
        start=(h-1)*(iteration+1)+1;
        p(h)=data(start,1);
       
        start=start+1;        
        EstT=data(start:start+iteration-1,4);
         
        
        ObMeanT(h)=mean(EstT);    %calculate average of closed wedges
        ObVarT(h)=var(EstT); %calculate variance of closed wedges
        ObRseT(h)=sqrt(ObVarT(h))/T;
        AreT=(ObMeanT(h)-T)/T;
        
        
       
       n= floor(p(h)*M);
       wg=floor(p(h)*W);
       q=n/wg;
       qq2=(n*(n-1))/(wg*(wg-1));
       
       a=T*(1-q*p(h))/(q*p(h));       
       b=2.0/5.0*K*(qq2-p(h)*q*q)/(p(h)*q*q);
       c=(T*T-2*K-T)*(qq2-q*q)/(q*q);
       EstRseT2(h)=sqrt(a+b+c)/T;
       
       
      
       t1=p(h)*M/W;
       t2=(p(h)*p(h)*M*M/(W*W));
       a=T*(t1-t2)/t2;
       t3=(p(h)*M*(p(h)*M-1)/(W*(p(h)*W-1)));
       b=(2.0/5.0)*K*(t3-t2)/t2;
       c=(T*T-2.0/5.0*K-T)*(p(h)*t3-t2)/t2;
       %d=(32/120)*K*(1-q)/q;
       EstRseT(h)=sqrt(a+b+c)/T;
       AproxEstRseT(h)=1/sqrt(T*t1);
       
       a=W/M;
       b=W*(M-1/p(h))/(M*(W-1/p(h)));
       b1=(M-W)/(p(h)*M*W-M);
       t1=(a-p(h))/(p(h)*T);
       t2=2*K*(b-p(h))/(5*p(h)*T*T);
       t3=(b1)*(T*T-2*K-T)/(T*T);
       EstRseT1(h)=sqrt(t1+t2+t3);
       
    end
    
    h1=subplot(4,6,i);
    s=1;e=5;
    line1=plot(p,ObRseT,'r-o','MarkerSize',5,'LineWidth',1.2);
    hold on
    %line2=plot(p,EstRseT,'k','LineWidth',0.6);
    %hold on
    line3=plot(p,AproxEstRseT,'b','LineWidth',1.2);
    %grid on
    set(h1,'Color',[0.96 0.98 0.99]);
%     if i==16
       ylim([0.07 0.48]);
%     else
%         ylim([0.08 0.48]);
%     end
    xlim([p(1)-p(1)/2 p(floor(point_number))+p(floor(point_number))/40]);
    box off;
    p = get(h1, 'position');
    if mod (i,6)==1
        p(1)=p(1)-0.05;
     elseif mod (i,6)==2
        p(1)=p(1)-0.08;
     elseif mod (i,6)==3
        p(1)=p(1)-0.11;
     elseif mod (i,6)==4
        p(1)=p(1)-0.14;
     elseif mod (i,6)==5
        p(1)=p(1)-0.17;
     elseif mod (i,6)==0
        p(1)=p(1)-0.2; 
     end 
     if i<=6
      p(2)=p(2)+0.042;
      
    elseif i<=12
      p(2)=p(2)+0.047;
      
    elseif i<=18
      p(2)=p(2)+0.058;
    else
         p(2)=p(2)+0.067;
    end
     p(3)=p(3)-0.01;  %+0.036;
     p(4)=p(4)+0.01;
    set(h1, 'position', p);

    %h=refline([0 0]);
    j = findobj(gca,'Type','text');
    set(j,'Position', [0 0 500 500]);
    %name=strcat('(',num2str(i),')',{' '},str1); 
    t=title(str1,'FontSize',12,'FontWeight','bold','Interpreter','Latex');
    p1=get(t,'Position');
    set(t,'Position',[p1(1) p1(2)*0.85  p1(3)]);
    
    if i==1||i==7 || i==13 || i==19 ||i==33||i==41||i==49||i==57
        ylabel('RSE','FontSize',9,'Interpreter','Latex');
    else
      set(gca,'YTickLabel',[]);  
    end
    if i>=19
        xlabel('$p$','Interpreter','Latex','FontSize',12);
    end
    
    
    set(gca,'FontSize',3,'XTickMode','auto');
    set(gca,'XTickMode','manual','YTickMode','manual');
    set(gca,'FontSize',8);
    NumTicks = 4;
    L = get(gca,'XLim');
    set(gca,'XTick',linspace(L(1),L(2),NumTicks));
    n=get(gca,'Xtick');
    b=size(n,2);
    z=n(1,1:b-1);
    format shortG;
    set(gca,'Xtick',[z'],'Xticklabel',num2str(get(gca,'XTick')','%1.0e'));
end
hL = legend([line1 line3],{'Obs. ','Est. '},'FontSize',11, 'LineWidth', 4,'interpreter','latex','Orientation','horizontal');
newPosition = [0.42 0 0.005 0.22];
newUnits = 'normalized';
set(hL,'Position', newPosition,'Units', newUnits);
legend('boxoff');
set(gcf,'color','w');
export_fig Fig5 -painters -eps

