The R Project SVN R-packages

#### Methods for the virtual class 'CsparseMatrix' of sparse matrices stored in
####  "column compressed" format.
#### -- many more specific things are e.g. in ./dgCMatrix.R

setAs("CsparseMatrix", "TsparseMatrix",
      function(from)
          ## |-> cholmod_C -> cholmod_T -> chm_triplet_to_SEXP
          ## modified to support triangular (../src/Csparse.c)
          .Call(Csparse_to_Tsparse, from, is(from, "triangularMatrix")))

## special cases (when a specific "to" class is specified)
setAs("dgCMatrix", "dgTMatrix",
      function(from) .Call(Csparse_to_Tsparse, from, FALSE))

setAs("dsCMatrix", "dsTMatrix",
      function(from) .Call(Csparse_to_Tsparse, from, FALSE))

setAs("dsCMatrix", "dgCMatrix",
      function(from) .Call(Csparse_symmetric_to_general, from))

setAs("dtCMatrix", "dtTMatrix",
      function(from) .Call(Csparse_to_Tsparse, from, TRUE))

## Current code loses symmetry and triangularity properties.  With suitable
## changes to chm_dense_to_SEXP (../src/chm_common.c) we can avoid this.
setAs("CsparseMatrix", "denseMatrix",
      function(from) {
          ## |-> cholmod_C -> cholmod_dense -> chm_dense_to_dense
          if (is(from, "triangularMatrix") && from@diag == "U")
              from <- .Call(Csparse_diagU2N, from)
          .Call(Csparse_to_dense, from)
      })

## special cases (when a specific "to" class is specified)
setAs("dgCMatrix", "dgeMatrix",
      function(from) .Call(Csparse_to_dense, from))

## cholmod_sparse_to_dense converts symmetric storage to general
## storage so symmetric classes are ok for conversion to matrix.
## unit triangular needs special handling
setAs("CsparseMatrix", "matrix",
      function(from) {
          ## |-> cholmod_C -> cholmod_dense -> chm_dense_to_matrix
          if (is(from, "triangularMatrix") && from@diag == "U")
              from <- .Call(Csparse_diagU2N, from)
          .Call(Csparse_to_matrix, from)
      })

### Some group methods:

setMethod("Arith",
      signature(e1 = "CsparseMatrix", e2 = "numeric"),
      function(e1, e2) {
          if(length(e2) == 1) { ## e.g.,  Mat ^ a
          f0 <- callGeneric(0, e2)
          if(!is.na(f0) && f0 == 0.) { # remain sparse, symm., tri.,...
              e1@x <- callGeneric(e1@x, e2)
              return(e1)
          }
          }
          ## all other (potentially non-sparse) cases: give up symm, tri,..
          callGeneric(as(e1, paste(.M.kind(e1), "gCMatrix", sep='')), e2)
      })

## The same,  e1 <-> e2 :
setMethod("Arith",
      signature(e1 = "numeric", e2 = "CsparseMatrix"),
      function(e1, e2) {
          if(length(e1) == 1) {
          f0 <- callGeneric(e1, 0)
          if(!is.na(f0) && f0 == 0.) {
              e2@x <- callGeneric(e1, e2@x)
              return(e2)
          }
          }
          callGeneric(e1, as(e2, paste(.M.kind(e2), "gCMatrix", sep='')))
      })


setMethod("Math",
      signature(x = "CsparseMatrix"),
      function(x) {
          f0 <- callGeneric(0.)
          if(!is.na(f0) && f0 == 0.) {
          ## sparseness, symm., triang.,... preserved
          x@x <- callGeneric(x@x)
          x
          } else { ## no sparseness
          callGeneric(as_dense(x))
          }
      })



