##### resample
resample <- function(x, ...) x[sample.int(length(x), ...)]
#####


##### sample from a Wishart dist
rwish<-function (S0, nu) 
{
    sS0 <- chol(S0)
    Z <- matrix(rnorm(nu * dim(S0)[1]), nu, dim(S0)[1]) %*% sS0
    t(Z) %*% Z
}
#####


##### sample from an array normal distribution
ranorm<-function(M,Sig)
{
  M+atrans( array(rnorm(prod(dim(M))),dim(M)) , lapply(Sig,mhalf)  )
}
#####

#####
ryi.fc<-function(idx,Y,M=array(0,dim=dim(Y)),S=NULL,P=lapply(S,solve))
{
  #samples from the full conditional distribution of an element 
  #of a gaussian array with separable covariance structure

  Z<-Y-M 
  qidx<-1 ; lidx<-1
  for(k in 1:length(P)) 
  {
    qidx<-qidx*P[[k]][idx[k],idx[k]] 
    lidx<-kronecker(P[[k]][,idx[k]],lidx)
  }
    
 lidx<- sum(lidx*c(Z)) - qidx*Z[matrix(idx,1,length(P))]   
 
 v<-1/qidx ; m<- M[matrix(idx,1,length(P))]  - lidx/qidx
 rnorm(1,m,sqrt(v))
}
#####

##### array and matrix functions
tr<-function(A) { sum(diag(A)) }
vec<-function(A){c(A)}
mat<-function(A,k){ K<-length(dim(A)) ; t(apply(A,k,"c"))  }  

amprod<-function(A,M,k)
{ #multiply k-vectors by M  and reform array
  K<-length(dim(A))
  AM<-M%*%mat(A,k) 
  AMA<-array(AM, dim=c(dim(M)[1],dim(A)[-k]) )
  aperm(AMA,  match(1:K,c(k,(1:K)[-k]) ) )
}

atrans<-function(A,B)
{
  #transform along each dimension
  X<-A
  for(k in 1:length(B)) { X<-amprod(X,B[[k]],k) }
  X
}

mhalf<-function(M) { 
  #symmetric square  root of a pos def matrix
  tmp<-eigen(M)
  tmp$vec%*%sqrt(diag(tmp$val))%*%t(tmp$vec)
                    }

mstrans<-function(Z,S)
{
  X<-Z
  K<-length(S)
  for(k in 1:K)
  {
    tmp<-eigen(S[[k]]) ; A<-tmp$vec%*%sqrt(diag(tmp$val))%*%t(tmp$vec)
    X<-aperm( apply( X , c(1:K)[-k], function(x) { A%*%x }) ,
              match(1:K,c(k,(1:K)[-k]) ))
 }
  X
}
#####





##### 
gof.stats<-function(Y)
{
   E<-Y -  outer( apply(Y,c(1,2,3),mean,na.rm=TRUE) ,rep(1,t))  
   E[is.na(E)]<-0
   SS<-list()
   for(k in 1:K)
   { 
     Ek<-mat(E,k)  
     SSk<-Ek%*%t(Ek) 
     SS[[k]]<-SSk/tr(SSk) 
   }

   lev<-NULL
   for(k in 1:(K-1))
   { 
     ev<-eigen(SS[[k]] )$val ; ev<-ev[ zapsmall(ev)>0 ]
     lev<-c(lev,sum(log(ev)) ) 
    }
   lev 
}
#####