# 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 1295, Fri Jun 9 21:47:22 2006 UTC
# Line 82  Line 82
82      da[2]      da[2]
83  }  }
84
85    ## Note: !isPacked(.)  i.e. `full' still contains
86    ## ----  "*sy" and "*tr" which have "undefined" lower or upper part
87    isPacked <- function(x)
88    {
89        ## Is 'x' a packed (dense) matrix ?
90        is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
91    }
92
93  emptyColnames <- function(x)  emptyColnames <- function(x)
94  {  {
95      ## Useful for compact printing of (parts) of sparse matrices      ## Useful for compact printing of (parts) of sparse matrices
# Line 134  Line 142
142  non0ind <- function(x) {  non0ind <- function(x) {
143      if(is.numeric(x))      if(is.numeric(x))
144          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))
145        ## else
## else return a (i,j) matrix of non-zero-indices

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

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