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 1592, Thu Sep 28 15:31:17 2006 UTC
# Line 1  Line 1
1  #### "Namespace private" Auxiliaries  such as method functions  #### "Namespace private" Auxiliaries  such as method functions
2  #### (called from more than one place --> need to be defined early)  #### (called from more than one place --> need to be defined early)
3
4    .isR_24 <- (paste(R.version$major, R.version$minor, sep=".") >= "2.4")
5
6    ## Need to consider NAs ;  "== 0" even works for logical & complex:
7    is0  <- function(x) !is.na(x) & x == 0
8    isN0 <- function(x)  is.na(x) | x != 0
9    all0 <- function(x) !any(is.na(x)) && all(x == 0)
10
11    allTrue  <- function(x) !any(is.na(x)) && all(x)
12    allFalse <- function(x) !any(is.na(x)) && !any(x)
13
14
15  ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):  ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):
16  .M.v <- function(x, y) callGeneric(x, as.matrix(y))  .M.v <- function(x, y) callGeneric(x, as.matrix(y))
17  .v.M <- function(x, y) callGeneric(rbind(x), y)  .v.M <- function(x, y) callGeneric(rbind(x), y)
18
19    .M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL)
20
21  .has.DN <- ## has non-trivial Dimnames slot?  .has.DN <- ## has non-trivial Dimnames slot?
22      function(x) !identical(list(NULL,NULL), x@Dimnames)      function(x) !identical(list(NULL,NULL), x@Dimnames)
23
# Line 17  Line 30
30                    fun, cl1, cl2), call. = FALSE)                    fun, cl1, cl2), call. = FALSE)
31  }  }
32
33    ## This should be done in C and be exported by 'methods':  [FIXME - ask JMC ]
34    copyClass <- function(x, newCl, sNames =
35                          intersect(slotNames(newCl), slotNames(x))) {
36        r <- new(newCl)
37        for(n in sNames)
38            slot(r, n) <- slot(x, n)
39        r
40    }
41
42  ## chol() via "dpoMatrix"  ## chol() via "dpoMatrix"
43  cholMat <- function(x, pivot, LINPACK) {  cholMat <- function(x, pivot, ...) {
44      px <- as(x, "dpoMatrix")      px <- as(x, "dpoMatrix")
45      if (isTRUE(validObject(px, test=TRUE))) chol(px)      if (isTRUE(validObject(px, test=TRUE))) chol(px)
46      else stop("'x' is not positive definite -- chol() undefined.")      else stop("'x' is not positive definite -- chol() undefined.")
# Line 82  Line 104
104      da[2]      da[2]
105  }  }
106
107    ## Note: !isPacked(.)  i.e. `full' still contains
108    ## ----  "*sy" and "*tr" which have "undefined" lower or upper part
109    isPacked <- function(x)
110    {
111        ## Is 'x' a packed (dense) matrix ?
112        is(x, "denseMatrix") &&
113        any("x" == slotNames(x)) && length(x@x) < prod(dim(x))
114    }
115
116  emptyColnames <- function(x)  emptyColnames <- function(x)
117  {  {
118      ## Useful for compact printing of (parts) of sparse matrices      ## Useful for compact printing of (parts) of sparse matrices
# Line 90  Line 121
121      x      x
122  }  }
123
124    ### TODO:  write in C and port to base (or 'utils') R
125    indTri <- function(n, upper = TRUE) {
126        ## == which(upper.tri(diag(n)) or
127        ##    which(lower.tri(diag(n)) -- but much more efficiently for largish 'n'
128        stopifnot(length(n) == 1, n == (n. <- as.integer(n)), (n <- n.) >= 0)
129        if(n <= 2)
130            return(if(n == 2) as.integer(if(upper) n+1 else n) else integer(0))
131        ## First, compute the 'diff(.)'  fast.  Use integers
132        one <- 1:1 ; two <- 2:2
133        n1 <- n - one
134        n2 <- n1 - one
135        r <- rep.int(one, n*n1/two - one)
136        r[cumsum(if(upper) 1:n2 else c(n1, if(n >= 4) n2:two))] <- if(upper) n:3 else 3:n
137        ## now have "dliu" difference; revert to "liu":
138        cumsum(c(if(upper) n+one else two, r))
139    }
140
141
142  prTriang <- function(x, digits = getOption("digits"),  prTriang <- function(x, digits = getOption("digits"),
143                         maxp = getOption("max.print"),
144                       justify = "none", right = TRUE)                       justify = "none", right = TRUE)
145  {  {
146      ## modeled along stats:::print.dist      ## modeled along stats:::print.dist
diag <- TRUE
147      upper <- x@uplo == "U"      upper <- x@uplo == "U"
148
149      m <- as(x, "matrix")      m <- as(x, "matrix")
# Line 103  Line 152
152          cf[row(cf) > col(cf)] <- "."          cf[row(cf) > col(cf)] <- "."
153      else      else
154          cf[row(cf) < col(cf)] <- "."          cf[row(cf) < col(cf)] <- "."
155      print(cf, quote = FALSE, right = right)      if(.isR_24)
156             print(cf, quote = FALSE, right = right, max = maxp)
157        else print(cf, quote = FALSE, right = right)
158      invisible(x)      invisible(x)
159  }  }
160
161  prMatrix <- function(x, digits = getOption("digits")) {  prMatrix <- function(x, digits = getOption("digits"),
162                         maxp = getOption("max.print")) {
163      d <- dim(x)      d <- dim(x)
164      cl <- class(x)      cl <- class(x)
165      cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))      cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
maxp <- getOption("max.print")
166      if(prod(d) <= maxp) {      if(prod(d) <= maxp) {
167          if(is(x, "triangularMatrix"))          if(is(x, "triangularMatrix"))
168              prTriang(x, digits = digits)              prTriang(x, digits = digits, maxp = maxp)
169          else          else {
170              print(as(x, "matrix"), digits = digits)              if(.isR_24)
171                     print(as(x, "matrix"), digits = digits, max = maxp)
172                else print(as(x, "matrix"), digits = digits)
173            }
174      }      }
175      else { ## d[1] > maxp / d[2] >= nr :      else { ## d[1] > maxp / d[2] >= nr :
176          m <- as(x, "matrix")          m <- as(x, "matrix")
# Line 130  Line 184
184      invisible(x)# as print() S3 methods do      invisible(x)# as print() S3 methods do
185  }  }
186
187    nonFALSE <- function(x) {
188        ## typically used for lMatrices:  (TRUE,NA,FALSE) |-> (TRUE,FALSE)
189        if(any(ix <- is.na(x))) x[ix] <- TRUE
190        x
191    }
192
193    nz.NA <- function(x, na.value) {
194        ## Non-Zeros of x
195        ## na.value: TRUE: NA's give TRUE, they are not 0
196        ##             NA: NA's are not known ==> result := NA
197        ##          FALSE: NA's give FALSE, could be 0
198        stopifnot(is.logical(na.value) && length(na.value) == 1)
199        if(is.na(na.value)) x != 0
200        else  if(na.value)  isN0(x)
201        else                x != 0 & !is.na(x)
202    }
203
204    ## Number of non-zeros :
205    ## FIXME? -- make this into a generic function (?)
206    nnzero <- function(x, na.counted = NA) {
207        ## na.counted: TRUE: NA's are counted, they are not 0
208        ##               NA: NA's are not known (0 or not) ==>  result := NA
209        ##            FALSE: NA's are omitted before counting
210        cl <- class(x)
211        if(!extends(cl, "Matrix"))
212            sum(nz.NA(x, na.counted))
213        else if(extends(cl, "sparseMatrix"))
214            ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!
215           switch(.sp.class(cl),
216                   "CsparseMatrix" = length(x@i),
217                   "TsparseMatrix" = length(x@i),
218                   "RsparseMatrix" = length(x@j))
219        else ## denseMatrix
220            sum(nz.NA(as_geClass(x, cl)@x, na.counted))
221    }
222
223  ## For sparseness handling  ## For sparseness handling
224    ## return a 2-column (i,j) matrix of
225    ## 0-based indices of non-zero entries  :
226  non0ind <- function(x) {  non0ind <- function(x) {
if(is.numeric(x))
return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))

