# SCM Repository

[matrix] Diff of /pkg/R/Auxiliaries.R
 [matrix] / pkg / R / Auxiliaries.R

# Diff of /pkg/R/Auxiliaries.R

revision 1198, Mon Jan 23 15:01:02 2006 UTC revision 1270, Thu May 18 06:44:43 2006 UTC
# Line 82  Line 82
82      da[2]      da[2]
83  }  }
84
85    isPacked <- function(x)
86    {
87        ## Is 'x' a packed (dense) matrix ?
88        is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
89    }
90
91  emptyColnames <- function(x)  emptyColnames <- function(x)
92  {  {
93      ## Useful for compact printing of (parts) of sparse matrices      ## Useful for compact printing of (parts) of sparse matrices
# Line 134  Line 140
140  non0ind <- function(x) {  non0ind <- function(x) {
141      if(is.numeric(x))      if(is.numeric(x))
142          return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))          return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))
143        ## else
## else return a (i,j) matrix of non-zero-indices

144      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
145        ## return a 2-column (i,j) matrix of
146        ## 0-based indices of non-zero entries  :
147      if(is(x, "TsparseMatrix"))      if(is(x, "TsparseMatrix"))
148          return(unique(cbind(x@i,x@j)))          return(unique(cbind(x@i,x@j)))
149        ## else:
150        isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
151        .Call("compressed_non_0_ij", x, isC, PACKAGE = "Matrix")
152    }
153
154    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
155    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
156    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
157    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
158
159      isCol <- function(M) any("i" == slotNames(M))  complementInd <- function(ij, dim)
160      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  {
161        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
162        ##          but as 1-based indices
163        n <- prod(dim)
164        if(n == 0) return(integer(0))
165        ii <- 1:n
166        ii[-(1 + encodeInd(ij, nr = dim[1]))]
167    }
168
169    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
170
171    intersectInd <- function(ij1, ij2, nrow) {
172        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
173        ## return only the *common* entries
174        decodeInd(intersect(encodeInd(ij1, nrow),
175                            encodeInd(ij2, nrow)), nrow)
176    }
177
178    WhichintersectInd <- function(ij1, ij2, nrow) {
179        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
180        ## find *where*  common entries are in ij1 & ij2
181        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
182        ni <- !is.na(m1)
183        list(which(ni), m1[ni])
184  }  }
185
186
187  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
188  uniq <- function(x) {  uniq <- function(x) {
189
190    ### Note: maybe, using
191    ### ----    xj <- .Call("Matrix_expand_pointers", x@p)
192    ### would be slightly more efficient than as( <dgC> , "dgTMatrix")
193    ### but really efficient would be to use only one .Call(.) for uniq(.) !
194    ### Try to do it particularly fast for the case where 'x' is already a 'uniq' <dgT>
195
196      if(is(x, "TsparseMatrix")) {      if(is(x, "TsparseMatrix")) {
197          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
198          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
199          ## -----------------------------------------------------------          ## -----------------------------------------------------------
200          ## The following is *not* efficient {but easy to program}:          ## The following is not quite efficient {but easy to program,
201            ## and both as() are based on C code
202          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
203          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
204          else stop("not implemented for class", class(x))          else stop("not implemented for class", class(x))
205
206      } else x      # not 'gT' ; i.e. "uniquely" represented in any case      } else x  ## not 'gT' ; i.e. "uniquely" represented in any case
207  }  }
208
209  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 211
211  {  {
212      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
213      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
214      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
215      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
216      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
217      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 273
273      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
274  }  }
275
276    dClass2 <- function(dClass, kind = "l") {
277        ## Find "corresponding" class for a dMatrix;
278        #  since pos.def. matrices have no pendant:
279        if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
280        else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
281        else sub("^d", kind, dClass)
282    }
283
284    geClass <- function(x) {
285        if(is(x, "dMatrix")) "dgeMatrix"
286        else if(is(x, "lMatrix")) "lgeMatrix"
287        else stop("general Matrix class not yet implemented for",
288                  class(x))
289    }
290
291  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
292  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
293      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
294      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
295                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
296                 from)                 from)
297      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 260  Line 322
322                           useInherited = c(from = TRUE, to = FALSE)))                           useInherited = c(from = TRUE, to = FALSE)))
323  }  }
324
325  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
326    isTriMat <- function(object, upper = NA) {
327      ## pretest: is it square?      ## pretest: is it square?
328      d <- dim(object)      d <- dim(object)
329      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 331
331      if(!is.matrix(object))      if(!is.matrix(object))
332          object <- as(object,"matrix")          object <- as(object,"matrix")
333      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
334      if(upper)      if(is.na(upper)) {
335            if(all(object[lower.tri(object)] == 0))
336                structure(TRUE, kind = "U")
337            else if(all(object[upper.tri(object)] == 0))
338                structure(TRUE, kind = "L")
339            else FALSE
340        } else if(upper)
341          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
342      else      else ## upper is FALSE
343          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
344  }  }
345
346    ## For Csparse matrices
347    isTriC <- function(x, upper = NA) {
348        ## pretest: is it square?
349        d <- dim(x)
350        if(d[1] != d[2]) return(FALSE)
351        ## else
352        if(d[1] == 0) return(TRUE)
353        ni <- 1:d[1]
354        ## the row indices split according to column:
355        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
356        lil <- unlist(lapply(ilist, length), use.names = FALSE)
357        if(any(lil == 0)) {
358            pos <- lil > 0
359            if(!any(pos)) ## matrix of all 0's
360                return(TRUE)
361            ilist <- ilist[pos]
362            ni <- ni[pos]
363        }
364        if(is.na(upper)) {
365            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni))
366                structure(TRUE, kind = "U")
367            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni))
368                structure(TRUE, kind = "L")
369            else FALSE
370        } else if(upper) {
371            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni)
372        } else { ## 'lower'
373            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni)
374        }
375    }
376
377  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
378      d <- dim(object)      d <- dim(object)
379      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
380      else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)      else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)
381  }  }
382
383    diagU2N <- function(x)
384    {
385        ## Purpose: Transform a *unit diagonal* triangular matrix
386        ##  into one with explicit diagonal entries '1'
387        xT <- as(x, "dgTMatrix")
388        ## leave it as  T* - the caller can always coerce to C* if needed:
389        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
390            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
391    }

Legend:
 Removed from v.1198 changed lines Added in v.1270