# 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 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: