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 1329, Fri Jul 21 06:47:06 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      if(is(x, "TsparseMatrix"))      ## return a 2-column (i,j) matrix of
148          return(unique(cbind(x@i,x@j)))      ## 0-based indices of non-zero entries  :
149        non0.i <- function(M) {
150            if(is(M, "TsparseMatrix"))
151                return(unique(cbind(M@i,M@j)))
152            if(is(M, "pMatrix"))
153                return(cbind(seq(length=nrow(M)), M@perm) - 1:1)
154            ## else:
155            isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
156            .Call(compressed_non_0_ij, M, isC)
157        }
158    
159        if(is(x, "symmetricMatrix")) { # also get "other" triangle
160            ij <- non0.i(x)
161            notdiag <- ij[,1] != ij[,2]# but not the diagonals again
162            rbind(ij, ij[notdiag, 2:1])
163        }
164        else
165            non0.i(x)
166    }
167    
168      isCol <- function(M) any("i" == slotNames(M))  ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
169      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
170    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
171    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
172    
173    complementInd <- function(ij, dim)
174    {
175        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
176        ##          but as 1-based indices
177        n <- prod(dim)
178        if(n == 0) return(integer(0))
179        ii <- 1:n
180        ii[-(1 + encodeInd(ij, nr = dim[1]))]
181    }
182    
183    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
184    
185    intersectInd <- function(ij1, ij2, nrow) {
186        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
187        ## return only the *common* entries
188        decodeInd(intersect(encodeInd(ij1, nrow),
189                            encodeInd(ij2, nrow)), nrow)
190    }
191    
192    WhichintersectInd <- function(ij1, ij2, nrow) {
193        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
194        ## find *where*  common entries are in ij1 & ij2
195        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
196        ni <- !is.na(m1)
197        list(which(ni), m1[ni])
198  }  }
199    
200    
201  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
202  uniq <- function(x) {  uniq <- function(x) {
203    
204    ### Note: maybe, using
205    ### ----    xj <- .Call(Matrix_expand_pointers, x@p)
206    ### would be slightly more efficient than as( <dgC> , "dgTMatrix")
207    ### but really efficient would be to use only one .Call(.) for uniq(.) !
208    ### Try to do it particularly fast for the case where 'x' is already a 'uniq' <dgT>
209    
210      if(is(x, "TsparseMatrix")) {      if(is(x, "TsparseMatrix")) {
211          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
212          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
213          ## -----------------------------------------------------------          ## -----------------------------------------------------------
214          ## The following is *not* efficient {but easy to program}:          ## The following is not quite efficient {but easy to program,
215            ## and both as() are based on C code
216          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
217          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
218          else stop("not implemented for class", class(x))          else stop("not implemented for class", class(x))
219    
220      } else x      # not 'gT' ; i.e. "uniquely" represented in any case      } else x  ## not 'gT' ; i.e. "uniquely" represented in any case
221  }  }
222    
223  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 225 
225  {  {
226      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
227      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
228      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
229      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
230      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
231      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 287 
287      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
288  }  }
289    
290    class2 <- function(cl, kind = "l", do.sub = TRUE) {
291        ## Find "corresponding" class; since pos.def. matrices have no pendant:
292        if     (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
293        else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
294        else if(do.sub) sub("^d", kind, cl)
295        else cl
296    }
297    
298    geClass <- function(x) {
299        if(is(x, "dMatrix")) "dgeMatrix"
300        else if(is(x, "lMatrix")) "lgeMatrix"
301        else stop("general Matrix class not yet implemented for ",
302                  class(x))
303    }
304    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
305    as_geClass <- function(x, cl) {
306        clx <- class(x)
307        kind <-
308            if(is.matrix(x)) { # 'old style matrix'
309                if(is.numeric(x)) "d" else if(is.logical(x)) "l"
310                else stop("general Matrix class not implemented for matrix w/ typeof ",
311                          typeof(x))
312            }
313            else if(extends(clx, "dMatrix")) "d"
314            else if(extends(clx, "lMatrix")) "l"
315            else
316                stop("general Matrix class not yet implemented for ", clx)
317    
318        clDia <- extends(cl, "diagonalMatrix")
319        clSym <- extends(cl, "symmetricMatrix")
320        clTri <- extends(cl, "triangularMatrix")
321        if     (clDia && .is.diagonal(x)) as(x, cl)
322        else if(clSym &&  isSymmetric(x)) as(x, class2(cl, kind, do.sub= kind != "d"))
323        else if(clTri && isTriangular(x)) as(x, cl)
324        else as(x, paste(kind, "geMatrix", sep=''))
325    }
326    
327  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
328  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
329      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
330      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
331                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
332                 from)                 from)
333      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 345 
345      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
346  }  }
347    
 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)))  
 }  
348    
349  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
350    isTriMat <- function(object, upper = NA) {
351      ## pretest: is it square?      ## pretest: is it square?
352      d <- dim(object)      d <- dim(object)
353      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 355 
355      if(!is.matrix(object))      if(!is.matrix(object))
356          object <- as(object,"matrix")          object <- as(object,"matrix")
357      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
358      if(upper)      if(is.na(upper)) {
359            if(all(object[lower.tri(object)] == 0))
360                structure(TRUE, kind = "U")
361            else if(all(object[upper.tri(object)] == 0))
362                structure(TRUE, kind = "L")
363            else FALSE
364        } else if(upper)
365          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
366      else      else ## upper is FALSE
367          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
368  }  }
369    
370    ## For Csparse matrices
371    isTriC <- function(x, upper = NA) {
372        ## pretest: is it square?
373        d <- dim(x)
374        if(d[1] != d[2]) return(FALSE)
375        ## else
376        if(d[1] == 0) return(TRUE)
377        ni <- 1:d[2]
378        ## the row indices split according to column:
379        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
380        lil <- unlist(lapply(ilist, length), use.names = FALSE)
381        if(any(lil == 0)) {
382            pos <- lil > 0
383            if(!any(pos)) ## matrix of all 0's
384                return(TRUE)
385            ilist <- ilist[pos]
386            ni <- ni[pos]
387        }
388        ni0 <- ni - 1:1 # '0-based ni'
389        if(is.na(upper)) {
390            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
391                structure(TRUE, kind = "U")
392            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
393                structure(TRUE, kind = "L")
394            else FALSE
395        } else if(upper) {
396            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
397        } else { ## 'lower'
398            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
399        }
400    }
401    
402  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
403      d <- dim(object)      d <- dim(object)
404      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
405      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)
406  }  }
407    
408    diagU2N <- function(x)
409    {
410        ## Purpose: Transform a *unit diagonal* triangular matrix
411        ##  into one with explicit diagonal entries '1'
412        xT <- as(x, "dgTMatrix")
413        ## leave it as  T* - the caller can always coerce to C* if needed:
414        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
415            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
416    }
417    
418    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
419        x <- as(x, "dgCMatrix")
420        callGeneric()
421    }
422    
423    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
424        x <- as(x, "dgTMatrix")
425        callGeneric()
426    }
427    
428    
429    ### Fast much simplified version of tapply()
430    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
431        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
432    }
433    
434    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
435    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
436    ## }
437    

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

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