library(lda)
Y<-as.matrix(read.table("CaltechY.dat"))
# replace the diagonals with 1s
diag(Y)<-1
#net<-network(Y,directed=FALSE)

link.predict<- function(fit, Y){
  N<- nrow(fit$net.assignments.left)
  theta<- thetahat(fit,Y)
  EY<- matrix(0,N,N)
  
  for(i in 1:N){
    for(j in (1:N)){
      if(i==j) EY[i,j]<- 0 else EY[i,j]<- theta[fit$net.assignments.right[i,j]+1 , fit$net.assignments.left[i,j]+1]
   }
 }

return(EY)
}

thetahat<- function(fit, Y){
   G<- nrow(fit$blocks.pos)
   theta.pos<- theta.neg<- matrix(0, G, G)
  
   for(g1 in 1:G){
      for(g2 in 1:G){
	  theta.pos[g1, g2]<- sum((fit$net.assignments.left==(g1-1))*(fit$net.assignments.right==(g2-1))*(Y))
	  theta.neg[g1, g2]<- sum((fit$net.assignments.left==(g1-1))*(fit$net.assignments.right==(g2-1))*(1-Y))
      }
  }
  theta.pos/(theta.pos + theta.neg)  
}

fit2<- mmsb.collapsed.gibbs.sampler(network=Y,K=2, num.iterations=5100, alpha=1/2, beta.prior=list(matrix(1,2,2), matrix(1,2,2)), burnin=as.integer(50000))
fit3<- mmsb.collapsed.gibbs.sampler(network=Y,K=3, num.iterations=5100, alpha=1/3, beta.prior=list(matrix(1,3,3), matrix(1,3,3)), burnin=as.integer(50000))
fit4<- mmsb.collapsed.gibbs.sampler(network=Y,K=4, num.iterations=5100, alpha=1/4, beta.prior=list(matrix(1,4,4), matrix(1,4,4)), burnin=as.integer(50000))
fit5<- mmsb.collapsed.gibbs.sampler(network=Y,K=5, num.iterations=5100, alpha=1/5, beta.prior=list(matrix(1,5,5), matrix(1,5,5)), burnin=as.integer(50000))

## Deliberately Overfit
fit9<- mmsb.collapsed.gibbs.sampler(network=Y,K=9, num.iterations=5100, alpha=1/9, beta.prior=list(matrix(1,9,9), matrix(1,9,9)), burnin=as.integer(50000))
fit10<- mmsb.collapsed.gibbs.sampler(network=Y,K=10, num.iterations=5100, alpha=1/10, beta.prior=list(matrix(1,10,10), matrix(1,10,10)), burnin=as.integer(50000))
fit11<- mmsb.collapsed.gibbs.sampler(network=Y,K=11, num.iterations=5100, alpha=1/11, beta.prior=list(matrix(1,11,11), matrix(1,11,11)), burnin=as.integer(50000))

l2<- link.predict(fit2, Y)
l3<- link.predict(fit3, Y)
l4<- link.predict(fit4, Y)
l5<- link.predict(fit5, Y)
l9<- link.predict(fit9, Y)
l10<- link.predict(fit10, Y)
l11<- link.predict(fit11, Y)