## else return a (i,j) matrix of non-zero-indices
227
228        if(is.numeric(x))
229            return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0))
230        ## else
231      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
232      if(is(x, "TsparseMatrix"))      non0.i <- function(M) {
233          return(unique(cbind(x@i,x@j)))          if(is(M, "TsparseMatrix"))
234                return(unique(cbind(M@i,M@j)))
235            if(is(M, "pMatrix"))
236                return(cbind(seq(length=nrow(M)), M@perm) - 1:1)
237            ## else:
238            isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
239            .Call(compressed_non_0_ij, M, isC)
240        }
241
242        if(is(x, "symmetricMatrix")) { # also get "other" triangle
243            ij <- non0.i(x)
244            notdiag <- ij[,1] != ij[,2]# but not the diagonals again
245            rbind(ij, ij[notdiag, 2:1])
246        }
247        else if(is(x, "triangularMatrix")) { # check for "U" diag
248            if(x@diag == "U") {
249                i <- seq(length = dim(x)[1]) - 1:1
250                rbind(non0.i(x), cbind(i,i))
251            } else non0.i(x)
252        }
253        else
254            non0.i(x)
255    }
256
257    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
258    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
259    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
260    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
261
262      isCol <- function(M) any("i" == slotNames(M))  complementInd <- function(ij, dim)
263      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  {
264        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
265        ##          but as 1-based indices
266        n <- prod(dim)
267        if(n == 0) return(integer(0))
268        ii <- 1:n
269        ii[-(1 + encodeInd(ij, nr = dim[1]))]
270    }
271
272    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
273
274    intersectInd <- function(ij1, ij2, nrow) {
275        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
276        ## return only the *common* entries
277        decodeInd(intersect(encodeInd(ij1, nrow),
278                            encodeInd(ij2, nrow)), nrow)
279    }
280
281    WhichintersectInd <- function(ij1, ij2, nrow) {
282        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
283        ## find *where*  common entries are in ij1 & ij2
284        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
285        ni <- !is.na(m1)
286        list(which(ni), m1[ni])
287  }  }
288
289
290  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
291  uniq <- function(x) {
292      if(is(x, "TsparseMatrix")) {  uniqTsparse <- function(x, class.x = c(class(x))) {
293          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
294          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
295          ## -----------------------------------------------------------          ## -----------------------------------------------------------
296          ## The following is *not* efficient {but easy to program}:      ## The following is not quite efficient {but easy to program,
297          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")      ## and as() are based on C code  (all of them?)
298          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")      ##
299          else stop("not implemented for class", class(x))      ## FIXME: Do it fast for the case where 'x' is already 'uniq'
300
301        switch(class.x,
302               "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
303               "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
304               "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
305               ## do we need this for "logical" ones, there's no sum() there!
306               "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
307               "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
308               "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
309               ## do we need this for "logical" ones, there's no sum() there!
310               "ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"),
311               "nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"),
312               "ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"),
313               ## otherwise:
314               stop("not yet implemented for class ", class.x))
315    }
316
317    ## Note: maybe, using
318    ## ----    xj <- .Call(Matrix_expand_pointers, x@p)
319    ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
320    ## but really efficient would be to use only one .Call(.) for uniq(.) !
321
322      } else x      # not 'gT' ; i.e. "uniquely" represented in any case  uniq <- function(x) {
323        if(is(x, "TsparseMatrix")) uniqTsparse(x) else x
324        ## else:  not 'Tsparse', i.e. "uniquely" represented in any case
325    }
326
327    asTuniq <- function(x) {
328        if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix")
329  }  }
330
331  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 333
333  {  {
334      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
335      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
336      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
337      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
338      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
339      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 220  Line 389
389      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
390  }  }
391
392    ## -> ./ndenseMatrix.R :
393    n2d_Matrix <- function(from) {
394        stopifnot(is(from, "nMatrix"))
395        fixupDense(new(sub("^n", "d", class(from)),
396                       x = as.double(from@x),
397                       Dim = from@Dim, Dimnames = from@Dimnames),
398                   from)
399        ## FIXME: treat 'factors' smartly {not for triangular!}
400    }
401    n2l_spMatrix <- function(from) {
402        stopifnot(is(from, "nMatrix"))
403        new(sub("^n", "l", class(from)),
404            ##x = as.double(from@x),
405            Dim = from@Dim, Dimnames = from@Dimnames)
406    }
407
408  if(FALSE)# unused  if(FALSE)# unused
409  l2d_meth <- function(x) {  l2d_meth <- function(x) {
410      cl <- class(x)      cl <- class(x)
411      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
412  }  }
413
414    ## return "d" or "l" or "n" or "z"
415    .M.kind <- function(x, clx = class(x)) {
416        if(is.matrix(x)) { ## 'old style matrix'
417            if     (is.numeric(x)) "d"
418            else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ??
419            else if(is.complex(x)) "z"
420            else stop("not yet implemented for matrix w/ typeof ", typeof(x))
421        }
422        else if(extends(clx, "dMatrix")) "d"
423        else if(extends(clx, "nMatrix")) "n"
424        else if(extends(clx, "lMatrix")) "l"
425        else if(extends(clx, "zMatrix")) "z"
426        else if(extends(clx, "pMatrix")) "n" # permutation -> pattern
427        else stop(" not yet be implemented for ", clx)
428    }
429
430    .M.shape <- function(x, clx = class(x)) {
431        if(is.matrix(x)) { ## 'old style matrix'
432            if     (isDiagonal  (x)) "d"
433            else if(isTriangular(x)) "t"
434            else if(isSymmetric (x)) "s"
435            else "g" # general
436        }
437        else if(extends(clx, "diagonalMatrix"))  "d"
438        else if(extends(clx, "triangularMatrix"))"t"
439        else if(extends(clx, "symmetricMatrix")) "s"
440        else "g"
441    }
442
443
444    class2 <- function(cl, kind = "l", do.sub = TRUE) {
445        ## Find "corresponding" class; since pos.def. matrices have no pendant:
446        if     (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
447        else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
448        else if(do.sub) sub("^d", kind, cl)
449        else cl
450    }
451
452    geClass <- function(x) {
453        if     (is(x, "dMatrix")) "dgeMatrix"
454        else if(is(x, "lMatrix")) "lgeMatrix"
455        else if(is(x, "nMatrix")) "ngeMatrix"
456        else if(is(x, "zMatrix")) "zgeMatrix"
457        else stop("general Matrix class not yet implemented for ",
458                  class(x))
459    }
460
461    .dense.prefixes <- c("d" = "di",
462                         "t" = "tr",
463                         "s" = "sy",
464                         "g" = "ge")
465
466    .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
467                          "t" = "t",
468                          "s" = "s",
469                          "g" = "g")
470
471    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
472    as_dense <- function(x) {
473        as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
474    }
475
476    .sp.class <- function(x) { ## find and return the "sparseness class"
477        if(!is.character(x)) x <- class(x)
478        for(cl in paste(c("C","T","R"), "sparseMatrix", sep=''))
479            if(extends(x, cl))
480                return(cl)
481        ## else (should rarely happen)
482        as.character(NA)
483    }
484
485    as_Csparse <- function(x) {
486        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
487    }
488    as_Rsparse <- function(x) {
489        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
490    }
491    as_Tsparse <- function(x) {
492        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
493    }
494
495    as_geClass <- function(x, cl) {
496        if     (extends(cl, "diagonalMatrix")  && isDiagonal(x))
497            as(x, cl)
498        else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x)) {
499            kind <- .M.kind(x)
500            as(x, class2(cl, kind, do.sub= kind != "d"))
501        } else if(extends(cl, "triangularMatrix") && isTriangular(x))
502            as(x, cl)
503        else
504            as(x, paste(.M.kind(x), "geMatrix", sep=''))
505    }
506
507    as_CspClass <- function(x, cl) {
508        if ((extends(cl, "diagonalMatrix")  && isDiagonal(x)) ||
509            (extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||
510            (extends(cl, "triangularMatrix")&& isTriangular(x)))
511            as(x, cl)
512        else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
513    }
514
515
516  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
517  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
518      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
519      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
520                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
521                 from)                 from)
522      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 534
534      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
535  }  }
536
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)))
}
537
538  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
539    isTriMat <- function(object, upper = NA) {
540      ## pretest: is it square?      ## pretest: is it square?
541      d <- dim(object)      d <- dim(object)
542      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
543      ## else slower test      ## else slower test
544      if(!is.matrix(object))      if(!is.matrix(object))
545          object <- as(object,"matrix")          object <- as(object,"matrix")
546      ## == 0 even works for logical & complex:      if(is.na(upper)) {
547      if(upper)          if(all0(object[lower.tri(object)]))
548          all(object[lower.tri(object)] == 0)              structure(TRUE, kind = "U")
549      else          else if(all0(object[upper.tri(object)]))
550          all(object[upper.tri(object)] == 0)              structure(TRUE, kind = "L")
551            else FALSE
552        } else if(upper)
553            all0(object[lower.tri(object)])
554        else ## upper is FALSE
555            all0(object[upper.tri(object)])
556    }
557
558    ## For Csparse matrices
559    isTriC <- function(x, upper = NA) {
560        ## pretest: is it square?
561        d <- dim(x)
562        if(d[1] != d[2]) return(FALSE)
563        ## else
564        if(d[1] == 0) return(TRUE)
565        ni <- 1:d[2]
566        ## the row indices split according to column:
567        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
568        lil <- unlist(lapply(ilist, length), use.names = FALSE)
569        if(any(lil == 0)) {
570            pos <- lil > 0
571            if(!any(pos)) ## matrix of all 0's
572                return(TRUE)
573            ilist <- ilist[pos]
574            ni <- ni[pos]
575        }
576        ni0 <- ni - 1:1 # '0-based ni'
577        if(is.na(upper)) {
578            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
579                structure(TRUE, kind = "U")
580            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
581                structure(TRUE, kind = "L")
582            else FALSE
583        } else if(upper) {
584            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
585        } else { ## 'lower'
586            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
587        }
588  }  }
589
590  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
591        ## "matrix" or "denseMatrix" (but not "diagonalMatrix")
592      d <- dim(object)      d <- dim(object)
593      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
594      else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)      else if(is.matrix(object))
595            ## requires that "vector-indexing" works for 'object' :
596            all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
597        else ## "denseMatrix" -- packed or unpacked
598            if(is(object, "generalMatrix")) # "dge", "lge", ...
599                all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
600            else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:
601                ## -> has 'uplo'  differentiate between packed and unpacked
602
603    ### .......... FIXME ...............
604
605                packed <- isPacked(object)
606                if(object@uplo == "U") {
607                } else { ## uplo == "L"
608                }
609
610    ### very cheap workaround
611                all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
612            }
613    }
614
615
616    ## FIXME? -- this should also work for "ltT", "ntT", ... :
617    diagU2N <- function(x)
618    {
619        ## Purpose: Transform a *unit diagonal* sparse triangular matrix
620        ##  into one with explicit diagonal entries '1'
621        xT <- as(x, "dgTMatrix")
622        ## leave it as  T* - the caller can always coerce to C* if needed:
623        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
624            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
625  }  }
626
627    ## FIXME: this should probably be dropped / replaced by as_Csparse
628    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
629        x <- as(x, "dgCMatrix")
630        callGeneric()
631    }
632
633    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
634        ## used e.g. inside colSums() etc methods
635        x <- as(x, "dgTMatrix")
636        callGeneric()
637    }
638
639
640    ### Fast much simplified version of tapply()
641    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
642        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
643    }
644
645    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
646    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
647    ## }
648

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

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: