vec.norm <- function(vec){
	return(vec/matrix(sqrt(rowSums(vec^2)),nrow(vec),ncol(vec)))
}

Hotelling.test <- function(SSef,SSer,dfef,dfer,exact=F){
	#D'apres Claude 2008, p.252 
	if (exact==T) warning("exact F-test not implemented! returns approx F")

	p <- qr(SSef+SSer)$rank
	k <- dfef
	w <- dfer
	s <- min(k,p)
	m <- (w-p-1)/2
	t1 <- (abs(p-k)-1)/2
	Ht <- sum(diag(SSef%*%ginv(SSer)))
	Fapprox <- Ht*(2*(s*m+1))/(s^2*(2*t1+s+1))
	ddfnum <- s*(2*t1+s+1)
	ddfden <- 2*(s*m+1)
	pval <- 1-pf(Fapprox,ddfnum,ddfden)
	
	unlist(list("dfeffect"=dfef,"dferror"=dfer,"T2"=Ht,"F_approx"=Fapprox,
	"df1"=ddfnum,"df2"=ddfden,"p"=pval))
}

meas.error <- function(SSeff,SSerr,dfeff,dferr){
	# Calcul des MS de l'erreur et de l'effet
	s2_within <- MS_within <- SSerr/dferr
	MS_among <- SSeff/dfeff
	# variance components: E[MS_among] = n.session*s2_among + s2_within
	s2_among <- (MS_among-MS_within)/n.session
	#Calcul de l'erreur de mesure selon Bayles and Burnes 1990
	me <- s2_within/(s2_within+s2_among)*100
	list(s2_among=s2_among,s2_within=s2_within,me=me)

}

add.mshape.links <-function(msh,myLinks,lwd=2,col="gray"){

for (i in 1:nrow(myLinks)){
	segments(msh[myLinks[,1],1],msh[myLinks[,1],2],
			msh[myLinks[,2],1],msh[myLinks[,2],2],lwd=lwd,col=col)
			}

}

back.trans <- function(b,eigenvec,n.land,n.dim){
	back.b <- b%*%t(eigenvec)
	return(matrix(back.b,n.land,n.dim,byrow=T))
}

reg.proj <- function(x,b){
	return(x %*% b %*% sqrt(solve(t(b) %*% b)))
	}

all.test.Manova <- function(mAoV){
	pillai <- summary(manova(mAoV),test="Pillai")
	Hoteling <- summary(manova(mAoV),test="Hotelling")
	Roy <- summary(manova(mAoV),test="Roy")
	stats <- rbind(pillai$stats[1,],Hoteling$stats[1,],Roy$stats[1,])
	stats <- rbind(stats,pillai$stats[2,])
	row.names(stats) <- c("Pillai","Hotelling-Lawley","Roy","Residuals")
	cnms <- colnames(stats)
	cnms[2] <- "statistics.value"
	colnames(stats) <- cnms
	return(stats)
	
}

myEigenDec.Real <- function(Sigma,x,dim){
	eig.dek <- eigen(Sigma)
	eig.vec.real <- matrix(as.double(eig.dek$vectors[,1:dim]),nrow(Sigma),dim)
	eig.val.real <- as.double(eig.dek$values[1:dim]) 
	scr <- x%*%(eig.vec.real*matrix(eig.val.real,nrow(Sigma),dim,byrow=T))
	list(values=eig.val.real,vectors=eig.vec.real,scores=scr)
}

angle <- function(vec1,vec2,dim){
	#! vec1 et vec2 doivent être des matrices [n x 2k-4]
	if (!is.matrix(vec1)) vec1 <- matrix(vec1,1,length(vec1))
	if (!is.matrix(vec2)) vec2 <- matrix(vec2,1,length(vec2))
	if (ncol(vec1)!=dim) vec1 <- t(vec1)
	if (ncol(vec2)!=dim) vec2 <- t(vec2)
	
	vec1 <- vec1/matrix(sqrt(diag(vec1%*%t(vec1))),1,dim)
	vec2 <- vec2/matrix(sqrt(diag(vec2%*%t(vec2))),1,dim)
	a <- abs(diag(vec1%*%t(vec2)))
	return(acos(a)*180/pi)
}

rangl <- function(value,dim,thousands){
	thousands <- 100;
	cum.a <- rep(1000,1)
	count <- 0
	for (r in 1:thousands){
		a <- replicate(1000,angle(rnorm(dim),rnorm(dim),dim))
		cum.a <- cum.a+sort(a)
		count <- count + sum(a<=value)
		}
	return(list(p.val=(count + 1)/(thousands*1000+1),alpha=cum.a/thousands))
}

Hotellingsp<-function(SSef, SSer, dfef, dfer, exact=F){library(MASS)
  #  p corresponds to the number of shape space dimensions.
  p <- qr(SSef+SSer)$rank
  k<-dfef; w<-dfer
  s<-min(k,p)
  m<-(w-p-1)/2
  t1<-(abs(p-k)-1)/2
  Ht<-sum(diag(SSef%*%ginv(SSer)))
  Fapprox<-Ht*(2 * (s*m+1))/(s^2*(2*t1+s+1))
  ddfnum<-s*(2*t1+s+1)
  ddfden<-2*(s*m+1)
  pval= 1-pf(Fapprox, ddfnum, ddfden)
  
  if (exact)
  {b<-(p+2*m)*(k+2*m)/((2*m+1)*(2*m-2))
  c1<-(2+(p*k+2)/(b-1))/(2*m)
  Fapprox<-((4+(p*k+2)/(b-1))/(p*k))*(Ht/c1)
  ddfnum<-p*k
  ddfden<-4+(p*k+2)/(b-1)}
  
  unlist(list("dfeffect"=dfef,"dferror"=dfer,"T2"=Ht,
              "Approx_F"=Fapprox,"df1"=ddfnum,"df2"=ddfden,"p"=pval))}

