# SCM Repository

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

# Annotation of /pkg/R/Auxiliaries.R

 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 1472 .isR_24 <- (paste(R.version$major, R.version$minor, sep=".") >= "2.4") 5 : 6 : maechler 1467 ## Need to consider NAs ; "== 0" even works for logical & complex: 7 : is0 <- function(x) !is.na(x) & x == 0 8 : maechler 1575 isN0 <- function(x) is.na(x) | x != 0 9 : maechler 1467 all0 <- function(x) !any(is.na(x)) && all(x == 0) 10 : 11 : maechler 1571 allTrue <- function(x) !any(is.na(x)) && all(x) 12 : allFalse <- function(x) !any(is.na(x)) && !any(x) 13 : 14 : 15 : maechler 656 ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}): 16 : .M.v <- function(x, y) callGeneric(x, as.matrix(y)) 17 : .v.M <- function(x, y) callGeneric(rbind(x), y) 18 : maechler 632 19 : maechler 1332 .M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL) 20 : 21 : maechler 656 .has.DN <- ## has non-trivial Dimnames slot? 22 : function(x) !identical(list(NULL,NULL), x@Dimnames) 23 : 24 : maechler 949 .bail.out.1 <- function(fun, cl) { 25 : maechler 1654 stop(gettextf('not-yet-implemented method for %s(<%s>).\n ->> Ask the package authors to implement the missing feature.', fun, cl), 26 : maechler 1238 call. = FALSE) 27 : maechler 949 } 28 : .bail.out.2 <- function(fun, cl1, cl2) { 29 : maechler 1654 stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>).\n ->> Ask the package authors to implement the missing feature.', 30 : maechler 1238 fun, cl1, cl2), call. = FALSE) 31 : maechler 949 } 32 : 33 : maechler 1571 ## This should be done in C and be exported by 'methods': [FIXME - ask JMC ] 34 : copyClass <- function(x, newCl, sNames = 35 : intersect(slotNames(newCl), slotNames(x))) { 36 : r <- new(newCl) 37 : for(n in sNames) 38 : slot(r, n) <- slot(x, n) 39 : r 40 : } 41 : 42 : maechler 632 ## chol() via "dpoMatrix" 43 : maechler 1571 cholMat <- function(x, pivot, ...) { 44 : maechler 632 px <- as(x, "dpoMatrix") 45 : bates 703 if (isTRUE(validObject(px, test=TRUE))) chol(px) 46 : maechler 632 else stop("'x' is not positive definite -- chol() undefined.") 47 : } 48 : maechler 908 49 : maechler 954 dimCheck <- function(a, b) { 50 : da <- dim(a) 51 : db <- dim(b) 52 : if(any(da != db)) 53 : stop(gettextf("Matrices must have same dimensions in %s", 54 : deparse(sys.call(sys.parent()))), 55 : call. = FALSE) 56 : da 57 : } 58 : 59 : maechler 956 dimNamesCheck <- function(a, b) { 60 : ## assume dimCheck() has happened before 61 : nullDN <- list(NULL,NULL) 62 : h.a <- !identical(nullDN, dna <- dimnames(a)) 63 : h.b <- !identical(nullDN, dnb <- dimnames(b)) 64 : if(h.a || h.b) { 65 : maechler 1084 if (!h.b) dna 66 : else if(!h.a) dnb 67 : maechler 956 else { ## both have non-trivial dimnames 68 : r <- dna # "default" result 69 : for(j in 1:2) { 70 : dn <- dnb[[j]] 71 : if(is.null(r[[j]])) 72 : r[[j]] <- dn 73 : else if (!is.null(dn) && any(r[[j]] != dn)) 74 : warning(gettextf("dimnames [%d] mismatch in %s", j, 75 : deparse(sys.call(sys.parent()))), 76 : call. = FALSE) 77 : } 78 : r 79 : } 80 : } 81 : else 82 : nullDN 83 : } 84 : 85 : maechler 908 rowCheck <- function(a, b) { 86 : da <- dim(a) 87 : db <- dim(b) 88 : if(da[1] != db[1]) 89 : stop(gettextf("Matrices must have same number of rows in %s", 90 : deparse(sys.call(sys.parent()))), 91 : call. = FALSE) 92 : ## return the common nrow() 93 : da[1] 94 : } 95 : 96 : colCheck <- function(a, b) { 97 : da <- dim(a) 98 : db <- dim(b) 99 : if(da[2] != db[2]) 100 : stop(gettextf("Matrices must have same number of columns in %s", 101 : deparse(sys.call(sys.parent()))), 102 : call. = FALSE) 103 : ## return the common ncol() 104 : da[2] 105 : } 106 : 107 : maechler 1285 ## Note: !isPacked(.) i.e. `full' still contains 108 : ## ---- "*sy" and "*tr" which have "undefined" lower or upper part 109 : maechler 1227 isPacked <- function(x) 110 : { 111 : maechler 1472 ## Is 'x' a packed (dense) matrix ? 112 : is(x, "denseMatrix") && 113 : any("x" == slotNames(x)) && length(x@x) < prod(dim(x)) 114 : maechler 1227 } 115 : 116 : maechler 956 emptyColnames <- function(x) 117 : { 118 : ## Useful for compact printing of (parts) of sparse matrices 119 : maechler 1238 ## possibly dimnames(x) "==" NULL : 120 : maechler 956 dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2])) 121 : x 122 : } 123 : maechler 908 124 : maechler 1592 ### TODO: write in C and port to base (or 'utils') R 125 : indTri <- function(n, upper = TRUE) { 126 : ## == which(upper.tri(diag(n)) or 127 : ## which(lower.tri(diag(n)) -- but much more efficiently for largish 'n' 128 : stopifnot(length(n) == 1, n == (n. <- as.integer(n)), (n <- n.) >= 0) 129 : if(n <= 2) 130 : return(if(n == 2) as.integer(if(upper) n+1 else n) else integer(0)) 131 : ## First, compute the 'diff(.)' fast. Use integers 132 : one <- 1:1 ; two <- 2:2 133 : n1 <- n - one 134 : n2 <- n1 - one 135 : r <- rep.int(one, n*n1/two - one) 136 : r[cumsum(if(upper) 1:n2 else c(n1, if(n >= 4) n2:two))] <- if(upper) n:3 else 3:n 137 : ## now have "dliu" difference; revert to "liu": 138 : cumsum(c(if(upper) n+one else two, r)) 139 : } 140 : 141 : 142 : maechler 919 prTriang <- function(x, digits = getOption("digits"), 143 : maechler 1389 maxp = getOption("max.print"), 144 : maechler 1238 justify = "none", right = TRUE) 145 : maechler 919 { 146 : ## modeled along stats:::print.dist 147 : upper <- x@uplo == "U" 148 : 149 : m <- as(x, "matrix") 150 : cf <- format(m, digits = digits, justify = justify) 151 : if(upper) 152 : maechler 1238 cf[row(cf) > col(cf)] <- "." 153 : maechler 919 else 154 : maechler 1238 cf[row(cf) < col(cf)] <- "." 155 : maechler 1472 if(.isR_24) 156 : print(cf, quote = FALSE, right = right, max = maxp) 157 : else print(cf, quote = FALSE, right = right) 158 : maechler 919 invisible(x) 159 : } 160 : 161 : maechler 1389 prMatrix <- function(x, digits = getOption("digits"), 162 : maxp = getOption("max.print")) { 163 : maechler 919 d <- dim(x) 164 : cl <- class(x) 165 : cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl)) 166 : if(prod(d) <= maxp) { 167 : maechler 1238 if(is(x, "triangularMatrix")) 168 : maechler 1472 prTriang(x, digits = digits, maxp = maxp) 169 : else { 170 : if(.isR_24) 171 : print(as(x, "matrix"), digits = digits, max = maxp) 172 : else print(as(x, "matrix"), digits = digits) 173 : } 174 : maechler 919 } 175 : else { ## d[1] > maxp / d[2] >= nr : 176 : maechler 1238 m <- as(x, "matrix") 177 : maechler 919 nr <- maxp %/% d[2] 178 : n2 <- ceiling(nr / 2) 179 : print(head(m, max(1, n2))) 180 : cat("\n ..........\n\n") 181 : print(tail(m, max(1, nr - n2))) 182 : } 183 : ## DEBUG: cat("str(.):\n") ; str(x) 184 : invisible(x)# as print() S3 methods do 185 : } 186 : 187 : maechler 1575 nonFALSE <- function(x) { 188 : ## typically used for lMatrices: (TRUE,NA,FALSE) |-> (TRUE,FALSE) 189 : if(any(ix <- is.na(x))) x[ix] <- TRUE 190 : x 191 : } 192 : 193 : nz.NA <- function(x, na.value) { 194 : ## Non-Zeros of x 195 : ## na.value: TRUE: NA's give TRUE, they are not 0 196 : ## NA: NA's are not known ==> result := NA 197 : ## FALSE: NA's give FALSE, could be 0 198 : stopifnot(is.logical(na.value) && length(na.value) == 1) 199 : if(is.na(na.value)) x != 0 200 : else if(na.value) isN0(x) 201 : else x != 0 & !is.na(x) 202 : } 203 : 204 : maechler 1592 ## Number of non-zeros : 205 : ## FIXME? -- make this into a generic function (?) 206 : maechler 1575 nnzero <- function(x, na.counted = NA) { 207 : ## na.counted: TRUE: NA's are counted, they are not 0 208 : ## NA: NA's are not known (0 or not) ==> result := NA 209 : ## FALSE: NA's are omitted before counting 210 : maechler 1548 cl <- class(x) 211 : if(!extends(cl, "Matrix")) 212 : maechler 1575 sum(nz.NA(x, na.counted)) 213 : maechler 1548 else if(extends(cl, "sparseMatrix")) 214 : ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}! 215 : switch(.sp.class(cl), 216 : "CsparseMatrix" = length(x@i), 217 : "TsparseMatrix" = length(x@i), 218 : "RsparseMatrix" = length(x@j)) 219 : else ## denseMatrix 220 : maechler 1575 sum(nz.NA(as_geClass(x, cl)@x, na.counted)) 221 : maechler 1548 } 222 : 223 : maechler 919 ## For sparseness handling 224 : maechler 1467 ## return a 2-column (i,j) matrix of 225 : ## 0-based indices of non-zero entries : 226 : maechler 919 non0ind <- function(x) { 227 : maechler 1467 228 : maechler 919 if(is.numeric(x)) 229 : maechler 1575 return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0)) 230 : maechler 1226 ## else 231 : maechler 919 stopifnot(is(x, "sparseMatrix")) 232 : maechler 1319 non0.i <- function(M) { 233 : if(is(M, "TsparseMatrix")) 234 : return(unique(cbind(M@i,M@j))) 235 : if(is(M, "pMatrix")) 236 : maechler 1654 return(cbind(seq_len(nrow(M)), M@perm) - 1:1) 237 : maechler 1319 ## else: 238 : isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse) 239 : .Call(compressed_non_0_ij, M, isC) 240 : } 241 : 242 : if(is(x, "symmetricMatrix")) { # also get "other" triangle 243 : ij <- non0.i(x) 244 : notdiag <- ij[,1] != ij[,2]# but not the diagonals again 245 : rbind(ij, ij[notdiag, 2:1]) 246 : } 247 : maechler 1389 else if(is(x, "triangularMatrix")) { # check for "U" diag 248 : if(x@diag == "U") { 249 : maechler 1654 i <- seq_len(dim(x)[1]) - 1:1 250 : maechler 1389 rbind(non0.i(x), cbind(i,i)) 251 : } else non0.i(x) 252 : } 253 : maechler 1319 else 254 : non0.i(x) 255 : maechler 1226 } 256 : maechler 954 257 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings: 258 : ## Further, these map to and from the usual "Fortran-indexing" (but 0-based) 259 : encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr 260 : decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr) 261 : 262 : complementInd <- function(ij, dim) 263 : { 264 : ## Purpose: Compute the complement of the 2-column 0-based ij-matrix 265 : maechler 1238 ## but as 1-based indices 266 : maechler 1226 n <- prod(dim) 267 : if(n == 0) return(integer(0)) 268 : ii <- 1:n 269 : ii[-(1 + encodeInd(ij, nr = dim[1]))] 270 : maechler 919 } 271 : 272 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2)) 273 : 274 : intersectInd <- function(ij1, ij2, nrow) { 275 : ## from 2-column (i,j) matrices where i in {0,.., nrow-1}, 276 : ## return only the *common* entries 277 : decodeInd(intersect(encodeInd(ij1, nrow), 278 : encodeInd(ij2, nrow)), nrow) 279 : } 280 : 281 : WhichintersectInd <- function(ij1, ij2, nrow) { 282 : ## from 2-column (i,j) matrices where i \in {0,.., nrow-1}, 283 : ## find *where* common entries are in ij1 & ij2 284 : m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow)) 285 : ni <- !is.na(m1) 286 : list(which(ni), m1[ni]) 287 : } 288 : 289 : 290 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R ! 291 : maechler 1270 292 : maechler 1331 uniqTsparse <- function(x, class.x = c(class(x))) { 293 : ## Purpose: produce a *unique* triplet representation: 294 : ## by having (i,j) sorted and unique 295 : ## ----------------------------------------------------------- 296 : ## The following is not quite efficient {but easy to program, 297 : ## and as() are based on C code (all of them?) 298 : ## 299 : ## FIXME: Do it fast for the case where 'x' is already 'uniq' 300 : maechler 1270 301 : maechler 1331 switch(class.x, 302 : "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"), 303 : "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"), 304 : "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"), 305 : ## do we need this for "logical" ones, there's no sum() there! 306 : "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"), 307 : "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"), 308 : "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"), 309 : maechler 1548 ## do we need this for "logical" ones, there's no sum() there! 310 : "ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"), 311 : "nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"), 312 : "ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"), 313 : maechler 1331 ## otherwise: 314 : maechler 1472 stop("not yet implemented for class ", class.x)) 315 : maechler 1331 } 316 : maechler 919 317 : maechler 1331 ## Note: maybe, using 318 : ## ---- xj <- .Call(Matrix_expand_pointers, x@p) 319 : ## would be slightly more efficient than as( , "dgTMatrix") 320 : ## but really efficient would be to use only one .Call(.) for uniq(.) ! 321 : 322 : maechler 1654 drop0 <- function(x, clx = c(class(x))) { 323 : ## FIXME: Csparse_drop should do this (not losing symm./triang.): 324 : ## Careful: 'Csparse_drop' also drops triangularity,... 325 : ## .Call(Csparse_drop, as_CspClass(x, clx), 0) 326 : as_CspClass(.Call(Csparse_drop, as_CspClass(x, clx), 0.), 327 : clx) 328 : } 329 : 330 : maechler 1331 uniq <- function(x) { 331 : maechler 1654 if(is(x, "TsparseMatrix")) uniqTsparse(x) else 332 : if(is(x, "sparseMatrix")) drop0(x) else x 333 : maechler 919 } 334 : 335 : maechler 1472 asTuniq <- function(x) { 336 : if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix") 337 : } 338 : 339 : maechler 919 if(FALSE) ## try an "efficient" version 340 : uniq_gT <- function(x) 341 : { 342 : ## Purpose: produce a *unique* triplet representation: 343 : ## by having (i,j) sorted and unique 344 : maechler 1226 ## ------------------------------------------------------------------ 345 : maechler 919 ## Arguments: a "gT" Matrix 346 : stopifnot(is(x, "gTMatrix")) 347 : if((n <- length(x@i)) == 0) return(x) 348 : ii <- order(x@i, x@j) 349 : if(any(ii != 1:n)) { 350 : maechler 1238 x@i <- x@i[ii] 351 : x@j <- x@j[ii] 352 : x@x <- x@x[ii] 353 : maechler 919 } 354 : ij <- x@i + nrow(x) * x@j 355 : if(any(dup <- duplicated(ij))) { 356 : 357 : } 358 : ### We should use a .Call() based utility for this! 359 : 360 : } 361 : 362 : maechler 946 t_geMatrix <- function(x) { 363 : x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here 364 : x@Dim <- x@Dim[2:1] 365 : x@Dimnames <- x@Dimnames[2:1] 366 : ## FIXME: how to set factors? 367 : x 368 : } 369 : 370 : ## t( [dl]trMatrix ) and t( [dl]syMatrix ) : 371 : t_trMatrix <- function(x) { 372 : x@x <- as.vector(t(as(x, "matrix"))) 373 : x@Dim <- x@Dim[2:1] 374 : x@Dimnames <- x@Dimnames[2:1] 375 : x@uplo <- if (x@uplo == "U") "L" else "U" 376 : # and keep x@diag 377 : x 378 : } 379 : maechler 956 380 : fixupDense <- function(m, from) { 381 : if(is(m, "triangularMatrix")) { 382 : maechler 1238 m@uplo <- from@uplo 383 : m@diag <- from@diag 384 : maechler 956 } else if(is(m, "symmetricMatrix")) { 385 : maechler 1238 m@uplo <- from@uplo 386 : maechler 956 } 387 : m 388 : } 389 : 390 : ## -> ./ldenseMatrix.R : 391 : l2d_Matrix <- function(from) { 392 : stopifnot(is(from, "lMatrix")) 393 : fixupDense(new(sub("^l", "d", class(from)), 394 : maechler 1238 x = as.double(from@x), 395 : Dim = from@Dim, Dimnames = from@Dimnames), 396 : from) 397 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!} 398 : maechler 956 } 399 : 400 : maechler 1548 ## -> ./ndenseMatrix.R : 401 : n2d_Matrix <- function(from) { 402 : stopifnot(is(from, "nMatrix")) 403 : fixupDense(new(sub("^n", "d", class(from)), 404 : x = as.double(from@x), 405 : Dim = from@Dim, Dimnames = from@Dimnames), 406 : from) 407 : ## FIXME: treat 'factors' smartly {not for triangular!} 408 : } 409 : n2l_spMatrix <- function(from) { 410 : stopifnot(is(from, "nMatrix")) 411 : new(sub("^n", "l", class(from)), 412 : ##x = as.double(from@x), 413 : Dim = from@Dim, Dimnames = from@Dimnames) 414 : } 415 : 416 : maechler 956 if(FALSE)# unused 417 : l2d_meth <- function(x) { 418 : cl <- class(x) 419 : as(callGeneric(as(x, sub("^l", "d", cl))), cl) 420 : } 421 : 422 : maechler 1548 ## return "d" or "l" or "n" or "z" 423 : maechler 1331 .M.kind <- function(x, clx = class(x)) { 424 : if(is.matrix(x)) { ## 'old style matrix' 425 : if (is.numeric(x)) "d" 426 : maechler 1548 else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ?? 427 : maechler 1331 else if(is.complex(x)) "z" 428 : else stop("not yet implemented for matrix w/ typeof ", typeof(x)) 429 : } 430 : else if(extends(clx, "dMatrix")) "d" 431 : maechler 1548 else if(extends(clx, "nMatrix")) "n" 432 : maechler 1331 else if(extends(clx, "lMatrix")) "l" 433 : else if(extends(clx, "zMatrix")) "z" 434 : maechler 1548 else if(extends(clx, "pMatrix")) "n" # permutation -> pattern 435 : maechler 1331 else stop(" not yet be implemented for ", clx) 436 : } 437 : 438 : .M.shape <- function(x, clx = class(x)) { 439 : if(is.matrix(x)) { ## 'old style matrix' 440 : if (isDiagonal (x)) "d" 441 : else if(isTriangular(x)) "t" 442 : else if(isSymmetric (x)) "s" 443 : else "g" # general 444 : } 445 : else if(extends(clx, "diagonalMatrix")) "d" 446 : else if(extends(clx, "triangularMatrix"))"t" 447 : else if(extends(clx, "symmetricMatrix")) "s" 448 : else "g" 449 : } 450 : 451 : 452 : maechler 1329 class2 <- function(cl, kind = "l", do.sub = TRUE) { 453 : ## Find "corresponding" class; since pos.def. matrices have no pendant: 454 : if (cl == "dpoMatrix") paste(kind, "syMatrix", sep='') 455 : else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='') 456 : else if(do.sub) sub("^d", kind, cl) 457 : else cl 458 : maechler 1226 } 459 : 460 : geClass <- function(x) { 461 : maechler 1331 if (is(x, "dMatrix")) "dgeMatrix" 462 : maechler 1226 else if(is(x, "lMatrix")) "lgeMatrix" 463 : maechler 1548 else if(is(x, "nMatrix")) "ngeMatrix" 464 : maechler 1331 else if(is(x, "zMatrix")) "zgeMatrix" 465 : maechler 1329 else stop("general Matrix class not yet implemented for ", 466 : maechler 1226 class(x)) 467 : } 468 : maechler 1331 469 : .dense.prefixes <- c("d" = "di", 470 : "t" = "tr", 471 : "s" = "sy", 472 : "g" = "ge") 473 : 474 : maechler 1349 .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular 475 : "t" = "t", 476 : "s" = "s", 477 : "g" = "g") 478 : maechler 1331 479 : maechler 1329 ## Used, e.g. after subsetting: Try to use specific class -- if feasible : 480 : maechler 1331 as_dense <- function(x) { 481 : as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep='')) 482 : } 483 : 484 : maechler 1548 .sp.class <- function(x) { ## find and return the "sparseness class" 485 : if(!is.character(x)) x <- class(x) 486 : for(cl in paste(c("C","T","R"), "sparseMatrix", sep='')) 487 : if(extends(x, cl)) 488 : return(cl) 489 : ## else (should rarely happen) 490 : as.character(NA) 491 : } 492 : 493 : maechler 1331 as_Csparse <- function(x) { 494 : maechler 1349 as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep='')) 495 : maechler 1331 } 496 : maechler 1349 as_Rsparse <- function(x) { 497 : as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep='')) 498 : } 499 : as_Tsparse <- function(x) { 500 : as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep='')) 501 : } 502 : maechler 1331 503 : maechler 1329 as_geClass <- function(x, cl) { 504 : maechler 1331 if (extends(cl, "diagonalMatrix") && isDiagonal(x)) 505 : as(x, cl) 506 : maechler 1357 else if(extends(cl, "symmetricMatrix") && isSymmetric(x)) { 507 : kind <- .M.kind(x) 508 : maechler 1331 as(x, class2(cl, kind, do.sub= kind != "d")) 509 : maechler 1357 } else if(extends(cl, "triangularMatrix") && isTriangular(x)) 510 : maechler 1331 as(x, cl) 511 : else 512 : as(x, paste(.M.kind(x), "geMatrix", sep='')) 513 : } 514 : maechler 1226 515 : maechler 1331 as_CspClass <- function(x, cl) { 516 : maechler 1654 if (## diagonal is *not* sparse: 517 : ##(extends(cl, "diagonalMatrix") && isDiagonal(x)) || 518 : (extends(cl, "symmetricMatrix") && isSymmetric(x)) || 519 : maechler 1331 (extends(cl, "triangularMatrix")&& isTriangular(x))) 520 : as(x, cl) 521 : maechler 1654 else if(is(x, "CsparseMatrix")) x 522 : maechler 1331 else as(x, paste(.M.kind(x), "gCMatrix", sep='')) 523 : maechler 1329 } 524 : 525 : maechler 1331 526 : maechler 956 ## -> ./ddenseMatrix.R : 527 : d2l_Matrix <- function(from) { 528 : stopifnot(is(from, "dMatrix")) 529 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here 530 : maechler 1238 Dim = from@Dim, Dimnames = from@Dimnames), 531 : from) 532 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!} 533 : maechler 956 } 534 : maechler 973 535 : 536 : try_as <- function(x, classes, tryAnyway = FALSE) { 537 : if(!tryAnyway && !is(x, "Matrix")) 538 : maechler 1238 return(x) 539 : maechler 973 ## else 540 : ok <- canCoerce(x, classes[1]) 541 : while(!ok && length(classes <- classes[-1])) { 542 : maechler 1238 ok <- canCoerce(x, classes[1]) 543 : maechler 973 } 544 : if(ok) as(x, classes[1]) else x 545 : } 546 : 547 : 548 : maechler 1238 ## For *dense* matrices 549 : isTriMat <- function(object, upper = NA) { 550 : maechler 1174 ## pretest: is it square? 551 : d <- dim(object) 552 : if(d[1] != d[2]) return(FALSE) 553 : ## else slower test 554 : if(!is.matrix(object)) 555 : maechler 1226 object <- as(object,"matrix") 556 : if(is.na(upper)) { 557 : maechler 1467 if(all0(object[lower.tri(object)])) 558 : maechler 1226 structure(TRUE, kind = "U") 559 : maechler 1467 else if(all0(object[upper.tri(object)])) 560 : maechler 1226 structure(TRUE, kind = "L") 561 : else FALSE 562 : } else if(upper) 563 : maechler 1467 all0(object[lower.tri(object)]) 564 : maechler 1226 else ## upper is FALSE 565 : maechler 1467 all0(object[upper.tri(object)]) 566 : maechler 1174 } 567 : 568 : maechler 1238 ## For Csparse matrices 569 : isTriC <- function(x, upper = NA) { 570 : ## pretest: is it square? 571 : d <- dim(x) 572 : if(d[1] != d[2]) return(FALSE) 573 : ## else 574 : if(d[1] == 0) return(TRUE) 575 : maechler 1315 ni <- 1:d[2] 576 : maechler 1238 ## the row indices split according to column: 577 : ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni)) 578 : lil <- unlist(lapply(ilist, length), use.names = FALSE) 579 : if(any(lil == 0)) { 580 : pos <- lil > 0 581 : if(!any(pos)) ## matrix of all 0's 582 : return(TRUE) 583 : ilist <- ilist[pos] 584 : ni <- ni[pos] 585 : } 586 : maechler 1315 ni0 <- ni - 1:1 # '0-based ni' 587 : maechler 1238 if(is.na(upper)) { 588 : maechler 1315 if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)) 589 : maechler 1238 structure(TRUE, kind = "U") 590 : maechler 1315 else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)) 591 : maechler 1238 structure(TRUE, kind = "L") 592 : else FALSE 593 : } else if(upper) { 594 : maechler 1315 all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0) 595 : maechler 1238 } else { ## 'lower' 596 : maechler 1315 all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0) 597 : maechler 1238 } 598 : } 599 : 600 : maechler 1174 .is.diagonal <- function(object) { 601 : maechler 1357 ## "matrix" or "denseMatrix" (but not "diagonalMatrix") 602 : maechler 1174 d <- dim(object) 603 : if(d[1] != (n <- d[2])) FALSE 604 : maechler 1357 else if(is.matrix(object)) 605 : ## requires that "vector-indexing" works for 'object' : 606 : maechler 1467 all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]) 607 : maechler 1357 else ## "denseMatrix" -- packed or unpacked 608 : if(is(object, "generalMatrix")) # "dge", "lge", ... 609 : maechler 1467 all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]) 610 : maechler 1357 else { ## "dense" but not {diag, general}, i.e. triangular or symmetric: 611 : ## -> has 'uplo' differentiate between packed and unpacked 612 : 613 : ### .......... FIXME ............... 614 : 615 : packed <- isPacked(object) 616 : if(object@uplo == "U") { 617 : } else { ## uplo == "L" 618 : } 619 : 620 : ### very cheap workaround 621 : maechler 1467 all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]) 622 : maechler 1357 } 623 : maechler 1174 } 624 : maechler 1253 625 : maechler 1357 626 : maechler 1592 ## FIXME? -- this should also work for "ltT", "ntT", ... : 627 : maechler 1253 diagU2N <- function(x) 628 : { 629 : maechler 1592 ## Purpose: Transform a *unit diagonal* sparse triangular matrix 630 : maechler 1253 ## into one with explicit diagonal entries '1' 631 : xT <- as(x, "dgTMatrix") 632 : ## leave it as T* - the caller can always coerce to C* if needed: 633 : new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim, 634 : Dimnames = x@Dimnames, uplo = x@uplo, diag = "N") 635 : } 636 : maechler 1290 637 : maechler 1659 ## Needed, e.g., in ./Csparse.R for colSums() etc: 638 : maechler 1295 .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) { 639 : x <- as(x, "dgCMatrix") 640 : callGeneric() 641 : } 642 : 643 : maechler 1659 .as.dgC.0.factors <- function(x) { 644 : if(!is(x, "dgCMatrix")) 645 : as(x, "dgCMatrix") # will not have 'factors' 646 : else ## dgCMatrix 647 : if(!length(x@factors)) x else { x@factors <- list() ; x } 648 : } 649 : 650 : maechler 1295 .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) { 651 : maechler 1349 ## used e.g. inside colSums() etc methods 652 : maechler 1295 x <- as(x, "dgTMatrix") 653 : callGeneric() 654 : } 655 : 656 : 657 : maechler 1290 ### Fast much simplified version of tapply() 658 : tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) { 659 : sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE) 660 : } 661 : 662 : ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) { 663 : ## tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify) 664 : ## } 665 :