# SCM Repository

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

# Diff of /pkg/R/Auxiliaries.R

revision 1225, Mon Mar 13 14:06:17 2006 UTC revision 1226, Mon Mar 13 14:07:58 2006 UTC
# Line 134  Line 134
134  non0ind <- function(x) {  non0ind <- function(x) {
135      if(is.numeric(x))      if(is.numeric(x))
136          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))
137        ## else
## else return a (i,j) matrix of non-zero-indices

138      stopifnot(is(x, "sparseMatrix"))      stopifnot(is(x, "sparseMatrix"))
139        ## return a 2-column (i,j) matrix of
140        ## 0-based indices of non-zero entries  :
141      if(is(x, "TsparseMatrix"))      if(is(x, "TsparseMatrix"))
142          return(unique(cbind(x@i,x@j)))          return(unique(cbind(x@i,x@j)))
143        ## else:
144        isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
145        .Call("compressed_non_0_ij", x, isC, PACKAGE = "Matrix")
146    }
147
148    ## nr= nrow: since  i in {0,1,.., nrow-1}  these are 1:1 "decimal" encodings:
149    ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
150    encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
151    decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
152
153      isCol <- function(M) any("i" == slotNames(M))  complementInd <- function(ij, dim)
154      .Call("compressed_non_0_ij", x, isCol(x), PACKAGE = "Matrix")  {
155        ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
156        ##          but as 1-based indices
157        n <- prod(dim)
158        if(n == 0) return(integer(0))
159        ii <- 1:n
160        ii[-(1 + encodeInd(ij, nr = dim[1]))]
161    }
162
163    unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
164
165    intersectInd <- function(ij1, ij2, nrow) {
166        ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
167        ## return only the *common* entries
168        decodeInd(intersect(encodeInd(ij1, nrow),
169                            encodeInd(ij2, nrow)), nrow)
170    }
171
172    WhichintersectInd <- function(ij1, ij2, nrow) {
173        ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
174        ## find *where*  common entries are in ij1 & ij2
175        m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
176        ni <- !is.na(m1)
177        list(which(ni), m1[ni])
178  }  }
179
180
181  ### There is a test on this in ../tests/dgTMatrix.R !  ### There is a test on this in ../tests/dgTMatrix.R !
182  uniq <- function(x) {  uniq <- function(x) {
183      if(is(x, "TsparseMatrix")) {      if(is(x, "TsparseMatrix")) {
184          ## Purpose: produce a *unique* triplet representation:          ## Purpose: produce a *unique* triplet representation:
185          ##              by having (i,j) sorted and unique          ##              by having (i,j) sorted and unique
186          ## -----------------------------------------------------------          ## -----------------------------------------------------------
187          ## The following is *not* efficient {but easy to program}:          ## The following is not quite efficient {but easy to program,
188            ## and both as() are based on C code
189          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")          if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
190          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")          else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
191          else stop("not implemented for class", class(x))          else stop("not implemented for class", class(x))
192
193      } else x      # not 'gT' ; i.e. "uniquely" represented in any case      } else x  ## not 'gT' ; i.e. "uniquely" represented in any case
194  }  }
195
196  if(FALSE) ## try an "efficient" version  if(FALSE) ## try an "efficient" version
# Line 164  Line 198
198  {  {
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      ## Arguments: a "gT" Matrix      ## Arguments: a "gT" Matrix
203      stopifnot(is(x, "gTMatrix"))      stopifnot(is(x, "gTMatrix"))
204      if((n <- length(x@i)) == 0) return(x)      if((n <- length(x@i)) == 0) return(x)
# Line 226  Line 260
260      as(callGeneric(as(x, sub("^l", "d", cl))), cl)      as(callGeneric(as(x, sub("^l", "d", cl))), cl)
261  }  }
262
263    dClass2 <- function(dClass, kind = "l") {
264        ## Find "corresponding" class for a dMatrix;
265        #  since pos.def. matrices have no pendant:
266        if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
267        else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
268        else sub("^d", kind, dClass)
269    }
270
271    geClass <- function(x) {
272        if(is(x, "dMatrix")) "dgeMatrix"
273        else if(is(x, "lMatrix")) "lgeMatrix"
274        else stop("general Matrix class not yet implemented for",
275                  class(x))
276    }
277
278  ## -> ./ddenseMatrix.R :  ## -> ./ddenseMatrix.R :
279  d2l_Matrix <- function(from) {  d2l_Matrix <- function(from) {
280      stopifnot(is(from, "dMatrix"))      stopifnot(is(from, "dMatrix"))
281      fixupDense(new(sub("^d", "l", class(from)),      fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
282                     Dim = from@Dim, Dimnames = from@Dimnames),                     Dim = from@Dim, Dimnames = from@Dimnames),
283                 from)                 from)
284      ## FIXME: treat 'factors' smartly {not for triangular!}      ## FIXME: treat 'factors' smartly {not for triangular!}
# Line 260  Line 309
309                           useInherited = c(from = TRUE, to = FALSE)))                           useInherited = c(from = TRUE, to = FALSE)))
310  }  }
311
312  .is.triangular <- function(object, upper = TRUE) {  .is.triangular <- function(object, upper = NA) {
313      ## pretest: is it square?      ## pretest: is it square?
314      d <- dim(object)      d <- dim(object)
315      if(d[1] != d[2]) return(FALSE)      if(d[1] != d[2]) return(FALSE)
# Line 268  Line 317
317      if(!is.matrix(object))      if(!is.matrix(object))
318          object <- as(object,"matrix")          object <- as(object,"matrix")
319      ## == 0 even works for logical & complex:      ## == 0 even works for logical & complex:
320      if(upper)      if(is.na(upper)) {
321            if(all(object[lower.tri(object)] == 0))
322                structure(TRUE, kind = "U")
323            else if(all(object[upper.tri(object)] == 0))
324                structure(TRUE, kind = "L")
325            else FALSE
326        } else if(upper)
327          all(object[lower.tri(object)] == 0)          all(object[lower.tri(object)] == 0)
328      else      else ## upper is FALSE
329          all(object[upper.tri(object)] == 0)          all(object[upper.tri(object)] == 0)
330  }  }
331

Legend:
 Removed from v.1225 changed lines Added in v.1226