### workhorse for "[<-" -- both for d* and l*  C-sparse matrices :
replCmat <- function (x, i, j, value)
{
    di <- dim(x)
    dn <- dimnames(x)
    i1 <- if(missing(i)) 0:(di[1] - 1:1) else .ind.prep2(i, 1, di, dn)
    i2 <- if(missing(j)) 0:(di[2] - 1:1) else .ind.prep2(j, 2, di, dn)
    dind <- c(length(i1), length(i2)) # dimension of replacement region
    lenRepl <- prod(dind)
    lenV <- length(value)
    if(lenV == 0) {
        if(lenRepl != 0)
            stop("nothing to replace with")
        else return(x)
    }
    ## else: lenV := length(value)   is > 0
    if(lenRepl %% lenV != 0)
        stop("number of items to replace is not a multiple of replacement length")
    if(lenV > lenRepl)
        stop("too many replacement values")

    if(is(x, "symmetricMatrix")) ## only half the indices are there..
    x <- .Call(Csparse_symmetric_to_general, x)
    clx <- c(class(x))

    xj <- .Call(Matrix_expand_pointers, x@p)
    sel <- (!is.na(match(x@i, i1)) &
            !is.na(match( xj, i2)))
    has.x <- any("x" == slotNames(x)) # i.e. *not* logical
    if(has.x && sum(sel) == lenRepl) { ## all entries to be replaced are non-zero:
        value <- rep(value, length = lenRepl)
        ## Ideally we only replace them where value != 0 and drop the value==0
        ## ones; but that would have to (?) go through dgT*
        ## v0 <- 0 == value
        ## if (lenRepl == 1) and v0 is TRUE, the following is not doing anything
        ##-  --> ./dgTMatrix.R  and its  replTmat()
        ## x@x[sel[!v0]] <- value[!v0]
        x@x[sel] <- value
        return(x)
    }
    ## else go via Tsparse.. {FIXME "waste": we already have 'xj' ..}
    x <- as(x, "TsparseMatrix")
    x[i,j] <- value
    as_CspClass(x, clx)
}

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
                                value = "replValue"),
                 function (x, i, value) replCmat(x, i=i, value=value))

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
                                value = "replValue"),
                 function (x, j, value) replCmat(x, j=j, value=value))

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "index",
                value = "replValue"),
                 replCmat)


setMethod("crossprod", signature(x = "CsparseMatrix", y = "missing"),
      function(x, y = NULL) {
          .Call(Csparse_crossprod, x, trans = FALSE, triplet = FALSE)
      })

## FIXME: Generalize the class of y.  This specific method is to replace one
##        in dgCMatrix.R
setMethod("crossprod", signature(x = "CsparseMatrix", y = "ddenseMatrix"),
      function(x, y = NULL) .Call(Csparse_dense_crossprod, x, y))

setMethod("crossprod", signature(x = "CsparseMatrix", y = "matrix"),
      function(x, y = NULL) .Call(Csparse_dense_crossprod, x, y))

setMethod("crossprod", signature(x = "CsparseMatrix", y = "numeric"),
      function(x, y = NULL) .Call(Csparse_dense_crossprod, x, y))
          
setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "missing"),
      function(x, y = NULL) {
              .Call(Csparse_crossprod, x, trans = TRUE, triplet = FALSE)
      })

setMethod("t", signature(x = "CsparseMatrix"),
      function(x) .Call(Csparse_transpose, x, is(x, "triangularMatrix")))

## FIXME (TODO):
## setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "CsparseMatrix"),
##    function(x, y)
##    .Call(Csparse_crossprod_2, x, y, trans = TRUE, triplet = FALSE)


setMethod("%*%", signature(x = "CsparseMatrix", y = "CsparseMatrix"),
          function(x, y) .Call(Csparse_Csparse_prod, x, y))

setMethod("%*%", signature(x = "CsparseMatrix", y = "ddenseMatrix"),
          function(x, y) .Call(Csparse_dense_prod, x, y))

setMethod("%*%", signature(x = "CsparseMatrix", y = "matrix"),
          function(x, y) .Call(Csparse_dense_prod, x, y))

