pcoa <- function(D, correction="none", rn=NULL)
#
# Principal coordinate analysis (PCoA) of a square distance matrix D
# with correction for negative eigenvalues.
# 
# References:
# Gower, J. C. 1966. Some distance properties of latent root and vector methods
#    used in multivariate analysis. Biometrika. 53: 325-338.
# Gower, J. C. and P. Legendre. 1986. Metric and Euclidean properties of 
#    dissimilarity coefficients. J. Classif. 3: 5-48.
# Legendre, P. and L. Legendre. 1998. Numerical ecology, 2nd English edition. 
#    Elsevier Science BV, Amsterdam. [PCoA: Section 9.2]
#
#      Pierre Legendre, October 2007
{
centre <- function(D,n)
# Centre a square matrix D by matrix algebra
# mat.cen = (I - 11'/n) D (I - 11'/n)
{	One <- matrix(1,n,n)
	mat <- diag(n) - One/n
	mat.cen <- mat %*% D %*% mat
}

bstick.def <- function (n, tot.var = 1, ...)   # 'bstick.default' from vegan
{
    res <- rev(cumsum(tot.var/n:1)/n)
    names(res) <- paste("Stick", seq(len = n), sep = "")
    return(res)
}

# ===== The PCoA function begins here =====

# Preliminary actions
	D <- as.matrix(D)
	n <- nrow(D)
	epsilon <- sqrt(.Machine$double.eps)
	if(length(rn)!=0) {
		names <- rn
		} else {
		names <- rownames(D)
		}
	CORRECTIONS <- c("none","lingoes","cailliez")
	correct <- pmatch(correction, CORRECTIONS)
	if(is.na(correct)) stop("Invalid correction method")
	# cat("Correction method =",correct,'\n')

# Gower centring of matrix D
# delta1 = (I - 11'/n) [-0.5 d^2] (I - 11'/n)
	delta1 <- centre((-0.5*D^2),n)
	trace <- sum(diag(delta1))

# Eigenvalue decomposition
	D.eig <- eigen(delta1)
	
# Negative eigenvalues?
	min.eig <- min(D.eig$values)
	zero.eig <- which(abs(D.eig$values) < epsilon)
	D.eig$values[zero.eig] <- 0

# No negative eigenvalue
	if(min.eig > -epsilon) {   # Curly 1
		correct <- 1
		eig <- D.eig$values
		k <- length(which(eig > epsilon))
		rel.eig <- eig[1:k]/trace
		cum.eig <- cumsum(rel.eig) 
		vectors <- sweep(D.eig$vectors[,1:k], 2, sqrt(eig[1:k]), FUN="*")
		bs <- bstick.def(k)
		cum.bs <- cumsum(bs)

		res <- data.frame(eig[1:k], rel.eig, bs, cum.eig, cum.bs)
		colnames(res) <- c("Eigenvalues","Relative_eig","Broken_stick","Cumul_eig","Cumul_br_stick")
		rownames(res) <- 1:nrow(res)

		rownames(vectors) <- names
		colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
		note <- paste("There were no negative eigenvalues. No correction was applied")
		out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace))

