.packageName <- "amap"
#-------------------------------------------------------
#
#  Created       : 29/10/02
#  Last Modified : Time-stamp: <2003-03-05 15:22:26 lucas>
#
#  Description   : Principal component analysis
#                  
#  Author        : Antoine Lucas
#                  lucas@toulouse.inra.fr
#
#  Licence       : GPL 
#
#-------------------------------------------------------



acp <- function(x,center=TRUE,reduce=TRUE)
{   
    x    <- as.matrix(x)
    x    <- scale(x ,center = center, scale = FALSE)
    if (reduce == TRUE)
         {
          x    <- apply(x,2,function(u) { u/sd(u)}) 
         }
    EIG  <- eigen( t(x) %*% x,symmetric=TRUE) 
    V    <- EIG$vector    # ou bien: V=svd(x)$v

    val  <- sqrt(EIG$values)

    scores <- x %*% V

    V      <- data.frame(V)
    scores <- data.frame(scores)
    dimnames(V)[[2]] <- paste("Comp",1:dim(x)[2])
    dimnames(V)[[1]] <- dimnames(x)[[2]]
    dimnames(scores)[[1]] <- dimnames(x)[[1]]
    dimnames(scores)[[2]] <- paste("Comp",1:dim(x)[2])
    sdev   <- apply(scores,2,sd)    
    res  <- list(eig=val,sdev=sdev,scores=scores,loadings=V)
    class(res) <- "acp"
    res
}
pca <- acp

print.acp <- function(x, ...)
{
    #cat("Call:\n"); dput(x$call)
    cat("\nStandard deviations:\n")
    print(x$sdev, ...)
    cat("\nEigen values:\n")
    print(x$eig, ...)
    invisible(x)
}


# 
#   SECTION GRAPHIQUES
#

plot.acp <- function(x,i=1,j=2,text=TRUE,label='Composante ',col='darkblue',main='ACP des individus',...)
{
    U    <- x$scores
    XLAB <- paste(label,i)
    YLAB <- paste(label,j)
    plot.new()
    plot.window(range(U[,i]),range(U[,j]))
    axis(1,label=TRUE,tick=TRUE)
    axis(2,label=TRUE,tick=TRUE)
    box()
    title(xlab=XLAB,ylab=YLAB,main=main)
    if(text){
        text(U[,i],U[,j],col=col,...)   
    }
    else{
        points(U[,i],U[,j],col=col,...) 
    }
}

biplot.acp <- function(x,i=1,j=2,label='Composante ',col='darkblue',length=0.1,main='ACP des variables',...)
{
    U    <- x$loadings
    LIM  <- c(-1.3,1.3)
    XLAB <- paste(label,i)
    YLAB <- paste(label,j)

    # PLOT DES AXES
    plot.new()
    plot.window(LIM,LIM)
    axis(1,label=TRUE,tick=TRUE)
    axis(2,label=TRUE,tick=TRUE)
    box()
    title(xlab=XLAB,ylab=YLAB,main=main)


    # PLOT DU NOM DES FLECHES
    text(U[,i]*1.3,U[,j]*1.3,labels=dimnames(U)[[1]],col=col)   

    # PLOT DES FLECHES
    arrows(0,0,U[,i],U[,j],length = length,col=col)

    # CERCLE
    t2p <- 2 * pi * seq(0,1, length = 200)
    xc <- cos(t2p)
    yc <- sin(t2p)
    lines(xc,yc,col='darkblue')
}

# Graphique: Eboulis des valeurs propres
plot2.acp <- function(x,pourcent=FALSE,eigen=TRUE,label='Comp.',col='lightgrey',main='Eboulis des valeurs propres',ylab='Valeurs propres')
{
    if(eigen){ U <- x$eig }
    else { U <- x$sdev }

    if(pourcent){U <- U/sum(U) }
    n     <- length(U)
    names <- paste(label,1:n)
    barplot(U,main=main,ylab=ylab,col=col,names.arg=names)
}


#-------------------------------------------------------
#
#  Created       : 30/10/02
#  Last Modified : Time-stamp: <2003-04-02 09:52:47 lucas>
#
#  Description   : Robust principal component analysis
#                  
#  Author        : Antoine Lucas
#                  lucas@toulouse.inra.fr
#
#  Licence       : GPL 
#
#-------------------------------------------------------

