# Function that calculates pagerank by the power method, and in doing so 
# calculates the state probabilities at each step along the way.
# 
# "g" is a graph object in the igraph library.  
#
# Examples:
# g2 = graph.formula(A++B,A++C,B-+C)
# g2 = add.edges(g2,c(1,1))
# 
# start is an initial vector of probabilities (default is 100% in 1st vertex)
# verbose=TRUE will print intermediate values and the transition matrix.
# display=TRUE will plot the intermediate values
# damp is the probability of following edges.  damp=1 will always follow edges.
# damp=0.85 gives google pagerank (apparently).

my_pagerank = function(g,start=NULL,n.step=100,verbose=FALSE,display=TRUE,damp=1,...){
  A = as.matrix(get.adjacency(g))
  for(i in 1:nrow(A)) {
    if (sum(A[i,])==0 & damp==1) A[i,i]=1  # a vertex w/o edges is assumed to have self edge.
    if (sum(A[i,])==0 & damp<1) A[i,]=1/nrow(A)  # ... unless we have damping
    A[i,] = A[i,]/sum(A[i,])
  }
  A = damp*A + (1-damp)*matrix(1,nrow=nrow(A),ncol=ncol(A))/ncol(A)
  if (is.null(start)) start = c(1,rep(0,nrow(A)-1)) 
  y = x = rbind(start)
  printlist=2^unique(round(log(1:n.step,2)))
  printlist=printlist[printlist<n.step]
  printlist=c(printlist,n.step)
  for (i in 1:n.step) {
    x = x%*%A
    y = rbind(y,x)
    if (verbose & i%in%printlist) cat('step ',i,':',x,'\n')
  }
  if (verbose) {
    cat('Transition Matrix\n')
    print(A)
  }
  if(display) 
    matplot(y,type='b',pch=names(V(g)),ylab='probability',xlab='step',log='x',...)
  y[nrow(y),]
}
