R/ltMatrices.R

In mvtnorm: Multivariate Normal and t Distributions

# R Header

###    Copyright (C) 2022- Torsten Hothorn
###
###    This file is part of the 'mvtnorm' R add-on package.
###
###    'mvtnorm' is free software: you can redistribute it and/or modify
###    it under the terms of the GNU General Public License as published by
###    the Free Software Foundation, version 2.
###
###    'mvtnorm' is distributed in the hope that it will be useful,
###    but WITHOUT ANY WARRANTY; without even the implied warranty of
###    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
###    GNU General Public License for more details.
###
###    You should have received a copy of the GNU General Public License
###    along with 'mvtnorm'.  If not, see <http://www.gnu.org/licenses/>.
###
###
###    DO NOT EDIT THIS FILE
###
###    Edit 'lmvnorm_src.w' and run 'nuweb -r lmvnorm_src.w'

# .lt

.lt <- function(J, diag = FALSE) {
    j <- seq_len(J)
    j1 <- rep(c(TRUE, FALSE), J)[- 2 * J]
    j2 <- c(rbind(rev(j), j))[- 2 * J]
    ret <- matrix(rep(j1, j2), nrow = J)
    if (!diag) diag(ret) <- FALSE
    ret
}
.ut <- function(J, diag = FALSE) {
    !.lt(J = J, diag = !diag)
}

# ltMatrices