# Negative eigenvalues present
	} else {   # Curly 1
		k <- n 
		eig <- D.eig$values
		rel.eig <- eig/trace
		rel.eig.cor <- (eig - min.eig)/(trace - (n-1)*min.eig) # Eq. 9.27 for a single dimension
		rel.eig.cor = c(rel.eig.cor[1:(zero.eig[1]-1)], rel.eig.cor[(zero.eig[1]+1):n], 0)
		cum.eig.cor <- cumsum(rel.eig.cor) 
		k2 <- length(which(eig > epsilon))
		k3 <- length(which(rel.eig.cor > epsilon))
		vectors <- sweep(D.eig$vectors[,1:k2], 2, sqrt(eig[1:k2]), FUN="*")
		# Only the eigenvectors with positive eigenvalues are shown

# Negative eigenvalues: three ways of handling the situation
	if((correct==2) | (correct==3)) {   # Curly 2

		if(correct == 2) {   # Curly 3
			# Lingoes correction: compute c1, then the corrected D
			c1 <- -min.eig
			note <- paste("Lingoes correction applied to negative eigenvalues: D' = -0.5*D^2 -",c1,", except diagonal elements")
			D <- -0.5*(D^2 + 2*c1)

			# Cailliez correction: compute c2, then the corrected D
			} else if(correct == 3) {
				delta2 <- centre((-0.5*D),n)
				upper <- cbind(matrix(0,n,n), 2*delta1)
				lower <- cbind(-diag(n), -4*delta2)
				sp.matrix <- rbind(upper, lower)
				c2 <- max(Re(eigen(sp.matrix, symmetric=FALSE, only.values=TRUE)$values))
				note <- paste("Cailliez correction applied to negative eigenvalues: D' = -0.5*(D +",c2,")^2, except diagonal elements")
				D <- -0.5*(D + c2)^2
			}   # End curly 3

	diag(D) <- 0
	mat.cor <- centre(D,n)
	toto.cor <- eigen(mat.cor)
	trace.cor <- sum(diag(mat.cor))

	# Negative eigenvalues present?
	min.eig.cor <- min(toto.cor$values)
	zero.eig.cor <- which((toto.cor$values < epsilon) & (toto.cor$values > -epsilon))
	toto.cor$values[zero.eig.cor] <- 0

	# No negative eigenvalue after correction: result OK
	if(min.eig.cor > -epsilon) {   # Curly 4
		eig.cor <- toto.cor$values
		rel.eig.cor <- eig.cor[1:k]/trace.cor
		cum.eig.cor <- cumsum(rel.eig.cor) 
		k2 <- length(which(eig.cor > epsilon))
		vectors.cor <- sweep(toto.cor$vectors[,1:k2], 2, sqrt(eig.cor[1:k2]), FUN="*")
		# bs <- broken.stick(k2)[,2]
		bs <- bstick.def(k2)
		bs <- c(bs, rep(0,(k-k2)))
		cum.bs <- cumsum(bs)

		# Negative eigenvalues still present after correction: incorrect result
		} else { 
		if(correct == 2) cat("Problem! Negative eigenvalues are still present after Lingoes",'\n') 
		if(correct == 3) cat("Problem! Negative eigenvalues are still present after Cailliez",'\n') 
		rel.eig.cor <- cum.eig.cor <- bs <- cum.bs <- rep(NA,n)
		vectors.cor <- matrix(NA,n,2)
		}   # End curly 4

	res <- data.frame(eig[1:k], eig.cor[1:k], rel.eig.cor, bs, cum.eig.cor, cum.bs)
	colnames(res) <- c("Eigenvalues", "Corr_eig", "Rel_corr_eig", "Broken_stick", "Cum_corr_eig", "Cum_br_stick")
	rownames(res) <- 1:nrow(res)

	rownames(vectors) <- names
	colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
	out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace, vectors.cor=vectors.cor, trace.cor=trace.cor))

	} else {   # Curly 2
			
	note <- "No correction was applied to the negative eigenvalues"
	bs <- bstick.def(k3)
	bs <- c(bs, rep(0,(k-k3)))
	cum.bs <- cumsum(bs)

	res <- data.frame(eig[1:k], rel.eig, rel.eig.cor, bs, cum.eig.cor, cum.bs)
	colnames(res) <- c("Eigenvalues","Relative_eig","Rel_corr_eig","Broken_stick","Cum_corr_eig","Cumul_br_stick")
	rownames(res) <- 1:nrow(res)

	rownames(vectors) <- names
	colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
	out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace))

			}   # End curly 2: three ways of handling the situation
		}	# End curly 1
	class(out) <- "pcoa"
	out
}	# End of PCoA