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',   'BerkStan', 'EmailEu','Stanford','skitter', 'youtube', ...
 'NotreDame',  'Gowalla_edges','soc-Epinions1','Google','Slashdot0902','indexedFacebook',...
   'flickrEdges','IMDB', 'dblp',   'Amazon0505',  'Gnutella',  'cit-Patents' };
           
fileprintname = {
    'WikiTalk',   'BerkStan', 'EmailEu','Stanford','Skitter', 'Youtube', ...
 'NotreDame',  'Gowalla','Epinions','Google','Slashdot','Facebook',...
   'Flickr','IMDB', 'DBLP',   'Amazon',  'Gnutella',  'Citation' };
           
%nFor20Percent=[ 624,  983,        1043, ... %re
%        673        1056 2164, ... 
%        314         709        3067, ...
%        489, 1937,   3278, ...
%        1176,         740,        1348, ...
%        1187,   3907,    9737];
nFor20Percent=[ 624         673         314         489        1176        1187 ...
         983        1056         709        1937         740        3907 ...
        1043        2164        3067        3278        1348        9737];
%nFor10Percent=[1839,  2336, 2246, ...        %re
%        1696, 2873, 4249, ...
%       669 , 1579, 5485,...        
%        1378, 3854, 6486,...        
%        2652, 1610,  2397, ...
%        2774,  7585,  17904];
nFor10Percent=[   1839        1696         669        1378        2652        2774 ...
        2336        2873        1579        3854        1610        7585 ...
        2246        4249        5485        6486        2397       17904];

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


metadata=[26.324065552126513	3.8985128631664163	10547.029128660128	4.822099379647787E8	2388953
14.508681997469608	20.104007135199197	85482.93731348755	3.5044769251372123E9	654782
13.664741385085675	3.023813336179903	1716.4550419869058	5861127.760763592	224832
11.509563810729752	15.214980510449925	30897.735678608504	5.922941542963822E8	255265
10.458082364509036	13.093478404547106	18921.996821698838	3.419314958030763E8	1694616
9.639867035793332	5.265045951590022	2603.72403316621	3.037149425649887E7	1134890
6.397981647914972	6.693343239318574	1878.6861593533274	8379023.923488544	325729
5.541705376283719	9.668062118815206	2964.025484381279	2.404552462091347E7	196591
4.024120270730777	10.694650552868458	1966.5210538107726	1384771.4308947376	75877
3.998273189542118	10.028843120254452	1708.4316746163247	4614806.502085763	855802
3.3463564232412297	12.273147697400448	1837.401883945088	1257929.9909210398	82168
3.1400908403932015	14.269994812112696	2211.486402560992	900960.337020682	2937612
2.6440420511284892	43.51835035944003	15133.65107832009	1.0264008488157397E7	105720
2.05352865142152	80.18370087927991	33542.18576223395	3.709788464797029E7	374511
1.6097380979814548	7.320835036964725	192.4720607712217	14993.062972085225	511163
1.2651700626532878	11.892847044140446	367.83619185054454	182617.38326719255	410236
1.2061037274836033	4.727481977589872	54.859992647176355	938.8709259762472	62561
1.1966123185313842	8.773234200743495	187.18100420364192	11388.604800541534	3764117
]; % gamma, average, first moment, second moment, thid moment \sum d_i^3/N
MUD3=metadata(:,4);
MUD =metadata(:,2);
gammas=metadata(:,1);
Ns=metadata(:,5);
filesize = size(filename,2);

for p = 1:filesize   
    gamma=gammas(p);
    Gamma =gamma^2+1;
    N=Ns(p);
    
    degrees = textread(char(strcat(samplepath,filename(p),'_redegree_10000000.txt')));
    ids = textread(char(strcat(samplepath,filename(p),'_renode_10000000.txt')));
    repeats=100;
    n20=nFor20Percent(p);
    for i=1:10
        n=n20*i;
        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;
        hat_var(i,p)=2*N^3/n^2/Gamma*(1+2*n*MUD3(p)/(N*Gamma*MUD(p)^3));
        hat_rse(i,p)=(2*N/Gamma)^(0.5)/n;
        var(i,p)=rse^2*N^2;
    end;    
end

vars=figure;
pp=1:18;
for p=1:filesize
    subplot(3,6,p);
    plot(ns(:,p), RSE(:,p),'-o', ns(:,p), (hat_var(:,p)).^(0.5)./Ns(p),'-r');
    title(fileprintname(p));
    legend('Obs.', 'Est.')
    legend('boxoff');
    xlabel('n'); ylabel('RSE');
    box off;
end;
saveas(vars,'c:\size\plots\var_18_ns.eps', 'psc2');

	