ltMatrices <- function(object, diag = FALSE, byrow = FALSE, names = TRUE) {

    if (!is.matrix(object)) 
        object <- matrix(object, ncol = 1L)

    # ltMatrices input
    
    if (is.ltMatrices(object)) {
        cls <- class(object)                ### keep inheriting classes
        ret <- .reorder(object, byrow = byrow)
        class(ret) <- class(object)
        return(ret)
    }
    
    # ltMatrices dim
    
    J <- floor((1 + sqrt(1 + 4 * 2 * nrow(object))) / 2 - diag)
    if (nrow(object) != J * (J - 1) / 2 + diag * J)
        stop("Dimension of object does not correspond to lower 
              triangular part of a square matrix")
    
    # ltMatrices names
    
    nonames <- FALSE
    if (!isTRUE(names)) {
        if (is.character(names))
            stopifnot(is.character(names) &&
                      length(unique(names)) == J)
        else
            nonames <- TRUE
    } else {
        names <- as.character(seq_len(J))
    }

    if (!nonames) {
        idx <- matrix(1:J, nrow = J, ncol = J)
        idx <- idx[.lt(J, diag  = diag)]
        j <- seq_len(J - !(diag + 0L))
        idx2 <- rep(j, rev(j))
        rn <- paste(names[idx], names[idx2], sep = ".")
        if (byrow)
           rn <- rn[.c2r(J, diag = diag)]
        rownames(object) <- rn
    } else {        ### add later
        warning("ltMatrices objects should be properly named")
    }
    

    attr(object, "J")       <- J
    attr(object, "diag")    <- diag
    attr(object, "byrow")   <- byrow
    attr(object, "rcnames") <- names

    class(object) <- c("ltMatrices", class(object))
    object
}

# syMatrices

as.syMatrices <- function(x) {
    if (is.syMatrices(x))
        return(x)
    if (is.ltMatrices(x)) {
        class(x) <- gsub("ltMatrices", "syMatrices", class(x))
        return(x)
    }
    as.syMatrices(as.ltMatrices(x * .lt(nrow(x), diag = TRUE)))
}
syMatrices <- function(object, diag = FALSE, byrow = FALSE, names = TRUE) {
    if (inherits(object, "syMatrices"))
        class(object) <- gsub("syMatrices", "ltMatrices", class(object))
    as.syMatrices(ltMatrices(object = object, diag = diag, byrow = byrow, 
                             names = names))
}

# dim ltMatrices

dim.ltMatrices <- function(x) {
    J <- attr(x, "J")
    return(c(attr(x, "dim")[2L], J, J)) ### ncol(unclass(x)) may trigger gc
}
dim.syMatrices <- dim.ltMatrices

# dimnames ltMatrices

dimnames.ltMatrices <- function(x)
    return(list(attr(x, "dimnames")[[2L]], attr(x, "rcnames"), attr(x, "rcnames")))
dimnames.syMatrices <- dimnames.ltMatrices

# names ltMatrices

names.ltMatrices <- function(x) {
    return(attr(x, "dimnames")[[1L]])
}
names.syMatrices <- names.ltMatrices

# is.ltMatrices

is.ltMatrices <- function(x) inherits(x, "ltMatrices")
is.syMatrices <- function(x) inherits(x, "syMatrices")
as.ltMatrices <- function(x) UseMethod("as.ltMatrices")
as.ltMatrices.syMatrices <- function(x) {
    cls <- class(x)
    class(x) <- cls[which(cls == "syMatrices"):length(cls)]
    class(x)[1L] <- "ltMatrices"
    return(x)
}
as.ltMatrices.ltMatrices <- function(x) {
    cls <- class(x)
    class(x) <- cls[which(cls == "ltMatrices"):length(cls)]
    return(x)
}

# as.ltMatrices

as.ltMatrices.default <- function(x) {
    stopifnot(is.numeric(x))
    if (!is.matrix(x)) x <- matrix(x)    
    DIAG <- max(abs(diag(x) - 1)) > .Machine$double.eps
    DIAG <- DIAG & (nrow(x) > 1)
    lt <- x[.lt(nrow(x), diag = DIAG)]
    up <- x[.ut(nrow(x), diag = FALSE)]
    stopifnot(max(abs(up)) < .Machine$double.eps)
    nm <- rownames(x)
    if (!is.null(nm))
        return(ltMatrices(lt, diag = DIAG, names = nm))
    return(ltMatrices(lt, diag = DIAG))
}

# print ltMatrices

as.array.ltMatrices <- function(x, symmetric = FALSE, ...) {

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    x <- unclass(x)

    L <- matrix(1L, nrow = J, ncol = J)
    diag(L) <- 2L
    if (byrow) {
        L[.ut(J, diag = diag)] <- floor(2L + 1:(J * (J - 1) / 2L + diag * J))
        L <- t(L)
    } else {
        L[.lt(J, diag = diag)] <- floor(2L + 1:(J * (J - 1) / 2L + diag * J))
    }
    if (symmetric) {
        L[.ut(J)] <- 0L
        dg <- diag(L)
        L <- L + t(L)
        diag(L) <- dg
    }
    ret <- rbind(0, 1, x)[c(L), , drop = FALSE]
    class(ret) <- "array"
    dim(ret) <- d[3:1]
    dimnames(ret) <- dn[3:1]
    return(ret)
}

as.array.syMatrices <- function(x, ...)
    return(as.array.ltMatrices(x, symmetric = TRUE))

print.ltMatrices <- function(x, zero.print = ".", ...) {
    if (is.null(attr(x, "dimnames")[[2L]]))
        attr(x, "dimnames")[[2L]] <- as.character(seq_len(dim(x)[1L]))
    print(as.table(as.array(x)), zero.print = zero.print, ...)
}

print.syMatrices <- function(x, ...)
    print(as.array(x))

# reorder ltMatrices

.r2c <- function(J, diag = FALSE) {

    rL <- cL <- diag(0L, nrow = J)
    lt <- .lt(J, diag = diag)
    ut <- !lt
    diag(ut) <- !diag(ut)
    x <- seq_len(J * (J + c(-1, 1)[diag + 1L]) / 2)
    cL[ut] <- seq_len(length(x))
    x[t(cL)[lt]]
}

.c2r <- function(J, diag = FALSE) {

    rL <- cL <- diag(0L, nrow = J)
    lt <- .lt(J, diag = diag)
    ut <- !lt
    diag(ut) <- !diag(ut)
    x <- seq_len(J * (J + c(-1, 1)[diag + 1L]) / 2)
    rL[lt] <- seq_len(length(x))
    x[t(rL)[ut]]
}

.reorder <- function(x, byrow = FALSE) {

    ### only we call this function
    ### stopifnot(is.ltMatrices(x))
    if (attr(x, "byrow") == byrow) return(x)

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    x <- unclass(x)

    if (byrow) { ### row -> col order
        return(ltMatrices(x[.r2c(J, diag = diag), , drop = FALSE], 
                          diag = diag, byrow = FALSE, names = dn[[2L]]))
    }
    return(ltMatrices(x[.c2r(J, diag = diag), , drop = FALSE], 
                      diag = diag, byrow = TRUE, names = dn[[2L]]))
}

# subset ltMatrices

# .subset ltMatrices

.subset_ltMatrices <- function(x, i, j, ..., drop = FALSE) {

    if (drop) warning("argument drop is ignored")
    if (missing(i) && missing(j)) return(x)

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    x <- unclass(x) 

    if (!missing(j)) {

        if (is.character(j)) {
            stopifnot(all(j %in% dn[[2L]]))
            j <- match(j, dn[[2L]])
        }

        j <- (1:J)[j] ### get rid of negative indices

        if (length(j) == 1L && !diag) {
            return(ltMatrices(matrix(1, ncol = ncol(x), nrow = 1), diag = TRUE, 
                              byrow = byrow, names = dn[[2L]][j]))
        }
        L <- diag(0L, nrow = J)
        Jp <- sum(.ut(J, diag = diag))
        if (byrow) {
            L[.ut(J, diag = diag)] <- 1:Jp
            L <- L + t(L)
            diag(L) <- diag(L) / 2
            L <- L[j, j, drop = FALSE]
            L <- L[.ut(nrow(L), diag = diag)]
        } else {
            L[.lt(J, diag = diag)] <- 1:Jp
            L <- L + t(L)
            diag(L) <- diag(L) / 2
            L <- L[j, j, drop = FALSE]
            L <- L[.lt(nrow(L), diag = diag)]
        }
        if (missing(i)) {
            return(ltMatrices(x[c(L), , drop = FALSE], diag = diag, 
                              byrow = byrow, names = dn[[2L]][j]))
        }
        return(ltMatrices(x[c(L), i, drop = FALSE], diag = diag, 
                          byrow = byrow, names = dn[[2L]][j]))
    }
    return(ltMatrices(x[, i, drop = FALSE], diag = diag, 
                      byrow = byrow, names = dn[[2L]]))
}

### if j is not ordered, result is not a lower triangular matrix
"[.ltMatrices" <- function(x, i, j, ..., drop = FALSE) {
    if (!missing(j)) {
        if (is.character(j)) {
            stopifnot(all(j %in% dimnames(x)[[2L]]))
            j <- match(j, dimnames(x)[[2L]])
        }
        if (all(j > 0)) {
            if (any(diff(j) < 0)) stop("invalid subset argument j")
        }
    }

    return(.subset_ltMatrices(x = x, i = i, j = j, ..., drop = drop))
}

"[.syMatrices" <- function(x, i, j, ..., drop = FALSE) {
    x <- as.syMatrices(x)
    ret <- .subset_ltMatrices(x = x, i = i, j = j, ..., drop = drop)
    class(ret)[1L] <- "syMatrices"
    ret
}

# lower triangular elements

Lower_tri <- function(x, diag = FALSE, byrow = attr(x, "byrow")) {

    if (is.syMatrices(x))
        x <- as.ltMatrices(x)
    adiag <- diag
    x <- ltMatrices(x, byrow = byrow)

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    if (diag == adiag)
        return(unclass(x)[,,drop = FALSE]) ### remove attributes

    if (!diag && adiag) {
        diagonals(x) <- 1
        return(unclass(x)[,,drop = FALSE]) ### remove attributes
    }

    x <- unclass(x)
    if (J == 1) {
        idx <- 1L
    } else {
      if (byrow)
          idx <- cumsum(c(1, 2:J))
      else
          idx <- cumsum(c(1, J:2))
    }
    return(x[-idx,,drop = FALSE])
}

# diagonals ltMatrices

diagonals <- function(x, ...)
    UseMethod("diagonals")

diagonals.ltMatrices <- function(x, ...) {

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    x <- unclass(x)

    if (!diag) {
        ret <- matrix(1, nrow = J, ncol = ncol(x))
        colnames(ret) <- dn[[1L]]
        rownames(ret) <- dn[[2L]]
        return(ret)
    } else {
        if (J == 1L) return(x)
        if (byrow)
            idx <- cumsum(c(1, 2:J))
        else
            idx <- cumsum(c(1, J:2))
        ret <- x[idx, , drop = FALSE]
        rownames(ret) <- dn[[2L]]
        return(ret)
    }
}

diagonals.syMatrices <- diagonals.ltMatrices

diagonals.matrix <- function(x, ...) diag(x)

# diagonal matrix

diagonals.integer <- function(x, ...)
    ltMatrices(rep(0, x * (x - 1) / 2), diag = FALSE, ...)

# mult ltMatrices

### C %*% y
Mult <- function(x, y, ...)
    UseMethod("Mult")
Mult.default <- function(x, y, transpose = FALSE, ...) {
    if (!transpose) return(x %*% y)
    return(crossprod(x, y))
}
Mult.ltMatrices <- function(x, y, transpose = FALSE, ...) {

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    stopifnot(is.numeric(y))
    if (!is.matrix(y)) y <- matrix(y, nrow = d[2L], ncol = d[1L])
    N <- ifelse(d[1L] == 1, ncol(y), d[1L])
    stopifnot(nrow(y) == d[2L])
    if (ncol(y) != N)
        return(sapply(1:ncol(y), function(i) Mult(x, y[,i], transpose = transpose)))

    # mult ltMatrices transpose
    
    if (transpose) {
        x <- ltMatrices(x, byrow = FALSE)
        if (!is.double(x)) storage.mode(x) <- "double"
        if (!is.double(y)) storage.mode(y) <- "double"

        ret <- .Call(mvtnorm_R_ltMatrices_Mult_transpose, x, y, as.integer(N), 
                     as.integer(d[2L]), as.logical(diag))

        rownames(ret) <- dn[[2L]]
        if (length(dn[[1L]]) == N)
            colnames(ret) <- dn[[1L]]
        return(ret)
    }
    

    x <- ltMatrices(x, byrow = TRUE)
    if (!is.double(x)) storage.mode(x) <- "double"
    if (!is.double(y)) storage.mode(y) <- "double"

    ret <- .Call(mvtnorm_R_ltMatrices_Mult, x, y, as.integer(N), 
                 as.integer(d[2L]), as.logical(diag))
    
    rownames(ret) <- dn[[2L]]
    if (length(dn[[1L]]) == N)
        colnames(ret) <- dn[[1L]]
    return(ret)
}

# mult syMatrices

Mult.syMatrices <- function(x, y, ...) {

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    x <- as.ltMatrices(x)
    stopifnot(is.numeric(y))
    if (!is.matrix(y)) y <- matrix(y, nrow = d[2L], ncol = d[1L])
    N <- ifelse(d[1L] == 1, ncol(y), d[1L])
    stopifnot(nrow(y) == d[2L])
    stopifnot(ncol(y) == N)

    ret <- Mult(x, y) + Mult(x, y, transpose = TRUE) - y * c(diagonals(x))
    return(ret)
}

# solve ltMatrices

solve.ltMatrices <- function(a, b, transpose = FALSE, ...) {

    byrow_orig <- attr(a, "byrow")

    x <- ltMatrices(a, byrow = FALSE)
    diag <- attr(x, "diag")
    ### dtptri and dtpsv require diagonal elements being present
    if (!diag) diagonals(x) <- diagonals(x)
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    if (!is.double(x)) storage.mode(x) <- "double"

    if (!missing(b)) {
        if (!is.matrix(b)) b <- matrix(b, nrow = J, ncol = d[1L])
        stopifnot(nrow(b) == J)
        N <- ifelse(d[1L] == 1, ncol(b), d[1L])
        stopifnot(ncol(b) == N)
        if (!is.double(b)) storage.mode(b) <- "double"
        ret <- .Call(mvtnorm_R_ltMatrices_solve, x, b, 
                     as.integer(N), as.integer(J), as.logical(diag),
                     as.logical(transpose))
        if (d[1L] == N) {
            colnames(ret) <- dn[[1L]]
        } else {
            colnames(ret) <- colnames(b)
        }
        rownames(ret) <- dn[[2L]]
        return(ret)
    }

    if (transpose) stop("cannot compute inverse of t(a)")
    ret <- .Call(mvtnorm_R_ltMatrices_solve_C, x, 
                 as.integer(d[1L]), as.integer(J), as.logical(diag),
                 as.logical(FALSE))
    ### chol = solve(invchol); invchol = solve(chol)
    if (is.chol(ret)) {
        class(ret) <- gsub("chol", "invchol", class(ret))
    } else if (is.invchol(ret)) {
        class(ret) <- gsub("invchol", "chol", class(ret))
    }

    if (!diag) {
        ### ret always includes diagonal elements, remove here
        ret <- Lower_tri(ret, diag = FALSE)
        ret <- ltMatrices(ret, diag = FALSE, byrow = FALSE, names = dn[[2L]])
    }

    ret <- ltMatrices(ret, byrow = byrow_orig)
    return(ret)
}

# logdet ltMatrices

logdet <- function(x) {

    if (!is.ltMatrices(x))
        stop("x is not an ltMatrices object")

    byrow <- attr(x, "byrow")
    diag <- attr(x, "diag")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_ltMatrices_logdet, x, 
                 as.integer(d[1L]), as.integer(J), as.logical(diag), 
                 as.logical(byrow))
    names(ret) <- dn[[1L]]
    return(ret) 
}