K <- function(u,kernel="gaussien") {
    switch(kernel,
        gaussien = (2*pi)^(-1/2) * exp(-u^2/2),
        quartic   = 15/16 * (1-u^2)^2 * (abs(u)<1),
        triweight = 35/32 * (1-u^2)^3 * (abs(u)<1),
        epanechikov = 3/4 * (1-u^2) *   (abs(u)<1),
        cosinus = pi/4 * cos (u*pi/2) * (abs(u)<1),
        uniform = 1/2 * (abs(u)<1),
    )
}

# Variance locale
W <- function(x,h,D=NULL,kernel="gaussien")
{
    x   <- as.matrix(x)
    n   <- dim(x)[1]
    p   <- dim(x)[2]
    if (is.null(D)) {
        D <- diag(1,p)
    }
    x <- as.vector(x)
    D <- as.vector(D)
    kernel <- substr(kernel,1,1)

    VarLoc <- .C(
                 "W",
                 as.double(x),
                 as.double(h),
                 as.double(D),
                 as.integer(n),
                 as.integer(p),
                 as.character(kernel),
                 res=double(p*p),
                 result = as.integer(1),
                 PACKAGE= "amap"
                 )

    if(VarLoc$result == 2)
      stop("Cannot allocate memory")
    if(VarLoc$result == 1)
      stop("Error")

    matrix(VarLoc$res,p)
}


WsansC <- function(x,h,D=NULL,kernel="gaussien")
{
    x   <- as.matrix(x)
    n   <- dim(x)[1]
    p   <- dim(x)[2]
    if (is.null(D)) {
        D <- diag(1,p)
    }

    som    <- 0
    result <- 0
    for(i in 1:(n-1) )
        {
        for(j in (i+1):n)
            {
            Delta <- x[i,]-x[j,]
            norm <- sqrt(t(Delta) %*%D %*%Delta)
            k <- K ( norm /h,kernel)   # K ( |Delta|/h )
            k <- as.numeric(k)
            som <- som + k
            result <- result + k * Delta %*% t(Delta)
            }
        }
    result /   som
}

varrob <- function(x,h,D=NULL,kernel="gaussien")
{
    x   <- as.matrix(x)
    x   <- scale(x, center = TRUE, scale = FALSE)
    n   <- dim(x)[1]
    p   <- dim(x)[2]
    if (is.null(D)) {
        D <- diag(1,p)
    }
    x <- as.vector(x)
    D <- as.vector(D)
    kernel <- substr(kernel,1,1)

    Calcul <- .C(
                 "VarRob",
                 as.double(x),
                 as.double(h),
                 as.double(D),
                 as.integer(n),
                 as.integer(p),
                 as.character(kernel),
                 res=double(p*p),
                 result = as.integer(1),
                 PACKAGE= "amap")

    if(Calcul$result == 2)
      stop("Cannot allocate memory")
    if(Calcul$result == 1)
      stop("Error")

    S <- matrix(Calcul$res,p)
    Sinv <- solve(S)
    solve ( Sinv - D / h)
}


varrobsansC <- function(x,h,D=NULL,kernel="gaussien")
{
    n   <- dim(x)[1]
    p   <- dim(x)[2]
    if (is.null(D)) {
        D <- diag(1,p)
    }
    x   <- as.matrix(x)
    x   <- scale(x ,center = TRUE, scale = FALSE)
    som <- 0
    res <- 0
    for (i in 1:n )
    {
        k <- K (  sqrt(t(x[i,]) %*%D %*% x[i,]) /h,kernel)
        k <- as.numeric(k)
        res <- res + k * x[i,] %*% t( x[i,])
        som <- som + k
    }
    S <- res / som
    Sinv <- solve(S)
    solve ( Sinv -  D / h )

}



