source("http://www.stat.duke.edu/~pdh10/Teaching/594/Code/functions.r")

m<-c(8,7,6) 
r<-c(4,3,2) 
K<-length(m) 

set.seed(1)

## generate a random rank-r tensor of dimension m
A<-list()
for(k in 1:K){  A[[k]]<-matrix(rnorm(m[k]*r[k]),m[k],r[k]) } 
Z<-array(rnorm(prod(r)),r) 
M<-atrans(Z,A) 

## now let M be the "mean matrix" 
s2<-4
Y<- M + array(rnorm(prod(m)),m)*sqrt(s2) 
plot(M,Y) ; abline(0,1) 


hYo<-US.als(Y,r) 
names(hYa)
plot(hYa$IMP)
Mo<-atrans(hYo$S,hYo$U) 

hYt<-hosvd(Y,r)
Mt<-atrans(hYt$S,hYt$U)

plot(Mt,Mo)

mean( (Y-Mo)^2 ) 
mean( (Y-Mt)^2 )




## ALS versus truncated TSVD for 
## (a) approximation of Y 
## (b) estimation of M

for(r1 in 1:min(m))
{
  r0<-rep(r1,3)
  hYt<-hosvd(Y,r0) 
  Mt<- atrans(hYt$S,hYt$U) 

  hYo<-US.als(Y,r0)  
  Mo<- atrans(hYo$S,hYo$U) 

  plot(M,Mt)  
  points(M,Mo,col="green") 
  abline(0,1)  

  cat(r1, mean( (Y-Mt)^2),  mean( (Y-Mo)^2)," ", 
          mean( (M-Mt)^2),  mean( (M-Mo)^2),"\n")

}


# Is the best OLS M better than best TSVD? 
r0<-c(2,2,2)  

hYt<-hosvd(Y,r0)
Mt<- atrans(hYt$S,hYt$U)
mean( (M-Mt)^2)

hYo<-US.als(Y,r0)
Mo<- atrans(hYo$S,hYo$U)
mean( (M-Mo)^2)