# tcrossprod ltMatrices

### C %*% t(C) => returns object of class syMatrices
### diag(C %*% t(C)) => returns matrix of diagonal elements
.Tcrossprod <- function(x, diag_only = FALSE, transpose = FALSE) {

    if (!is.ltMatrices(x)) {
        ret <- tcrossprod(x)
        if (diag_only) ret <- diag(ret)
        return(ret)
    }

    byrow_orig <- attr(x, "byrow")
    diag <- attr(x, "diag")
    d <- dim(x)
    N <- d[1L]
    J <- d[2L]
    dn <- dimnames(x)

    x <- ltMatrices(x, byrow = FALSE)
    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_ltMatrices_tcrossprod, x, as.integer(N), as.integer(J), 
                 as.logical(diag), as.logical(diag_only), as.logical(transpose))
    colnames(ret) <- dn[[1L]]
    if (diag_only) {
        rownames(ret) <- dn[[2L]]
    } else {
        ret <- ltMatrices(ret, diag = TRUE, byrow = FALSE, names = dn[[2L]])
        ret <- as.syMatrices(ltMatrices(ret, byrow = byrow_orig))
    }
    return(ret)
}
Tcrossprod <- function(x, diag_only = FALSE)
    .Tcrossprod(x = x, diag_only = diag_only, transpose = FALSE)

