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 (inherits(object, "ltMatrices")) {
        ret <- .reorder(object, byrow = byrow)
        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)]
    }
    

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

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

# dim ltMatrices

dim.ltMatrices <- function(x) {
    J <- attr(x, "J")
    class(x) <- class(x)[-1L]
    return(c(ncol(x), J, J))
}
dim.syMatrices <- dim.ltMatrices

# dimnames ltMatrices

dimnames.ltMatrices <- function(x)
    return(list(colnames(unclass(x)), attr(x, "rcnames"), attr(x, "rcnames")))
dimnames.syMatrices <- dimnames.ltMatrices

# names ltMatrices

names.ltMatrices <- function(x) {
    return(rownames(unclass(x)))
}
names.syMatrices <- names.ltMatrices

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

    class(x) <- class(x)[-1L]

    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(inherits(x, "ltMatrices"))
    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)
    

    class(x) <- class(x)[-1L]

    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)
    

    class(x) <- class(x)[-1L]

    if (!missing(j)) {

        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 (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) {
    class(x)[1L] <- "ltMatrices"
    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 (inherits(x, "syMatrices"))
        class(x)[1L] <- "ltMatrices"
    stopifnot(inherits(x, "ltMatrices"))
    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))

    if (!diag && adiag) {
        diagonals(x) <- 1
        return(unclass(x))
    }

    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)
    

    class(x) <- class(x)[-1L]

    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, transpose = FALSE) {

    if (!inherits(x, "ltMatrices")) {
        if (!transpose) return(x %*% y)
        return(crossprod(x, y))
    }

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

    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) {
        J <- dim(x)[2L]
        if (dim(x)[1L] == 1L) x <- x[rep(1, N),]
        ax <- as.array(x)
        ay <- array(y[rep(1:J, J),,drop = FALSE], dim = dim(ax), 
                    dimnames = dimnames(ax))
        ret <- ay * ax
        ### was: return(margin.table(ret, 2:3))
        ret <- matrix(colSums(matrix(ret, nrow = dim(ret)[1L])), 
                      nrow = dim(ret)[2L], ncol = dim(ret)[3L],
                      dimnames = dimnames(ret)[-1L])
        return(ret)
    }
    

    x <- ltMatrices(x, byrow = TRUE)

    class(x) <- class(x)[-1L]
    storage.mode(x) <- "double"
    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)
}

# solve ltMatrices

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

    byrow_orig <- attr(a, "byrow")

    x <- ltMatrices(a, byrow = FALSE)
    diag <- attr(x, "diag")
    d <- dim(x)
    J <- d[2L]
    dn <- dimnames(x)
    class(x) <- class(x)[-1L]
    storage.mode(x) <- "double"

    if (!missing(b)) {
        if (!is.matrix(b)) b <- matrix(b, nrow = J, ncol = ncol(x))
        stopifnot(nrow(b) == J)
        N <- ifelse(d[1L] == 1, ncol(b), d[1L])
        stopifnot(ncol(b) == N)
        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 <- try(.Call(mvtnorm_R_ltMatrices_solve, x, NULL,
                     as.integer(ncol(x)), 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)
}

# 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 (!inherits(x, "ltMatrices")) {
        ret <- tcrossprod(x)
        if (diag_only) ret <- diag(ret)
        return(ret)
    }

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

    x <- ltMatrices(x, byrow = FALSE)
    class(x) <- class(x)[-1L]
    N <- d[1L]
    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 <- ltMatrices(ret, byrow = byrow_orig)
        class(ret)[1L] <- "syMatrices"
    }
    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]
    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(inherits(x, "ltMatrices")) 

    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) <- colnames(unclass(x))
    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) {

    class(x)[1L] <- "ltMatrices"
    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
    
    stopifnot(inherits(C, "ltMatrices"))
    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]
    storage.mode(C) <- "double"
    
    # check S argument
    
    SltM <- inherits(S, "ltMatrices")
    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)
        }
    }
    storage.mode(S) <- "double"
    
    # check A argument
    
    if (missing(A)) {
        A <- C
    } else {
        stopifnot(inherits(A, "ltMatrices"))
        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]
        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

# D times C

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

    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]]

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

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

# L times D

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

    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]]

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

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


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

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

### C -> L
chol2invchol <- function(x)
    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.ltMatrices <- function(a, perm, is_chol = FALSE, ...) {

    if (is_chol) { ### a is Cholesky of covariance
        Sperm <- chol2cov(a)[,perm]
        return(chol(Sperm))
    }

    Sperm <- invchol2cov(a)[,perm]
    chol2invchol(chol(Sperm))
}

