SCM

SCM Repository

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

Diff of /pkg/R/Auxiliaries.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1198, Mon Jan 23 15:01:02 2006 UTC revision 1673, Mon Nov 6 20:54:26 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    .isR_25 <- (paste(R.version$major, R.version$minor, sep=".") >= "2.5")
6    
7    ## Need to consider NAs ;  "== 0" even works for logical & complex:
8    is0  <- function(x) !is.na(x) & x == 0
9    isN0 <- function(x)  is.na(x) | x != 0
10    all0 <- function(x) !any(is.na(x)) && all(x == 0)
11    
12    allTrue  <- function(x) !any(is.na(x)) && all(x)
13    allFalse <- function(x) !any(is.na(x)) && !any(x)
14    
15    
16  ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):  ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):
17  .M.v <- function(x, y) callGeneric(x, as.matrix(y))  .M.v <- function(x, y) callGeneric(x, as.matrix(y))
18  .v.M <- function(x, y) callGeneric(rbind(x), y)  .v.M <- function(x, y) callGeneric(rbind(x), y)
19    
20    .M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL)
21    
22  .has.DN <- ## has non-trivial Dimnames slot?  .has.DN <- ## has non-trivial Dimnames slot?
23      function(x) !identical(list(NULL,NULL), x@Dimnames)      function(x) !identical(list(NULL,NULL), x@Dimnames)
24    
25  .bail.out.1 <- function(fun, cl) {  .bail.out.1 <- function(fun, cl) {
26      stop(gettextf('not-yet-implemented method for %s(<%s>)', fun, cl),      stop(gettextf('not-yet-implemented method for %s(<%s>).\n ->>  Ask the package authors to implement the missing feature.', fun, cl),
27           call. = FALSE)           call. = FALSE)
28  }  }
29  .bail.out.2 <- function(fun, cl1, cl2) {  .bail.out.2 <- function(fun, cl1, cl2) {
30      stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)',      stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>).\n ->>  Ask the package authors to implement the missing feature.',
31                    fun, cl1, cl2), call. = FALSE)                    fun, cl1, cl2), call. = FALSE)
32  }  }
33    
34    ## This should be done in C and be exported by 'methods':  [FIXME - ask JMC ]
35    copyClass <- function(x, newCl, sNames =
36                          intersect(slotNames(newCl), slotNames(x))) {
37        r <- new(newCl)
38        for(n in sNames)
39            slot(r, n) <- slot(x, n)
40        r
41    }
42    
43  ## chol() via "dpoMatrix"  ## chol() via "dpoMatrix"
44  cholMat <- function(x, pivot, LINPACK) {  cholMat <- function(x, pivot, ...) {
45      px <- as(x, "dpoMatrix")      px <- as(x, "dpoMatrix")
46      if (isTRUE(validObject(px, test=TRUE))) chol(px)      if (isTRUE(validObject(px, test=TRUE))) chol(px)
47      else stop("'x' is not positive definite -- chol() undefined.")      else stop("'x' is not positive definite -- chol() undefined.")
# Line 82  Line 105 
105      da[2]      da[2]
106  }  }
107    
108    ## Note: !isPacked(.)  i.e. `full' still contains
109    ## ----  "*sy" and "*tr" which have "undefined" lower or upper part
110    isPacked <- function(x)
111    {
112        ## Is 'x' a packed (dense) matrix ?
113        is(x, "denseMatrix") &&
114        any("x" == slotNames(x)) && length(x@x) < prod(dim(x))
115    }
116    
117  emptyColnames <- function(x)  emptyColnames <- function(x)
118  {  {
119      ## Useful for compact printing of (parts) of sparse matrices      ## Useful for compact printing of (parts) of sparse matrices
# Line 90  Line 122 
122      x      x
123  }  }
124    
125    ### TODO:  write in C and port to base (or 'utils') R
126    indTri <- function(n, upper = TRUE) {
127        ## == which(upper.tri(diag(n)) or
128        ##    which(lower.tri(diag(n)) -- but much more efficiently for largish 'n'
129        stopifnot(length(n) == 1, n == (n. <- as.integer(n)), (n <- n.) >= 0)
130        if(n <= 2)
131            return(if(n == 2) as.integer(if(upper) n+1 else n) else integer(0))
132        ## First, compute the 'diff(.)'  fast.  Use integers
133        one <- 1:1 ; two <- 2:2
134        n1 <- n - one
135        n2 <- n1 - one
136        r <- rep.int(one, n*n1/two - one)
137        r[cumsum(if(upper) 1:n2 else c(n1, if(n >= 4) n2:two))] <- if(upper) n:3 else 3:n
138        ## now have "dliu" difference; revert to "liu":
139        cumsum(c(if(upper) n+one else two, r))
140    }
141    
142    
143  prTriang <- function(x, digits = getOption("digits"),  prTriang <- function(x, digits = getOption("digits"),
144                         maxp = getOption("max.print"),
145                       justify = "none", right = TRUE)                       justify = "none", right = TRUE)
146  {  {
147      ## modeled along stats:::print.dist      ## modeled along stats:::print.dist
     diag <- TRUE  
148      upper <- x@uplo == "U"      upper <- x@uplo == "U"
149    
150      m <- as(x, "matrix")      m <- as(x, "matrix")
# Line 103  Line 153 
153          cf[row(cf) > col(cf)] <- "."          cf[row(cf) > col(cf)] <- "."
154      else      else
155          cf[row(cf) < col(cf)] <- "."          cf[row(cf) < col(cf)] <- "."
156      print(cf, quote = FALSE, right = right)      if(.isR_24)
157             print(cf, quote = FALSE, right = right, max = maxp)
158        else print(cf, quote = FALSE, right = right)
159      invisible(x)      invisible(x)
160  }  }
161    
162  prMatrix <- function(x, digits = getOption("digits")) {  prMatrix <- function(x, digits = getOption("digits"),
163                         maxp = getOption("max.print")) {
164      d <- dim(x)      d <- dim(x)
165      cl <- class(x)      cl <- class(x)
166      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")  
167      if(prod(d) <= maxp) {      if(prod(d) <= maxp) {
168          if(is(x, "triangularMatrix"))          if(is(x, "triangularMatrix"))
169              prTriang(x, digits = digits)              prTriang(x, digits = digits, maxp = maxp)
170          else          else {
171              print(as(x, "matrix"), digits = digits)              if(.isR_24)
172                     print(as(x, "matrix"), digits = digits, max = maxp)
173                else print(as(x, "matrix"), digits = digits)
174            }
175      }      }
176      else { ## d[1] > maxp / d[2] >= nr :      else { ## d[1] > maxp / d[2] >= nr :
177          m <- as(x, "matrix")          m <- as(x, "matrix")
# Line 130  Line 185 
185      invisible(x)# as print() S3 methods do      invisible(x)# as print() S3 methods do
186  }  }
187    
188    nonFALSE <- function(x) {
189        ## typically used for lMatrices:  (TRUE,NA,FALSE) |-> (TRUE,FALSE)
190        if(any(ix <- is.na(x))) x[ix] <- TRUE
191        x
192    }
193    
194    nz.NA <- function(x, na.value) {
195        ## Non-Zeros of x
196        ## na.value: TRUE: NA's give TRUE, they are not 0
197        ##             NA: NA's are not known ==> result := NA
198        ##          FALSE: NA's give FALSE, could be 0
199        stopifnot(is.logical(na.value) && length(na.value) == 1)
200        if(is.na(na.value)) x != 0
201        else  if(na.value)  isN0(x)
202        else                x != 0 & !is.na(x)
203    }
204    
205    ## Number of non-zeros :
206    ## FIXME? -- make this into a generic function (?)
207    nnzero <- function(x, na.counted = NA) {
208        ## na.counted: TRUE: NA's are counted, they are not 0
209        ##               NA: NA's are not known (0 or not) ==>  result := NA
210        ##            FALSE: NA's are omitted before counting
211        cl <- class(x)
212        if(!extends(cl, "Matrix"))
213            sum(nz.NA(x, na.counted))
214        else if(extends(cl, "sparseMatrix"))
215            ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!
216           switch(.sp.class(cl),
217                   "CsparseMatrix" = length(x@i),
218                   "TsparseMatrix" = length(x@i),
219                   "RsparseMatrix" = length(x@j))
220        else ## denseMatrix
221            sum(nz.NA(as_geClass(x, cl)@x, na.counted))
222    }
223    
224  ## For sparseness handling  ## For sparseness handling
225    ## return a 2-column (i,j) matrix of
226    ## 0-based indices of non-zero entries  :
227  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  
228    
229        if(is.numeric(x))
230            return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0))
231        ## else
232      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
233      if(is(x, "TsparseMatrix"))      non0.i <- function(M) {
234          return(unique(cbind(x@i,x@j)))          if(is(M, "TsparseMatrix"))
235                return(unique(cbind(M@i,M@j)))
236            if(is(M, "pMatrix"))
237                return(cbind(seq_len(nrow(M)), M@perm) - 1:1)
238            ## else:
239            isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
240            .Call(compressed_non_0_ij, M, isC)
241        }
242    
243        if(is(x, "symmetricMatrix")) { # also get "other" triangle
244            ij <- non0.i(x)
245            notdiag <- ij[,1] != ij[,2]# but not the diagonals again
246            rbind(ij, ij[notdiag, 2:1])
247        }
248        else if(is(x, "triangularMatrix")) { # check for "U" diag
249            if(x@diag == "U") {
250                i <- seq_len(dim(x)[1]) - 1:1
251                rbind(non0.i(x), cbind(i,i))
252            } else non0.i(x)
253        }
254        else
255            non0.i(x)
256    }
257    
258    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
259    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
260    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
261    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
262    
263      isCol <- function(M) any("i" == slotNames(M))  complementInd <- function(ij, dim)
264      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  {
265        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
266        ##          but as 1-based indices
267        n <- prod(dim)
268        if(n == 0) return(integer(0))
269        ii <- 1:n
270        ii[-(1 + encodeInd(ij, nr = dim[1]))]
271    }
272    
273    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
274    
275    intersectInd <- function(ij1, ij2, nrow) {
276        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
277        ## return only the *common* entries
278        decodeInd(intersect(encodeInd(ij1, nrow),
279                            encodeInd(ij2, nrow)), nrow)
280    }
281    
282    WhichintersectInd <- function(ij1, ij2, nrow) {
283        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
284        ## find *where*  common entries are in ij1 & ij2
285        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
286        ni <- !is.na(m1)
287        list(which(ni), m1[ni])
288  }  }
289    
290    
291  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
292  uniq <- function(x) {  
293      if(is(x, "TsparseMatrix")) {  uniqTsparse <- function(x, class.x = c(class(x))) {
294          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
295          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
296          ## -----------------------------------------------------------          ## -----------------------------------------------------------
297          ## The following is *not* efficient {but easy to program}:      ## The following is not quite efficient {but easy to program,
298          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")      ## and as() are based on C code  (all of them?)
299          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")      ##
300          else stop("not implemented for class", class(x))      ## FIXME: Do it fast for the case where 'x' is already 'uniq'
301    
302        switch(class.x,
303               "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
304               "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
305               "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
306               ## do we need this for "logical" ones, there's no sum() there!
307               "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
308               "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
309               "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
310               ## do we need this for "logical" ones, there's no sum() there!
311               "ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"),
312               "nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"),
313               "ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"),
314               ## otherwise:
315               stop("not yet implemented for class ", class.x))
316    }
317    
318    ## Note: maybe, using
319    ## ----    xj <- .Call(Matrix_expand_pointers, x@p)
320    ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
321    ## but really efficient would be to use only one .Call(.) for uniq(.) !
322    
323    drop0 <- function(x, clx = c(class(x))) {
324        ## FIXME: Csparse_drop should do this (not losing symm./triang.):
325        ## Careful: 'Csparse_drop' also drops triangularity,...
326        ## .Call(Csparse_drop, as_CspClass(x, clx), 0)
327        as_CspClass(.Call(Csparse_drop, as_CspClass(x, clx), 0.),
328                    clx)
329    }
330    
331    uniq <- function(x) {
332        if(is(x, "TsparseMatrix")) uniqTsparse(x) else
333        if(is(x, "sparseMatrix")) drop0(x) else x
334    }
335    
336      } else x      # not 'gT' ; i.e. "uniquely" represented in any case  asTuniq <- function(x) {
337        if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix")
338  }  }
339    
340  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 342 
342  {  {
343      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
344      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
345      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
346      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
347      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
348      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 220  Line 398 
398      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
399  }  }
400    
401    ## -> ./ndenseMatrix.R :
402    n2d_Matrix <- function(from) {
403        stopifnot(is(from, "nMatrix"))
404        fixupDense(new(sub("^n", "d", class(from)),
405                       x = as.double(from@x),
406                       Dim = from@Dim, Dimnames = from@Dimnames),
407                   from)
408        ## FIXME: treat 'factors' smartly {not for triangular!}
409    }
410    n2l_spMatrix <- function(from) {
411        stopifnot(is(from, "nMatrix"))
412        new(sub("^n", "l", class(from)),
413            ##x = as.double(from@x),
414            Dim = from@Dim, Dimnames = from@Dimnames)
415    }
416    
417  if(FALSE)# unused  if(FALSE)# unused
418  l2d_meth <- function(x) {  l2d_meth <- function(x) {
419      cl <- class(x)      cl <- class(x)
420      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
421  }  }
422    
423    ## return "d" or "l" or "n" or "z"
424    .M.kind <- function(x, clx = class(x)) {
425        if(is.matrix(x)) { ## 'old style matrix'
426            if     (is.numeric(x)) "d"
427            else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ??
428            else if(is.complex(x)) "z"
429            else stop("not yet implemented for matrix w/ typeof ", typeof(x))
430        }
431        else if(extends(clx, "dMatrix")) "d"
432        else if(extends(clx, "nMatrix")) "n"
433        else if(extends(clx, "lMatrix")) "l"
434        else if(extends(clx, "zMatrix")) "z"
435        else if(extends(clx, "pMatrix")) "n" # permutation -> pattern
436        else stop(" not yet be implemented for ", clx)
437    }
438    
439    .type.kind <- c("d" = "double",
440                    "l" = "logical",
441                    "n" = "logical",
442                    "z" = "complex")
443    
444    .M.shape <- function(x, clx = class(x)) {
445        if(is.matrix(x)) { ## 'old style matrix'
446            if     (isDiagonal  (x)) "d"
447            else if(isTriangular(x)) "t"
448            else if(isSymmetric (x)) "s"
449            else "g" # general
450        }
451        else if(extends(clx, "diagonalMatrix"))  "d"
452        else if(extends(clx, "triangularMatrix"))"t"
453        else if(extends(clx, "symmetricMatrix")) "s"
454        else "g"
455    }
456    
457    
458    class2 <- function(cl, kind = "l", do.sub = TRUE) {
459        ## Find "corresponding" class; since pos.def. matrices have no pendant:
460        if     (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
461        else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
462        else if(do.sub) sub("^d", kind, cl)
463        else cl
464    }
465    
466    geClass <- function(x) {
467        if     (is(x, "dMatrix")) "dgeMatrix"
468        else if(is(x, "lMatrix")) "lgeMatrix"
469        else if(is(x, "nMatrix")) "ngeMatrix"
470        else if(is(x, "zMatrix")) "zgeMatrix"
471        else stop("general Matrix class not yet implemented for ",
472                  class(x))
473    }
474    
475    .dense.prefixes <- c("d" = "di",
476                         "t" = "tr",
477                         "s" = "sy",
478                         "g" = "ge")
479    
480    .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
481                          "t" = "t",
482                          "s" = "s",
483                          "g" = "g")
484    
485    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
486    as_dense <- function(x) {
487        as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
488    }
489    
490    .sp.class <- function(x) { ## find and return the "sparseness class"
491        if(!is.character(x)) x <- class(x)
492        for(cl in paste(c("C","T","R"), "sparseMatrix", sep=''))
493            if(extends(x, cl))
494                return(cl)
495        ## else (should rarely happen)
496        as.character(NA)
497    }
498    
499    as_Csparse <- function(x) {
500        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
501    }
502    as_Rsparse <- function(x) {
503        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
504    }
505    as_Tsparse <- function(x) {
506        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
507    }
508    
509    as_geClass <- function(x, cl) {
510        if     (extends(cl, "diagonalMatrix")  && isDiagonal(x))
511            as(x, cl)
512        else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x)) {
513            kind <- .M.kind(x)
514            as(x, class2(cl, kind, do.sub= kind != "d"))
515        } else if(extends(cl, "triangularMatrix") && isTriangular(x))
516            as(x, cl)
517        else
518            as(x, paste(.M.kind(x), "geMatrix", sep=''))
519    }
520    
521    as_CspClass <- function(x, cl) {
522        if (## diagonal is *not* sparse:
523            ##(extends(cl, "diagonalMatrix") && isDiagonal(x)) ||
524            (extends(cl, "symmetricMatrix") && isSymmetric(x)) ||
525            (extends(cl, "triangularMatrix")&& isTriangular(x)))
526            as(x, cl)
527        else if(is(x, "CsparseMatrix")) x
528        else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
529    }
530    
531    
532  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
533  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
534      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
535      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
536                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
537                 from)                 from)
538      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 550 
550      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
551  }  }
552    
 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)))  
 }  
