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

population = [2388953,325729,105720,654782,196591,374511,224832,75877,511163,255265,855802,410236,1694616,82168,62561,1134890,2937612,3764117];


%data from http://cs.uwindsor.ca/~wang115o/graph/Size/Size/sizeSimpleRWT40/
%A=textread('c:/osn/Amazon0505WCCRECOLLISION_18115.dat');%5434.dat');%2717.txt');%7246.dat'); %18115.dat');%7246.dat');%2717.txt');
%A=textread('c:/osn/youtubeWCCRECOLLISION_13559.dat');%2717.txt');
%A=textread('c:/osn/WikiTalkWCCRECOLLISION_6557.dat'); %43716.dat');%6557.dat');%27118.dat');%2717.txt');
%A=textread('c:/wanghao/RE/WikiTalkWCCRECOLLISION_20.dat');
%A=textread('c:/wanghao/RE/WikiTalkWCCRECOLLISION_10.dat');
%A=textread('c:/wanghao/RE/Amazon0505WCCRECOLLISION_20.dat');
A=textread('c:/wanghao/RE/youtubeWCCRECOLLISION_15.dat');

N=410000;%population(1);
EC=20^2;
mud=A(:,3);
mude=A(:,4);
C=A(:,1);
Gamma=mean(mude./mud)
y=mude./mud;
n=sqrt(2*EC*N/Gamma)
x=n^2./(2.*C);

N=mean(x.*y)

mu_x=mean(x);
mu_y=mean(y);

v_x=var(x);
v_y=var(y);

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=mu_x^2*v_y
t1=mu_y^2*v_x
t3=mean(delta_x.^2.*delta_y.^2);
[t3,v_x*v_y]
t4= 2*mu_x*mu_y*E11
t5= 2*mu_x*E12
t6= 2*mu_y*E21
t7= E11^2

mean(x.*y)
mean(x)*mean(y)
cov=mean(x.*y)-mean(x)*mean(y)
2*mu_x*mu_y*cov
cor=(mean(x.*y)-mean(x)*mean(y))/(v_x*v_y)^(0.5)

%cooreltion is negative
subplot(3,2,1);
plot(x,y,'o')

% gamma^2 is correct
subplot(3,2,2);
plot(mude./mud, A(:,2).^2+1,'o');

grid on;
gy=v_y^(0.5)/mu_y
gx=v_x^(0.5)/mu_x

v_xy_observed=var(x.*y)
v_xy_goodman=mu_x^2*v_y+mu_y^2*v_x+ t3+   2*mu_x*mu_y*E11  +2*mu_x*E12+2*mu_y*E21-E11^2
v_xy_independent=mu_x^2*v_y+mu_y^2*v_x+v_x*v_y
v_xy_t4=mu_x^2*v_y+mu_y^2*v_x+t4
v_xy_gamma=mu_y^2*v_x;%mu_x^2*mu_y^2*(gx^2);%+gy^2)
gx/gy
subplot(3,2,3);
bar([(v_xy_observed)^(0.5)/N,(v_xy_goodman)^(0.5)/N,v_xy_independent^(0.5)/N,(v_xy_t4)^(0.5)/N, (v_xy_gamma)^(0.5)/N, t1^(0.5)/N]);
grid on;

%seven terms of V(xy)
subplot(3,2,4);
bar([v_xy_observed,t1,t2,t3,t4,t5,t6,t7, N^2/mean(C)]);
%xlabel('observed', 't1','t2','t3','t4','t5','t6','t7', 'N^2/mean(C)');

subplot(3,2,5);
plot(C,A(:,2),'o')
hist(C)

subplot(3,2,6);
bar([N^2/EC, N^2/mean(C),t1])
bar([n^4/(4*EC^3),v_x ])
v_x/(n^4/(4*EC^3))
var(C)/EC
