# initial conditions
# Y, an m x n matrix 
# X, an m x n x p array of regressors (can be absent) 
# R, the integer rank for the latent factors
# S, the number of iterations of the Markov chain

source("http://www.stat.washington.edu/~hoff/Code/Misc/SVDMODEL/svdmodel.f.r")

##### starting values and setup
Y0<-Y; Y0[is.na(Y0)]<-mean(Y,na.rm=TRUE)
Z<-matrix(qnorm(rank(Y0,ties="random")/(1+prod(dim(Y)))), dim(Y)[1],dim(Y)[2])
rm(Y0)
if(exists("X") && dim(X)[3]>0) 
{
  xt<-NULL; for(j in seq(length=dim(X)[3])) {xt<-rbind(xt,c(X[,,j])) }
  ixtx<-solve(xt%*%t(xt)) 
  B<-ixtx%*%xt%*%c(Z)
}
if(!exists("X") || dim(X)[3]==0)
{ 
X<-array(dim=c(dim(Z)[1],dim(Z)[2],0)) ; B<-rep(0,0) 
}
U<-svd(Z-XB(X,B))$u[,1:R]%*%diag(sqrt(svd(Z-XB(X,B))$d[1:R]))
V<-svd(Z-XB(X,B))$v[,1:R]%*%diag(sqrt(svd(Z-XB(X,B))$d[1:R]))
C.U<-cov(U) ; C.V<-cov(V) ; E.U<-apply(U,2,mean) ; E.V<-apply(V,2,mean)
s2<-1
B.ps<-D.ps<-NULL ; Z.psum<-M.psum<-matrix(0,dim(Y)[1],dim(Y)[2])
odens<-round(S/1000)
#####


##### Gibbs sampler
for(s in 1:S) 
{

  if(dim(X)[3]>0) { B<-rB.fc(Z-U%*%t(V),xt,ixtx,s2) }

  Z<-rZ.fc(XB(X,B)+U%*%t(V),Y)
 
  U<-rF.fc(Z-XB(X,B),V,E.U,C.U,s2)
 
  V<-rF.fc(t(Z-XB(X,B)),U,E.V,C.V,s2)
 
  C.U<-rC.fc(U,E.U)

  C.V<-rC.fc(V,E.V)

  E.U<-rE.fc(U,C.U)

  E.V<-rE.fc(V,C.V)

  ### a little output
  if(s%%odens==0)
  {  
    D.ps<-rbind(D.ps, svd(U%*%t(V))$d[1:R] ) 
    B.ps<-rbind(B.ps,B)
    M.psum<-M.psum+U%*%t(V)  
    Z.psum<-Z.psum+Z
    cat(s,"  ",round(svd(U%*%t(V))$d[1:R],2),"  ",round(B,2),"  ",date(),"\n")  
  }
  ### 
}
#####

  M.pm<-M.psum/dim(D.ps)[1]   
  Z.pm<-Z.psum/dim(D.ps)[1]