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 : stop(gettextf('not-yet-implemented method for %s(<%s>)', fun, cl), 26 : maechler 1238 call. = FALSE) 27 : maechler 949 } 28 : .bail.out.2 <- function(fun, cl1, cl2) { 29 : stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)', 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 919 prTriang <- function(x, digits = getOption("digits"), 125 : maechler 1389 maxp = getOption("max.print"), 126 : maechler 1238 justify = "none", right = TRUE) 127 : maechler 919 { 128 : ## modeled along stats:::print.dist 129 : upper <- x@uplo == "U" 130 : 131 : m <- as(x, "matrix") 132 : cf <- format(m, digits = digits, justify = justify) 133 : if(upper) 134 : maechler 1238 cf[row(cf) > col(cf)] <- "." 135 : maechler 919 else 136 : maechler 1238 cf[row(cf) < col(cf)] <- "." 137 : maechler 1472 if(.isR_24) 138 : print(cf, quote = FALSE, right = right, max = maxp) 139 : else print(cf, quote = FALSE, right = right) 140 : maechler 919 invisible(x) 141 : } 142 : 143 : maechler 1389 prMatrix <- function(x, digits = getOption("digits"), 144 : maxp = getOption("max.print")) { 145 : maechler 919 d <- dim(x) 146 : cl <- class(x) 147 : cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl)) 148 : if(prod(d) <= maxp) { 149 : maechler 1238 if(is(x, "triangularMatrix")) 150 : maechler 1472 prTriang(x, digits = digits, maxp = maxp) 151 : else { 152 : if(.isR_24) 153 : print(as(x, "matrix"), digits = digits, max = maxp) 154 : else print(as(x, "matrix"), digits = digits) 155 : } 156 : maechler 919 } 157 : else { ## d[1] > maxp / d[2] >= nr : 158 : maechler 1238 m <- as(x, "matrix") 159 : maechler 919 nr <- maxp %/% d[2] 160 : n2 <- ceiling(nr / 2) 161 : print(head(m, max(1, n2))) 162 : cat("\n ..........\n\n") 163 : print(tail(m, max(1, nr - n2))) 164 : } 165 : ## DEBUG: cat("str(.):\n") ; str(x) 166 : invisible(x)# as print() S3 methods do 167 : } 168 : 169 : maechler 1575 nonFALSE <- function(x) { 170 : ## typically used for lMatrices: (TRUE,NA,FALSE) |-> (TRUE,FALSE) 171 : if(any(ix <- is.na(x))) x[ix] <- TRUE 172 : x 173 : } 174 : 175 : nz.NA <- function(x, na.value) { 176 : ## Non-Zeros of x 177 : ## na.value: TRUE: NA's give TRUE, they are not 0 178 : ## NA: NA's are not known ==> result := NA 179 : ## FALSE: NA's give FALSE, could be 0 180 : stopifnot(is.logical(na.value) && length(na.value) == 1) 181 : if(is.na(na.value)) x != 0 182 : else if(na.value) isN0(x) 183 : else x != 0 & !is.na(x) 184 : } 185 : 186 : maechler 1548 ### FIXME? -- make this into a generic function (?) 187 : maechler 1575 nnzero <- function(x, na.counted = NA) { 188 : ## na.counted: TRUE: NA's are counted, they are not 0 189 : ## NA: NA's are not known (0 or not) ==> result := NA 190 : ## FALSE: NA's are omitted before counting 191 : maechler 1548 cl <- class(x) 192 : if(!extends(cl, "Matrix")) 193 : maechler 1575 sum(nz.NA(x, na.counted)) 194 : maechler 1548 else if(extends(cl, "sparseMatrix")) 195 : ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}! 196 : switch(.sp.class(cl), 197 : "CsparseMatrix" = length(x@i), 198 : "TsparseMatrix" = length(x@i), 199 : "RsparseMatrix" = length(x@j)) 200 : else ## denseMatrix 201 : maechler 1575 sum(nz.NA(as_geClass(x, cl)@x, na.counted)) 202 : maechler 1548 } 203 : 204 : maechler 919 ## For sparseness handling 205 : maechler 1467 ## return a 2-column (i,j) matrix of 206 : ## 0-based indices of non-zero entries : 207 : maechler 919 non0ind <- function(x) { 208 : maechler 1467 209 : maechler 919 if(is.numeric(x)) 210 : maechler 1575 return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0)) 211 : maechler 1226 ## else 212 : maechler 919 stopifnot(is(x, "sparseMatrix")) 213 : maechler 1319 non0.i <- function(M) { 214 : if(is(M, "TsparseMatrix")) 215 : return(unique(cbind(M@i,M@j))) 216 : if(is(M, "pMatrix")) 217 : return(cbind(seq(length=nrow(M)), M@perm) - 1:1) 218 : ## else: 219 : isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse) 220 : .Call(compressed_non_0_ij, M, isC) 221 : } 222 : 223 : if(is(x, "symmetricMatrix")) { # also get "other" triangle 224 : ij <- non0.i(x) 225 : notdiag <- ij[,1] != ij[,2]# but not the diagonals again 226 : rbind(ij, ij[notdiag, 2:1]) 227 : } 228 : maechler 1389 else if(is(x, "triangularMatrix")) { # check for "U" diag 229 : if(x@diag == "U") { 230 : maechler 1390 i <- seq(length = dim(x)[1]) - 1:1 231 : maechler 1389 rbind(non0.i(x), cbind(i,i)) 232 : } else non0.i(x) 233 : } 234 : maechler 1319 else 235 : non0.i(x) 236 : maechler 1226 } 237 : maechler 954 238 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings: 239 : ## Further, these map to and from the usual "Fortran-indexing" (but 0-based) 240 : encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr 241 : decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr) 242 : 243 : complementInd <- function(ij, dim) 244 : { 245 : ## Purpose: Compute the complement of the 2-column 0-based ij-matrix 246 : maechler 1238 ## but as 1-based indices 247 : maechler 1226 n <- prod(dim) 248 : if(n == 0) return(integer(0)) 249 : ii <- 1:n 250 : ii[-(1 + encodeInd(ij, nr = dim[1]))] 251 : maechler 919 } 252 : 253 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2)) 254 : 255 : intersectInd <- function(ij1, ij2, nrow) { 256 : ## from 2-column (i,j) matrices where i in {0,.., nrow-1}, 257 : ## return only the *common* entries 258 : decodeInd(intersect(encodeInd(ij1, nrow), 259 : encodeInd(ij2, nrow)), nrow) 260 : } 261 : 262 : WhichintersectInd <- function(ij1, ij2, nrow) { 263 : ## from 2-column (i,j) matrices where i \in {0,.., nrow-1}, 264 : ## find *where* common entries are in ij1 & ij2 265 : m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow)) 266 : ni <- !is.na(m1) 267 : list(which(ni), m1[ni]) 268 : } 269 : 270 : 271 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R ! 272 : maechler 1270 273 : maechler 1331 uniqTsparse <- function(x, class.x = c(class(x))) { 274 : ## Purpose: produce a *unique* triplet representation: 275 : ## by having (i,j) sorted and unique 276 : ## ----------------------------------------------------------- 277 : ## The following is not quite efficient {but easy to program, 278 : ## and as() are based on C code (all of them?) 279 : ## 280 : ## FIXME: Do it fast for the case where 'x' is already 'uniq' 281 : maechler 1270 282 : maechler 1331 switch(class.x, 283 : "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"), 284 : "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"), 285 : "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"), 286 : ## do we need this for "logical" ones, there's no sum() there! 287 : "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"), 288 : "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"), 289 : "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"), 290 : maechler 1548 ## do we need this for "logical" ones, there's no sum() there! 291 : "ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"), 292 : "nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"), 293 : "ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"), 294 : maechler 1331 ## otherwise: 295 : maechler 1472 stop("not yet implemented for class ", class.x)) 296 : maechler 1331 } 297 : maechler 919 298 : maechler 1331 ## Note: maybe, using 299 : ## ---- xj <- .Call(Matrix_expand_pointers, x@p) 300 : ## would be slightly more efficient than as( , "dgTMatrix") 301 : ## but really efficient would be to use only one .Call(.) for uniq(.) ! 302 : 303 : uniq <- function(x) { 304 : if(is(x, "TsparseMatrix")) uniqTsparse(x) else x 305 : ## else: not 'Tsparse', i.e. "uniquely" represented in any case 306 : maechler 919 } 307 : 308 : maechler 1472 asTuniq <- function(x) { 309 : if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix") 310 : } 311 : 312 : maechler 919 if(FALSE) ## try an "efficient" version 313 : uniq_gT <- function(x) 314 : { 315 : ## Purpose: produce a *unique* triplet representation: 316 : ## by having (i,j) sorted and unique 317 : maechler 1226 ## ------------------------------------------------------------------ 318 : maechler 919 ## Arguments: a "gT" Matrix 319 : stopifnot(is(x, "gTMatrix")) 320 : if((n <- length(x@i)) == 0) return(x) 321 : ii <- order(x@i, x@j) 322 : if(any(ii != 1:n)) { 323 : maechler 1238 x@i <- x@i[ii] 324 : x@j <- x@j[ii] 325 : x@x <- x@x[ii] 326 : maechler 919 } 327 : ij <- x@i + nrow(x) * x@j 328 : if(any(dup <- duplicated(ij))) { 329 : 330 : } 331 : ### We should use a .Call() based utility for this! 332 : 333 : } 334 : 335 : maechler 946 t_geMatrix <- function(x) { 336 : x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here 337 : x@Dim <- x@Dim[2:1] 338 : x@Dimnames <- x@Dimnames[2:1] 339 : ## FIXME: how to set factors? 340 : x 341 : } 342 : 343 : ## t( [dl]trMatrix ) and t( [dl]syMatrix ) : 344 : t_trMatrix <- function(x) { 345 : x@x <- as.vector(t(as(x, "matrix"))) 346 : x@Dim <- x@Dim[2:1] 347 : x@Dimnames <- x@Dimnames[2:1] 348 : x@uplo <- if (x@uplo == "U") "L" else "U" 349 : # and keep x@diag 350 : x 351 : } 352 : maechler 956 353 : fixupDense <- function(m, from) { 354 : if(is(m, "triangularMatrix")) { 355 : maechler 1238 m@uplo <- from@uplo 356 : m@diag <- from@diag 357 : maechler 956 } else if(is(m, "symmetricMatrix")) { 358 : maechler 1238 m@uplo <- from@uplo 359 : maechler 956 } 360 : m 361 : } 362 : 363 : ## -> ./ldenseMatrix.R : 364 : l2d_Matrix <- function(from) { 365 : stopifnot(is(from, "lMatrix")) 366 : fixupDense(new(sub("^l", "d", class(from)), 367 : maechler 1238 x = as.double(from@x), 368 : Dim = from@Dim, Dimnames = from@Dimnames), 369 : from) 370 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!} 371 : maechler 956 } 372 : 373 : maechler 1548 ## -> ./ndenseMatrix.R : 374 : n2d_Matrix <- function(from) { 375 : stopifnot(is(from, "nMatrix")) 376 : fixupDense(new(sub("^n", "d", class(from)), 377 : x = as.double(from@x), 378 : Dim = from@Dim, Dimnames = from@Dimnames), 379 : from) 380 : ## FIXME: treat 'factors' smartly {not for triangular!} 381 : } 382 : n2l_spMatrix <- function(from) { 383 : stopifnot(is(from, "nMatrix")) 384 : new(sub("^n", "l", class(from)), 385 : ##x = as.double(from@x), 386 : Dim = from@Dim, Dimnames = from@Dimnames) 387 : } 388 : 389 : maechler 956 if(FALSE)# unused 390 : l2d_meth <- function(x) { 391 : cl <- class(x) 392 : as(callGeneric(as(x, sub("^l", "d", cl))), cl) 393 : } 394 : 395 : maechler 1548 ## return "d" or "l" or "n" or "z" 396 : maechler 1331 .M.kind <- function(x, clx = class(x)) { 397 : if(is.matrix(x)) { ## 'old style matrix' 398 : if (is.numeric(x)) "d" 399 : maechler 1548 else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ?? 400 : maechler 1331 else if(is.complex(x)) "z" 401 : else stop("not yet implemented for matrix w/ typeof ", typeof(x)) 402 : } 403 : else if(extends(clx, "dMatrix")) "d" 404 : maechler 1548 else if(extends(clx, "nMatrix")) "n" 405 : maechler 1331 else if(extends(clx, "lMatrix")) "l" 406 : else if(extends(clx, "zMatrix")) "z" 407 : maechler 1548 else if(extends(clx, "pMatrix")) "n" # permutation -> pattern 408 : maechler 1331 else stop(" not yet be implemented for ", clx) 409 : } 410 : 411 : .M.shape <- function(x, clx = class(x)) { 412 : if(is.matrix(x)) { ## 'old style matrix' 413 : if (isDiagonal (x)) "d" 414 : else if(isTriangular(x)) "t" 415 : else if(isSymmetric (x)) "s" 416 : else "g" # general 417 : } 418 : else if(extends(clx, "diagonalMatrix")) "d" 419 : else if(extends(clx, "triangularMatrix"))"t" 420 : else if(extends(clx, "symmetricMatrix")) "s" 421 : else "g" 422 : } 423 : 424 : 425 : maechler 1329 class2 <- function(cl, kind = "l", do.sub = TRUE) { 426 : ## Find "corresponding" class; since pos.def. matrices have no pendant: 427 : if (cl == "dpoMatrix") paste(kind, "syMatrix", sep='') 428 : else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='') 429 : else if(do.sub) sub("^d", kind, cl) 430 : else cl 431 : maechler 1226 } 432 : 433 : geClass <- function(x) { 434 : maechler 1331 if (is(x, "dMatrix")) "dgeMatrix" 435 : maechler 1226 else if(is(x, "lMatrix")) "lgeMatrix" 436 : maechler 1548 else if(is(x, "nMatrix")) "ngeMatrix" 437 : maechler 1331 else if(is(x, "zMatrix")) "zgeMatrix" 438 : maechler 1329 else stop("general Matrix class not yet implemented for ", 439 : maechler 1226 class(x)) 440 : } 441 : maechler 1331 442 : .dense.prefixes <- c("d" = "di", 443 : "t" = "tr", 444 : "s" = "sy", 445 : "g" = "ge") 446 : 447 : maechler 1349 .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular 448 : "t" = "t", 449 : "s" = "s", 450 : "g" = "g") 451 : maechler 1331 452 : maechler 1329 ## Used, e.g. after subsetting: Try to use specific class -- if feasible : 453 : maechler 1331 as_dense <- function(x) { 454 : as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep='')) 455 : } 456 : 457 : maechler 1548 .sp.class <- function(x) { ## find and return the "sparseness class" 458 : if(!is.character(x)) x <- class(x) 459 : for(cl in paste(c("C","T","R"), "sparseMatrix", sep='')) 460 : if(extends(x, cl)) 461 : return(cl) 462 : ## else (should rarely happen) 463 : as.character(NA) 464 : } 465 : 466 : maechler 1331 as_Csparse <- function(x) { 467 : maechler 1349 as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep='')) 468 : maechler 1331 } 469 : maechler 1349 as_Rsparse <- function(x) { 470 : as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep='')) 471 : } 472 : as_Tsparse <- function(x) { 473 : as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep='')) 474 : } 475 : maechler 1331 476 : maechler 1329 as_geClass <- function(x, cl) { 477 : maechler 1331 if (extends(cl, "diagonalMatrix") && isDiagonal(x)) 478 : as(x, cl) 479 : maechler 1357 else if(extends(cl, "symmetricMatrix") && isSymmetric(x)) { 480 : kind <- .M.kind(x) 481 : maechler 1331 as(x, class2(cl, kind, do.sub= kind != "d")) 482 : maechler 1357 } else if(extends(cl, "triangularMatrix") && isTriangular(x)) 483 : maechler 1331 as(x, cl) 484 : else 485 : as(x, paste(.M.kind(x), "geMatrix", sep='')) 486 : } 487 : maechler 1226 488 : maechler 1331 as_CspClass <- function(x, cl) { 489 : if ((extends(cl, "diagonalMatrix") && isDiagonal(x)) || 490 : (extends(cl, "symmetricMatrix") && isSymmetric(x)) || 491 : (extends(cl, "triangularMatrix")&& isTriangular(x))) 492 : as(x, cl) 493 : else as(x, paste(.M.kind(x), "gCMatrix", sep='')) 494 : maechler 1329 } 495 : 496 : maechler 1331 497 : maechler 956 ## -> ./ddenseMatrix.R : 498 : d2l_Matrix <- function(from) { 499 : stopifnot(is(from, "dMatrix")) 500 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here 501 : maechler 1238 Dim = from@Dim, Dimnames = from@Dimnames), 502 : from) 503 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!} 504 : maechler 956 } 505 : maechler 973 506 : 507 : try_as <- function(x, classes, tryAnyway = FALSE) { 508 : if(!tryAnyway && !is(x, "Matrix")) 509 : maechler 1238 return(x) 510 : maechler 973 ## else 511 : ok <- canCoerce(x, classes[1]) 512 : while(!ok && length(classes <- classes[-1])) { 513 : maechler 1238 ok <- canCoerce(x, classes[1]) 514 : maechler 973 } 515 : if(ok) as(x, classes[1]) else x 516 : } 517 : 518 : 519 : maechler 1238 ## For *dense* matrices 520 : isTriMat <- function(object, upper = NA) { 521 : maechler 1174 ## pretest: is it square? 522 : d <- dim(object) 523 : if(d[1] != d[2]) return(FALSE) 524 : ## else slower test 525 : if(!is.matrix(object)) 526 : maechler 1226 object <- as(object,"matrix") 527 : if(is.na(upper)) { 528 : maechler 1467 if(all0(object[lower.tri(object)])) 529 : maechler 1226 structure(TRUE, kind = "U") 530 : maechler 1467 else if(all0(object[upper.tri(object)])) 531 : maechler 1226 structure(TRUE, kind = "L") 532 : else FALSE 533 : } else if(upper) 534 : maechler 1467 all0(object[lower.tri(object)]) 535 : maechler 1226 else ## upper is FALSE 536 : maechler 1467 all0(object[upper.tri(object)]) 537 : maechler 1174 } 538 : 539 : maechler 1238 ## For Csparse matrices 540 : isTriC <- function(x, upper = NA) { 541 : ## pretest: is it square? 542 : d <- dim(x) 543 : if(d[1] != d[2]) return(FALSE) 544 : ## else 545 : if(d[1] == 0) return(TRUE) 546 : maechler 1315 ni <- 1:d[2] 547 : maechler 1238 ## the row indices split according to column: 548 : ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni)) 549 : lil <- unlist(lapply(ilist, length), use.names = FALSE) 550 : if(any(lil == 0)) { 551 : pos <- lil > 0 552 : if(!any(pos)) ## matrix of all 0's 553 : return(TRUE) 554 : ilist <- ilist[pos] 555 : ni <- ni[pos] 556 : } 557 : maechler 1315 ni0 <- ni - 1:1 # '0-based ni' 558 : maechler 1238 if(is.na(upper)) { 559 : maechler 1315 if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)) 560 : maechler 1238 structure(TRUE, kind = "U") 561 : maechler 1315 else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)) 562 : maechler 1238 structure(TRUE, kind = "L") 563 : else FALSE 564 : } else if(upper) { 565 : maechler 1315 all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0) 566 : maechler 1238 } else { ## 'lower' 567 : maechler 1315 all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0) 568 : maechler 1238 } 569 : } 570 : 571 : maechler 1174 .is.diagonal <- function(object) { 572 : maechler 1357 ## "matrix" or "denseMatrix" (but not "diagonalMatrix") 573 : maechler 1174 d <- dim(object) 574 : if(d[1] != (n <- d[2])) FALSE 575 : maechler 1357 else if(is.matrix(object)) 576 : ## requires that "vector-indexing" works for 'object' : 577 : maechler 1467 all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]) 578 : maechler 1357 else ## "denseMatrix" -- packed or unpacked 579 : if(is(object, "generalMatrix")) # "dge", "lge", ... 580 : maechler 1467 all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]) 581 : maechler 1357 else { ## "dense" but not {diag, general}, i.e. triangular or symmetric: 582 : ## -> has 'uplo' differentiate between packed and unpacked 583 : 584 : ### .......... FIXME ............... 585 : 586 : packed <- isPacked(object) 587 : if(object@uplo == "U") { 588 : } else { ## uplo == "L" 589 : } 590 : 591 : ### very cheap workaround 592 : maechler 1467 all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)]) 593 : maechler 1357 } 594 : maechler 1174 } 595 : maechler 1253 596 : maechler 1357 597 : maechler 1253 diagU2N <- function(x) 598 : { 599 : ## Purpose: Transform a *unit diagonal* triangular matrix 600 : ## into one with explicit diagonal entries '1' 601 : xT <- as(x, "dgTMatrix") 602 : ## leave it as T* - the caller can always coerce to C* if needed: 603 : new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim, 604 : Dimnames = x@Dimnames, uplo = x@uplo, diag = "N") 605 : } 606 : maechler 1290 607 : maechler 1349 ## FIXME: this should probably be dropped / replaced by as_Csparse 608 : maechler 1295 .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) { 609 : x <- as(x, "dgCMatrix") 610 : callGeneric() 611 : } 612 : 613 : .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) { 614 : maechler 1349 ## used e.g. inside colSums() etc methods 615 : maechler 1295 x <- as(x, "dgTMatrix") 616 : callGeneric() 617 : } 618 : 619 : 620 : maechler 1290 ### Fast much simplified version of tapply() 621 : tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) { 622 : sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE) 623 : } 624 : 625 : ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) { 626 : ## tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify) 627 : ## } 628 :