# SCM Repository

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

# Diff of /pkg/R/Auxiliaries.R

revision 1390, Fri Aug 18 15:41:36 2006 UTC revision 1654, Fri Oct 27 16:58:15 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)
# Line 11  Line 22
22      function(x) !identical(list(NULL,NULL), x@Dimnames)      function(x) !identical(list(NULL,NULL), x@Dimnames)
23
24  .bail.out.1 <- function(fun, cl) {  .bail.out.1 <- function(fun, cl) {
25      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),
26           call. = FALSE)           call. = FALSE)
27  }  }
28  .bail.out.2 <- function(fun, cl1, cl2) {  .bail.out.2 <- function(fun, cl1, cl2) {
29      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.',
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 88  Line 108
108  ## ----  "*sy" and "*tr" which have "undefined" lower or upper part  ## ----  "*sy" and "*tr" which have "undefined" lower or upper part
109  isPacked <- function(x)  isPacked <- function(x)
110  {  {
111      ## Is 'x' a packed (dense) matrix ? -- gives also TRUE for sparse      ## Is 'x' a packed (dense) matrix ?
112      is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))      is(x, "denseMatrix") &&
113        any("x" == slotNames(x)) && length(x@x) < prod(dim(x))
114  }  }
115
116  emptyColnames <- function(x)  emptyColnames <- function(x)
# Line 100  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"),                       maxp = getOption("max.print"),
144                       justify = "none", right = TRUE)                       justify = "none", right = TRUE)
# Line 113  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  ## TODO (in R): print() should accept 'maxp'      if(.isR_24)
156      print(cf, quote = FALSE, right = right)           print(cf, quote = FALSE, right = right, max = maxp)
157        else print(cf, quote = FALSE, right = right)
158      invisible(x)      invisible(x)
159  }  }
160
# Line 125  Line 165
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))
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 141  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) {
227
228      if(is.numeric(x))      if(is.numeric(x))
229          return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))          return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0))
230      ## else      ## else
231      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
## return a 2-column (i,j) matrix of
## 0-based indices of non-zero entries  :
232      non0.i <- function(M) {      non0.i <- function(M) {
233          if(is(M, "TsparseMatrix"))          if(is(M, "TsparseMatrix"))
234              return(unique(cbind(M@i,M@j)))              return(unique(cbind(M@i,M@j)))
235          if(is(M, "pMatrix"))          if(is(M, "pMatrix"))
236              return(cbind(seq(length=nrow(M)), M@perm) - 1:1)              return(cbind(seq_len(nrow(M)), M@perm) - 1:1)
237          ## else:          ## else:
238          isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)          isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
239          .Call(compressed_non_0_ij, M, isC)          .Call(compressed_non_0_ij, M, isC)
# Line 166  Line 246
246      }      }
247      else if(is(x, "triangularMatrix")) { # check for "U" diag      else if(is(x, "triangularMatrix")) { # check for "U" diag
248          if(x@diag == "U") {          if(x@diag == "U") {
249              i <- seq(length = dim(x)[1]) - 1:1              i <- seq_len(dim(x)[1]) - 1:1
250              rbind(non0.i(x), cbind(i,i))              rbind(non0.i(x), cbind(i,i))
251          } else non0.i(x)          } else non0.i(x)
252      }      }
# Line 226  Line 306
306             "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),             "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
307             "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),             "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
308             "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),             "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:             ## otherwise:
314             stop("not yet implemented for class ", clx))             stop("not yet implemented for class ", class.x))
315  }  }
316
317  ## Note: maybe, using  ## Note: maybe, using
# Line 235  Line 319
319  ## would be slightly more efficient than as( <dgC> , "dgTMatrix")  ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
320  ## but really efficient would be to use only one .Call(.) for uniq(.) !  ## but really efficient would be to use only one .Call(.) for uniq(.) !
321
322    drop0 <- function(x, clx = c(class(x))) {
323        ## FIXME: Csparse_drop should do this (not losing symm./triang.):
324        ## Careful: 'Csparse_drop' also drops triangularity,...
325        ## .Call(Csparse_drop, as_CspClass(x, clx), 0)
326        as_CspClass(.Call(Csparse_drop, as_CspClass(x, clx), 0.),
327                    clx)
328    }
329
330  uniq <- function(x) {  uniq <- function(x) {
331      if(is(x, "TsparseMatrix")) uniqTsparse(x) else x      if(is(x, "TsparseMatrix")) uniqTsparse(x) else
332      ## else:  not 'Tsparse', i.e. "uniquely" represented in any case      if(is(x, "sparseMatrix")) drop0(x) else x
333    }
334
335    asTuniq <- function(x) {
336        if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix")
337  }  }
338
339  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 301  Line 397
397      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
398  }  }
399
400    ## -> ./ndenseMatrix.R :
401    n2d_Matrix <- function(from) {
402        stopifnot(is(from, "nMatrix"))
403        fixupDense(new(sub("^n", "d", class(from)),
404                       x = as.double(from@x),
405                       Dim = from@Dim, Dimnames = from@Dimnames),
406                   from)
407        ## FIXME: treat 'factors' smartly {not for triangular!}
408    }
409    n2l_spMatrix <- function(from) {
410        stopifnot(is(from, "nMatrix"))
411        new(sub("^n", "l", class(from)),
412            ##x = as.double(from@x),
413            Dim = from@Dim, Dimnames = from@Dimnames)
414    }
415
416  if(FALSE)# unused  if(FALSE)# unused
417  l2d_meth <- function(x) {  l2d_meth <- function(x) {
418      cl <- class(x)      cl <- class(x)
419      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
420  }  }
421
422  ## return "d" or "l" or "z"  ## return "d" or "l" or "n" or "z"
423  .M.kind <- function(x, clx = class(x)) {  .M.kind <- function(x, clx = class(x)) {
424      if(is.matrix(x)) { ## 'old style matrix'      if(is.matrix(x)) { ## 'old style matrix'
425          if     (is.numeric(x)) "d"          if     (is.numeric(x)) "d"
426          else if(is.logical(x)) "l"          else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ??
427          else if(is.complex(x)) "z"          else if(is.complex(x)) "z"
428          else stop("not yet implemented for matrix w/ typeof ", typeof(x))          else stop("not yet implemented for matrix w/ typeof ", typeof(x))
429      }      }
430      else if(extends(clx, "dMatrix")) "d"      else if(extends(clx, "dMatrix")) "d"
431        else if(extends(clx, "nMatrix")) "n"
432      else if(extends(clx, "lMatrix")) "l"      else if(extends(clx, "lMatrix")) "l"
433      else if(extends(clx, "zMatrix")) "z"      else if(extends(clx, "zMatrix")) "z"
434        else if(extends(clx, "pMatrix")) "n" # permutation -> pattern
435      else stop(" not yet be implemented for ", clx)      else stop(" not yet be implemented for ", clx)
436  }  }
437
# Line 346  Line 460
460  geClass <- function(x) {  geClass <- function(x) {
461      if     (is(x, "dMatrix")) "dgeMatrix"      if     (is(x, "dMatrix")) "dgeMatrix"
462      else if(is(x, "lMatrix")) "lgeMatrix"      else if(is(x, "lMatrix")) "lgeMatrix"
463        else if(is(x, "nMatrix")) "ngeMatrix"
464      else if(is(x, "zMatrix")) "zgeMatrix"      else if(is(x, "zMatrix")) "zgeMatrix"
465      else stop("general Matrix class not yet implemented for ",      else stop("general Matrix class not yet implemented for ",
466                class(x))                class(x))
# Line 366  Line 481
481      as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))      as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
482  }  }
483
484    .sp.class <- function(x) { ## find and return the "sparseness class"
485        if(!is.character(x)) x <- class(x)
486        for(cl in paste(c("C","T","R"), "sparseMatrix", sep=''))
487            if(extends(x, cl))
488                return(cl)
489        ## else (should rarely happen)
490        as.character(NA)
491    }
492
493  as_Csparse <- function(x) {  as_Csparse <- function(x) {
494      as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))      as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
495  }  }
# Line 389  Line 513
513  }  }
514
515  as_CspClass <- function(x, cl) {  as_CspClass <- function(x, cl) {
516      if ((extends(cl, "diagonalMatrix")  && isDiagonal(x)) ||      if (## diagonal is *not* sparse:
517            ##(extends(cl, "diagonalMatrix") && isDiagonal(x)) ||
518          (extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||          (extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||
519          (extends(cl, "triangularMatrix")&& isTriangular(x)))          (extends(cl, "triangularMatrix")&& isTriangular(x)))
520          as(x, cl)          as(x, cl)
521        else if(is(x, "CsparseMatrix")) x
522      else as(x, paste(.M.kind(x), "gCMatrix", sep=''))      else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
523  }  }
524
# Line 427  Line 553
553      ## else slower test      ## else slower test
554      if(!is.matrix(object))      if(!is.matrix(object))
555          object <- as(object,"matrix")          object <- as(object,"matrix")
## == 0 even works for logical & complex:
556      if(is.na(upper)) {      if(is.na(upper)) {
557          if(all(object[lower.tri(object)] == 0))          if(all0(object[lower.tri(object)]))
558              structure(TRUE, kind = "U")              structure(TRUE, kind = "U")
559          else if(all(object[upper.tri(object)] == 0))          else if(all0(object[upper.tri(object)]))
560              structure(TRUE, kind = "L")              structure(TRUE, kind = "L")
561          else FALSE          else FALSE
562      } else if(upper)      } else if(upper)
563          all(object[lower.tri(object)] == 0)          all0(object[lower.tri(object)])
564      else ## upper is FALSE      else ## upper is FALSE
565          all(object[upper.tri(object)] == 0)          all0(object[upper.tri(object)])
566  }  }
567
568  ## For Csparse matrices  ## For Csparse matrices
# Line 478  Line 603
603      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
604      else if(is.matrix(object))      else if(is.matrix(object))
605          ## requires that "vector-indexing" works for 'object' :          ## requires that "vector-indexing" works for 'object' :
606          all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)          all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
607      else ## "denseMatrix" -- packed or unpacked      else ## "denseMatrix" -- packed or unpacked
608          if(is(object, "generalMatrix")) # "dge", "lge", ...          if(is(object, "generalMatrix")) # "dge", "lge", ...
609              all(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)              all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
610          else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:          else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:
611              ## -> has 'uplo'  differentiate between packed and unpacked              ## -> has 'uplo'  differentiate between packed and unpacked
612
# Line 493  Line 618
618              }              }
619
620  ### very cheap workaround  ### very cheap workaround
621              all(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]              all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
== 0)
622          }          }
623  }  }
624
625
626    ## FIXME? -- this should also work for "ltT", "ntT", ... :
627  diagU2N <- function(x)  diagU2N <- function(x)
628  {  {
629      ## Purpose: Transform a *unit diagonal* triangular matrix      ## Purpose: Transform a *unit diagonal* sparse triangular matrix
630      ##  into one with explicit diagonal entries '1'      ##  into one with explicit diagonal entries '1'
631      xT <- as(x, "dgTMatrix")      xT <- as(x, "dgTMatrix")
632      ## leave it as  T* - the caller can always coerce to C* if needed:      ## leave it as  T* - the caller can always coerce to C* if needed:

Legend:
 Removed from v.1390 changed lines Added in v.1654