# 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(inherits(x, "ltMatrices"))

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    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 <- solve(x)
        tmp <- base::chol(Tcrossprod(x)[,which])
        if (missing(chol)) tmp <- solve(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(inherits(x, "ltMatrices"))

    N <- dim(x)[1L]
    J <- dim(x)[2L]
    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))[-which,,drop = FALSE]
        L <- L[,-which]
        if (missing(invchol)) {
            if (center)
                return(list(center = centerm, chol = solve(L)))
            return(list(mean = -solve(L, centerm), chol = solve(L)))
        }
        if (center)
            return(list(center = centerm, invchol = L))
        return(list(mean = -solve(L, centerm), invchol = L))
    }
    
    # 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)
}

# ldmvnorm

ldmvnorm <- function(obs, mean = 0, chol, invchol, logLik = TRUE) {

    stopifnot(xor(missing(chol), missing(invchol)))
    if (!is.matrix(obs)) obs <- matrix(obs, ncol = 1L)
    p <- ncol(obs)

    if (!missing(chol)) {
         # ldmvnorm chol
         
         if (missing(chol))
             stop("either chol or invchol must be given")
         ## chol is given
         if (!inherits(chol, "ltMatrices"))
             stop("chol is not an object of class ltMatrices")
         N <- dim(chol)[1L]
         N <- ifelse(N == 1, p, N)
         J <- dim(chol)[2L]
         obs <- .check_obs(obs = obs, mean = mean, J = J, N = N)
         logretval <- colSums(dnorm(solve(chol, obs), log = TRUE))
         if (attr(chol, "diag"))
             logretval <- logretval - colSums(log(diagonals(chol)))
         
    } else {
         # ldmvnorm invchol
         
         ## invchol is given
         if (!inherits(invchol, "ltMatrices"))
             stop("invchol is not an object of class ltMatrices")
         N <- dim(invchol)[1L]
         N <- ifelse(N == 1, p, N)
         J <- dim(invchol)[2L]
         obs <- .check_obs(obs = obs, mean = mean, J = J, N = N)
         ## use dnorm (gets the normalizing factors right)
         ## NOTE: obs is (J x N) 
         logretval <- colSums(dnorm(Mult(invchol, obs), log = TRUE))
         ## note that the second summand gets recycled the correct number
         ## of times in case dim(invchol)[1L] == 1 but ncol(obs) > 1
         if (attr(invchol, "diag"))
             logretval <- logretval + colSums(log(diagonals(invchol)))
         
    }

    names(logretval) <- colnames(obs)
    if (logLik) return(sum(logretval))
    return(logretval)
}

# sldmvnorm

sldmvnorm <- function(obs, mean = 0, chol, invchol, logLik = TRUE) {

    stopifnot(xor(missing(chol), missing(invchol)))
    if (!is.matrix(obs)) obs <- matrix(obs, ncol = 1L)

    if (!missing(invchol)) {

        N <- dim(invchol)[1L]
        N <- ifelse(N == 1, ncol(obs), N)
        J <- dim(invchol)[2L]
        obs <- .check_obs(obs = obs, mean = mean, J = J, N = N)

        Mix <- Mult(invchol, obs)
        sobs <- - Mult(invchol, Mix, transpose = TRUE)

        Y <- matrix(obs, byrow = TRUE, nrow = J, ncol = N * J)
        ret <- - matrix(Mix[, rep(1:N, each = J)] * Y, ncol = N)

        M <- matrix(1:(J^2), nrow = J, byrow = FALSE)
        ret <- ltMatrices(ret[M[lower.tri(M, diag = attr(invchol, "diag"))],,drop = FALSE], 
                          diag = attr(invchol, "diag"), byrow = FALSE)
        ret <- ltMatrices(ret, 
                          diag = attr(invchol, "diag"), byrow = attr(invchol, "byrow"))
        if (attr(invchol, "diag")) {
            ### recycle properly
            diagonals(ret) <- diagonals(ret) + c(1 / diagonals(invchol))
        } else {
            diagonals(ret) <- 0
        }
        ret <- list(obs = sobs, invchol = ret)
        if (logLik) 
            ret$logLik <- ldmvnorm(obs = obs, mean = mean, invchol = invchol, logLik = FALSE)
        return(ret)
    }

    invchol <- solve(chol)
    ret <- sldmvnorm(obs = obs, mean = mean, invchol = invchol)
    ### this means: ret$chol <- - vectrick(invchol, ret$invchol, invchol)
    ret$chol <- - vectrick(invchol, ret$invchol)
    ret$invchol <- NULL
    return(ret)
}

# ldpmvnorm

