# SCM Repository

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

# Diff of /pkg/R/Auxiliaries.R

revision 1198, Mon Jan 23 15:01:02 2006 UTC revision 1331, Sat Jul 22 17:59:53 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) {
203      if(is(x, "TsparseMatrix")) {  uniqTsparse <- function(x, class.x = c(class(x))) {
204          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
205          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
206          ## -----------------------------------------------------------          ## -----------------------------------------------------------
207          ## The following is *not* efficient {but easy to program}:      ## The following is not quite efficient {but easy to program,
208          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")      ## and as() are based on C code  (all of them?)
209          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")      ##
210          else stop("not implemented for class", class(x))      ## FIXME: Do it fast for the case where 'x' is already 'uniq'
211
212        switch(class.x,
213               "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
214               "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
215               "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
216               ## do we need this for "logical" ones, there's no sum() there!
217               "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
218               "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
219               "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
220               ## otherwise:
221               stop("not yet implemented for class ", clx))
222    }
223
224    ## Note: maybe, using
225    ## ----    xj <- .Call(Matrix_expand_pointers, x@p)
226    ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
227    ## but really efficient would be to use only one .Call(.) for uniq(.) !
228
229      } else x      # not 'gT' ; i.e. "uniquely" represented in any case  uniq <- function(x) {
230        if(is(x, "TsparseMatrix")) uniqTsparse(x) else x
231        ## else:  not 'Tsparse', i.e. "uniquely" represented in any case
232  }  }
233
234  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 236
236  {  {
237      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
238      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
239      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
240      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
241      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
242      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 298
298      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
299  }  }
300
301    ## return "d" or "l" or "z"
302    .M.kind <- function(x, clx = class(x)) {
303        if(is.matrix(x)) { ## 'old style matrix'
304            if     (is.numeric(x)) "d"
305            else if(is.logical(x)) "l"
306            else if(is.complex(x)) "z"
307            else stop("not yet implemented for matrix w/ typeof ", typeof(x))
308        }
309        else if(extends(clx, "dMatrix")) "d"
310        else if(extends(clx, "lMatrix")) "l"
311        else if(extends(clx, "zMatrix")) "z"
312        else stop(" not yet be implemented for ", clx)
313    }
314
315    .M.shape <- function(x, clx = class(x)) {
316        if(is.matrix(x)) { ## 'old style matrix'
317            if     (isDiagonal  (x)) "d"
318            else if(isTriangular(x)) "t"
319            else if(isSymmetric (x)) "s"
320            else "g" # general
321        }
322        else if(extends(clx, "diagonalMatrix"))  "d"
323        else if(extends(clx, "triangularMatrix"))"t"
324        else if(extends(clx, "symmetricMatrix")) "s"
325        else "g"
326    }
327
328
329    class2 <- function(cl, kind = "l", do.sub = TRUE) {
330        ## Find "corresponding" class; since pos.def. matrices have no pendant:
331        if     (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
332        else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
333        else if(do.sub) sub("^d", kind, cl)
334        else cl
335    }
336
337    geClass <- function(x) {
338        if     (is(x, "dMatrix")) "dgeMatrix"
339        else if(is(x, "lMatrix")) "lgeMatrix"
340        else if(is(x, "zMatrix")) "zgeMatrix"
341        else stop("general Matrix class not yet implemented for ",
342                  class(x))
343    }
344
345    .dense.prefixes <- c("d" = "di",
346                         "t" = "tr",
347                         "s" = "sy",
348                         "g" = "ge")
349
350    .Csparse.prefix <- function(ch) {
351        switch(ch,
352               "d" =, "t" = "tC",
353               "s" = "sC",
354               "g" = "gC",
355               stop("invalid Matrix shape: ", ch))
356    }
357
358    ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
359    as_dense <- function(x) {
360        as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
361    }
362
363    as_Csparse <- function(x) {
364        as(x, paste(.M.kind(x), .Csparse.prefix(.M.shape(x)), "Matrix", sep=''))
365    }
366
367    as_geClass <- function(x, cl) {
368        if     (extends(cl, "diagonalMatrix")  && isDiagonal(x))
369            as(x, cl)
370        else if(extends(cl, "symmetricMatrix") &&  isSymmetric(x))
371            as(x, class2(cl, kind, do.sub= kind != "d"))
372        else if(extends(cl, "triangularMatrix") && isTriangular(x))
373            as(x, cl)
374        else
375            as(x, paste(.M.kind(x), "geMatrix", sep=''))
376    }
377
378    as_CspClass <- function(x, cl) {
379        if ((extends(cl, "diagonalMatrix")  && isDiagonal(x)) ||
380            (extends(cl, "symmetricMatrix") &&  isSymmetric(x)) ||
381            (extends(cl, "triangularMatrix")&& isTriangular(x)))
382            as(x, cl)
383        else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
384    }
385
386
387  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
388  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
389      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
390      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
391                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
392                 from)                 from)
393      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 405
405      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
406  }  }
407
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)))
}
408
409  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
410    isTriMat <- function(object, upper = NA) {
411      ## pretest: is it square?      ## pretest: is it square?
412      d <- dim(object)      d <- dim(object)
413      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 415
415      if(!is.matrix(object))      if(!is.matrix(object))
416          object <- as(object,"matrix")          object <- as(object,"matrix")
417      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
418      if(upper)      if(is.na(upper)) {
419            if(all(object[lower.tri(object)] == 0))
420                structure(TRUE, kind = "U")
421            else if(all(object[upper.tri(object)] == 0))
422                structure(TRUE, kind = "L")
423            else FALSE
424        } else if(upper)
425          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
426      else      else ## upper is FALSE
427          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
428  }  }
429
430    ## For Csparse matrices
431    isTriC <- function(x, upper = NA) {
432        ## pretest: is it square?
433        d <- dim(x)
434        if(d[1] != d[2]) return(FALSE)
435        ## else
436        if(d[1] == 0) return(TRUE)
437        ni <- 1:d[2]
438        ## the row indices split according to column:
439        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
440        lil <- unlist(lapply(ilist, length), use.names = FALSE)
441        if(any(lil == 0)) {
442            pos <- lil > 0
443            if(!any(pos)) ## matrix of all 0's
444                return(TRUE)
445            ilist <- ilist[pos]
446            ni <- ni[pos]
447        }
448        ni0 <- ni - 1:1 # '0-based ni'
449        if(is.na(upper)) {
450            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
451                structure(TRUE, kind = "U")
452            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
453                structure(TRUE, kind = "L")
454            else FALSE
455        } else if(upper) {
456            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
457        } else { ## 'lower'
458            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
459        }
460    }
461
462  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
463      d <- dim(object)      d <- dim(object)
464      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
465      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)
466  }  }
467
468    diagU2N <- function(x)
469    {
470        ## Purpose: Transform a *unit diagonal* triangular matrix
471        ##  into one with explicit diagonal entries '1'
472        xT <- as(x, "dgTMatrix")
473        ## leave it as  T* - the caller can always coerce to C* if needed:
474        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
475            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
476    }
477
478    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
479        x <- as(x, "dgCMatrix")
480        callGeneric()
481    }
482
483    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
484        x <- as(x, "dgTMatrix")
485        callGeneric()
486    }
487
488
489    ### Fast much simplified version of tapply()
490    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
491        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
492    }
493
494    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
495    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
496    ## }
497

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