acpgen <- function(x,h1,h2,center=TRUE,reduce=TRUE,kernel="gaussien")
{
    # CENTRONS, ET REDUISONS
    x    <- as.matrix(x)
    x    <- scale(x ,center = center, scale = FALSE)
    if (reduce == TRUE)
         {
          x    <- apply(x,2,function(u) { u/sd(u)}) 
         }

    # ESTIMATION DE W et VarRob
    n <- dim(x)[1]
    VarInv   <- solve(var(x)*(n-1)/n) # solve= inverser
    leU    <- varrob(x,h1,D=VarInv,kernel=kernel)
    leW    <- W(x,h2,D=VarInv,kernel=kernel)
    Winv   <- solve(leW) 


    # anal. spec de Var.W^-1 :
    EIG    <- eigen(leU %*% Winv)  
    V      <- EIG$vector

    #EIG    <- eigen( x %*% Winv %*% t(x)  )
    #U      <- EIG$vector
    #n      <- dim(x)[1]
    #p      <- dim(x)[2]
    #S      <- diag(Re(EIG$values),n)   
    #S1     <- diag(Re(1/EIG$values),n)
    #S      <- sqrt(S[,1:p])
    #S1     <- sqrt(S1[,1:p])
    #V      <- t(x)%*% U%*% S1
    # X=U.S.V' -> V = X' U S^-1
    

    # AFFICHAGE DES RESULTATS


    scores <- x %*% Winv %*% V

    V      <- data.frame(V)
    scores <- data.frame(scores)
    dimnames(V)[[2]] <- paste("Comp",1:dim(x)[2])
    dimnames(V)[[1]] <- dimnames(x)[[2]]
    dimnames(scores)[[1]] <- dimnames(x)[[1]]
    dimnames(scores)[[2]] <- paste("Comp",1:dim(x)[2])
    eig    <- sqrt(EIG$values)
    sdev   <- apply(scores,2,sd)    
    res    <- list(eig=eig,sdev=sdev,scores=scores,loadings=V)
    class(res) <- "acp"
    res
}


acprob <- function(x,h=1,center=TRUE,reduce=TRUE,kernel="gaussien")
{   
    x    <- as.matrix(x)
    x    <- scale(x ,center = center, scale = FALSE)
    if (reduce == TRUE)
         {
          x    <- apply(x,2,function(u) { u/sd(u)}) 
         }
    EIG  <- eigen( varrob(x,h),symmetric=TRUE) 
    V    <- EIG$vector    # ou bien: V=svd(x)$v

    val  <- sqrt(EIG$values)

    scores <- x %*% V

    V      <- data.frame(V)
    scores <- data.frame(scores)
    dimnames(V)[[2]] <- paste("Comp",1:dim(x)[2])
    dimnames(V)[[1]] <- dimnames(x)[[2]]
    dimnames(scores)[[1]] <- dimnames(x)[[1]]
    dimnames(scores)[[2]] <- paste("Comp",1:dim(x)[2])
    sdev   <- apply(scores,2,sd)    
    res  <- list(eig=val,sdev=sdev,scores=scores,loadings=V)
    class(res) <- "acp"
    res
}
dist <- function(x, method="euclidean", diag=FALSE, upper=FALSE)
{
    ## account for possible spellings of euclid?an
    if(!is.na(pmatch(method, "euclidian")))
	method <- "euclidean"

    METHODS <- c("euclidean", "maximum",
                 "manhattan", "canberra", "binary","pearson","correlation")
    method <- pmatch(method, METHODS)
    if(is.na(method))
	stop("invalid distance method")
    if(method == -1)
	stop("ambiguous distance method")

    N <- nrow(x <- as.matrix(x))
    d <- .C("R_distance",
	    x = as.double(x),
	    nr= N,
	    nc= ncol(x),
	    d = double(N*(N - 1)/2),
	    diag  = as.integer(FALSE),
	    method= as.integer(method),
            ierr=as.integer(0),
	    DUP = FALSE,
            NAOK=TRUE,
            PACKAGE="amap"
            )$d
    attr(d, "Size") <- N
    attr(d, "Labels") <- dimnames(x)[[1]]
    attr(d, "Diag") <- diag
    attr(d, "Upper") <- upper
    attr(d, "method") <- METHODS[method]
    attr(d, "call") <- match.call()
    class(d) <- "dist"
    return(d)
}

names.dist <- function(x) attr(x, "Labels")

"names<-.dist" <- function(x, value)
{
    if(length(value) != attr(x, "Size"))
	stop("invalid names for dist object")
    attr(x, "Labels") <- value
    x
}

## Because names(d) != length(d) for "dist"-object d, we need
format.dist <- function(x, ...) format(as.vector(x), ...)

