%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  demo script to illlustrate the page rank algorithm on a small 6 node graph
%          demo taken from Page Rank notes from Amy Langville 
%              http://www.cofc.edu/~langvillea/UDelIntroWebLinkAnalysis.pdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
format short;  % make output legible
%
% Make the adjacency matrix 
%
nPages=6;
A=spalloc(nPages,nPages,3*nPages);  % allocate space for Adjacency matrix with max 3 links per page
A([ 2 3],1)=1;
A([ 1 2 5 ], 3)=1;
A([ 5 6 ],4)=1;
A([ 4 6 ],5)=1;
A( 4,6)=1;

%
% show the Adjacency Matrix  as a Full matrix
disp('Adjacency Matrix');
disp(full(A))

pause;

%
%  Calculate the Transitional Probability matrix H
%
nLinks=sum(A)'; % calculate the number of outlinks for each page
a=~nLinks;  % find dead-end nodes
noLinks=find(a)
invLinks=1./nLinks;
invLinks(noLinks)=0;
D=spdiags(invLinks,0,nPages,nPages);  % create sparse Scaling matrix to turn A into near markov matrix
H=A*D;

%
% Display Full Transitional Probability matrix H
%
disp('Probabalistic Surfer Matrix')
disp(full((H)))

pause;

%
% calculate stocastic fix for dead-end nodes (random-teleporter fix)
%

e=ones(nPages,1); % create a colume vector of ones
v=1/nPages*e;  % calculate the uniform probability to jump to any page

disp('Stochastic Fix')
S=H+v*a'

pause;

%
% Calculate full Google matrix by adding small random probability to make
%   the matrix irreducible
%
alpha=.85;  % Google weighting between S and fully connected graph
disp('Full Google Matrix')
G=alpha*S+(1-alpha)*v*e';
disp(G);

pause;

%
%  now calculate steady do power method starting with uniform probability vector
s=alpha*a'+(1-alpha)*e';  % calculate rank-1 update vector
xOld=v; % start with uniform probability vector
xPlot=v;
its=1:100;
for k=its
   x=(alpha*H*xOld+v*(s*xOld));
   xPlot(:,k+1)=x;
   error=x-xOld;
   xOld=x;
   norms(k,:)=[norm(x) norm(error)];
end

x
[pageRank,pageSorted]=sort(x);
disp('page   Rank');
disp([ pageSorted(length(pageSorted):-1:1), pageRank(length(pageRank):-1:1) ])

figure
h=semilogy(its,norms(:,1),'r',its,norms(:,2),'b');
set(h,'linewidth',[2]);
grid;
set(gca,'fontsize',14,'fontweight','bold');
xlabel('Iterations');
ylabel('Norms');
legend('||x||','||e||')
title('Convergence of Power method for PageRank Algorithm')
print -depsc2 PageRankExampleConvergence.eps

figure
plot([1:6],xPlot(:,[2:9]),'k--');
hold on;
h=plot(xPlot(:,[1 10 100]));
set(h,'linewidth',[2]);
legend(h,'x_0','x_{10}','x_{100}');
grid;
set(gca,'fontsize',14,'fontweight','bold');
xlabel('Page Number');
ylabel('Page Rank');
title('Evolution of PageRank probability')
print -depsc2 PageRankVectorEvolution.eps