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'

# 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(1:J)
    }

    if (!nonames) {
        L1 <- matrix(names, nrow = J, ncol = J)
        L2 <- matrix(names, nrow = J, ncol = J, byrow = TRUE)
        L <- matrix(paste(L1, L2, sep = "."), nrow = J, ncol = J)
        if (byrow)
            rownames(object) <- t(L)[upper.tri(L, diag = diag)]
        else
            rownames(object) <- L[lower.tri(L, diag = diag)]
    } # 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)
    x <- as.ltMatrices(x)       ### make sure "ltMatrices"
                                ### is first class
    class(x)[1L] <- "syMatrices"
    return(x)
}
syMatrices <- function(object, diag = FALSE, byrow = FALSE, names = TRUE)
    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[lower.tri(x, diag = DIAG)]
    up <- x[upper.tri(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[upper.tri(L, diag = diag)] <- floor(2L + 1:(J * (J - 1) / 2L + diag * J))
        L <- t(L)
    } else {
        L[lower.tri(L, diag = diag)] <- floor(2L + 1:(J * (J - 1) / 2L + diag * J))
    }
    if (symmetric) {
        L[upper.tri(L)] <- 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, ...)
    print(as.array(x))

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

# reorder ltMatrices

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

    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)

    rL <- cL <- diag(0, nrow = J)
    rL[lower.tri(rL, diag = diag)] <- cL[upper.tri(cL, diag = diag)] <- 1:nrow(x)
    cL <- t(cL)
    if (byrow) ### row -> col order
        return(ltMatrices(x[cL[lower.tri(cL, diag = diag)], , drop = FALSE], 
                          diag = diag, byrow = FALSE, names = dn[[2L]]))
    ### col -> row order
    return(ltMatrices(x[t(rL)[upper.tri(rL, 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(upper.tri(L, diag = diag))
        if (byrow) {
            L[upper.tri(L, diag = diag)] <- 1:Jp
            L <- L + t(L)
            diag(L) <- diag(L) / 2
            L <- L[j, j, drop = FALSE]
            L <- L[upper.tri(L, diag = diag)]
        } else {
            L[lower.tri(L, diag = diag)] <- 1:Jp
            L <- L + t(L)
            diag(L) <- diag(L) / 2
            L <- L[j, j, drop = FALSE]
            L <- L[lower.tri(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))
    colnames(ret) <- dn[[1L]]

    if (!diag)
        ### ret always includes diagonal elements, remove here
        ret <- ret[- cumsum(c(1, J:2)), , drop = FALSE]

    ret <- ltMatrices(ret, diag = diag, 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)

# chol syMatrices

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

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

    ### x is of class syMatrices, coerse to ltMatrices first and re-arrange
    ### second
    x <- ltMatrices(unclass(x), diag = TRUE, 
                    byrow = byrow_orig, names = dnm[[2L]])
    x <- ltMatrices(x, byrow = FALSE)
    # class(x) <- class(x)[-1]
    if (!is.double(x)) storage.mode(x) <- "double"

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

    ret <- ltMatrices(ret, diag = TRUE,
                      byrow = FALSE, names = dnm[[2L]])
    ret <- ltMatrices(ret, byrow = byrow_orig)

    return(ret)
}

# add diagonal elements

.adddiag <- function(x) {

    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 <- diag(J)
    L[lower.tri(L, diag = TRUE)] <- 1:(J * (J + 1) / 2)

    D <- diag(J)
    ret <- matrix(D[lower.tri(D, diag = TRUE)], 
                  nrow = J * (J + 1) / 2, ncol = N)
    colnames(ret) <- dimnames(x)[[1L]]
    ret[L[lower.tri(L, diag = FALSE)],] <- 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[lower.tri(L, 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

    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

    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)

### 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))

### 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(chol2invchol(chol(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)) x <- invchol2chol(x)
        ### note: aperm would work but computes
        ### Cholesky of J^2, here only length(which)^2
        ### is needed
        tmp <- base::chol(chol2cov(x)[,which])
        if (missing(chol)) tmp <- chol2invchol(tmp)
    }

    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(base::chol(Pw))
    Pa <- as.array(P)
    Sa <- as.array(S <- Crossprod(chol))
    if (dim(chol)[1L] == 1L) {
       Pa <- Pa[,,1]
       Sa <- Sa[,,1]
       mean <- -Sa %*% Pa[-which, which, drop = FALSE] %*% given
    } else {
       if (ncol(given) == N) {
           mean <- sapply(1:N, function(i) 
               -Sa[,,i] %*% Pa[-which,which,i] %*% given[,i,drop = FALSE])
       } else {  ### compare to Mult() with ncol(y) !%in% (1, N)
           mean <- sapply(1:N, function(i) 
               -Sa[,,i] %*% Pa[-which,which,i] %*% given)
       }
    }
    

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

    return(list(mean = mean, invchol = solve(chol)))
}

# check obs

.check_obs <- 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)
}

# 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)
}