## Not needed because of c("Matrix", "numeric") method
##setMethod("%*%", signature(x = "CsparseMatrix", y = "numeric"),
##          function(x, y) .Call(Csparse_dense_prod, x, y))

setMethod("%*%", signature(x = "ddenseMatrix", y = "CsparseMatrix"),
          function(x, y) t(.Call(Csparse_dense_crossprod, y, t(x))),
          valueClass = "dgeMatrix")

setMethod("%*%", signature(x = "matrix", y = "CsparseMatrix"),
          function(x, y) t(.Call(Csparse_dense_crossprod, y, t(x))),
          valueClass = "dgeMatrix")

## Not needed because of c("numeric", "Matrix") method
##setMethod("%*%", signature(x = "numeric", y = "CsparseMatrix"),
##          function(x, y) t(.Call(Csparse_dense_crossprod, y, x)),
##          valueClass = "dgeMatrix")
          
## NB: have extra tril(), triu() methods for symmetric ["dsC" and "lsC"]
setMethod("tril", "CsparseMatrix",
      function(x, k = 0, ...) {
          k <- as.integer(k[1])
          dd <- dim(x); sqr <- dd[1] == dd[2]
          stopifnot(-dd[1] <= k, k <= dd[1]) # had k <= 0
          r <- .Call(Csparse_band, x, -dd[1], k)
          ## return "lower triangular" if k <= 0
          if(sqr && k <= 0)
          as(r, paste(.M.kind(x), "tCMatrix", sep='')) else r
      })

setMethod("triu", "CsparseMatrix",
      function(x, k = 0, ...) {
          k <- as.integer(k[1])
          dd <- dim(x); sqr <- dd[1] == dd[2]
          stopifnot(-dd[1] <= k, k <= dd[1]) # had k >= 0
          r <- .Call(Csparse_band, x, k, dd[2])
          ## return "upper triangular" if k >= 0
          if(sqr && k >= 0)
          as(r, paste(.M.kind(x), "tCMatrix", sep='')) else r
      })

setMethod("band", "CsparseMatrix",
      function(x, k1, k2, ...) {
          k1 <- as.integer(k1[1])
          k2 <- as.integer(k2[1])
          dd <- dim(x); sqr <- dd[1] == dd[2]
          stopifnot(-dd[1] <= k1, k1 <= k2, k2 <= dd[1])
          r <- .Call(Csparse_band, x, k1, k2)
          if(sqr && k1 * k2 >= 0) ## triangular
          as(r, paste(.M.kind(x), "tCMatrix", sep=''))
          else if (k1 < 0  &&  k1 == -k2  && isSymmetric(x)) ## symmetric
          as(r, paste(.M.kind(x), "sCMatrix", sep=''))
          else
          r
      })

setMethod("diag", "CsparseMatrix",
          function(x, nrow, ncol = n) {
              dm <- .Call(Csparse_band, x, 0, 0)
              dlen <- min(dm@Dim)
              ind1 <- dm@i + 1:1        # 1-based index vector
              if (is(dm, "lMatrix")) {
                  val <- rep.int(FALSE, dlen)
                  val[ind1] <- TRUE
                  return(val)
              }
              val <- rep.int(0, dlen)
              ## cMatrix not yet active but for future expansion
              if (is(dm, "cMatrix")) val <- as.complex(val)
              val[ind1] <- dm@x
              val
          })


setMethod("colSums", signature(x = "CsparseMatrix"), .as.dgC.Fun,
      valueClass = "numeric")
setMethod("colMeans", signature(x = "CsparseMatrix"), .as.dgC.Fun,
      valueClass = "numeric")
setMethod("rowSums", signature(x = "CsparseMatrix"), .as.dgC.Fun,
      valueClass = "numeric")
setMethod("rowMeans", signature(x = "CsparseMatrix"), .as.dgC.Fun,
      valueClass = "numeric")