# SCM Repository

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

# Annotation of /pkg/R/Auxiliaries.R

Revision 1227 - (view) (download)

 1 : maechler 632 #### "Namespace private" Auxiliaries such as method functions 2 : #### (called from more than one place --> need to be defined early) 3 : 4 : maechler 656 ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}): 5 : .M.v <- function(x, y) callGeneric(x, as.matrix(y)) 6 : .v.M <- function(x, y) callGeneric(rbind(x), y) 7 : maechler 632 8 : maechler 656 .has.DN <- ## has non-trivial Dimnames slot? 9 : function(x) !identical(list(NULL,NULL), x@Dimnames) 10 : 11 : maechler 949 .bail.out.1 <- function(fun, cl) { 12 : stop(gettextf('not-yet-implemented method for %s(<%s>)', fun, cl), 13 : call. = FALSE) 14 : } 15 : .bail.out.2 <- function(fun, cl1, cl2) { 16 : stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)', 17 : fun, cl1, cl2), call. = FALSE) 18 : } 19 : 20 : maechler 632 ## chol() via "dpoMatrix" 21 : cholMat <- function(x, pivot, LINPACK) { 22 : px <- as(x, "dpoMatrix") 23 : bates 703 if (isTRUE(validObject(px, test=TRUE))) chol(px) 24 : maechler 632 else stop("'x' is not positive definite -- chol() undefined.") 25 : } 26 : maechler 908 27 : maechler 954 dimCheck <- function(a, b) { 28 : da <- dim(a) 29 : db <- dim(b) 30 : if(any(da != db)) 31 : stop(gettextf("Matrices must have same dimensions in %s", 32 : deparse(sys.call(sys.parent()))), 33 : call. = FALSE) 34 : da 35 : } 36 : 37 : maechler 956 dimNamesCheck <- function(a, b) { 38 : ## assume dimCheck() has happened before 39 : nullDN <- list(NULL,NULL) 40 : h.a <- !identical(nullDN, dna <- dimnames(a)) 41 : h.b <- !identical(nullDN, dnb <- dimnames(b)) 42 : if(h.a || h.b) { 43 : maechler 1084 if (!h.b) dna 44 : else if(!h.a) dnb 45 : maechler 956 else { ## both have non-trivial dimnames 46 : r <- dna # "default" result 47 : for(j in 1:2) { 48 : dn <- dnb[[j]] 49 : if(is.null(r[[j]])) 50 : r[[j]] <- dn 51 : else if (!is.null(dn) && any(r[[j]] != dn)) 52 : warning(gettextf("dimnames [%d] mismatch in %s", j, 53 : deparse(sys.call(sys.parent()))), 54 : call. = FALSE) 55 : } 56 : r 57 : } 58 : } 59 : else 60 : nullDN 61 : } 62 : 63 : maechler 908 rowCheck <- function(a, b) { 64 : da <- dim(a) 65 : db <- dim(b) 66 : if(da[1] != db[1]) 67 : stop(gettextf("Matrices must have same number of rows in %s", 68 : deparse(sys.call(sys.parent()))), 69 : call. = FALSE) 70 : ## return the common nrow() 71 : da[1] 72 : } 73 : 74 : colCheck <- function(a, b) { 75 : da <- dim(a) 76 : db <- dim(b) 77 : if(da[2] != db[2]) 78 : stop(gettextf("Matrices must have same number of columns in %s", 79 : deparse(sys.call(sys.parent()))), 80 : call. = FALSE) 81 : ## return the common ncol() 82 : da[2] 83 : } 84 : 85 : maechler 1227 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 : maechler 956 emptyColnames <- function(x) 92 : { 93 : ## Useful for compact printing of (parts) of sparse matrices 94 : ## possibly dimnames(x) "==" NULL : 95 : dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2])) 96 : x 97 : } 98 : maechler 908 99 : maechler 919 prTriang <- function(x, digits = getOption("digits"), 100 : justify = "none", right = TRUE) 101 : { 102 : ## modeled along stats:::print.dist 103 : diag <- TRUE 104 : upper <- x@uplo == "U" 105 : 106 : m <- as(x, "matrix") 107 : cf <- format(m, digits = digits, justify = justify) 108 : if(upper) 109 : cf[row(cf) > col(cf)] <- "." 110 : else 111 : cf[row(cf) < col(cf)] <- "." 112 : print(cf, quote = FALSE, right = right) 113 : invisible(x) 114 : } 115 : 116 : prMatrix <- function(x, digits = getOption("digits")) { 117 : d <- dim(x) 118 : cl <- class(x) 119 : cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl)) 120 : maxp <- getOption("max.print") 121 : if(prod(d) <= maxp) { 122 : if(is(x, "triangularMatrix")) 123 : prTriang(x, digits = digits) 124 : else 125 : print(as(x, "matrix"), digits = digits) 126 : } 127 : else { ## d[1] > maxp / d[2] >= nr : 128 : m <- as(x, "matrix") 129 : nr <- maxp %/% d[2] 130 : n2 <- ceiling(nr / 2) 131 : print(head(m, max(1, n2))) 132 : cat("\n ..........\n\n") 133 : print(tail(m, max(1, nr - n2))) 134 : } 135 : ## DEBUG: cat("str(.):\n") ; str(x) 136 : invisible(x)# as print() S3 methods do 137 : } 138 : 139 : ## For sparseness handling 140 : non0ind <- function(x) { 141 : if(is.numeric(x)) 142 : maechler 1226 return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0)) 143 : ## else 144 : maechler 919 stopifnot(is(x, "sparseMatrix")) 145 : maechler 1226 ## return a 2-column (i,j) matrix of 146 : ## 0-based indices of non-zero entries : 147 : maechler 925 if(is(x, "TsparseMatrix")) 148 : maechler 954 return(unique(cbind(x@i,x@j))) 149 : maechler 1226 ## else: 150 : isC <- any("i" == slotNames(x))# is Csparse (not Rsparse) 151 : .Call("compressed_non_0_ij", x, isC, PACKAGE = "Matrix") 152 : } 153 : maechler 954 154 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings: 155 : ## 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 : maechler 919 } 168 : 169 : maechler 1226 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 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R ! 188 : maechler 919 uniq <- function(x) { 189 : maechler 925 if(is(x, "TsparseMatrix")) { 190 : maechler 919 ## Purpose: produce a *unique* triplet representation: 191 : ## by having (i,j) sorted and unique 192 : ## ----------------------------------------------------------- 193 : maechler 1226 ## The following is not quite efficient {but easy to program, 194 : ## and both as() are based on C code 195 : maechler 919 if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix") 196 : else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix") 197 : else stop("not implemented for class", class(x)) 198 : 199 : maechler 1226 } else x ## not 'gT' ; i.e. "uniquely" represented in any case 200 : maechler 919 } 201 : 202 : if(FALSE) ## try an "efficient" version 203 : uniq_gT <- function(x) 204 : { 205 : ## Purpose: produce a *unique* triplet representation: 206 : ## by having (i,j) sorted and unique 207 : maechler 1226 ## ------------------------------------------------------------------ 208 : maechler 919 ## Arguments: a "gT" Matrix 209 : stopifnot(is(x, "gTMatrix")) 210 : if((n <- length(x@i)) == 0) return(x) 211 : ii <- order(x@i, x@j) 212 : if(any(ii != 1:n)) { 213 : x@i <- x@i[ii] 214 : x@j <- x@j[ii] 215 : x@x <- x@x[ii] 216 : } 217 : ij <- x@i + nrow(x) * x@j 218 : if(any(dup <- duplicated(ij))) { 219 : 220 : } 221 : ### We should use a .Call() based utility for this! 222 : 223 : } 224 : 225 : maechler 946 t_geMatrix <- function(x) { 226 : x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here 227 : x@Dim <- x@Dim[2:1] 228 : x@Dimnames <- x@Dimnames[2:1] 229 : ## FIXME: how to set factors? 230 : x 231 : } 232 : 233 : ## t( [dl]trMatrix ) and t( [dl]syMatrix ) : 234 : t_trMatrix <- function(x) { 235 : x@x <- as.vector(t(as(x, "matrix"))) 236 : x@Dim <- x@Dim[2:1] 237 : x@Dimnames <- x@Dimnames[2:1] 238 : x@uplo <- if (x@uplo == "U") "L" else "U" 239 : # and keep x@diag 240 : x 241 : } 242 : maechler 956 243 : fixupDense <- function(m, from) { 244 : if(is(m, "triangularMatrix")) { 245 : m@uplo <- from@uplo 246 : m@diag <- from@diag 247 : } else if(is(m, "symmetricMatrix")) { 248 : m@uplo <- from@uplo 249 : } 250 : m 251 : } 252 : 253 : ## -> ./ldenseMatrix.R : 254 : l2d_Matrix <- function(from) { 255 : stopifnot(is(from, "lMatrix")) 256 : fixupDense(new(sub("^l", "d", class(from)), 257 : x = as.double(from@x), 258 : maechler 1198 Dim = from@Dim, Dimnames = from@Dimnames), 259 : maechler 956 from) 260 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!} 261 : maechler 956 } 262 : 263 : if(FALSE)# unused 264 : l2d_meth <- function(x) { 265 : cl <- class(x) 266 : as(callGeneric(as(x, sub("^l", "d", cl))), cl) 267 : } 268 : 269 : maechler 1226 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 : maechler 956 ## -> ./ddenseMatrix.R : 285 : d2l_Matrix <- function(from) { 286 : stopifnot(is(from, "dMatrix")) 287 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here 288 : maechler 1198 Dim = from@Dim, Dimnames = from@Dimnames), 289 : maechler 956 from) 290 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!} 291 : maechler 956 } 292 : maechler 973 293 : 294 : try_as <- function(x, classes, tryAnyway = FALSE) { 295 : if(!tryAnyway && !is(x, "Matrix")) 296 : return(x) 297 : ## else 298 : ok <- canCoerce(x, classes[1]) 299 : while(!ok && length(classes <- classes[-1])) { 300 : ok <- canCoerce(x, classes[1]) 301 : } 302 : if(ok) as(x, classes[1]) else x 303 : } 304 : 305 : maechler 1174 if(paste(R.version$major, R.version$minor, sep=".") < "2.3") 306 : ## This will be in R 2.3.0 307 : maechler 973 canCoerce <- function(object, Class) { 308 : ## Purpose: test if 'object' is coercable to 'Class', i.e., 309 : ## as(object, Class) will {typically} work 310 : ## ---------------------------------------------------------------------- 311 : ## Author: John Chambers, Date: 6 Oct 2005 312 : is(object, Class) || 313 : !is.null(selectMethod("coerce", c(class(object), Class), 314 : optional = TRUE, 315 : useInherited = c(from = TRUE, to = FALSE))) 316 : } 317 : 318 : maechler 1226 .is.triangular <- function(object, upper = NA) { 319 : maechler 1174 ## pretest: is it square? 320 : d <- dim(object) 321 : if(d[1] != d[2]) return(FALSE) 322 : ## else slower test 323 : if(!is.matrix(object)) 324 : maechler 1226 object <- as(object,"matrix") 325 : maechler 1174 ## == 0 even works for logical & complex: 326 : maechler 1226 if(is.na(upper)) { 327 : if(all(object[lower.tri(object)] == 0)) 328 : structure(TRUE, kind = "U") 329 : else if(all(object[upper.tri(object)] == 0)) 330 : structure(TRUE, kind = "L") 331 : else FALSE 332 : } else if(upper) 333 : all(object[lower.tri(object)] == 0) 334 : else ## upper is FALSE 335 : all(object[upper.tri(object)] == 0) 336 : maechler 1174 } 337 : 338 : .is.diagonal <- function(object) { 339 : d <- dim(object) 340 : if(d[1] != (n <- d[2])) FALSE 341 : else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0) 342 : }

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: