%function []=sizeEst()
clear all;
close all;
% July 2013, Jianguo
% Size est on 18 datasets. RE sampling
% t1 is the dominant term. t4 is negative, reducing the variance.
%  overestimate, can be corrected by C+1
%   RW/RN sample size ratio is the same as gamma, better fit than RE/RN.  
% M2^2=M3

% t1 is the dominant term, t1> t2
% corr(t1,t4)<0, or corr(,x,y)<0, except flickr  and  gnutella. in the sample, if gamma happens to be high by chance, collissions will be higher, therefore n^2/C lower. 

datapath = 'c:\size\networks\';
plotpath = 'c:\size\plots\';
samplepath='c:\size\samples\');

filename = {'WikiTalk','NotreDame',    'flickrEdges', ... 
        'BerkStan','Gowalla_edges','IMDB',...
        'EmailEu','soc-Epinions1','dblp',...
        'Stanford','Google','Amazon0505'...
        'skitter','Slashdot0902','Gnutella',...
        'youtube','indexedFacebook','cit-Patents'};

nFor20Percent=[ 624,  983,        1043, ... %re
        673        1056 2164, ... 
        314         709        3067, ...
        489, 1937,   3278, ...
        1176,         740,        1348, ...
        1187,   3907,    9737];
nFor10Percent=[1839,  2336, 2246, ...        %re
        1696, 2873, 4249, ...
        669 , 1579, 5485,...        
        1378, 3854, 6486,...        
        2652, 1610,  2397, ...
        2774,  7585,  17904];

%nFor10PercentRW=[     1929        2446        2326        1786        3713
%5     3354        8255       18394] %rw restart=50, nForRSER step=100

population = [2388953,325729,105720,654782,196591,374511,224832,75877,511163,255265,855802,410236,1694616,82168,62561,1134890,2937612,3764117];
MU3=[8138403, 27942.4280, 124.5373,...  
        431295,26608, 71.9599, ...  
        211990, 1132, 38.2117,...  
        168160, 4575.1941, 108.5653, ...
        152325, 680.4380, 8.8862, ... % is it normalized?
        208094, 310.0517, 16.86519]; % third moment about the origin, sum di^3 /N/<d>^3. 


%GAMMA = [26.32,6.398,2.6441,... % gamma from wang hao. estimated gamma? 
%        14.5087,5.5417,2.0535,...
%        13.6648,4.0241,1.6097,...
%        11.5096,3.9983,1.2652,...
%        10.4581,3.3464,1.2061,...
%        9.6399,3.1401,1.1966];%,14.69,13.93,11.79,10.46,9.65,6.40,5.54,4.02,4.02,3.35,3.14,2.65,1.56,1.27,1.20,0.35];


%gammas=[9.851841416235846, 6.2053643297295595, 2.236119508793825, ...
%        20.605701982536015,6.679767495619156, 2.065359253363805, ...
%        3.7064345196346515, 2.828263336430937, 1.1633611015443646, ...
%        15.12981392615206, 4.915477877718687, 1.3193968313413063, ...
%        13.822156748337083,2.534635781518198,0.2833490280697166, ...
%        9.28640643250971,1.4671857363179228,0.681864913941838
%];  % true gamma before removing direction

gammas=[25.34305273108365,5.475659747628559, 1.8268283230744222, ...
        13.543103290003089, 4.6312075505642545, 1.2840709100658603,...
        12.70128304653521,3.146509852069128, 0.8950611452121939,...
        10.552924266749955, 3.1214303946823243, 0.6126547328656944,...
        9.505783490195052, 2.4925780322517976, 0.5667438212579112,...
        8.691596177502182, 2.295477277409341, 0.5594489542338521]; % true gamma after removing direction

filesize = size(filename,2);

