#'
#' @import MASS
#'
NULL

library(MASS)

#' Project a distance matrix in a euclidean space (NMDS).
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
nmds = function(distances,
                maxit = 50, trace = TRUE,
                tol = 0.001, p = 2) {
  if (inherits(distances,"matrix"))
    distances=as.dist(distances)

  stopifnot(inherits(distances,"dist"))

  k = attr(distances,"Size") - 1

  n = isoMDS(distances,
             k=k,
             maxit = maxit,
             trace = trace,
             tol = tol,
             p=p)

  p = n$points
  attr(p,"stress") = n$stress
  attr(x,"method")="nmds"

  return(p)
}

#' Project a distance matrix in a euclidean space (PCOA).
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
pcoa = function(distances) {
  if (inherits(distances,"matrix"))
    distances=as.dist(distances)

  stopifnot(inherits(distances,"dist"))

  k = attr(distances,"Size") - 1
  x = cmdscale(distances,k=k)
  attr(x,"method")="pcoa"

  return()
}

#' Project a set of points in a euclidean space (PCA).
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
pca = function(data,scale=FALSE) {
  k = min(nrow(data)-1,ncol(data))
  p = prcomp(data,
             retx = TRUE,
             center = TRUE,
             scale. = scale,
             tol=0)

  x = p$x
  attr(x,"method")="pca"

  return(x)
}

#' Double centering of a matrix.
#'
#' colSums and rowSums of the returned matrix are all equal to zero.
#'
#' Inspired from the algorithm described in stackoverflow
#' \url{https://stackoverflow.com/questions/43639063/double-centering-in-r}
#'
#' @export
bicentered = function(m) {

  # compute the row-wise and column-wise mean matrices
  R = m*0 + rowMeans(m)
  C = t(m*0 + colMeans(m))

  # substract them and add the grand mean
  return(m - R - C + mean(m[]))
}

#' Test if the distance matrix is euclidean.
#'
#' Actually a simplified version of the ADE4 implementation
#' (\code{\link[ade4]{is.euclid}}).
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
is.euclid = function (distances, tol = 1e-07)
{
  if (!inherits(distances, "dist"))
    stop("Object of class 'dist' expected")

  if (any(distances < tol))
    warning("Zero distance(s)")

  distances <- as.matrix(distances)
  n <- ncol(distances)
  delta <- -0.5 * ProcMod::bicenter(distances * distances)
  lambda <- eigen(delta, symmetric = TRUE, only.values = TRUE)$values
  w0 <- lambda[n]/lambda[1]
  return((w0 > -tol))
}

#' Project a distance matrix in a euclidean space.
#'
#' Project a set of points defined by a distance matrix in
#' an eucleadean space. If the used distance is a metric,
#' this is done using the \code{\link[ProcMod]{pcoa}} function,
#' otherwise the \code{\link[ProcMod]{nmds}} is used.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
dist2orthospace = function(distances, tol = 1e-07) {
  if (inherits(distances,"matrix"))
    distances=as.dist(distances)

  if (ProcMod::is.euclid(distances,tol = tol))
    return(ProcMod::pcoa(distances))
  else
    return(ProcMod::nmds(distances))
}