Update indMatrix->arma::SpMat and pMatrix->arma::SpMat for Matrix 1.5-5

[jaganmn]

In Matrix 1.5-5, class indMatrix will gain a slot margin with value 1L or 2L. If x is an indMatrix and x@margin = 1L (the "prototype"), then the remaining slots of x are interpreted as in older versions of Matrix: x has dimensions x@Dim with names x@Dimnames and is zero except for ones in positions cbind(i = seq_along(x@perm), j = x@perm). If x@margin = 2L, then the ones occur in positions cbind(i = x@perm, j = seq_along(x@perm)) instead.

In general, we try to avoid making such changes to class definitions, as they can invalidate serialized objects. But this particular change (a very useful generalization) is long overdue and, perhaps amazingly, RcppArmadillo is the only broken reverse dependency.

AFAICT, it should suffice to update the indMatrix and pMatrix branches in inst/include/RcppArmadillo/interface/RcppArmadilloAs.h to behave conditional on margin, checking first for the absence of a margin attribute (implying that an old Matrix is installed) and in that case behaving as if margin = 1L.

We have not discussed a release date for 1.5-5, so this is just an FYI (for now). In the mean time, I can advise on a patch if necessary.

Note that R-Forge is currently hosting an out-of-date Matrix tarball. Until that is fixed, you can use the R-universe tarball for testing (or install from sources if you have Subversion).

[eddelbuettel]

A patch, even a sketch, would be great!

[jaganmn]

For a sketch (written in R), we can modify the coercion method selectMethod("coerce", c("indMatrix", "CsparseMatrix")).

For indMatrix->dgCMatrix:

function (from) {
    perm <- from@perm
    margin <- attr(from, "margin")
    if (is.null(margin))
        margin <- 1L
    J <- new("dgCMatrix")
    J@Dim <- d <- from@Dim
    J@Dimnames <- from@Dimnames
    if (margin == 1L) {
        J@x <- rep.int(1, d[1L])
        J@p <- c(0L, cumsum(tabulate(perm, d[2L])))
        J@i <- sort.list(perm) - 1L
    } else {
        J@x <- rep.int(1, d[2L])
        J@p <- 0:d[2L]
        J@i <- perm - 1L
    }
    J
}

For pMatrix->dgCMatrix, a special case of the above:

function (from) {
    perm <- from@perm
    margin <- attr(from, "margin")
    if (is.null(margin))
        margin <- 1L
    J <- new("dgCMatrix")
    J@Dim <- d <- from@Dim
    J@Dimnames <- from@Dimnames
    J@x <- rep.int(1, d[2L])
    J@p <- 0:d[2L]
    if (margin == 1L)
        J@i <- sort.list(perm) - 1L
    else
        J@i <- perm - 1L
    J
}