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 1290, Thu Jun 8 09:30:21 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 260  Line 324 
324                           useInherited = c(from = TRUE, to = FALSE)))                           useInherited = c(from = TRUE, to = FALSE)))
325  }  }
326    
327  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
328    isTriMat <- function(object, upper = NA) {
329      ## pretest: is it square?      ## pretest: is it square?
330      d <- dim(object)      d <- dim(object)
331      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 333 
333      if(!is.matrix(object))      if(!is.matrix(object))
334          object <- as(object,"matrix")          object <- as(object,"matrix")
335      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
336      if(upper)      if(is.na(upper)) {
337            if(all(object[lower.tri(object)] == 0))
338                structure(TRUE, kind = "U")
339            else if(all(object[upper.tri(object)] == 0))
340                structure(TRUE, kind = "L")
341            else FALSE
342        } else if(upper)
343          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
344      else      else ## upper is FALSE
345          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
346  }  }
347    
348    ## For Csparse matrices
349    isTriC <- function(x, upper = NA) {
350        ## pretest: is it square?
351        d <- dim(x)
352        if(d[1] != d[2]) return(FALSE)
353        ## else
354        if(d[1] == 0) return(TRUE)
355        ni <- 1:d[1]
356        ## the row indices split according to column:
357        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
358        lil <- unlist(lapply(ilist, length), use.names = FALSE)
359        if(any(lil == 0)) {
360            pos <- lil > 0
361            if(!any(pos)) ## matrix of all 0's
362                return(TRUE)
363            ilist <- ilist[pos]
364            ni <- ni[pos]
365        }
366        if(is.na(upper)) {
367            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni))
368                structure(TRUE, kind = "U")
369            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni))
370                structure(TRUE, kind = "L")
371            else FALSE
372        } else if(upper) {
373            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni)
374        } else { ## 'lower'
375            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni)
376        }
377    }
378    
379  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
380      d <- dim(object)      d <- dim(object)
381      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
382      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)
383  }  }
384    
385    diagU2N <- function(x)
386    {
387        ## Purpose: Transform a *unit diagonal* triangular matrix
388        ##  into one with explicit diagonal entries '1'
389        xT <- as(x, "dgTMatrix")
390        ## leave it as  T* - the caller can always coerce to C* if needed:
391        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
392            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
393    }
394    
395    ### Fast much simplified version of tapply()
396    tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
397        sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
398    }
399    
400    ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
401    ##     tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
402    ## }
403    

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

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