"cc" =
function (X, Y)
{
    Xnames = dimnames(X)[[2]]
    Ynames = dimnames(Y)[[2]]
    ind.names = dimnames(X)[[1]]
    res = rcc(X, Y, 0, 0)
    return(res)
}


"comput" =
function (X, Y, res)
{
    X.aux = scale(X, center=TRUE, scale=FALSE)
    Y.aux = scale(Y, center=TRUE, scale=FALSE)
    X.aux[is.na(X.aux)] = 0
    Y.aux[is.na(Y.aux)] = 0

    xscores = X.aux%*%res$xcoef
    yscores = Y.aux%*%res$ycoef

    corr.X.xscores = cor(X, xscores, use = "pairwise")
    corr.Y.xscores = cor(Y, xscores, use = "pairwise")
    corr.X.yscores = cor(X, yscores, use = "pairwise")
    corr.Y.yscores = cor(Y, yscores, use = "pairwise")

    return(list(xscores = xscores, yscores = yscores,
        corr.X.xscores = corr.X.xscores,
        corr.Y.xscores = corr.Y.xscores,
        corr.X.yscores = corr.X.yscores,
        corr.Y.yscores = corr.Y.yscores))
}

"estim.regul" <-
function (X, Y, grid1 = NULL, grid2 = NULL, plt = TRUE)
{
    if (is.null(grid1)) {grid1 = seq(0.001, 1, length = 5)}
    if (is.null(grid2)) {grid2 = seq(0.001, 1, length = 5)}
    grid = expand.grid(grid1, grid2)
    res = apply(grid, 1, function(x) {loo(X, Y, x[1], x[2])})
    res.grid = cbind(grid, res)
    mat = matrix(res, nc = length(grid2))

    if (isTRUE(plt)) {
        img.estim.regul(list(grid1 = grid1, grid2 = grid2, mat = mat))
    }

    opt = res.grid[res.grid[, 3] == max(res.grid[, 3]), ]
    cat("  lambda1 = ", opt[[1]], "\n", " lambda2 = ", opt[[2]], "\n",
    "CV-score = ", opt[[3]], "\n")
    return(invisible(list(lambda1 = opt[[1]], lambda2 = opt[[2]],
    CVscore = opt[[3]], grid1 = grid1, grid2 = grid2, mat = mat)))
}

"img.estim.regul" <-
function (estim)
{
    grid1 = estim$grid1
    grid2 = estim$grid2
    mat = estim$mat
    nlevel = min(512, length(c(unique(mat))))

    par(mar = c(4.5, 4.5, 4.5, 6))
    image(grid1, grid2, mat, col = heat.colors(nlevel),
    xlab = expression(lambda[1]), ylab=expression(lambda[2]),
    main = expression(CV(lambda[1], lambda[2])), axes = FALSE,
    zlim = c(min(mat), max(mat)), oldstyle = TRUE)

    if (length(grid1) > 10) {
        grid1 = seq(min(grid1), max(grid1), length = 11)
    }
    if (length(grid2) > 10) {
        grid2 = seq(min(grid2), max(grid2), length = 11)
    }

    axis(1, at = grid1, labels = as.character(round(grid1, 4)))
    axis(2, at = grid2, labels = as.character(round(grid2, 4)))
    box()

    image.plot(legend.only = TRUE, nlevel = nlevel,
            zlim = c(min(mat), max(mat)), col = heat.colors(nlevel),
            legend.width = 1.8)
}
"img.matcor" <-
function (correl, type = 1)
{
    matcorX = correl$Xcor
    matcorY = correl$Ycor
    matcorXY = correl$XYcor
    lX = ncol(matcorX)
    lY = ncol(matcorY)
    def.par <- par(no.readonly = TRUE)
    if (type == 1) {
        par(mfrow = c(1, 1), pty = "s")
        image(1:(lX + lY), 1:(lX + lY), t(matcorXY[nrow(matcorXY):1, ]),
            zlim = c(-1, 1), main = "XY correlation",
            col = tim.colors(64),
            axes = FALSE, , xlab = "", ylab = "")
        box()
        abline(h = lY + 0.5, v = lX + 0.5, lwd = 2, lty = 2)
        image.plot(legend.only = TRUE, zlim = c(-1, 1),
            col = tim.colors(64), horizontal = TRUE)
    }
    else {
        layout(matrix(c(1, 2, 3, 3, 0, 0), ncol = 2, nrow = 3,
             byrow = TRUE), widths = 1, heights = c(0.8, 1, 0.06))

    #layout 1
        par(pty = "s", mar = c(2, 2, 2, 1))
        image(1:lX, 1:lX, t(matcorX[lX:1, ]), zlim = c(-1, 1),
            main = "X correlation", axes = FALSE, xlab = "", ylab = "",
            col = tim.colors(64))
        box()

        image(1:lY, 1:lY, t(matcorY[lY:1, ]), zlim = c(-1, 1),
            main = "Y correlation",
            col = tim.colors(64), axes = FALSE,
            xlab = "", ylab = "")
        box()

    #layout 2
        partXY = matcorXY[lX:1, (lX + 1):(lX + lY)]
        if (lX > lY) {
            partXY = matcorXY[(lX + lY):(lX + 1), 1:lX]
            lX = ncol(matcorY)
            lY = ncol(matcorX)
        }

        par(pty = "m", mar = c(5,4,2,3), mai = c(0.8, 0.5, 0.3, 0.4))
        image(1:lY, 1:lX, t(partXY), zlim = c(-1, 1),
            main = "Cross-correlation", axes = FALSE, xlab = "",
            ylab = "", col = tim.colors(64))
        box()

        image.plot(legend.only = TRUE, zlim = c(-1, 1),
            horizontal = TRUE, col = tim.colors(64),
            legend.width = 2.5)
    }
    par(def.par)
}