for p = 1:filesize   
    gamma=gammas(p);
    Gamma =gamma^2+1;
    N=population(p);
    
    %    degrees = textread(char(strcat(path,filename(p),'_rwdegree_0_10000000.txt')));
    %   ids = textread(char(strcat(path,filename(p),'_rwnode_0_10000000.txt')));
    %   degrees=degrees(1:(10000000-100),1);
    %  ids=ids(1:(10000000-100),1);
    
    % len=size(degrees,1);
    % thinning=100;
    % height=floor(len/thinning);
    % degrees=reshape(reshape(degrees,thinning,height)',len,1);     
    % ids=reshape(reshape(ids,thinning,height)',len,1);     
    degrees = textread(char(strcat(samplepath,filename(p),'_redegree_10000000.txt')));
    ids = textread(char(strcat(samplepath,filename(p),'_renode_10000000.txt')));
    repeats=500;
    rse=0.10;
    %    [required_n,est_Ns]=nForRSE(degrees, ids, N, Gamma, rse, repeats);
    %    n(p)=required_n;
    %    hatN(p,:)=est_Ns;
    
    n=nFor10Percent(p);
    ns(p)=n;
    n0(p)=(2*N/rse^2)^(0.5);
    %%%% analyse the variance
    for i = 1:repeats
        sample_ids=ids((i-1)*n+1:i*n);
        sample_degrees=degrees((i-1)*n+1:i*n);
        C(i)=getCollisions(sample_ids);
        
        if C==0
            C=0.1
        end
        mud(i)=n/sum(1./sample_degrees);
        mude(i)=sum(sample_degrees)/n; 
    end
    Cs(p,:)=C;
    muds(p,:)=mud;
    mudes(p,:)=mude;
    y=mude./mud;
    x=n^2./(2.*(C));
    hatNs(p,:)=x.*y;
    
    mu_x=mean(x);
    mu_y=mean(y);
    
    v_x=var(x);
    v_y=var(y);
    v_c=var(C);
    delta_x=x-mu_x;
    delta_y=y-mu_y;
    
    E11=mean(delta_x.*delta_y);
    E12=mean(delta_x.*delta_y.^2);
    E21=mean(delta_x.^2.*delta_y);
    E22=mean(delta_x.^2.*delta_y.^2);
    
    t2(p)=mu_x^2*v_y;
    t1(p)=mu_y^2*v_x;
    t4(p)=mean(delta_x.^2.*delta_y.^2);
    t3(p)= 2*mu_x*mu_y*E11;
    t5(p)= 2*mu_x*E12;
    t6(p)= 2*mu_y*E21;
    t7(p)= E11^2;
    
    RSE_C(p)=std(C)./mean(C);
    V_C(p)=v_c;
    E_C(p)=Gamma*n^2/(2*N);
    RSE_N(p)=std(hatNs(p,:))/N;
    RB_N(p)=mean(hatNs(p,:))/N-1;
    %for i=1:10
    %   C=100*i;
    %  n=floor((2*N*C)^(0.5));
    % ns(i, p)=n;
    %[rb, rse, rrmse,hatN]=rse_by_n(degrees, ids, N, Gamma,n, repeats); 
    %RSE(i,p)=rse;
    %RB(i,p)=rb;
    %RRMSE(i,p)=rrmse;
    %hatNs(i,:)=hatN;   
    % end;
    
end


%%%% fig 1: t1 is the dominant term, t3 is negative as expected. 
vars=figure;
pp=1:18;
reshaped=reshape(reshape(pp,6,3)',18,1);
for i=1:filesize
 subplot(3,6,reshaped(i));
 bar([t1(i),t2(i),t3(i),t4(i),t5(i),t6(i),t7(i)]);
 title(filename(i));
 %plot(RSE(:,i),'o');
end;
saveas(vars,'c:\size\plots\vars10percent.eps', 'psc2');

%%%% fig 2: n0/n grows linearly with 
rnre=figure;    
x=(gammas.^2+1).^(0.5); y=n0./ns;
loglog(x, y, 'o', x,x,'-r','LineWidth',2);
xlabel('\gamma', 'FontSize', 14);
ylabel('RN/RE of sample size','FontSize', 14);
legend('RN/RE', 'When x=y');
legend('boxoff');
ylim([1,20]);
%pf = polyfit(x,y,1);   % p returns 2 coefficients fitting r = a_1 * x + a_2
%r = pf(1) .* x; % compute a new vector r that has matching datapoints in x
%hold on;
%loglog(sortrows(x), sortrows(x), '-r');
%hold off;
corrcoef(x,y)
saveas(rnre,'c:\size\plots\rern10percent.eps', 'psc2');

%%%% fig3:  V(C) approx E(C)for small gamma
vc=figure;    
for i=1:filesize
 subplot(3,6,reshaped(i));
 bar([V_C(i),E_C(i)]);
 title(filename(i));
 ylabel('V_C vs. E_C');
end;
saveas(vc,'c:\wanghao\size\vc.eps', 'psc2');

%%%% fig4: t1 vs var(N)
t1std=figure;    
for i=1:filesize
 subplot(3,6,reshaped(i));
 bar([t1(i)^(0.5), std(hatNs(i,:))]);
 title(filename(i));
end;
saveas(t1std,'c:\wanghao\size\t1std.eps', 'psc2');
title('t1');

%%%% fig5: relative bias of N
rb=figure;    
bar(RB_N);
saveas(rb,'c:\wanghao\size\rb.eps', 'psc2');
title('RB');

%%%% fig6: t1+t2+t3 vs var(N)
t1t2t3=figure;    
for i=1:filesize
 subplot(3,6,reshaped(i));
 bar([(t1(i)+t2(i)+t3(i))^(0.5), std(hatNs(i,:))]);
 title(filename(i));
end;
saveas(t1t2t3,'c:\wanghao\size\t123std.eps', 'psc2');


%fig 7: est V(C) = E(C)+n^3 sum p_i^3.
% when MU3 is small, hatVC overestimate
% when MU3 is large, hat
vc_est=figure;

ratio=(E_C+ns.^3./population.^2.*MU3)./V_C;
plot(ratio(reshaped), 'o');
ylabel('observed/estimated var');
xlabel('Graphs');

RSE_RW=[    0.0319    0.1335    0.3103    0.2050    0.0737    0.0561
    0.2730    0.0526    0.0669    0.1915    0.0871    0.0831
    0.0461    0.0562    0.0648    0.0427    0.0417    0.0594
]; % rw without restart, repeat=300, n^2=2NC, C=100. RSE=0.1 for RN.   
%conclustion: RW better than RN in most cases. 

RSE_RE=[
    0.0300    0.0517    0.0526    0.0368    0.0698    0.0521
    0.0361    0.0576    0.0554    0.0450    0.0453    0.0740
    0.0296    0.0561    0.0713    0.0419    0.0366    0.0615];

%when n^2=2NC,C=100
RRMSE_RE=[0.0301    0.0516    0.0526    0.0367    0.0704    0.0521
    0.0363    0.0578    0.0554    0.0453    0.0455    0.0748
    0.0296    0.0560    0.0713    0.0420    0.0366    0.0616];


nForRSE10PercentRW=[1915         362       58951           
    330000        3180    3729    
    0        1394        1883     
    33318     997        3980        
    8081        1927        2769
    3580        8220       15486];

%rw repeats=100; thinning=100. 
nFor10PercentRW=[
    2367       12893       18992       69062        3171
    3465       67708        1236        4829       11138
    3833        5572        2634        1738        2414
    2904        8044       17395
];