# crossprod ltMatrices

Crossprod <- function(x, diag_only = FALSE)
    .Tcrossprod(x, diag_only = diag_only, transpose = TRUE)

# crossprod tcrossprod methods

crossprod.ltMatrices <- function(x, y = NULL, ...) {

    if (is.null(y))
        return(Crossprod(x = x))
    return(Mult(x, y, transpose = TRUE))
}

crossprod.syMatrices <- crossprod.ltMatrices

tcrossprod.ltMatrices <- function(x, y = NULL, ...) {

    if (is.null(y))
        return(Tcrossprod(x = x))
    return(Mult(x, y, transpose = FALSE))
}

tcrossprod.syMatrices <- tcrossprod.ltMatrices

'%*%.ltMatrices' <- function(x, y)
    Mult(x, y)

'%*%.syMatrices' <- function(x, y)
    Mult(x, y)

# chol syMatrices

chol.syMatrices <- function(x, ...) {

    byrow_orig <- attr(x, "byrow")
    stopifnot(attr(x, "diag"))
    d <- dim(x)

    ### x is of class syMatrices, coerse to ltMatrices first and re-arrange
    ### second
    class(x) <- gsub("syMatrices", "ltMatrices", class(x))    
    x <- ltMatrices(x, byrow = FALSE)

    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_syMatrices_chol, x, 
                 as.integer(d[1L]), as.integer(d[2L]))

    ret <- ltMatrices(ret, byrow = byrow_orig)

    return(ret)
}