553    
554  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
555    isTriMat <- function(object, upper = NA) {
556      ## pretest: is it square?      ## pretest: is it square?
557      d <- dim(object)      d <- dim(object)
558      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
559      ## else slower test      ## else slower test
560      if(!is.matrix(object))      if(!is.matrix(object))
561          object <- as(object,"matrix")          object <- as(object,"matrix")
562      ## == 0 even works for logical & complex:      if(is.na(upper)) {
563      if(upper)          if(all0(object[lower.tri(object)]))
564          all(object[lower.tri(object)] == 0)              structure(TRUE, kind = "U")
565      else          else if(all0(object[upper.tri(object)]))
566          all(object[upper.tri(object)] == 0)              structure(TRUE, kind = "L")
567            else FALSE
568        } else if(upper)
569            all0(object[lower.tri(object)])
570        else ## upper is FALSE
571            all0(object[upper.tri(object)])
572    }
573    
574    ## For Csparse matrices
575    isTriC <- function(x, upper = NA) {
576        ## pretest: is it square?
577        d <- dim(x)
578        if(d[1] != d[2]) return(FALSE)
579        ## else
580        if(d[1] == 0) return(TRUE)
581        ni <- 1:d[2]
582        ## the row indices split according to column:
583        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
584        lil <- unlist(lapply(ilist, length), use.names = FALSE)
585        if(any(lil == 0)) {
586            pos <- lil > 0
587            if(!any(pos)) ## matrix of all 0's
588                return(TRUE)
589            ilist <- ilist[pos]
590            ni <- ni[pos]
591        }
592        ni0 <- ni - 1:1 # '0-based ni'
593        if(is.na(upper)) {
594            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
595                structure(TRUE, kind = "U")
596            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
597                structure(TRUE, kind = "L")
598            else FALSE
599        } else if(upper) {
600            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
601        } else { ## 'lower'
602            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
603        }
604  }  }
605    
606  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
607        ## "matrix" or "denseMatrix" (but not "diagonalMatrix")
608      d <- dim(object)      d <- dim(object)
609      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
610      else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)      else if(is.matrix(object))
611            ## requires that "vector-indexing" works for 'object' :
612            all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
613        else ## "denseMatrix" -- packed or unpacked
614            if(is(object, "generalMatrix")) # "dge", "lge", ...
615                all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
616            else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:
617                ## -> has 'uplo'  differentiate between packed and unpacked
618    
619    ### .......... FIXME ...............
620    
621                packed <- isPacked(object)
622                if(object@uplo == "U") {
623                } else { ## uplo == "L"
624  }  }
625    
626    ### very cheap workaround
627                all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
628            }
629    }
630    
631    
632    diagU2N <- function(x)
633    {
634        ## Purpose: Transform a *unit diagonal* sparse triangular matrix
635        ##  into one with explicit diagonal entries '1'
636        if(is(x, "CsparseMatrix"))
637            return(.Call(Csparse_diagU2N, x))
638        ## else
639    
640        ## FIXME! -- for "ltT", "ntT", ... :
641        xT <- as(x, "dgTMatrix")
642        ## leave it as  T* - the caller can always coerce to C* if needed:
643        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
644            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
645    }
646    
647    ## Needed, e.g., in ./Csparse.R for colSums() etc:
648    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
649        x <- as(x, "dgCMatrix")
650        callGeneric()
651    }
652    
653    .as.dgC.0.factors <- function(x) {
654        if(!is(x, "dgCMatrix"))
655            as(x, "dgCMatrix") # will not have 'factors'
656        else ## dgCMatrix
657            if(!length(x@factors)) x else { x@factors <- list() ; x }
658    }
659    
660    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
661        ## used e.g. inside colSums() etc methods
662        x <- as(x, "dgTMatrix")
663        callGeneric()
664    }
665    
666    
667    ### Fast much simplified version of tapply()
668    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
669        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
670    }
671    
672    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
673    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
674    ## }
675    

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

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge