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 1332, Thu Jul 27 13:42:18 2006 UTC
# Line 5  Line 5 
5  .M.v <- function(x, y) callGeneric(x, as.matrix(y))  .M.v <- function(x, y) callGeneric(x, as.matrix(y))
6  .v.M <- function(x, y) callGeneric(rbind(x), y)  .v.M <- function(x, y) callGeneric(rbind(x), y)
7    
8    .M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL)
9    
10  .has.DN <- ## has non-trivial Dimnames slot?  .has.DN <- ## has non-trivial Dimnames slot?
11      function(x) !identical(list(NULL,NULL), x@Dimnames)      function(x) !identical(list(NULL,NULL), x@Dimnames)
12    
# Line 82  Line 84 
84      da[2]      da[2]
85  }  }
86    
87    ## Note: !isPacked(.)  i.e. `full' still contains
88    ## ----  "*sy" and "*tr" which have "undefined" lower or upper part
89    isPacked <- function(x)
90    {
91        ## Is 'x' a packed (dense) matrix ?
92        is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
93    }
94    
95  emptyColnames <- function(x)  emptyColnames <- function(x)
96  {  {
97      ## Useful for compact printing of (parts) of sparse matrices      ## Useful for compact printing of (parts) of sparse matrices
# Line 134  Line 144 
144  non0ind <- function(x) {  non0ind <- function(x) {
145      if(is.numeric(x))      if(is.numeric(x))
146          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))
147        ## else
     ## else return a (i,j) matrix of non-zero-indices  
   
148      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
149      if(is(x, "TsparseMatrix"))      ## return a 2-column (i,j) matrix of
150          return(unique(cbind(x@i,x@j)))      ## 0-based indices of non-zero entries  :
151        non0.i <- function(M) {
152            if(is(M, "TsparseMatrix"))
153                return(unique(cbind(M@i,M@j)))
154            if(is(M, "pMatrix"))
155                return(cbind(seq(length=nrow(M)), M@perm) - 1:1)
156            ## else:
157            isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
158            .Call(compressed_non_0_ij, M, isC)
159        }
160    
161        if(is(x, "symmetricMatrix")) { # also get "other" triangle
162            ij <- non0.i(x)
163            notdiag <- ij[,1] != ij[,2]# but not the diagonals again
164            rbind(ij, ij[notdiag, 2:1])
165        }
166        else
167            non0.i(x)
168    }
169    
170    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
171    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
172    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
173    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
174    
175      isCol <- function(M) any("i" == slotNames(M))  complementInd <- function(ij, dim)
176      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  {
177        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
178        ##          but as 1-based indices
179        n <- prod(dim)
180        if(n == 0) return(integer(0))
181        ii <- 1:n
182        ii[-(1 + encodeInd(ij, nr = dim[1]))]
183    }
184    
185    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
186    
187    intersectInd <- function(ij1, ij2, nrow) {
188        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
189        ## return only the *common* entries
190        decodeInd(intersect(encodeInd(ij1, nrow),
191                            encodeInd(ij2, nrow)), nrow)
192    }
193    
194    WhichintersectInd <- function(ij1, ij2, nrow) {
195        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
196        ## find *where*  common entries are in ij1 & ij2
197        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
198        ni <- !is.na(m1)
199        list(which(ni), m1[ni])
200  }  }
201    
202    
203  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
204  uniq <- function(x) {  
205      if(is(x, "TsparseMatrix")) {  uniqTsparse <- function(x, class.x = c(class(x))) {
206          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
207          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
208          ## -----------------------------------------------------------          ## -----------------------------------------------------------
209          ## The following is *not* efficient {but easy to program}:      ## The following is not quite efficient {but easy to program,
210          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")      ## and as() are based on C code  (all of them?)
211          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")      ##
212          else stop("not implemented for class", class(x))      ## FIXME: Do it fast for the case where 'x' is already 'uniq'
213    
214        switch(class.x,
215               "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
216               "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
217               "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
218               ## do we need this for "logical" ones, there's no sum() there!
219               "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
220               "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
221               "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
222               ## otherwise:
223               stop("not yet implemented for class ", clx))
224    }
225    
226    ## Note: maybe, using
227    ## ----    xj <- .Call(Matrix_expand_pointers, x@p)
228    ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
229    ## but really efficient would be to use only one .Call(.) for uniq(.) !
230    
231      } else x      # not 'gT' ; i.e. "uniquely" represented in any case  uniq <- function(x) {
232        if(is(x, "TsparseMatrix")) uniqTsparse(x) else x
233        ## else:  not 'Tsparse', i.e. "uniquely" represented in any case
234  }  }
235    
236  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 238 
238  {  {
239      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
240      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
241      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
242      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
243      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
244      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 300 
300      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
301  }  }
302    
303    ## return "d" or "l" or "z"
304    .M.kind <- function(x, clx = class(x)) {
305        if(is.matrix(x)) { ## 'old style matrix'
306            if     (is.numeric(x)) "d"
307            else if(is.logical(x)) "l"
308            else if(is.complex(x)) "z"
309            else stop("not yet implemented for matrix w/ typeof ", typeof(x))
310        }
311        else if(extends(clx, "dMatrix")) "d"
312        else if(extends(clx, "lMatrix")) "l"
313        else if(extends(clx, "zMatrix")) "z"
314        else stop(" not yet be implemented for ", clx)
315    }
316    
317    .M.shape <- function(x, clx = class(x)) {
318        if(is.matrix(x)) { ## 'old style matrix'
319            if     (isDiagonal  (x)) "d"
320            else if(isTriangular(x)) "t"
321            else if(isSymmetric (x)) "s"
322            else "g" # general
323        }
324        else if(extends(clx, "diagonalMatrix"))  "d"
325        else if(extends(clx, "triangularMatrix"))"t"
326        else if(extends(clx, "symmetricMatrix")) "s"
327        else "g"
328    }
329    
330    
331    class2 <- function(cl, kind = "l", do.sub = TRUE) {
332        ## Find "corresponding" class; since pos.def. matrices have no pendant:
333        if     (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
334        else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
335        else if(do.sub) sub("^d", kind, cl)
336        else cl
337    }
338    
339    geClass <- function(x) {
340        if     (is(x, "dMatrix")) "dgeMatrix"
341        else if(is(x, "lMatrix")) "lgeMatrix"
342        else if(is(x, "zMatrix")) "zgeMatrix"
343        else stop("general Matrix class not yet implemented for ",
344                  class(x))
345    }
346    
347    .dense.prefixes <- c("d" = "di",
348                         "t" = "tr",
349                         "s" = "sy",
350                         "g" = "ge")
351    
352    .Csparse.prefix <- function(ch) {
353        switch(ch,
354               "d" =, "t" = "tC",
355               "s" = "sC",
356               "g" = "gC",
357               stop("invalid Matrix shape: ", ch))
358    }
359    
360    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
361    as_dense <- function(x) {
362        as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
363    }
364    
365    as_Csparse <- function(x) {
366        as(x, paste(.M.kind(x), .Csparse.prefix(.M.shape(x)), "Matrix", sep=''))
367    }
368    
369    as_geClass <- function(x, cl) {
370        if     (extends(cl, "diagonalMatrix")  && isDiagonal(x))
371            as(x, cl)
372        else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x))
373            as(x, class2(cl, kind, do.sub= kind != "d"))
374        else if(extends(cl, "triangularMatrix") && isTriangular(x))
375            as(x, cl)
376        else
377            as(x, paste(.M.kind(x), "geMatrix", sep=''))
378    }
379    
380    as_CspClass <- function(x, cl) {
381        if ((extends(cl, "diagonalMatrix")  && isDiagonal(x)) ||
382            (extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||
383            (extends(cl, "triangularMatrix")&& isTriangular(x)))
384            as(x, cl)
385        else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
386    }
387    
388    
389  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
390  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
391      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
392      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
393                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
394                 from)                 from)
395      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 407 
407      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
408  }  }
409    
 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)))  
 }  
410    
411  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
412    isTriMat <- function(object, upper = NA) {
413      ## pretest: is it square?      ## pretest: is it square?
414      d <- dim(object)      d <- dim(object)
415      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 417 
417      if(!is.matrix(object))      if(!is.matrix(object))
418          object <- as(object,"matrix")          object <- as(object,"matrix")
419      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
420      if(upper)      if(is.na(upper)) {
421            if(all(object[lower.tri(object)] == 0))
422                structure(TRUE, kind = "U")
423            else if(all(object[upper.tri(object)] == 0))
424                structure(TRUE, kind = "L")
425            else FALSE
426        } else if(upper)
427          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
428      else      else ## upper is FALSE
429          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
430  }  }
431    
432    ## For Csparse matrices
433    isTriC <- function(x, upper = NA) {
434        ## pretest: is it square?
435        d <- dim(x)
436        if(d[1] != d[2]) return(FALSE)
437        ## else
438        if(d[1] == 0) return(TRUE)
439        ni <- 1:d[2]
440        ## the row indices split according to column:
441        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
442        lil <- unlist(lapply(ilist, length), use.names = FALSE)
443        if(any(lil == 0)) {
444            pos <- lil > 0
445            if(!any(pos)) ## matrix of all 0's
446                return(TRUE)
447            ilist <- ilist[pos]
448            ni <- ni[pos]
449        }
450        ni0 <- ni - 1:1 # '0-based ni'
451        if(is.na(upper)) {
452            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
453                structure(TRUE, kind = "U")
454            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
455                structure(TRUE, kind = "L")
456            else FALSE
457        } else if(upper) {
458            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
459        } else { ## 'lower'
460            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
461        }
462    }
463    
464  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
465      d <- dim(object)      d <- dim(object)
466      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
467      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)
468  }  }
469    
470    diagU2N <- function(x)
471    {
472        ## Purpose: Transform a *unit diagonal* triangular matrix
473        ##  into one with explicit diagonal entries '1'
474        xT <- as(x, "dgTMatrix")
475        ## leave it as  T* - the caller can always coerce to C* if needed:
476        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
477            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
478    }
479    
480    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
481        x <- as(x, "dgCMatrix")
482        callGeneric()
483    }
484    
485    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
486        x <- as(x, "dgTMatrix")
487        callGeneric()
488    }
489    
490    
491    ### Fast much simplified version of tapply()
492    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
493        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
494    }
495    
496    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
497    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
498    ## }
499    

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

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