# invchol syMatrices

invchol <- function(x, ...)
    UseMethod("invchol")

invchol.default <- function(x, ...)
    invchol(as.syMatrices(as.matrix(x)))

invchol.syMatrices <- function(x, ...) {

    byrow_orig <- attr(x, "byrow")
    stopifnot(attr(x, "diag"))
    d <- dim(x)

    ### x is of class syMatrices, coerse to ltMatrices first and re-arrange
    ### second
    class(x) <- gsub("syMatrices", "ltMatrices", class(x))
    x <- ltMatrices(x, byrow = TRUE)

    if (!is.double(x)) storage.mode(x) <- "double"

    ret <- .Call(mvtnorm_R_syMatrices_invchol, x, 
                 as.integer(d[1L]), as.integer(d[2L]))

    ret <- ltMatrices(ret, byrow = byrow_orig)

    return(ret)
}

# add diagonal elements

.adddiag <- function(x) {

    ### only we call this function
    # stopifnot(is.ltMatrices(x))

    if (attr(x, "diag")) return(x)

    byrow_orig <- attr(x, "byrow")

    x <- ltMatrices(x, byrow = FALSE)

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    nm <- dimnames(x)[[2L]]

    L <- D <- diag(J)
    lt <- .lt(J, diag = TRUE)
    L[lt] <- seq_len(J * (J + 1) / 2)

    ### by default new diagonal elements are 1
    ret <- matrix(D[lt], nrow = J * (J + 1) / 2, ncol = N)
    colnames(ret) <- dimnames(x)[[1L]]
    diag(lt) <- FALSE
    ret[L[lt],] <- unclass(x)

    ret <- ltMatrices(ret, diag = TRUE, byrow = FALSE, names = nm)
    ret <- ltMatrices(ret, byrow = byrow_orig)

    ret
}

# assign diagonal elements

"diagonals<-" <- function(x, value)
    UseMethod("diagonals<-")

"diagonals<-.ltMatrices" <- function(x, value) {

    # extract slots
    
    diag <- attr(x, "diag")
    byrow <- attr(x, "byrow")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    

    if (byrow)
        idx <- cumsum(c(1, 2:J))
    else
        idx <- cumsum(c(1, J:2))

    ### diagonals(x) <- NULL returns ltMatrices(..., diag = FALSE)
    if (is.null(value)) {
        if (!attr(x, "diag")) return(x)
        if (J == 1L) {
            x[] <- 1
            return(x)
        }
        return(ltMatrices(unclass(x)[-idx,,drop = FALSE], diag = FALSE, 
                          byrow = byrow, names = dn[[2L]]))
    }

    x <- .adddiag(x)

    if (!is.matrix(value))
        value <- matrix(value, nrow = J, ncol = d[1L])

    stopifnot(is.matrix(value) && nrow(value) == J 
                               && ncol(value) == d[1L])

    if (J == 1L) {
        x[] <- value
        return(x)
    }

    x[idx, ] <- value

    return(x)
}

"diagonals<-.syMatrices" <- function(x, value) {

    x <- as.ltMatrices(x)
    diagonals(x) <- value
    class(x)[1L] <- "syMatrices"

    return(x)
}

# kronecker vec trick