"loo" =
function (X, Y, lambda1, lambda2)
{
    n = nrow(X)
    xscore = vector(mode = "numeric", length = n)
    yscore = vector(mode = "numeric", length = n)
    for (i in 1:n) {
        Xcv = X[-i, ]
        Ycv = Y[-i, ]
        res = rcc(Xcv, Ycv, lambda1, lambda2)
        xscore[i] = t(X[i, ]) %*% res$xcoef[, 1]
        yscore[i] = t(Y[i, ]) %*% res$ycoef[, 1]
    }
    cv = cor(xscore, yscore, use = "pairwise")
    return(invisible(cv))
}

"plt.cc" <-
function (res, d1 = 1, d2 = 2, int = 0.5, type = "b", ind.names=NULL, var.label=FALSE, Xnames=NULL, Ynames=NULL)
{
    par(mfrow = c(1, 1), pty = "s")
    if (type == "v")
        plt.var(res, d1, d2, int, var.label, Xnames, Ynames)
    if (type == "i")
        plt.indiv(res, d1, d2, ind.names)
    if (type == "b") {
        def.par <- par(no.readonly = TRUE)
        layout(matrix(c(0, 0, 1, 2, 0, 0), ncol = 2, nrow = 3,
        byrow = TRUE), widths = 1, heights = c(0.1, 1, 0.1))
        par(pty = "s", mar = c(4, 4.5, 0, 1))
        plt.var(res, d1, d2, int, var.label, Xnames, Ynames)
        plt.indiv(res, d1, d2, ind.names)
        par(def.par)
    }
}

"plt.indiv" <-
function (res, d1, d2, ind.names = NULL)
{
    if (is.null(ind.names)) ind.names = res$names$ind.names
    if (is.null(ind.names)) ind.names = 1:nrow(res$scores$xscores)
    plot(res$scores$xscores[, d1], res$scores$xscores[, d2],
        type = "n", main = "", xlab = paste("Dimension ", d1),
        ylab = paste("Dimension ", d2))
    text(res$scores$xscores[, d1], res$scores$xscores[, d2],
        ind.names)
    abline(v = 0, h = 0, lty=2)
}

"plt.var" =
function (res, d1, d2, int = 0.5, var.label = FALSE,
                     Xnames = NULL, Ynames = NULL)
{
    if (!var.label) {
        plot(0, type = "n", xlim = c(-1, 1), ylim = c(-1, 1),
        xlab = paste("Dimension ", d1), ylab = paste("Dimension ", d2))
        points(res$scores$corr.X.xscores[, d1],
               res$scores$corr.X.xscores[, d2],
               pch = 20, cex = 1.2, col = "red")
        points(res$scores$corr.Y.xscores[, d1],
               res$scores$corr.Y.xscores[, d2],
               pch = 24, cex = 0.7, col = "blue")
    }
    else {
        if (is.null(Xnames)) Xnames = res$names$Xnames
        if (is.null(Ynames)) Ynames = res$names$Ynames
        plot(0, type = "n", xlim = c(-1, 1), ylim = c(-1, 1),
        xlab = paste("Dimension ", d1), ylab = paste("Dimension ", d2))
        text(res$scores$corr.X.xscores[, d1],
             res$scores$corr.X.xscores[, d2],
             Xnames, col = "red", font = 2)
        text(res$scores$corr.Y.xscores[, d1],
             res$scores$corr.Y.xscores[, d2],
             Ynames, col = "blue", font = 3)
    }

    abline(v = 0, h = 0)
    lines(cos(seq(0, 2 * pi, l = 100)), sin(seq(0, 2 * pi, l = 100)))
    lines(int * cos(seq(0, 2 * pi, l = 100)),
          int * sin(seq(0, 2 * pi, l = 100)))
}

