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 1349, Mon Aug 7 09:05:42 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    .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
353                          "t" = "t",
354                          "s" = "s",
355                          "g" = "g")
356    
357    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
358    as_dense <- function(x) {
359        as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
360    }
361    
362    as_Csparse <- function(x) {
363        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
364    }
365    as_Rsparse <- function(x) {
366        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
367    }
368    as_Tsparse <- function(x) {
369        as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
370    }
371    
372    as_geClass <- function(x, cl) {
373        if     (extends(cl, "diagonalMatrix")  && isDiagonal(x))
374            as(x, cl)
375        else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x))
376            as(x, class2(cl, kind, do.sub= kind != "d"))
377        else if(extends(cl, "triangularMatrix") && isTriangular(x))
378            as(x, cl)
379        else
380            as(x, paste(.M.kind(x), "geMatrix", sep=''))
381    }
382    
383    as_CspClass <- function(x, cl) {
384        if ((extends(cl, "diagonalMatrix")  && isDiagonal(x)) ||
385            (extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||
386            (extends(cl, "triangularMatrix")&& isTriangular(x)))
387            as(x, cl)
388        else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
389    }
390    
391    
392  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
393  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
394      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
395      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
396                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
397                 from)                 from)
398      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 410 
410      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
411  }  }
412    
 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)))  
 }  
413    
414  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
415    isTriMat <- function(object, upper = NA) {
416      ## pretest: is it square?      ## pretest: is it square?
417      d <- dim(object)      d <- dim(object)
418      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 420 
420      if(!is.matrix(object))      if(!is.matrix(object))
421          object <- as(object,"matrix")          object <- as(object,"matrix")
422      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
423      if(upper)      if(is.na(upper)) {
424            if(all(object[lower.tri(object)] == 0))
425                structure(TRUE, kind = "U")
426            else if(all(object[upper.tri(object)] == 0))
427                structure(TRUE, kind = "L")
428            else FALSE
429        } else if(upper)
430          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
431      else      else ## upper is FALSE
432          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
433  }  }
434    
435    ## For Csparse matrices
436    isTriC <- function(x, upper = NA) {
437        ## pretest: is it square?
438        d <- dim(x)
439        if(d[1] != d[2]) return(FALSE)
440        ## else
441        if(d[1] == 0) return(TRUE)
442        ni <- 1:d[2]
443        ## the row indices split according to column:
444        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
445        lil <- unlist(lapply(ilist, length), use.names = FALSE)
446        if(any(lil == 0)) {
447            pos <- lil > 0
448            if(!any(pos)) ## matrix of all 0's
449                return(TRUE)
450            ilist <- ilist[pos]
451            ni <- ni[pos]
452        }
453        ni0 <- ni - 1:1 # '0-based ni'
454        if(is.na(upper)) {
455            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
456                structure(TRUE, kind = "U")
457            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
458                structure(TRUE, kind = "L")
459            else FALSE
460        } else if(upper) {
461            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
462        } else { ## 'lower'
463            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
464        }
465    }
466    
467  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
468      d <- dim(object)      d <- dim(object)
469      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
470      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)
471  }  }
472    
473    diagU2N <- function(x)
474    {
475        ## Purpose: Transform a *unit diagonal* triangular matrix
476        ##  into one with explicit diagonal entries '1'
477        xT <- as(x, "dgTMatrix")
478        ## leave it as  T* - the caller can always coerce to C* if needed:
479        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
480            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
481    }
482    
483    ## FIXME: this should probably be dropped / replaced by as_Csparse
484    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
485        x <- as(x, "dgCMatrix")
486        callGeneric()
487    }
488    
489    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
490        ## used e.g. inside colSums() etc methods
491        x <- as(x, "dgTMatrix")
492        callGeneric()
493    }
494    
495    
496    ### Fast much simplified version of tapply()
497    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
498        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
499    }
500    
501    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
502    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
503    ## }
504    

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

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