vectrick <- function(C, S, A, transpose = c(TRUE, TRUE)) {

    stopifnot(all(is.logical(transpose)))
    stopifnot(length(transpose) == 2L)

    # check C argument
    
    C <- as.ltMatrices(C)
    if (!attr(C, "diag")) diagonals(C) <- 1
    C_byrow_orig <- attr(C, "byrow")
    C <- ltMatrices(C, byrow = FALSE)
    dC <- dim(C)
    nm <- attr(C, "rcnames")
    N <- dC[1L]
    J <- dC[2L]
    class(C) <- class(C)[-1L]   ### works because of as.ltMatrices(C)
    if (!is.double(C)) storage.mode(C) <- "double"
    
    # check S argument
    
    SltM <- is.ltMatrices(S)
    if (SltM) {
        if (!attr(S, "diag")) diagonals(S) <- 1
        S_byrow_orig <- attr(S, "byrow")
        stopifnot(S_byrow_orig == C_byrow_orig)
        S <- ltMatrices(S, byrow = FALSE)
        dS <- dim(S)
        stopifnot(dC[2L] == dS[2L])
        if (dC[1] != 1L) {
            stopifnot(dC[1L] == dS[1L])
        } else {
            N <- dS[1L]
        }
        ## argument A in dtrmm is not in packed form, so expand in J x J
        ## matrix
        S <- matrix(as.array(S), ncol = dS[1L])
    } else {
        stopifnot(is.matrix(S))
        stopifnot(nrow(S) == J^2)
        if (dC[1] != 1L) {
            stopifnot(dC[1L] == ncol(S))
        } else {
            N <- ncol(S)
        }
    }
    if (!is.double(S)) storage.mode(S) <- "double"
    
    # check A argument
    
    if (missing(A)) {
        A <- C
    } else {
        A <- as.ltMatrices(A)
        if (!attr(A, "diag")) diagonals(A) <- 1
        A_byrow_orig <- attr(A, "byrow")
        stopifnot(C_byrow_orig == A_byrow_orig)
        A <- ltMatrices(A, byrow = FALSE)
        dA <- dim(A)
        stopifnot(dC[2L] == dA[2L])
        class(A) <- class(A)[-1L]
        if (!is.double(A)) storage.mode(A) <- "double"
        if (dC[1L] != dA[1L]) {
            if (dC[1L] == 1L)
                C <- C[, rep(1, N), drop = FALSE]
            if (dA[1L] == 1L)
                A <- A[, rep(1, N), drop = FALSE]
            stopifnot(ncol(A) == ncol(C))
        }
    }
    

    ret <- .Call(mvtnorm_R_vectrick, C, as.integer(N), as.integer(J), S, A, 
                 as.logical(TRUE), as.logical(transpose))

    if (!SltM) return(matrix(c(ret), ncol = N))

    L <- matrix(1:(J^2), nrow = J)
    ret <- ltMatrices(ret[L[.lt(J, diag = TRUE)],,drop = FALSE], 
                      diag = TRUE, byrow = FALSE, names = nm)
    ret <- ltMatrices(ret, byrow = C_byrow_orig)
    return(ret)
}

# convenience functions

# chol classes

is.chol <- function(x) inherits(x, "chol")
as.chol <- function(x) {
    stopifnot(is.ltMatrices(x))
    if (is.chol(x)) return(x)
    if (is.invchol(x))
        return(invchol2chol(x))
    class(x) <- c("chol", class(x))
    return(x)
}
is.invchol <- function(x) inherits(x, "invchol")
as.invchol <- function(x) {
    stopifnot(is.ltMatrices(x))
    if (is.invchol(x)) return(x)
    if (is.chol(x))
        return(chol2invchol(x))
    class(x) <- c("invchol", class(x))
    return(x)
}

# D times C

Dchol <- function(x, D = 1 / sqrt(Tcrossprod(x, diag_only = TRUE))) {

    if (is.invchol(x)) stop("Dchol cannot work with invchol objects")

    x <- .adddiag(x)

    byrow_orig <- attr(x, "byrow")

    x <- ltMatrices(x, byrow = TRUE)

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    nm <- dimnames(x)[[2L]]

    ### for some parameter configurations logdet(ret) would
    ### be -Inf; make sure this does't happen
    if (any(D < .Machine$double.eps))
        D[D < .Machine$double.eps] <- 2 * .Machine$double.eps
    if (any(D > 1 / .Machine$double.eps))
        D[D > 1 / .Machine$double.eps] <-  (1 / .Machine$double.eps) / 2

    x <- unclass(x) * D[rep(1:J, 1:J),,drop = FALSE]

    ret <- ltMatrices(x, diag = TRUE, byrow = TRUE, names = nm)
    ret <- as.chol(ltMatrices(ret, byrow = byrow_orig))
    return(ret)
}

# L times D

### invcholD = solve(Dchol)
invcholD <- function(x, D = sqrt(Tcrossprod(solve(x), diag_only = TRUE))) {

    if (is.chol(x)) stop("invcholD cannot work with chol objects")

    x <- .adddiag(x)

    byrow_orig <- attr(x, "byrow")

    x <- ltMatrices(x, byrow = FALSE)

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    nm <- dimnames(x)[[2L]]

    ### for some parameter configurations logdet(ret) would
    ### be -Inf; make sure this does't happen
    if (any(D < .Machine$double.eps))
        D[D < .Machine$double.eps] <- 2 * .Machine$double.eps
    if (any(D > 1 / .Machine$double.eps))
        D[D > 1 / .Machine$double.eps] <-  (1 / .Machine$double.eps) / 2

    x <- unclass(x) * D[rep(1:J, J:1),,drop = FALSE]

    ret <- ltMatrices(x, diag = TRUE, byrow = FALSE, names = nm)
    ret <- as.invchol(ltMatrices(ret, byrow = byrow_orig))
    return(ret)
}