"rcc" <-
function (X, Y, lambda1, lambda2)
{
    Xnames <- dimnames(X)[[2]]
    Ynames <- dimnames(Y)[[2]]
    ind.names <- dimnames(X)[[1]]
    Cxx <- var(X, na.rm = TRUE, use = "pairwise") + diag(lambda1,
        ncol(X))
    Cyy <- var(Y, na.rm = TRUE, use = "pairwise") + diag(lambda2,
        ncol(Y))
    Cxy <- cov(X, Y, use = "pairwise")
    res <- geigen(Cxy, Cxx, Cyy)
    names(res) <- c("cor", "xcoef", "ycoef")
    scores <- comput(X, Y, res)
    return(list(cor = res$cor, names = list(Xnames = Xnames,
        Ynames = Ynames, ind.names = ind.names), xcoef = res$xcoef,
        ycoef = res$ycoef, scores = scores))
}


"matcor" <-
function (X, Y) 
{
    matcorX = cor(X, use = "pairwise")
    matcorY = cor(Y, use = "pairwise")
    matcorXY = cor(cbind(X, Y), use = "pairwise")
    return(list(Xcor = matcorX, Ycor = matcorY, XYcor = matcorXY))
}

"geigen" <-
function (Amat, Bmat, Cmat) 
{
    Bdim <- dim(Bmat)
    Cdim <- dim(Cmat)
    if (Bdim[1] != Bdim[2]) 
	stop("BMAT is not square")
    if (Cdim[1] != Cdim[2]) 
	stop("CMAT is not square")
    p <- Bdim[1]
    q <- Cdim[1]
    s <- min(c(p, q))
    if (max(abs(Bmat - t(Bmat)))/max(abs(Bmat)) > 1e-10) 
	stop("BMAT not symmetric.")
    if (max(abs(Cmat - t(Cmat)))/max(abs(Cmat)) > 1e-10) 
	stop("CMAT not symmetric.")
    Bmat <- (Bmat + t(Bmat))/2
    Cmat <- (Cmat + t(Cmat))/2
    Bfac <- chol(Bmat)
    Cfac <- chol(Cmat)
    Bfacinv <- solve(Bfac)
    Cfacinv <- solve(Cfac)
    Dmat <- t(Bfacinv) %*% Amat %*% Cfacinv
    if (p >= q) {
        result <- svd2(Dmat)
        values <- result$d
        Lmat <- Bfacinv %*% result$u
        Mmat <- Cfacinv %*% result$v
    }
    else {
        result <- svd2(t(Dmat))
        values <- result$d
        Lmat <- Bfacinv %*% result$v
        Mmat <- Cfacinv %*% result$u
    }
    geigenlist <- list(values, Lmat, Mmat)
    names(geigenlist) <- c("values", "Lmat", "Mmat")
    return(geigenlist)
}

"svd2" <-
function (x, nu = min(n, p), nv = min(n, p), LINPACK = FALSE) 
{
    dx <- dim(x)
    n <- dx[1]
    p <- dx[2]
    svd.x <- try(svd(x, nu, nv, LINPACK))
    if (class(svd.x) == "try-error") {
        nNA <- sum(is.na(x))
        nInf <- sum(abs(x) == Inf)
        if ((nNA > 0) || (nInf > 0)) {
            msg <- paste("sum(is.na(x)) = ", nNA, "; sum(abs(x)==Inf) = ", 
						 nInf, ".  'x stored in .svd.x.NA.Inf'", sep = "")
            assign(".svd.x.NA.Inf", x, env = .GlobalEnv)
            stop(msg)
        }
        attr(x, "n") <- n
        attr(x, "p") <- p
        attr(x, "LINPACK") <- LINPACK
        .x2 <- c(".svd.LAPACK.error.matrix", ".svd.LINPACK.error.matrix")
        .x <- .x2[1 + LINPACK]
        assign(.x, x, env = .GlobalEnv)
        msg <- paste("svd failed using LINPACK = ", LINPACK, 
					 " with n = ", n, " and p = ", p, ";  x stored in '", 
					 .x, "'", sep = "")
        warning(msg)
        svd.x <- try(svd(x, nu, nv, !LINPACK))
        if (class(svd.x) == "try-error") {
            .xc <- .x2[1 + (!LINPACK)]
            assign(.xc, x, env = .GlobalEnv)
            stop("svd also failed using LINPACK = ", !LINPACK, 
				 ";  x stored in '", .xc, "'")
        }
    }
    svd.x
}