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 1295, Fri Jun 9 21:47:22 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        ## return a 2-column (i,j) matrix of
148        ## 0-based indices of non-zero entries  :
149      if(is(x, "TsparseMatrix"))      if(is(x, "TsparseMatrix"))
150          return(unique(cbind(x@i,x@j)))          return(unique(cbind(x@i,x@j)))
151        ## else:
152        isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
153        .Call(compressed_non_0_ij, x, isC)
154    }
155    
156    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
157    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
158    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
159    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
160    
161      isCol <- function(M) any("i" == slotNames(M))  complementInd <- function(ij, dim)
162      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  {
163        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
164        ##          but as 1-based indices
165        n <- prod(dim)
166        if(n == 0) return(integer(0))
167        ii <- 1:n
168        ii[-(1 + encodeInd(ij, nr = dim[1]))]
169    }
170    
171    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
172    
173    intersectInd <- function(ij1, ij2, nrow) {
174        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
175        ## return only the *common* entries
176        decodeInd(intersect(encodeInd(ij1, nrow),
177                            encodeInd(ij2, nrow)), nrow)
178    }
179    
180    WhichintersectInd <- function(ij1, ij2, nrow) {
181        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
182        ## find *where*  common entries are in ij1 & ij2
183        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
184        ni <- !is.na(m1)
185        list(which(ni), m1[ni])
186  }  }
187    
188    
189  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
190  uniq <- function(x) {  uniq <- function(x) {
191    
192    ### Note: maybe, using
193    ### ----    xj <- .Call(Matrix_expand_pointers, x@p)
194    ### would be slightly more efficient than as( <dgC> , "dgTMatrix")
195    ### but really efficient would be to use only one .Call(.) for uniq(.) !
196    ### Try to do it particularly fast for the case where 'x' is already a 'uniq' <dgT>
197    
198      if(is(x, "TsparseMatrix")) {      if(is(x, "TsparseMatrix")) {
199          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
200          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
201          ## -----------------------------------------------------------          ## -----------------------------------------------------------
202          ## The following is *not* efficient {but easy to program}:          ## The following is not quite efficient {but easy to program,
203            ## and both as() are based on C code
204          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
205          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
206          else stop("not implemented for class", class(x))          else stop("not implemented for class", class(x))
207    
208      } else x      # not 'gT' ; i.e. "uniquely" represented in any case      } else x  ## not 'gT' ; i.e. "uniquely" represented in any case
209  }  }
210    
211  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 213 
213  {  {
214      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
215      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
216      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
217      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
218      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
219      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 275 
275      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
276  }  }
277    
278    dClass2 <- function(dClass, kind = "l") {
279        ## Find "corresponding" class for a dMatrix;
280        #  since pos.def. matrices have no pendant:
281        if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
282        else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
283        else sub("^d", kind, dClass)
284    }
285    
286    geClass <- function(x) {
287        if(is(x, "dMatrix")) "dgeMatrix"
288        else if(is(x, "lMatrix")) "lgeMatrix"
289        else stop("general Matrix class not yet implemented for",
290                  class(x))
291    }
292    
293  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
294  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
295      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
296      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
297                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
298                 from)                 from)
299      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 247  Line 311 
311      if(ok) as(x, classes[1]) else x      if(ok) as(x, classes[1]) else x
312  }  }
313    
 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)))  
 }  
314    
315  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
316    isTriMat <- function(object, upper = NA) {
317      ## pretest: is it square?      ## pretest: is it square?
318      d <- dim(object)      d <- dim(object)
319      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 321 
321      if(!is.matrix(object))      if(!is.matrix(object))
322          object <- as(object,"matrix")          object <- as(object,"matrix")
323      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
324      if(upper)      if(is.na(upper)) {
325            if(all(object[lower.tri(object)] == 0))
326                structure(TRUE, kind = "U")
327            else if(all(object[upper.tri(object)] == 0))
328                structure(TRUE, kind = "L")
329            else FALSE
330        } else if(upper)
331          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
332      else      else ## upper is FALSE
333          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
334  }  }
335    
336    ## For Csparse matrices
337    isTriC <- function(x, upper = NA) {
338        ## pretest: is it square?
339        d <- dim(x)
340        if(d[1] != d[2]) return(FALSE)
341        ## else
342        if(d[1] == 0) return(TRUE)
343        ni <- 1:d[1]
344        ## the row indices split according to column:
345        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
346        lil <- unlist(lapply(ilist, length), use.names = FALSE)
347        if(any(lil == 0)) {
348            pos <- lil > 0
349            if(!any(pos)) ## matrix of all 0's
350                return(TRUE)
351            ilist <- ilist[pos]
352            ni <- ni[pos]
353        }
354        if(is.na(upper)) {
355            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni))
356                structure(TRUE, kind = "U")
357            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni))
358                structure(TRUE, kind = "L")
359            else FALSE
360        } else if(upper) {
361            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni)
362        } else { ## 'lower'
363            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni)
364        }
365    }
366    
367  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
368      d <- dim(object)      d <- dim(object)
369      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
370      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)
371  }  }
372    
373    diagU2N <- function(x)
374    {
375        ## Purpose: Transform a *unit diagonal* triangular matrix
376        ##  into one with explicit diagonal entries '1'
377        xT <- as(x, "dgTMatrix")
378        ## leave it as  T* - the caller can always coerce to C* if needed:
379        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
380            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
381    }
382    
383    .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
384        x <- as(x, "dgCMatrix")
385        callGeneric()
386    }
387    
388    .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
389        x <- as(x, "dgTMatrix")
390        callGeneric()
391    }
392    
393    
394    ### Fast much simplified version of tapply()
395    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
396        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
397    }
398    
399    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
400    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
401    ## }
402    

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

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