### C -> Sigma
chol2cov <- function(x)
    Tcrossprod(x)

### Sigma -> C
cov2chol <- function(x)
    as.chol(chol(x))

### L -> C
invchol2chol <- function(x)
    as.chol(solve(x))

### C -> L
chol2invchol <- function(x)
    as.invchol(solve(x))

### L -> Sigma
invchol2cov <- function(x)
    chol2cov(invchol2chol(x))

### Sigma -> L
cov2invchol <- function(x)
    as.invchol(invchol(x))

### L -> Precision
invchol2pre <- function(x)
    Crossprod(x)

### C -> Precision
chol2pre <- function(x)
    Crossprod(chol2invchol(x))

### C -> R
chol2cor <- function(x) {
    ret <- Tcrossprod(Dchol(x))
    diagonals(ret) <- NULL
    return(ret)
}

### L -> R
invchol2cor <- function(x) {
    ret <- chol2cor(invchol2chol(x))
    diagonals(ret) <- NULL
    return(ret)
}

### L -> A
invchol2pc <- function(x) {
    ret <- -Crossprod(invcholD(x, D = 1 / sqrt(Crossprod(x, diag_only = TRUE))))
    diagonals(ret) <- 0
    ret
}

### C -> A
chol2pc <- function(x)
    invchol2pc(solve(x))

# aperm

aperm.chol <- function(a, perm, ...) {

    # aperm checks
    
    J <- dim(a)[2L]
    if (missing(perm)) return(a)
    if (is.character(perm)) 
        perm <- match(perm, dimnames(a)[[2L]])
    stopifnot(all(perm %in% 1:J))

    args <- list(...)
    if (length(args) > 0L)
        warning("Additional arguments", names(args), "ignored")
    

    return(as.chol(chol(chol2cov(a)[,perm])))
}
aperm.invchol <- function(a, perm, ...) {

    # aperm checks
    
    J <- dim(a)[2L]
    if (missing(perm)) return(a)
    if (is.character(perm)) 
        perm <- match(perm, dimnames(a)[[2L]])
    stopifnot(all(perm %in% 1:J))

    args <- list(...)
    if (length(args) > 0L)
        warning("Additional arguments", names(args), "ignored")
    

    return(as.invchol(invchol(invchol2cov(a)[,perm])))
}

