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) 

## check rank
svd(mat(M,1))$d

S<-hosvd(M)$S
zapsmall(S) 

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

## truncated TSVD for estimation of M
Mhat<-Y
mean( (M-Mhat)^2) 

r0<-r 
hY<-hosvd(Y,r0)
Mhat<- atrans(hY$S,hY$U)
plot(M,Mhat) ; abline(0,1)
mean( (M-Mhat)^2) 

## what is the "best" rank for estimating M?
for(r1 in 1:min(m))
{
  r0<-rep(r1,3) 
  hY<-hosvd(Y,r0)
  Mhat<- atrans(hY$S,hY$U)
  plot(M,Mhat) ; abline(0,1)
  cat(r1, mean( (M-Mhat)^2),"\n") 

}

par(mfrow=c(1,3))
for(k in 1:K)
{
  plot(svd(mat(Y,k))$d,type="h") 
}


hY<-hosvd(Y,c(1,2,3))
Mhat<- atrans(hY$S,hY$U)
mean( (M-Mhat)^2)