ldpmvnorm <- function(obs, lower, upper, mean = 0, chol, invchol, 
                      logLik = TRUE, ...) {

    if (missing(obs) || is.null(obs))
        return(lpmvnorm(lower = lower, upper = upper, mean = mean,
                        chol = chol, invchol = invchol, logLik = logLik, ...))
    if (missing(lower) && missing(upper) || is.null(lower) && is.null(upper))
        return(ldmvnorm(obs = obs, mean = mean,
                        chol = chol, invchol = invchol, logLik = logLik))

    # dp input checks
    
    stopifnot(xor(missing(chol), missing(invchol)))
    cJ <- nrow(obs)
    dJ <- nrow(lower)
    N <- ncol(obs)
    stopifnot(N == ncol(lower))
    stopifnot(N == ncol(upper))
    if (all(mean == 0)) {
        cmean <- 0
        dmean <- 0
    } else {
        if (!is.matrix(mean)) 
            mean <- matrix(mean, nrow = cJ + dJ, ncol = N)
        stopifnot(nrow(mean) == cJ + dJ)
        stopifnot(ncol(mean) == N)
        cmean <- mean[1:cJ,, drop = FALSE]
        dmean <- mean[-(1:cJ),, drop = FALSE]
    }
        

    if (!missing(invchol)) {
        J <- dim(invchol)[2L]
        stopifnot(cJ + dJ == J)

        md <- marg_mvnorm(invchol = invchol, which = 1:cJ)
        ret <- ldmvnorm(obs = obs, mean = cmean, invchol = md$invchol, 
                        logLik = logLik)

        cd <- cond_mvnorm(invchol = invchol, which_given = 1:cJ, 
                          given = obs - cmean, center = TRUE)
        ret <- ret + lpmvnorm(lower = lower, upper = upper, mean = dmean, 
                              invchol = cd$invchol, center = cd$center, 
                              logLik = logLik, ...)
        return(ret)
    }

    J <- dim(chol)[2L]
    stopifnot(cJ + dJ == J)

    md <- marg_mvnorm(chol = chol, which = 1:cJ)
    ret <- ldmvnorm(obs = obs, mean = cmean, chol = md$chol, logLik = logLik)

    cd <- cond_mvnorm(chol = chol, which_given = 1:cJ, 
                      given = obs - cmean, center = TRUE)
    ret <- ret + lpmvnorm(lower = lower, upper = upper, mean = dmean, 
                          chol = cd$chol, center = cd$center, 
                          logLik = logLik, ...)
    return(ret)
}

# sldpmvnorm

