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 1571, Sat Sep 16 21:03:12 2006 UTC revision 1654, Fri Oct 27 16:58:15 2006 UTC
# Line 5  Line 5 
5    
6  ## Need to consider NAs ;  "== 0" even works for logical & complex:  ## Need to consider NAs ;  "== 0" even works for logical & complex:
7  is0  <- function(x) !is.na(x) & x == 0  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)  all0 <- function(x) !any(is.na(x)) && all(x == 0)
10    
11  allTrue  <- function(x) !any(is.na(x)) && all(x)  allTrue  <- function(x) !any(is.na(x)) && all(x)
# Line 21  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    
# Line 120  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 165  Line 184 
184      invisible(x)# as print() S3 methods do      invisible(x)# as print() S3 methods do
185  }  }
186    
187  ### FIXME? -- make this into a generic function (?)  nonFALSE <- function(x) {
188  nnzero <- function(x) {      ## 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)      cl <- class(x)
211      if(!extends(cl, "Matrix"))      if(!extends(cl, "Matrix"))
212          sum(x != 0)          sum(nz.NA(x, na.counted))
213      else if(extends(cl, "sparseMatrix"))      else if(extends(cl, "sparseMatrix"))
214          ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!          ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!
215         switch(.sp.class(cl),         switch(.sp.class(cl),
# Line 177  Line 217 
217                 "TsparseMatrix" = length(x@i),                 "TsparseMatrix" = length(x@i),
218                 "RsparseMatrix" = length(x@j))                 "RsparseMatrix" = length(x@j))
219      else ## denseMatrix      else ## denseMatrix
220          sum(as_geClass(x)@x != 0)          sum(nz.NA(as_geClass(x, cl)@x, na.counted))
221  }  }
222    
223  ## For sparseness handling  ## For sparseness handling
# Line 186  Line 226 
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"))
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 206  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 279  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) {  asTuniq <- function(x) {
# Line 465  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 573  Line 623 
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.1571  
changed lines
  Added in v.1654

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