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 1253, Wed Apr 19 09:17:14 2006 UTC
# Line 82  Line 82 
82      da[2]      da[2]
83  }  }
84    
85    isPacked <- function(x)
86    {
87        ## Is 'x' a packed (dense) matrix ?
88        is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
89    }
90    
91  emptyColnames <- function(x)  emptyColnames <- function(x)
92  {  {
93      ## Useful for compact printing of (parts) of sparse matrices      ## Useful for compact printing of (parts) of sparse matrices
# Line 134  Line 140 
140  non0ind <- function(x) {  non0ind <- function(x) {
141      if(is.numeric(x))      if(is.numeric(x))
142          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))
143        ## else
     ## else return a (i,j) matrix of non-zero-indices  
   
144      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
145        ## return a 2-column (i,j) matrix of
146        ## 0-based indices of non-zero entries  :
147      if(is(x, "TsparseMatrix"))      if(is(x, "TsparseMatrix"))
148          return(unique(cbind(x@i,x@j)))          return(unique(cbind(x@i,x@j)))
149        ## else:
150        isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
151        .Call("compressed_non_0_ij", x, isC, PACKAGE = "Matrix")
152    }
153    
154      isCol <- function(M) any("i" == slotNames(M))  ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
155      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
156    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
157    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
158    
159    complementInd <- function(ij, dim)
160    {
161        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
162        ##          but as 1-based indices
163        n <- prod(dim)
164        if(n == 0) return(integer(0))
165        ii <- 1:n
166        ii[-(1 + encodeInd(ij, nr = dim[1]))]
167    }
168    
169    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
170    
171    intersectInd <- function(ij1, ij2, nrow) {
172        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
173        ## return only the *common* entries
174        decodeInd(intersect(encodeInd(ij1, nrow),
175                            encodeInd(ij2, nrow)), nrow)
176    }
177    
178    WhichintersectInd <- function(ij1, ij2, nrow) {
179        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
180        ## find *where*  common entries are in ij1 & ij2
181        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
182        ni <- !is.na(m1)
183        list(which(ni), m1[ni])
184  }  }
185    
186    
187  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
188  uniq <- function(x) {  uniq <- function(x) {
189      if(is(x, "TsparseMatrix")) {      if(is(x, "TsparseMatrix")) {
190          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
191          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
192          ## -----------------------------------------------------------          ## -----------------------------------------------------------
193          ## The following is *not* efficient {but easy to program}:          ## The following is not quite efficient {but easy to program,
194            ## and both as() are based on C code
195          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
196          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
197          else stop("not implemented for class", class(x))          else stop("not implemented for class", class(x))
198    
199      } else x      # not 'gT' ; i.e. "uniquely" represented in any case      } else x  ## not 'gT' ; i.e. "uniquely" represented in any case
200  }  }
201    
202  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 204 
204  {  {
205      ## Purpose: produce a *unique* triplet representation:      ## Purpose: produce a *unique* triplet representation:
206      ##          by having (i,j) sorted and unique      ##          by having (i,j) sorted and unique
207      ## ----------------------------------------------------------------------      ## ------------------------------------------------------------------
208      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
209      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
210      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 266 
266      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
267  }  }
268    
269    dClass2 <- function(dClass, kind = "l") {
270        ## Find "corresponding" class for a dMatrix;
271        #  since pos.def. matrices have no pendant:
272        if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
273        else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
274        else sub("^d", kind, dClass)
275    }
276    
277    geClass <- function(x) {
278        if(is(x, "dMatrix")) "dgeMatrix"
279        else if(is(x, "lMatrix")) "lgeMatrix"
280        else stop("general Matrix class not yet implemented for",
281                  class(x))
282    }
283    
284  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
285  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
286      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
287      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
288                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
289                 from)                 from)
290      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 260  Line 315 
315                           useInherited = c(from = TRUE, to = FALSE)))                           useInherited = c(from = TRUE, to = FALSE)))
316  }  }
317    
318  .is.triangular <- function(object, upper = TRUE) {  ## For *dense* matrices
319    isTriMat <- function(object, upper = NA) {
320      ## pretest: is it square?      ## pretest: is it square?
321      d <- dim(object)      d <- dim(object)
322      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 324 
324      if(!is.matrix(object))      if(!is.matrix(object))
325          object <- as(object,"matrix")          object <- as(object,"matrix")
326      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
327      if(upper)      if(is.na(upper)) {
328            if(all(object[lower.tri(object)] == 0))
329                structure(TRUE, kind = "U")
330            else if(all(object[upper.tri(object)] == 0))
331                structure(TRUE, kind = "L")
332            else FALSE
333        } else if(upper)
334          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
335      else      else ## upper is FALSE
336          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
337  }  }
338    
339    ## For Csparse matrices
340    isTriC <- function(x, upper = NA) {
341        ## pretest: is it square?
342        d <- dim(x)
343        if(d[1] != d[2]) return(FALSE)
344        ## else
345        if(d[1] == 0) return(TRUE)
346        ni <- 1:d[1]
347        ## the row indices split according to column:
348        ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
349        lil <- unlist(lapply(ilist, length), use.names = FALSE)
350        if(any(lil == 0)) {
351            pos <- lil > 0
352            if(!any(pos)) ## matrix of all 0's
353                return(TRUE)
354            ilist <- ilist[pos]
355            ni <- ni[pos]
356        }
357        if(is.na(upper)) {
358            if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni))
359                structure(TRUE, kind = "U")
360            else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni))
361                structure(TRUE, kind = "L")
362            else FALSE
363        } else if(upper) {
364            all(sapply(ilist, max, USE.NAMES = FALSE) <= ni)
365        } else { ## 'lower'
366            all(sapply(ilist, min, USE.NAMES = FALSE) >= ni)
367        }
368    }
369    
370  .is.diagonal <- function(object) {  .is.diagonal <- function(object) {
371      d <- dim(object)      d <- dim(object)
372      if(d[1] != (n <- d[2])) FALSE      if(d[1] != (n <- d[2])) FALSE
373      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)
374  }  }
375    
376    diagU2N <- function(x)
377    {
378        ## Purpose: Transform a *unit diagonal* triangular matrix
379        ##  into one with explicit diagonal entries '1'
380        xT <- as(x, "dgTMatrix")
381        ## leave it as  T* - the caller can always coerce to C* if needed:
382        new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
383            Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
384    }

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

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