sldpmvnorm <- function(obs, lower, upper, mean = 0, chol, invchol, logLik = TRUE, ...) {

    if (missing(obs) || is.null(obs))
        return(slpmvnorm(lower = lower, upper = upper, mean = mean,
                         chol = chol, invchol = invchol, logLik = logLik, ...))
    if (missing(lower) && missing(upper) || is.null(lower) && is.null(upper))
        return(sldmvnorm(obs = obs, mean = mean,
                         chol = chol, invchol = invchol, logLik = logLik))

    # dp input checks
    
    stopifnot(xor(missing(chol), missing(invchol)))
    cJ <- nrow(obs)
    dJ <- nrow(lower)
    N <- ncol(obs)
    stopifnot(N == ncol(lower))
    stopifnot(N == ncol(upper))
    if (all(mean == 0)) {
        cmean <- 0
        dmean <- 0
    } else {
        if (!is.matrix(mean)) 
            mean <- matrix(mean, nrow = cJ + dJ, ncol = N)
        stopifnot(nrow(mean) == cJ + dJ)
        stopifnot(ncol(mean) == N)
        cmean <- mean[1:cJ,, drop = FALSE]
        dmean <- mean[-(1:cJ),, drop = FALSE]
    }
        

    if (!missing(invchol)) {
        # sldpmvnorm invchol
        
        byrow_orig <- attr(invchol, "byrow")
        invchol <- ltMatrices(invchol, byrow = TRUE)

        J <- dim(invchol)[2L]
        stopifnot(cJ + dJ == J)

        md <- marg_mvnorm(invchol = invchol, which = 1:cJ)
        cs <- sldmvnorm(obs = obs, mean = cmean, invchol = md$invchol)

        obs_cmean <- obs - cmean
        cd <- cond_mvnorm(invchol = invchol, which_given = 1:cJ, 
                          given = obs_cmean, center = TRUE)
        ds <- slpmvnorm(lower = lower, upper = upper, mean = dmean, 
                        center = cd$center, invchol = cd$invchol, 
                        logLik = logLik, ...)

        tmp0 <- solve(cd$invchol, ds$mean, transpose = TRUE)
        tmp <- - tmp0[rep(1:dJ, each = cJ),,drop = FALSE] * 
                 obs_cmean[rep(1:cJ, dJ),,drop = FALSE]

        Jp <- nrow(unclass(invchol))
        diag <- attr(invchol, "diag")
        M <- as.array(ltMatrices(1:Jp, diag = diag, byrow = TRUE))[,,1]
        ret <- matrix(0, nrow = Jp, ncol = ncol(obs))
        M1 <- M[1:cJ, 1:cJ]
        idx <- t(M1)[upper.tri(M1, diag = diag)]
        ret[idx,] <- Lower_tri(cs$invchol, diag = diag)

        idx <- c(t(M[-(1:cJ), 1:cJ]))
        ret[idx,] <- tmp

        M3 <- M[-(1:cJ), -(1:cJ)]
        idx <- t(M3)[upper.tri(M3, diag = diag)]
        ret[idx,] <- Lower_tri(ds$invchol, diag = diag)

        ret <- ltMatrices(ret, diag = diag, byrow = TRUE)
        if (!diag) diagonals(ret) <- 0
        ret <- ltMatrices(ret, byrow = byrow_orig)

        ### post differentiate mean 
        aL <- as.array(invchol)[-(1:cJ), 1:cJ,,drop = FALSE]
        lst <- tmp0[rep(1:dJ, cJ),,drop = FALSE]
        if (dim(aL)[3] == 1)
              aL <- aL[,,rep(1, ncol(lst)), drop = FALSE]
        dim <- dim(aL)
        dobs <- -margin.table(aL * array(lst, dim = dim), 2:3)

        ret <- c(list(invchol = ret, obs = cs$obs + dobs), 
                 ds[c("lower", "upper")])
        ret$mean <- rbind(-ret$obs, ds$mean)
        return(ret)
        
    }

    invchol <- solve(chol)
    ret <- sldpmvnorm(obs = obs, lower = lower, upper = upper, 
                      mean = mean, invchol = invchol, logLik = logLik, ...)
    ### this means: ret$chol <- - vectrick(invchol, ret$invchol, invchol)
    ret$chol <- - vectrick(invchol, ret$invchol)
    ret$invchol <- NULL
    return(ret)
}

# standardize

standardize <- function(chol, invchol) {
    stopifnot(xor(missing(chol), missing(invchol)))
    if (!missing(invchol)) {
        stopifnot(!attr(invchol, "diag"))
        return(invcholD(invchol))
    }
    stopifnot(!attr(chol, "diag"))
    return(Dchol(chol))
}

# destandardize

destandardize <- function(chol = solve(invchol), invchol, score_schol)
{
    stopifnot(inherits(chol, "ltMatrices"))
    J <- dim(chol)[2L]
    stopifnot(!attr(chol, "diag"))
    byrow_orig <- attr(chol, "byrow")
    chol <- ltMatrices(chol, byrow = FALSE)
    
    if (inherits(score_schol, "ltMatrices"))
        score_schol <- matrix(as.array(score_schol), 
                              nrow = dim(score_schol)[2L]^2)
    stopifnot(is.matrix(score_schol))
    N <- ncol(score_schol)
    stopifnot(J^2 == nrow(score_schol))

    CCt <- Tcrossprod(chol, diag_only = TRUE)
    DC <- Dchol(chol, D = Dinv <- 1 / sqrt(CCt))
    SDC <- solve(DC)

    IDX <- t(M <- matrix(1:J^2, nrow = J, ncol = J))
    i <- cumsum(c(1, rep(J + 1, J - 1)))
    ID <- diagonals(as.integer(J), byrow = FALSE)
    if (dim(ID)[1L] != dim(chol)[1L])
        ID <- ID[rep(1, dim(chol)[1L]),]

    B <- vectrick(ID, score_schol, chol)
    B[i,] <- B[i,] * (-.5) * c(CCt)^(-3/2)
    B[-i,] <- 0

    Dtmp <- Dchol(ID, D = Dinv)

    ret <- vectrick(ID, B, chol, transpose = c(TRUE, FALSE)) +
           vectrick(chol, B, ID)[IDX,] +
           vectrick(Dtmp, score_schol, ID)

    if (!missing(invchol)) {
        ### this means: ret <- - vectrick(chol, ret, chol)
        ret <- - vectrick(chol, ret)
    }
    ret <- ltMatrices(ret[M[lower.tri(M)],,drop = FALSE],
                      diag = FALSE, byrow = FALSE)
    ret <- ltMatrices(ret, byrow = byrow_orig)
    diagonals(ret) <- 0
    return(ret)
}