aperm.ltMatrices <- function(a, perm, ...)
    stop("Cannot permute objects of class ltMatrices, 
          consider calling as.chol() or as.invchol() first")

aperm.syMatrices <- function(a, perm, ...)
    return(a[,perm])

# marginal

marg_mvnorm <- function(chol, invchol, which = 1L) {

    # mc input checks
    
    stopifnot(xor(missing(chol), missing(invchol)))
    x <- if (missing(chol)) invchol else chol

    stopifnot(is.ltMatrices(x))

    N <- dim(x)[1L]
    J <- dim(x)[2L]

    if (missing(which)) return(x)

    if (is.character(which)) which <- match(which, dimnames(x)[[2L]])
    stopifnot(all(which %in% 1:J))
    

    if (which[1] == 1L && (length(which) == 1L || 
                           all(diff(which) == 1L))) {
        ### which is 1:j
        tmp <- x[,which]
    } else {
        if (missing(chol)) { 
            cv <- invchol2cov(x)
            tmp <- cov2invchol(cv[,which])
        } else {
            cv <- chol2cov(x)
            tmp <- cov2chol(cv[,which])
        }
    }

    if (missing(chol))
        ret <- list(invchol = tmp)
    else
        ret <- list(chol = tmp)

    ret
}

# conditional

cond_mvnorm <- function(chol, invchol, which_given = 1L, given, center = FALSE) {

    which <- which_given
    # mc input checks
    
    stopifnot(xor(missing(chol), missing(invchol)))
    x <- if (missing(chol)) invchol else chol

    stopifnot(is.ltMatrices(x))

    N <- dim(x)[1L]
    J <- dim(x)[2L]

    if (missing(which)) return(x)

    if (is.character(which)) which <- match(which, dimnames(x)[[2L]])
    stopifnot(all(which %in% 1:J))
    

    if (N == 1) N <- NCOL(given)
    stopifnot(is.matrix(given) && nrow(given) == length(which))

    # cond simple
    
    if (which[1] == 1L && (length(which) == 1L || 
                           all(diff(which) == 1L))) {
        ### which is 1:j
        L <- if (missing(invchol)) solve(chol) else invchol
        tmp <- matrix(0, ncol = ncol(given), nrow = J - length(which))
        centerm <- Mult(L, rbind(given, tmp))               
        ### if ncol(given) is not N = dim(L)[1L] > 1, then
        ### solve() below won't work and we loop over
        ### columns of centerm
        if (dim(L)[1L] > 1 && ncol(given) != N) {
            centerm <- lapply(1:ncol(centerm), function(j)
                matrix(centerm[,j], nrow = J, ncol = N)[-which,,drop = FALSE]
            )
        } else {
            centerm <- centerm[-which,,drop = FALSE]
        }
        L <- L[,-which]
        ct <- centerm
        if (!is.matrix(ct)) ct <- do.call("rbind", ct)
        if (is.matrix(centerm)) {
            m <- -solve(L, centerm)
        } else {
            m <- do.call("rbind", lapply(centerm, function(cm) -solve(L, cm)))
        }
        if (missing(invchol)) {
            if (center)
                return(list(center = ct, chol = solve(L)))
            return(list(mean = m, chol = solve(L)))
        }
        if (center)
            return(list(center = ct, invchol = L))
        return(list(mean = m, invchol = L))
    }
    

    ### general with center = TRUE => permute first and go simple
    if (center) {
        perm <- c(which, (1:J)[!(1:J) %in% which])
        if (!missing(chol))
        return(cond_mvnorm(chol = aperm(as.chol(chol), perm = perm),
                           which_given = 1:length(which), given = given,
                           center = center))
        return(cond_mvnorm(invchol = aperm(as.invchol(invchol), perm = perm),
                           which_given = 1:length(which), given = given,
                           center = center))
    }

    # cond general
    
    stopifnot(!center)

    if (!missing(chol)) ### chol is C = Cholesky of covariance
        P <- Crossprod(solve(chol)) ### P = t(L) %*% L with L = C^-1
    else                ### invcol is L = Cholesky of precision
        P <- Crossprod(invchol)

    Pw <- P[, -which]
    # chol <- solve(A <- base::chol(Pw)) ### Pw = A A^\top
    ### arg invchol is missing, masking mvtnorm::invchol
    ### which we can not access bc R CMD check is not happy about it
    chol <- getFromNamespace("invchol", "mvtnorm")(Pw) ### Pw = A A^\top
    g0 <- matrix(0, nrow = J, ncol = NCOL(given))
    g0[which,] <- given
    S <- Crossprod(chol) ### P^{-1}_jj = A^-top A^-1
    if (ncol(g0) %in% c(1L, N)) {
        mean <- - S %*% (P %*% g0)[-which,,drop = FALSE]
    } else { ### compare to Mult() with ncol(y) !%in% (1, N)
        ### <FIXME> check dimension and if t() is needed
        mean <- t(sapply(seq_len(ncol(g0)), function(i)
            - S %*% (P %*% g0[,i])[-which,,drop = FALSE]))
    }
    

    if (missing(invchol)) 
        return(list(mean = mean, chol = base::chol(S)))

    return(list(mean = mean, invchol = invchol(S)))
}

# check obs

.check_obs_mean <- function(obs, mean, J, N) {

    nr <- nrow(obs)
    nc <- ncol(obs)
    if (nc != N)
        stop("obs and (inv)chol have non-conforming size")
    if (nr != J)
        stop("obs and (inv)chol have non-conforming size")
    if (identical(unique(mean), 0)) return(obs)
    if (length(mean) == J) 
        return(obs - c(mean))
    if (!is.matrix(mean))
        stop("obs and mean have non-conforming size")
    if (nrow(mean) != nr)
        stop("obs and mean have non-conforming size")
    if (ncol(mean) != nc)
        stop("obs and mean have non-conforming size")
    return(obs - mean)
}
.check_obs_invcholmean <- function(obs, invcholmean, J, N) {

    nr <- nrow(obs)
    nc <- ncol(obs)
    if (nc != N)
        stop("obs and (inv)chol have non-conforming size")
    if (nr != J)
        stop("obs and (inv)chol have non-conforming size")

    if (is.null(invcholmean)) 
        return(matrix(0, nrow = J, ncol = N))
    if (identical(unique(invcholmean), 0)) 
        return(matrix(0, nrow = J, ncol = N))

    if (!is.matrix(invcholmean))
        invcholmean <- matrix(invcholmean, nrow = J)
    nr <- nrow(invcholmean)
    nc <- ncol(invcholmean)
    if (!(nc %in% c(1L, N)))
        stop("obs and (inv)chol have non-conforming size")
    if (nr != J)
        stop("invcholmean and (inv)chol have non-conforming size")

    if (ncol(invcholmean) == N) 
        return(invcholmean)
    return(matrix(invcholmean, nrow = J, ncol = N))
}

# colSumsdnorm ltMatrices

.colSumsdnorm <- function(z) {
    stopifnot(is.numeric(z))
    if (!is.matrix(z))
        z <- matrix(z, nrow = 1, ncol = length(z))
    ret <- .Call(mvtnorm_R_ltMatrices_colSumsdnorm, z, ncol(z), nrow(z))
    names(ret) <- colnames(z)
    return(ret)
}