as.matrix.dist <- function(x)
{
    size <- attr(x, "Size")
    df <- matrix(0, size, size)
    df[row(df) > col(df)] <- x
    df <- df + t(df)
    labels <- attr(x, "Labels")
    dimnames(df) <-
	if(is.null(labels)) list(1:size,1:size) else list(labels,labels)
    df
}


as.dist <- function(m, diag = FALSE, upper=FALSE)
{
    m <- as.matrix(m)

    retval <-  m[row(m) > col(m)]

    attributes(retval) <- NULL

    if(!is.null(rownames(m)))
        attr(retval,"Labels") <- rownames(m)
    else if(!is.null(colnames(m)))
        attr(retval,"Labels") <- colnames(m)

    attr(retval,"Size") <- nrow(m)
    attr(retval,"Diag") <- diag
    attr(retval,"Upper") <- upper
    attr(retval, "call") <- match.call()
    class(retval) <- "dist"
    retval
}


print.dist <- function(x, diag = NULL, upper = NULL, ...)
{
    if(is.null(diag))
	diag <-  if(is.null(a <- attr(x, "Diag"))) FALSE else a
    if(is.null(upper))
	upper <- if(is.null(a <- attr(x,"Upper"))) FALSE else a

    size <- attr(x, "Size")
    df <- as.matrix.dist(x)
    if(!upper)
	df[row(df) < col(df)] <- NA
    if(!diag)
	df[row(df) == col(df)] <- NA
    print(if(diag || upper) df else df[-1, -size], na = "", ...)
    invisible(x)
}
## Hierarchical clustering
##
## Created       : 18/11/02
## Last Modified : Time-stamp: <2003-05-28 10:33:31 lucas>
##
## This function is a "mix" of function dist and function hclust.
##
## Author : Antoine Lucas
##



hcluster <- function (x, method = "euclidean", diag = FALSE, upper = FALSE, link = "complete", members = NULL)
{
                                        # take from dist
  if (!is.na(pmatch(method, "euclidian"))) 
    method <- "euclidean"
  METHODS <- c("euclidean", "maximum", "manhattan", "canberra", 
               "binary","pearson","correlation")
  method <- pmatch(method, METHODS)
  if (is.na(method)) 
    stop("invalid distance method")
  if (method == -1) 
    stop("ambiguous distance method")
  N <- nrow(x <- as.matrix(x))
  

  
                                        #take from hclust
  METHODSLINKS <- c("ward", "single", "complete", "average", "mcquitty", 
                    "median", "centroid")
  
  link <- pmatch(link, METHODSLINKS)
  if (is.na(link)) 
    stop("invalid clustering method")
  if (link == -1) 
    stop("ambiguous clustering method")
    if (N < 2) 
        stop("Must have n >= 2 objects to cluster")
  if (is.null(members)) 
    members <- rep(1, N)
  if (length(members) != N) 
    stop("Invalid length of members")
  n <- N
  
  hcl <- .C("hcluster",
            x = as.double(x),
            nr = as.integer(n),
            nc = as.integer(ncol(x)),
            diag = as.integer(FALSE),
            method = as.integer(method), 
            iopt = as.integer(link),
            ia = integer(n),
            ib = integer(n),
            order = integer(n),
            crit = double(n),
            members = as.double(members),
            res  = as.integer (1),
            DUP = FALSE,
            NAOK=TRUE,
            PACKAGE= "amap")

  if(hcl$res == 2)
    stop("Cannot allocate memory")
  if(hcl$res == 3)
    stop("Missing values in distance Matrix")
  if(hcl$res == 1)
    stop("Error")
  
  tree <- list(merge = cbind(hcl$ia[1:(N - 1)],
                 hcl$ib[1:(N -  1)]),
               height = hcl$crit[1:(N - 1)],
               order = hcl$order, 
               labels = dimnames(x)[[1]],
               dist.method = METHODS[method],
               method = METHODSLINKS[link],
               call = match.call())
  class(tree) <- "hclust"
  tree
}
.noGenerics <- TRUE
.conflicts.OK <- TRUE

.onLoad <- .First.lib <- function(lib, pkg)
{
    library.dynam("amap", pkg, lib)
    have.mva <- "package:mva" %in% search()
    if(!have.mva) require("mva")
}

