SCM Repository

[matrix] Annotation of /pkg/Matrix/R/diagMatrix.R
 [matrix] / pkg / Matrix / R / diagMatrix.R

Annotation of /pkg/Matrix/R/diagMatrix.R

Revision 2840 - (view) (download)

 1 : maechler 1617 #### All methods for "diagonalMatrix" and its subclasses, 2 : #### currently "ddiMatrix", "ldiMatrix" 3 : 4 : maechler 1109 ## Purpose: Constructor of diagonal matrices -- ~= diag() , 5 : ## but *not* diag() extractor! 6 : Diagonal <- function(n, x = NULL) 7 : { 8 : mmaechler 2820 ## Allow Diagonal(4), Diagonal(x=1:5), and Diagonal(4, TRUE) 9 : maechler 1109 if(missing(n)) 10 : maechler 1575 n <- length(x) 11 : maechler 1109 else { 12 : maechler 1575 stopifnot(length(n) == 1, n == as.integer(n), n >= 0) 13 : n <- as.integer(n) 14 : maechler 1109 } 15 : 16 : maechler 1654 if(missing(x)) ## unit diagonal matrix 17 : maechler 1575 new("ddiMatrix", Dim = c(n,n), diag = "U") 18 : maechler 1109 else { 19 : maechler 2128 lx <- length(x) 20 : mmaechler 2820 lx.1 <- lx == 1L 21 : stopifnot(lx.1 || lx == n) # but keep 'x' short for now 22 : maechler 1575 if(is.logical(x)) 23 : cl <- "ldiMatrix" 24 : else if(is.numeric(x)) { 25 : cl <- "ddiMatrix" 26 : x <- as.numeric(x) 27 : } 28 : else if(is.complex(x)) { 29 : cl <- "zdiMatrix" # will not yet work 30 : } else stop("'x' has invalid data type") 31 : mmaechler 2820 if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal.. 32 : new(cl, Dim = c(n,n), diag = "U") 33 : else 34 : new(cl, Dim = c(n,n), diag = "N", 35 : x = if(lx.1) rep.int(x,n) else x) 36 : maechler 1109 } 37 : } 38 : 39 : mmaechler 2840 .sparseDiagonal <- function(n, x = 1, uplo = "U", 40 : mmaechler 2508 shape = if(missing(cols)) "t" else "g", 41 : mmaechler 2840 unitri, kind, 42 : cols = if(n) 0:(n - 1L) else integer(0)) 43 : mmaechler 2508 { 44 : maechler 2144 stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0) 45 : mmaechler 2840 if(!(mcols <- missing(cols))) 46 : mmaechler 2508 stopifnot(0 <= (cols <- as.integer(cols)), cols < n) 47 : m <- length(cols) 48 : if(missing(kind)) 49 : kind <- 50 : if(is.double(x)) "d" 51 : else if(is.logical(x)) "l" 52 : else { ## for now 53 : storage.mode(x) <- "double" 54 : "d" 55 : } 56 : else stopifnot(any(kind == c("d","l","n"))) 57 : mmaechler 2226 stopifnot(is.character(shape), nchar(shape) == 1, 58 : any(shape == c("t","s","g"))) # triangular / symmetric / general 59 : mmaechler 2840 if((missing(unitri) || unitri) && shape == "t" && 60 : (mcols || cols == 0:(n-1L)) && 61 : ((any(kind == c("l", "n")) && allTrue(x)) || 62 : ( kind == "d" && allTrue(x == 1)))) { ## uni-triangular 63 : new(paste0(kind,"tCMatrix"), Dim = c(n,n), 64 : uplo = uplo, diag = "U", p = rep.int(0L, n+1L)) 65 : } 66 : else if(kind == "n") { 67 : mmaechler 2508 if(shape == "g") 68 : new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m) 69 : else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo, 70 : i = cols, p = 0:m) 71 : } 72 : mmaechler 2840 else { ## kind != "n" -- have x slot : 73 : if((lx <- length(x)) == 1) x <- rep.int(x, m) 74 : else if(lx != m) stop("length(x) must be either 1 or #{cols}") 75 : if(shape == "g") 76 : new(paste0(kind, "gCMatrix"), Dim = c(n,m), 77 : x = x, i = cols, p = 0:m) 78 : else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo, 79 : x = x, i = cols, p = 0:m) 80 : } 81 : maechler 2144 } 82 : 83 : mmaechler 2226 ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I() 84 : .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") 85 : .sparseDiagonal(n, x, uplo, shape = "s") 86 : 87 : mmaechler 2608 # instead of diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case: 88 : mmaechler 2840 .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE) 89 : .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri) 90 : mmaechler 2226 91 : 92 : mmaechler 2764 ## This is modified from a post of Bert Gunter to R-help on 1 Sep 2005. 93 : ## Bert's code built on a post by Andy Liaw who most probably was influenced 94 : ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002 95 : ## who posted his bdiag() function written in December 1995. 96 : mmaechler 2209 if(FALSE)##--- no longer used: 97 : mmaechler 2206 .bdiag <- function(lst) { 98 : mmaechler 2764 ## block-diagonal matrix [a dgTMatrix] from list of matrices 99 : mmaechler 2206 stopifnot(is.list(lst), length(lst) >= 1) 100 : mmaechler 2813 dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE) 101 : maechler 1617 ## make sure we had all matrices: 102 : if(!(is.matrix(dims) && nrow(dims) == 2)) 103 : stop("some arguments are not matrices") 104 : maechler 1845 csdim <- rbind(rep.int(0L, 2), 105 : mmaechler 2206 apply(dims, 1, cumsum)) 106 : r <- new("dgTMatrix") 107 : r@Dim <- as.integer(csdim[nrow(csdim),]) 108 : maechler 1617 add1 <- matrix(1:0, 2,2) 109 : mmaechler 2206 for(i in seq_along(lst)) { 110 : maechler 1617 indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2]) 111 : if(is.null(dim(indx))) ## non-square matrix 112 : mmaechler 2206 r[indx[[1]],indx[[2]]] <- lst[[i]] 113 : maechler 1617 else ## square matrix 114 : mmaechler 2206 r[indx[,1], indx[,2]] <- lst[[i]] 115 : maechler 1617 } 116 : mmaechler 2206 r 117 : maechler 1617 } 118 : mmaechler 2209 ## expand() needed something like bdiag() for lower-triangular 119 : ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient 120 : ## implementation for those; now extended and generalized: 121 : mmaechler 2206 .bdiag <- function(lst) { 122 : mmaechler 2209 ## block-diagonal matrix [a dgTMatrix] from list of matrices 123 : mmaechler 2206 stopifnot(is.list(lst), (nl <- length(lst)) >= 1) 124 : maechler 1617 125 : mmaechler 2206 Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N" 126 : mmaechler 2209 as, "TsparseMatrix") 127 : mmaechler 2206 if(nl == 1) return(Tlst[[1]]) 128 : ## else 129 : mmaechler 2813 i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L))) 130 : j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L))) 131 : mmaechler 2209 132 : mmaechler 2813 clss <- vapply(Tlst, class, "") 133 : mmaechler 2741 typ <- substr(clss, 2, 2) 134 : knd <- substr(clss, 1, 1) 135 : sym <- typ == "s" # symmetric ones 136 : tri <- typ == "t" # triangular ones 137 : use.n <- any(is.n <- knd == "n") 138 : if(use.n && !(use.n <- all(is.n))) { 139 : mmaechler 2209 Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix") 140 : mmaechler 2741 knd [is.n] <- "l" 141 : } 142 : use.l <- !use.n && all(knd == "l") 143 : mmaechler 2209 if(all(sym)) { ## result should be *symmetric* 144 : mmaechler 2813 uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L" 145 : mmaechler 2209 tLU <- table(uplos)# of length 1 or 2 .. 146 : if(length(tLU) == 1) { ## all "U" or all "L" 147 : useU <- uplos[1] == "U" 148 : } else { ## length(tLU) == 2, counting "L" and "U" 149 : useU <- diff(tLU) >= 0 150 : if(useU && (hasL <- tLU[1] > 0)) 151 : Tlst[hasL] <- lapply(Tlst[hasL], t) 152 : else if(!useU && (hasU <- tLU[2] > 0)) 153 : Tlst[hasU] <- lapply(Tlst[hasU], t) 154 : } 155 : if(use.n) { ## return nsparseMatrix : 156 : r <- new("nsTMatrix") 157 : } else { 158 : mmaechler 2741 r <- new(paste0(if(use.l) "l" else "d", "sTMatrix")) 159 : mmaechler 2209 r@x <- unlist(lapply(Tlst, slot, "x")) 160 : } 161 : r@uplo <- if(useU) "U" else "L" 162 : } 163 : mmaechler 2813 else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")## "U" or "L" 164 : mmaechler 2209 all(ULs[1L] == ULs[-1L]) } ## all upper or all lower 165 : ){ ## *triangular* result 166 : 167 : if(use.n) { ## return nsparseMatrix : 168 : r <- new("ntTMatrix") 169 : } else { 170 : mmaechler 2741 r <- new(paste0(if(use.l) "l" else "d", "tTMatrix")) 171 : mmaechler 2209 r@x <- unlist(lapply(Tlst, slot, "x")) 172 : } 173 : r@uplo <- ULs[1L] 174 : } 175 : else { 176 : if(any(sym)) 177 : Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix") 178 : if(use.n) { ## return nsparseMatrix : 179 : r <- new("ngTMatrix") 180 : } else { 181 : mmaechler 2741 r <- new(paste0(if(use.l) "l" else "d", "gTMatrix")) 182 : mmaechler 2209 r@x <- unlist(lapply(Tlst, slot, "x")) 183 : } 184 : } 185 : r@Dim <- c(i_off[nl+1], j_off[nl + 1]) 186 : r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k])) 187 : r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k])) 188 : r 189 : mmaechler 2206 } 190 : 191 : bdiag <- function(...) { 192 : if((nA <- nargs()) == 0) return(new("dgCMatrix")) 193 : if(nA == 1 && !is.list(...)) 194 : return(as(..., "CsparseMatrix")) 195 : alis <- if(nA == 1 && is.list(..1)) ..1 else list(...) 196 : if(length(alis) == 1) 197 : return(as(alis[[1]], "CsparseMatrix")) 198 : 199 : ## else : two or more arguments 200 : as(.bdiag(alis), "CsparseMatrix") 201 : } 202 : 203 : 204 : maechler 2136 .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) { 205 : ## to triangular Tsparse 206 : maechler 1845 i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L 207 : mmaechler 2778 new(paste0(kind, "tTMatrix"), 208 : maechler 1654 diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames, 209 : maechler 2128 uplo = uplo, 210 : maechler 1654 x = from@x, # <- ok for diag = "U" and "N" (!) 211 : i = i, j = i) 212 : } 213 : maechler 1109 214 : maechler 2136 .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) { 215 : ## to symmetric Tsparse 216 : maechler 2128 n <- from@Dim[1] 217 : i <- seq_len(n) - 1L 218 : mmaechler 2778 new(paste0(kind, "sTMatrix"), 219 : maechler 1845 Dim = from@Dim, Dimnames = from@Dimnames, 220 : maechler 2128 i = i, j = i, uplo = uplo, 221 : x = if(from@diag == "N") from@x else ## "U"-diag 222 : rep.int(switch(kind, 223 : "d" = 1., 224 : "l" =, 225 : "n" = TRUE, 226 : ## otherwise 227 : stop("'", kind,"' kind not yet implemented")), n)) 228 : maechler 1845 } 229 : 230 : maechler 2128 ## diagonal -> triangular, upper / lower depending on "partner": 231 : maechler 2144 diag2tT.u <- function(d, x, kind = .M.kind(d)) 232 : .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind) 233 : maechler 2128 234 : mmaechler 2175 ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner": 235 : diag2Tsmart <- function(d, x, kind = .M.kind(d)) { 236 : clx <- getClassDef(class(x)) 237 : if(extends(clx, "symmetricMatrix")) 238 : .diag2sT(d, uplo = x@uplo, kind) 239 : else 240 : .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind) 241 : } 242 : maechler 2136 243 : mmaechler 2814 ## FIXME: should not be needed {when ddi* is dsparse* etc}: 244 : setMethod("is.finite", signature(x = "diagonalMatrix"), 245 : function(x) is.finite(.diag2tT(x))) 246 : setMethod("is.infinite", signature(x = "diagonalMatrix"), 247 : function(x) is.infinite(.diag2tT(x))) 248 : mmaechler 2175 249 : maechler 2136 ## In order to evade method dispatch ambiguity warnings, 250 : ## and because we can save a .M.kind() call, we use this explicit 251 : ## "hack" instead of signature x = "diagonalMatrix" : 252 : ## 253 : ## ddi*: 254 : diag2tT <- function(from) .diag2tT(from, "U", "d") 255 : setAs("ddiMatrix", "triangularMatrix", diag2tT) 256 : mmaechler 2239 ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT) 257 : maechler 1805 ## needed too (otherwise -> Tsparse is taken): 258 : maechler 2136 setAs("ddiMatrix", "TsparseMatrix", diag2tT) 259 : mmaechler 2482 setAs("ddiMatrix", "dsparseMatrix", diag2tT) 260 : maechler 2136 setAs("ddiMatrix", "CsparseMatrix", 261 : function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix")) 262 : setAs("ddiMatrix", "symmetricMatrix", 263 : function(from) .diag2sT(from, "U", "d")) 264 : ## 265 : ## ldi*: 266 : diag2tT <- function(from) .diag2tT(from, "U", "l") 267 : setAs("ldiMatrix", "triangularMatrix", diag2tT) 268 : mmaechler 2239 ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT) 269 : maechler 2136 ## needed too (otherwise -> Tsparse is taken): 270 : setAs("ldiMatrix", "TsparseMatrix", diag2tT) 271 : mmaechler 2482 setAs("ldiMatrix", "lsparseMatrix", diag2tT) 272 : maechler 2136 setAs("ldiMatrix", "CsparseMatrix", 273 : function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix")) 274 : setAs("ldiMatrix", "symmetricMatrix", 275 : function(from) .diag2sT(from, "U", "l")) 276 : maechler 1654 277 : maechler 1845 278 : maechler 2128 setAs("diagonalMatrix", "nMatrix", 279 : function(from) { 280 : n <- from@Dim[1] 281 : i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L 282 : new("ntTMatrix", i = i, j = i, diag = from@diag, 283 : Dim = from@Dim, Dimnames = from@Dimnames) 284 : }) 285 : 286 : mmaechler 2237 setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix")) 287 : maechler 2128 288 : mmaechler 2175 ## Cheap fast substitute for diag() which *does* preserve the mode of x : 289 : mkDiag <- function(x, n) { 290 : y <- matrix(as0(mod=mode(x)), n,n) 291 : mmaechler 2820 if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x 292 : mmaechler 2175 y 293 : } 294 : 295 : maechler 1109 setAs("diagonalMatrix", "matrix", 296 : function(from) { 297 : mmaechler 2175 ## want "ldiMatrix" -> "matrix" : 298 : mkDiag(if(from@diag == "U") as1(from@x) else from@x, 299 : n = from@Dim[1]) 300 : maechler 1109 }) 301 : 302 : maechler 2098 setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"), 303 : function(x, mode) { 304 : n <- x@Dim[1] 305 : mmaechler 2175 mod.x <- mode(x@x) 306 : r <- vector(mod.x, length = n^2) 307 : maechler 2098 if(n) 308 : mmaechler 2820 r[1 + 0:(n - 1L) * (n + 1)] <- 309 : mmaechler 2175 if(x@diag == "U") as1(mod=mod.x) else x@x 310 : maechler 2098 r 311 : }) 312 : 313 : maechler 1654 setAs("diagonalMatrix", "generalMatrix", # prefer sparse: 314 : maechler 2128 function(from) as(as(from, "CsparseMatrix"), "generalMatrix")) 315 : maechler 1174 316 : mmaechler 2268 setAs("diagonalMatrix", "denseMatrix", 317 : function(from) as(as(from, "CsparseMatrix"), "denseMatrix")) 318 : 319 : mmaechler 2185 .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x 320 : maechler 1655 321 : .diag.2N <- function(m) { 322 : if(m@diag == "U") m@diag <- "N" 323 : m 324 : } 325 : 326 : maechler 1174 setAs("ddiMatrix", "dgeMatrix", 327 : mmaechler 2183 function(from) .Call(dup_mMatrix_as_dgeMatrix, from)) 328 : setAs("ddiMatrix", "ddenseMatrix", 329 : function(from) as(as(from, "triangularMatrix"),"denseMatrix")) 330 : setAs("ldiMatrix", "ldenseMatrix", 331 : function(from) as(as(from, "triangularMatrix"),"denseMatrix")) 332 : maechler 1174 333 : mmaechler 2183 334 : maechler 1109 setAs("matrix", "diagonalMatrix", 335 : function(from) { 336 : maechler 1295 d <- dim(from) 337 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix") 338 : if(any(from[row(from) != col(from)] != 0)) 339 : mmaechler 2341 stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix") 340 : maechler 1295 x <- diag(from) 341 : if(is.logical(x)) { 342 : cl <- "ldiMatrix" 343 : mmaechler 2341 uni <- allTrue(x) ## uni := {is it unit-diagonal ?} 344 : maechler 1295 } else { 345 : cl <- "ddiMatrix" 346 : mmaechler 2341 uni <- allTrue(x == 1) 347 : maechler 1295 storage.mode(x) <- "double" 348 : maechler 1575 } ## TODO: complex 349 : maechler 1295 new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 350 : x = if(uni) x[FALSE] else x) 351 : maechler 1109 }) 352 : 353 : ## ``generic'' coercion to diagonalMatrix : build on isDiagonal() and diag() 354 : setAs("Matrix", "diagonalMatrix", 355 : function(from) { 356 : mmaechler 2341 d <- dim(from) 357 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix") 358 : mmaechler 2341 if(!isDiagonal(from)) stop("matrix is not diagonal") 359 : ## else: 360 : x <- diag(from) 361 : if(is.logical(x)) { 362 : cl <- "ldiMatrix" 363 : uni <- allTrue(x) 364 : } else { 365 : cl <- "ddiMatrix" 366 : uni <- allTrue(x == 1) 367 : storage.mode(x) <- "double" 368 : } ## TODO: complex 369 : new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 370 : x = if(uni) x[FALSE] else x) 371 : maechler 1109 }) 372 : 373 : maechler 1708 374 : mmaechler 2185 setMethod("diag", signature(x = "diagonalMatrix"), 375 : function(x = 1, nrow, ncol) .diag.x(x)) 376 : maechler 1708 377 : maechler 2098 subDiag <- function(x, i, j, ..., drop) { 378 : mmaechler 2508 x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now) 379 : maechler 1799 x <- if(missing(i)) 380 : x[, j, drop=drop] 381 : else if(missing(j)) 382 : mmaechler 2363 if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop] 383 : maechler 1799 else 384 : x[i,j, drop=drop] 385 : maechler 2120 if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x 386 : maechler 1799 } 387 : 388 : setMethod("[", signature(x = "diagonalMatrix", i = "index", 389 : j = "index", drop = "logical"), subDiag) 390 : setMethod("[", signature(x = "diagonalMatrix", i = "index", 391 : mmaechler 2363 j = "missing", drop = "logical"), 392 : function(x, i, j, ..., drop) { 393 : na <- nargs() 394 : Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2) 395 : if(na == 4) 396 : subDiag(x, i=i, , drop=drop) 397 : else subDiag(x, i=i, drop=drop) 398 : }) 399 : maechler 1799 setMethod("[", signature(x = "diagonalMatrix", i = "missing", 400 : j = "index", drop = "logical"), 401 : maechler 2098 function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop)) 402 : maechler 1799 403 : maechler 1617 ## When you assign to a diagonalMatrix, the result should be 404 : maechler 1708 ## diagonal or sparse --- 405 : ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch 406 : maechler 2098 ## Only(?) current bug: x[i] <- value is wrong when i is *vector* 407 : replDiag <- function(x, i, j, ..., value) { 408 : mmaechler 2811 x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07 409 : maechler 1710 if(missing(i)) 410 : x[, j] <- value 411 : maechler 2098 else if(missing(j)) { ## x[i , ] <- v *OR* x[i] <- v 412 : na <- nargs() 413 : ## message("diagnosing replDiag() -- nargs()= ", na) 414 : if(na == 4) 415 : x[i, ] <- value 416 : else if(na == 3) 417 : x[i] <- value 418 : else stop("Internal bug: nargs()=",na,"; please report") 419 : } else 420 : maechler 1710 x[i,j] <- value 421 : if(isDiagonal(x)) as(x, "diagonalMatrix") else x 422 : } 423 : maechler 1617 424 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", 425 : j = "index", value = "replValue"), replDiag) 426 : maechler 2098 427 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", 428 : j = "missing", value = "replValue"), 429 : maechler 2098 function(x,i,j, ..., value) { 430 : ## message("before replDiag() -- nargs()= ", nargs()) 431 : if(nargs() == 3) 432 : replDiag(x, i=i, value=value) 433 : else ## nargs() == 4 : 434 : replDiag(x, i=i, , value=value) 435 : }) 436 : 437 : mmaechler 2192 setReplaceMethod("[", signature(x = "diagonalMatrix", 438 : i = "matrix", # 2-col.matrix 439 : maechler 2096 j = "missing", value = "replValue"), 440 : maechler 2098 function(x,i,j, ..., value) { 441 : maechler 2096 if(ncol(i) == 2) { 442 : if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only 443 : mmaechler 2192 if(x@diag == "U") { 444 : one <- as1(x@x) 445 : if(any(value != one | is.na(value))) { 446 : x@diag <- "N" 447 : x@x <- rep.int(one, x@Dim[1]) 448 : mmaechler 2811 } else return(x) 449 : mmaechler 2192 } 450 : maechler 2096 x@x[ii] <- value 451 : x 452 : } else { ## no longer diagonal, but remain sparse: 453 : mmaechler 2239 x <- as(x, "TsparseMatrix") 454 : maechler 2096 x[i] <- value 455 : x 456 : } 457 : } 458 : else { # behave as "base R": use as if vector 459 : x <- as(x, "matrix") 460 : x[i] <- value 461 : Matrix(x) 462 : } 463 : }) 464 : 465 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", 466 : j = "index", value = "replValue"), 467 : maechler 2098 function(x,i,j, ..., value) replDiag(x, j=j, value=value)) 468 : maechler 1710 469 : mmaechler 2207 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index", 470 : mmaechler 2239 value = "sparseMatrix"), 471 : mmaechler 2207 function (x, i, j, ..., value) 472 : callGeneric(x=x, , j=j, value = as(value, "sparseVector"))) 473 : setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing", 474 : mmaechler 2239 value = "sparseMatrix"), 475 : mmaechler 2207 function (x, i, j, ..., value) 476 : callGeneric(x=x, i=i, , value = as(value, "sparseVector"))) 477 : setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index", 478 : mmaechler 2239 value = "sparseMatrix"), 479 : mmaechler 2207 function (x, i, j, ..., value) 480 : callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector"))) 481 : maechler 1710 482 : mmaechler 2207 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index", 483 : value = "sparseVector"), 484 : replDiag) 485 : setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing", 486 : value = "sparseVector"), 487 : replDiag) 488 : setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index", 489 : value = "sparseVector"), 490 : replDiag) 491 : 492 : 493 : maechler 1109 setMethod("t", signature(x = "diagonalMatrix"), 494 : function(x) { x@Dimnames <- x@Dimnames[2:1] ; x }) 495 : 496 : maechler 1331 setMethod("isDiagonal", signature(object = "diagonalMatrix"), 497 : function(object) TRUE) 498 : setMethod("isTriangular", signature(object = "diagonalMatrix"), 499 : function(object) TRUE) 500 : maechler 1109 setMethod("isSymmetric", signature(object = "diagonalMatrix"), 501 : maechler 2112 function(object, ...) TRUE) 502 : maechler 1109 503 : maechler 2112 setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x) 504 : setMethod("skewpart", signature(x = "diagonalMatrix"), setZero) 505 : 506 : maechler 2133 setMethod("chol", signature(x = "ddiMatrix"), 507 : function(x, pivot, ...) { 508 : mmaechler 2183 if(x@diag == "U") return(x) 509 : ## else 510 : maechler 2133 if(any(x@x < 0)) 511 : stop("chol() is undefined for diagonal matrix with negative entries") 512 : maechler 1654 x@x <- sqrt(x@x) 513 : x 514 : }) 515 : ## chol(L) is L for logical diagonal: 516 : maechler 2133 setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x) 517 : maechler 1654 518 : mmaechler 2175 setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"), 519 : mmaechler 2185 function(x, logarithm, ...) mkDet(.diag.x(x), logarithm)) 520 : mmaechler 2175 521 : mmaechler 2183 setMethod("norm", signature(x = "diagonalMatrix", type = "character"), 522 : function(x, type, ...) { 523 : if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix" 524 : type <- toupper(substr(type[1], 1, 1)) 525 : isU <- (x@diag == "U") # unit-diagonal 526 : if(type == "F") sqrt(if(isU) n else sum(x@x^2)) 527 : else { ## norm == "I","1","O","M" : 528 : if(isU) 1 else max(abs(x@x)) 529 : } 530 : }) 531 : 532 : 533 : 534 : maechler 1109 ## Basic Matrix Multiplication {many more to add} 535 : maechler 1654 ## --------------------- 536 : ## Note that "ldi" logical are treated as numeric 537 : maechler 1109 diagdiagprod <- function(x, y) { 538 : mmaechler 2633 dimCheck(x,y) 539 : maechler 1109 if(x@diag != "U") { 540 : maechler 1654 if(y@diag != "U") { 541 : nx <- x@x * y@x 542 : if(is.numeric(nx) && !is.numeric(x@x)) 543 : x <- as(x, "dMatrix") 544 : x@x <- as.numeric(nx) 545 : } 546 : return(x) 547 : maechler 1109 } else ## x is unit diagonal 548 : return(y) 549 : } 550 : 551 : maechler 1654 setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 552 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 553 : 554 : maechler 1654 formals(diagdiagprod) <- alist(x=, y=x) 555 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 556 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 557 : maechler 1654 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 558 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 559 : maechler 1654 setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"), 560 : diagdiagprod, valueClass = "ddiMatrix") 561 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), 562 : diagdiagprod, valueClass = "ddiMatrix") 563 : maechler 1109 564 : 565 : diagmatprod <- function(x, y) { 566 : mmaechler 2183 ## x is diagonalMatrix 567 : maechler 1109 dx <- dim(x) 568 : dy <- dim(y) 569 : if(dx[2] != dy[1]) stop("non-matching dimensions") 570 : as(if(x@diag == "U") y else x@x * y, "Matrix") 571 : } 572 : setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), 573 : maechler 1654 diagmatprod) 574 : mmaechler 2183 ## sneaky .. : 575 : maechler 1109 formals(diagmatprod) <- alist(x=, y=NULL) 576 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"), 577 : maechler 1654 diagmatprod) 578 : maechler 1109 579 : mmaechler 2183 diagGeprod <- function(x, y) { 580 : maechler 1109 dx <- dim(x) 581 : dy <- dim(y) 582 : if(dx[2] != dy[1]) stop("non-matching dimensions") 583 : if(x@diag != "U") 584 : y@x <- x@x * y@x 585 : y 586 : } 587 : mmaechler 2183 setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod) 588 : setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod) 589 : formals(diagGeprod) <- alist(x=, y=NULL) 590 : maechler 1109 setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), 591 : mmaechler 2183 diagGeprod, valueClass = "dgeMatrix") 592 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"), 593 : diagGeprod) 594 : maechler 1109 595 : mmaechler 2183 matdiagprod <- function(x, y) { 596 : dx <- dim(x) 597 : dy <- dim(y) 598 : if(dx[2] != dy[1]) stop("non-matching dimensions") 599 : Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1])) 600 : } 601 : maechler 1109 setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), 602 : mmaechler 2183 matdiagprod) 603 : formals(matdiagprod) <- alist(x=, y=NULL) 604 : setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"), 605 : matdiagprod) 606 : maechler 1109 607 : mmaechler 2183 gediagprod <- function(x, y) { 608 : dx <- dim(x) 609 : dy <- dim(y) 610 : if(dx[2] != dy[1]) stop("non-matching dimensions") 611 : if(y@diag == "N") 612 : x@x <- x@x * rep(y@x, each = dx[1]) 613 : x 614 : } 615 : setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod) 616 : setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod) 617 : formals(gediagprod) <- alist(x=, y=NULL) 618 : setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"), 619 : gediagprod) 620 : setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"), 621 : gediagprod) 622 : maechler 1109 623 : maechler 1295 ## crossprod {more of these} 624 : maechler 1109 625 : maechler 1295 ## tcrossprod --- all are not yet there: do the dense ones here: 626 : maechler 1109 627 : mmaechler 2183 setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"), 628 : function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix")) 629 : setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"), 630 : function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y) 631 : 632 : 633 : maechler 1295 ## FIXME: 634 : ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"), 635 : ## function(x, y = NULL) { 636 : ## }) 637 : maechler 1109 638 : dmbates 2417 Cspdiagprod <- function(x, y) { 639 : dmbates 2458 dx <- dim(x <- .Call(Csparse_diagU2N, x)) 640 : dmbates 2417 dy <- dim(y) 641 : if(dx[2] != dy[1]) stop("non-matching dimensions") 642 : mmaechler 2674 if(y@diag == "N") { 643 : if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix")) 644 : x <- as(x, "generalMatrix") 645 : ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L]) 646 : x@x <- x@x * y@x[ind] 647 : } 648 : mmaechler 2608 if(is(x, "compMatrix") && length(xf <- x@factors)) { 649 : ## instead of dropping all factors, be smart about some 650 : ## TODO ...... 651 : x@factors <- list() 652 : } 653 : dmbates 2417 x 654 : } 655 : maechler 1109 656 : dmbates 2417 diagCspprod <- function(x, y) { 657 : dx <- dim(x) 658 : dmbates 2458 dy <- dim(y <- .Call(Csparse_diagU2N, y)) 659 : dmbates 2417 if(dx[2] != dy[1]) stop("non-matching dimensions") 660 : mmaechler 2674 if(x@diag == "N") { 661 : if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix")) 662 : y <- as(y, "generalMatrix") 663 : y@x <- y@x * x@x[y@i + 1L] 664 : } 665 : mmaechler 2608 if(is(y, "compMatrix") && length(yf <- y@factors)) { 666 : ## instead of dropping all factors, be smart about some 667 : ## TODO 668 : keep <- character() 669 : if(iLU <- names(yf) == "LU") { 670 : ## TODO keep <- "LU" 671 : } 672 : y@factors <- yf[keep] 673 : } 674 : dmbates 2417 y 675 : } 676 : 677 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"), 678 : function(x, y = NULL) diagCspprod(x, y)) 679 : 680 : maechler 1799 setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 681 : dmbates 2417 function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix"))) 682 : maechler 1295 683 : dmbates 2417 ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway 684 : ## x'y == (y'x)' 685 : setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"), 686 : function(x, y = NULL) t(diagCspprod(y, x))) 687 : 688 : maechler 1799 setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 689 : dmbates 2417 function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix")))) 690 : maechler 1799 691 : dmbates 2417 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"), 692 : function(x, y = NULL) diagCspprod(x, t(y))) 693 : 694 : maechler 1799 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 695 : dmbates 2417 function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix")))) 696 : maechler 1799 697 : dmbates 2417 setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"), 698 : function(x, y = NULL) Cspdiagprod(x, y)) 699 : 700 : maechler 1799 setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 701 : dmbates 2417 function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y)) 702 : mmaechler 2472 703 : dmbates 2417 setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"), 704 : function(x, y) diagCspprod(x, y)) 705 : maechler 1799 706 : setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"), 707 : dmbates 2417 function(x, y) diagCspprod(as(x, "CsparseMatrix"), y)) 708 : mmaechler 2472 709 : mmaechler 2183 setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"), 710 : dmbates 2417 function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y)) 711 : 712 : maechler 1799 setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"), 713 : dmbates 2417 function(x, y) Cspdiagprod(x, y)) 714 : 715 : maechler 1799 ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal* 716 : ## do indeed work by going through sparse (and *not* ddense)! 717 : 718 : setMethod("solve", signature(a = "diagonalMatrix", b = "missing"), 719 : function(a, b, ...) { 720 : a@x <- 1/ a@x 721 : a@Dimnames <- a@Dimnames[2:1] 722 : a 723 : }) 724 : 725 : solveDiag <- function(a, b, ...) { 726 : mmaechler 2633 if(a@Dim[1] != nrow(b)) 727 : maechler 1799 stop("incompatible matrix dimensions") 728 : ## trivially invert a 'in place' and multiply: 729 : a@x <- 1/ a@x 730 : a@Dimnames <- a@Dimnames[2:1] 731 : a %*% b 732 : } 733 : setMethod("solve", signature(a = "diagonalMatrix", b = "matrix"), 734 : solveDiag) 735 : setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"), 736 : solveDiag) 737 : 738 : maechler 2106 ## Schur() ---> ./eigen.R 739 : maechler 1799 740 : 741 : 742 : mmaechler 2175 ###---------------- (, , ) ---------------------- 743 : maechler 1295 744 : maechler 1655 ## Use function for several signatures, in order to evade 745 : mmaechler 2175 diagOdiag <- function(e1,e2) { 746 : ## result should also be diagonal _ if possible _ 747 : maechler 1655 r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible" 748 : mmaechler 2175 ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...: 749 : r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE, 750 : if(is.numeric(e2@x)) 0 else FALSE) 751 : if(is0(r00)) { ## r00 == 0 or FALSE --- result *is* diagonal 752 : mmaechler 2814 if(is.numeric(r)) { # "double" *or* "integer" 753 : mmaechler 2175 if(is.numeric(e2@x)) { 754 : e2@x <- r; return(.diag.2N(e2)) } 755 : if(!is.numeric(e1@x)) 756 : ## e.g. e1, e2 are logical; 757 : e1 <- as(e1, "dMatrix") 758 : mmaechler 2814 if(!is.double(r)) r <- as.double(r) 759 : mmaechler 2175 } 760 : else if(is.logical(r)) 761 : e1 <- as(e1, "lMatrix") 762 : else stop("intermediate 'r' is of type", typeof(r)) 763 : e1@x <- r 764 : .diag.2N(e1) 765 : maechler 1655 } 766 : mmaechler 2175 else { ## result not diagonal, but at least symmetric: 767 : mmaechler 2496 ## e.g., m == m 768 : mmaechler 2175 isNum <- (is.numeric(r) || is.numeric(r00)) 769 : isLog <- (is.logical(r) || is.logical(r00)) 770 : mmaechler 2496 Matrix.msg("exploding o into dense matrix", .M.level = 2) 771 : mmaechler 2175 d <- e1@Dim 772 : n <- d[1] 773 : stopifnot(length(r) == n) 774 : mmaechler 2814 if(isNum && !is.double(r)) r <- as.double(r) 775 : mmaechler 2175 xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n)) 776 : newcl <- 777 : mmaechler 2811 paste0(if(isNum) "d" else if(isLog) { 778 : mmaechler 2175 if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l" 779 : mmaechler 2811 } else stop("not yet implemented .. please report"), "syMatrix") 780 : mmaechler 2175 781 : new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx) 782 : } 783 : maechler 1655 } 784 : 785 : mmaechler 2175 ### This would be *the* way, but we get tons of "ambiguous method dispatch" 786 : mmaechler 2185 ## we use this hack instead of signature x = "diagonalMatrix" : 787 : diCls <- names(getClass("diagonalMatrix")@subclasses) 788 : mmaechler 2175 if(FALSE) { 789 : maechler 1655 setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"), 790 : diagOdiag) 791 : mmaechler 2175 } else { ## These are just for method disambiguation: 792 : for(c1 in diCls) 793 : for(c2 in diCls) 794 : setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag) 795 : } 796 : maechler 1655 797 : mmaechler 2493 ## diagonal o triangular |--> triangular 798 : ## diagonal o symmetric |--> symmetric 799 : maechler 1845 ## {also when other is sparse: do these "here" -- 800 : ## before conversion to sparse, since that loses "diagonality"} 801 : mmaechler 2811 diagOtri <- function(e1,e2) { 802 : ## result must be triangular 803 : r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible" 804 : ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...: 805 : e1.0 <- if(.n1 <- is.numeric(d1 )) 0 else FALSE 806 : r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE) 807 : if(is0(r00)) { ## r00 == 0 or FALSE --- result *is* triangular 808 : diag(e2) <- r 809 : ## check what happens "in the triangle" 810 : e2.2 <- if(.n2) 2 else TRUE 811 : if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change: 812 : n <- dim(e2)[1L] 813 : it <- indTri(n, upper = (e2@uplo == "U")) 814 : e2[it] <- callGeneric(e1.0, e2[it]) 815 : } 816 : e2 817 : } 818 : else { ## result not triangular ---> general 819 : rr <- as(e2, "generalMatrix") 820 : diag(rr) <- r 821 : rr 822 : } 823 : } 824 : maechler 1845 825 : mmaechler 2811 826 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"), 827 : diagOtri) 828 : ## For the reverse, Ops == "Arith" | "Compare" | "Logic" 829 : ## 'Arith' := '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"' 830 : setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"), 831 : function(e1,e2) 832 : { ## this must only trigger for *dense* e1 833 : switch(.Generic, 834 : "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), .diag.x(e2)), 835 : "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)), 836 : "*" = { 837 : n <- e2@Dim[1L] 838 : d2 <- if(e2@diag == "U") { # unit-diagonal 839 : d <- rep.int(as1(e2@x), n) 840 : e2@x <- d 841 : e2@diag <- "N" 842 : d 843 : } else e2@x 844 : e2@x <- diag(e1) * d2 845 : e2 846 : }, 847 : "^" = { ## will be dense ( as ^ 0 == 1 ): 848 : e1 ^ as(e2, "denseMatrix") 849 : }, 850 : ## otherwise: 851 : callGeneric(e1, diag2Tsmart(e2,e1))) 852 : }) 853 : 854 : ## Compare --> 'swap' (e.g. e1 < e2 <==> e2 > e1 ): 855 : setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"), 856 : .Cmp.swap) 857 : ## '&' and "|' are commutative: 858 : setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"), 859 : function(e1,e2) callGeneric(e2, e1)) 860 : 861 : maechler 1655 ## For almost everything else, diag* shall be treated "as sparse" : 862 : maechler 1295 ## These are cheap implementations via coercion 863 : 864 : maechler 2144 ## For disambiguation --- define this for "sparseMatrix" , then for "ANY"; 865 : ## and because we can save an .M.kind() call, we use this explicit 866 : ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" : 867 : ## 868 : ## ddi*: 869 : setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"), 870 : mmaechler 2175 function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)) 871 : maechler 2144 setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"), 872 : mmaechler 2175 function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d"))) 873 : maechler 2144 ## ldi* 874 : setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"), 875 : mmaechler 2175 function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)) 876 : maechler 2144 setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"), 877 : mmaechler 2175 function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l"))) 878 : maechler 1654 879 : mmaechler 2820 ## Ops: Arith --> numeric : "dMatrix" 880 : ## Compare --> logical 881 : ## Logic --> logical: "lMatrix" 882 : mmaechler 2175 883 : mmaechler 2820 ## Other = "numeric" : stay diagonal if possible 884 : mmaechler 2175 ## ddi*: Arith: result numeric, potentially ddiMatrix 885 : mmaechler 2813 for(arg2 in c("numeric","logical")) 886 : setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2), 887 : maechler 2144 function(e1,e2) { 888 : mmaechler 2820 n <- e1@Dim[1] 889 : maechler 2144 f0 <- callGeneric(0, e2) 890 : if(all(is0(f0))) { # remain diagonal 891 : maechler 2157 L1 <- (le <- length(e2)) == 1L 892 : mmaechler 2820 if(e1@diag == "U") { 893 : if(any((r <- callGeneric(1, e2)) != 1)) { 894 : e1@diag <- "N" 895 : e1@x[seq_len(n)] <- r # possibly recycling r 896 : } ## else: result = e1 (is "U" diag) 897 : } else { 898 : maechler 2144 r <- callGeneric(e1@x, e2) 899 : mmaechler 2820 ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix 900 : e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le] 901 : } 902 : e1 903 : } else 904 : callGeneric(diag2tT.u(e1,e2, "d"), e2) 905 : maechler 2144 }) 906 : maechler 1655 907 : mmaechler 2813 for(arg1 in c("numeric","logical")) 908 : setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"), 909 : maechler 2144 function(e1,e2) { 910 : mmaechler 2820 n <- e2@Dim[1] 911 : maechler 2144 f0 <- callGeneric(e1, 0) 912 : if(all(is0(f0))) { # remain diagonal 913 : maechler 2157 L1 <- (le <- length(e1)) == 1L 914 : mmaechler 2820 if(e2@diag == "U") { 915 : if(any((r <- callGeneric(e1, 1)) != 1)) { 916 : e2@diag <- "N" 917 : e2@x[seq_len(n)] <- r # possibly recycling r 918 : } ## else: result = e2 (is "U" diag) 919 : } else { 920 : maechler 2144 r <- callGeneric(e1, e2@x) 921 : mmaechler 2820 ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix 922 : e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le] 923 : } 924 : e2 925 : } else 926 : callGeneric(e1, diag2tT.u(e2,e1, "d")) 927 : maechler 2144 }) 928 : mmaechler 2175 929 : ## ldi* Arith --> result numeric, potentially ddiMatrix 930 : mmaechler 2813 for(arg2 in c("numeric","logical")) 931 : setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2), 932 : mmaechler 2175 function(e1,e2) { 933 : mmaechler 2820 n <- e1@Dim[1] 934 : mmaechler 2175 f0 <- callGeneric(0, e2) 935 : if(all(is0(f0))) { # remain diagonal 936 : L1 <- (le <- length(e2)) == 1L 937 : mmaechler 2820 E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE 938 : ## E <- copyClass(e1, "ddiMatrix", check=FALSE) 939 : ## storage.mode(E@x) <- "double" 940 : if(e1@diag == "U") { 941 : if(any((r <- callGeneric(1, e2)) != 1)) { 942 : E@diag <- "N" 943 : E@x[seq_len(n)] <- r # possibly recycling r 944 : } ## else: result = E (is "U" diag) 945 : } else { 946 : mmaechler 2175 r <- callGeneric(e1@x, e2) 947 : mmaechler 2820 ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix 948 : E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le] 949 : mmaechler 2814 } 950 : mmaechler 2820 E 951 : } else 952 : callGeneric(diag2tT.u(e1,e2, "l"), e2) 953 : mmaechler 2175 }) 954 : 955 : mmaechler 2813 for(arg1 in c("numeric","logical")) 956 : setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"), 957 : mmaechler 2175 function(e1,e2) { 958 : mmaechler 2820 n <- e2@Dim[1] 959 : mmaechler 2175 f0 <- callGeneric(e1, 0) 960 : if(all(is0(f0))) { # remain diagonal 961 : L1 <- (le <- length(e1)) == 1L 962 : mmaechler 2820 E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE 963 : ## E <- copyClass(e2, "ddiMatrix", check=FALSE) 964 : ## storage.mode(E@x) <- "double" 965 : if(e2@diag == "U") { 966 : if(any((r <- callGeneric(e1, 1)) != 1)) { 967 : E@diag <- "N" 968 : E@x[seq_len(n)] <- r # possibly recycling r 969 : } ## else: result = E (is "U" diag) 970 : } else { 971 : mmaechler 2175 r <- callGeneric(e1, e2@x) 972 : mmaechler 2820 ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix 973 : E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le] 974 : mmaechler 2814 } 975 : mmaechler 2820 E 976 : } else 977 : callGeneric(e1, diag2tT.u(e2,e1, "l")) 978 : mmaechler 2175 }) 979 : 980 : mmaechler 2820 ## ddi*: for "Ops" without "Arith": or --> result logical, potentially ldi 981 : ## 982 : ## Note that ("numeric", "ddiMatrix") is simply swapped, e.g., 983 : if(FALSE) { 984 : selectMethod("<", c("numeric","lMatrix"))# Compare 985 : selectMethod("&", c("numeric","lMatrix"))# Logic 986 : } ## so we don't need to define a method here : 987 : 988 : mmaechler 2813 for(arg2 in c("numeric","logical")) 989 : setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2), 990 : mmaechler 2175 function(e1,e2) { 991 : mmaechler 2820 n <- e1@Dim[1] 992 : mmaechler 2175 f0 <- callGeneric(0, e2) 993 : if(all(is0(f0))) { # remain diagonal 994 : L1 <- (le <- length(e2)) == 1L 995 : mmaechler 2820 E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE 996 : ## E <- copyClass(e1, "ldiMatrix", check=FALSE) 997 : ## storage.mode(E@x) <- "logical" 998 : if(e1@diag == "U") { 999 : if(any((r <- callGeneric(1, e2)) != 1)) { 1000 : E@diag <- "N" 1001 : E@x[seq_len(n)] <- r # possibly recycling r 1002 : } ## else: result = E (is "U" diag) 1003 : } else { 1004 : mmaechler 2175 r <- callGeneric(e1@x, e2) 1005 : mmaechler 2820 ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix 1006 : E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le] 1007 : mmaechler 2814 } 1008 : mmaechler 2820 E 1009 : } else 1010 : callGeneric(diag2tT.u(e1,e2, "d"), e2) 1011 : mmaechler 2175 }) 1012 : 1013 : mmaechler 2820 ## ldi*: for "Ops" without "Arith": or --> result logical, potentially ldi 1014 : mmaechler 2813 for(arg2 in c("numeric","logical")) 1015 : setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2), 1016 : maechler 2144 function(e1,e2) { 1017 : mmaechler 2820 n <- e1@Dim[1] 1018 : maechler 2144 f0 <- callGeneric(FALSE, e2) 1019 : if(all(is0(f0))) { # remain diagonal 1020 : maechler 2157 L1 <- (le <- length(e2)) == 1L 1021 : mmaechler 2820 1022 : if(e1@diag == "U") { 1023 : if(any((r <- callGeneric(TRUE, e2)) != 1)) { 1024 : e1@diag <- "N" 1025 : e1@x[seq_len(n)] <- r # possibly recycling r 1026 : } ## else: result = e1 (is "U" diag) 1027 : } else { 1028 : maechler 2144 r <- callGeneric(e1@x, e2) 1029 : mmaechler 2820 ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix 1030 : e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le] 1031 : } 1032 : e1 1033 : } else 1034 : callGeneric(diag2tT.u(e1,e2, "l"), e2) 1035 : maechler 2144 }) 1036 : maechler 1655 1037 : maechler 2144 1038 : ## Not {"sparseMatrix", "numeric} : {"denseMatrix", "matrix", ... } 1039 : mmaechler 2196 for(other in c("ANY", "Matrix", "dMatrix")) { 1040 : ## ddi*: 1041 : setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other), 1042 : mmaechler 2493 function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)) 1043 : mmaechler 2196 setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"), 1044 : mmaechler 2493 function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d"))) 1045 : mmaechler 2196 ## ldi*: 1046 : setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other), 1047 : mmaechler 2493 function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)) 1048 : mmaechler 2196 setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"), 1049 : mmaechler 2493 function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l"))) 1050 : mmaechler 2196 } 1051 : mmaechler 2331 1052 : ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... : 1053 : dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses 1054 : mmaechler 2813 names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] }) 1055 : mmaechler 2268 for(DI in diCls) { 1056 : mmaechler 2814 dMeth <- if(extends(DI, "dMatrix")) 1057 : function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2) 1058 : else # "lMatrix", the only other kind for now 1059 : function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2) 1060 : mmaechler 2331 for(c2 in c(dense.subCl, "Matrix")) { 1061 : mmaechler 2814 for(Fun in c("*", "&")) { 1062 : mmaechler 2268 setMethod(Fun, signature(e1 = DI, e2 = c2), 1063 : function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2)))) 1064 : setMethod(Fun, signature(e1 = c2, e2 = DI), 1065 : function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2)) 1066 : } 1067 : mmaechler 2814 setMethod("^", signature(e1 = c2, e2 = DI), 1068 : function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2)) 1069 : for(Fun in c("%%", "%/%", "/")) ## 0 0 |--> NaN for these. 1070 : setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth) 1071 : mmaechler 2175 } 1072 : } 1073 : maechler 2144 1074 : 1075 : mmaechler 2175 ### "Summary" : "max" "min" "range" "prod" "sum" "any" "all" 1076 : mmaechler 2345 ### ---------- the last 4: separately here 1077 : mmaechler 2175 for(cl in diCls) { 1078 : setMethod("any", cl, 1079 : function (x, ..., na.rm) { 1080 : if(any(x@Dim == 0)) FALSE 1081 : else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm) 1082 : }) 1083 : mmaechler 2341 setMethod("all", cl, function (x, ..., na.rm) { 1084 : n <- x@Dim[1] 1085 : if(n >= 2) FALSE 1086 : else if(n == 0 || x@diag == "U") TRUE 1087 : else all(x@x, ..., na.rm = na.rm) 1088 : }) 1089 : setMethod("prod", cl, function (x, ..., na.rm) { 1090 : n <- x@Dim[1] 1091 : if(n >= 2) 0 1092 : else if(n == 0 || x@diag == "U") 1 1093 : else ## n == 1, diag = "N" : 1094 : prod(x@x, ..., na.rm = na.rm) 1095 : }) 1096 : mmaechler 2175 1097 : setMethod("sum", cl, 1098 : function(x, ..., na.rm) { 1099 : r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly 1100 : if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r 1101 : }) 1102 : } 1103 : 1104 : ## The remaining ones are max, min, range : 1105 : 1106 : setMethod("Summary", "ddiMatrix", 1107 : function(x, ..., na.rm) { 1108 : if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm) 1109 : else if(x@diag == "U") 1110 : callGeneric(x@x, 0, 1, ..., na.rm=na.rm) 1111 : else callGeneric(x@x, 0, ..., na.rm=na.rm) 1112 : }) 1113 : setMethod("Summary", "ldiMatrix", 1114 : function(x, ..., na.rm) { 1115 : if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm) 1116 : else if(x@diag == "U") 1117 : callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm) 1118 : else callGeneric(x@x, FALSE, ..., na.rm=na.rm) 1119 : }) 1120 : 1121 : 1122 : 1123 : maechler 1109 ## similar to prTriang() in ./Auxiliaries.R : 1124 : prDiag <- 1125 : function(x, digits = getOption("digits"), justify = "none", right = TRUE) 1126 : { 1127 : cf <- array(".", dim = x@Dim, dimnames = x@Dimnames) 1128 : cf[row(cf) == col(cf)] <- 1129 : sapply(diag(x), format, digits = digits, justify = justify) 1130 : print(cf, quote = FALSE, right = right) 1131 : invisible(x) 1132 : } 1133 : 1134 : mmaechler 2237 ## somewhat consistent with "print" for sparseMatrix : 1135 : setMethod("print", signature(x = "diagonalMatrix"), prDiag) 1136 : 1137 : maechler 1109 setMethod("show", signature(object = "diagonalMatrix"), 1138 : maechler 1592 function(object) { 1139 : d <- dim(object) 1140 : cl <- class(object) 1141 : mmaechler 2203 cat(sprintf('%d x %d diagonal matrix of class "%s"', 1142 : maechler 1592 d[1], d[2], cl)) 1143 : mmaechler 2203 if(d[1] < 50) { 1144 : cat("\n") 1145 : prDiag(object) 1146 : } else { 1147 : cat(", with diagonal entries\n") 1148 : show(diag(object)) 1149 : invisible(object) 1150 : } 1151 : maechler 1592 }) 1152 : mmaechler 2813 1153 : rm(dense.subCl, diCls)# not used elsewhere 1154 : mmaechler 2820 1155 : setMethod("summary", signature(object = "diagonalMatrix"), 1156 : function(object, ...) { 1157 : d <- dim(object) 1158 : r <- summary(object@x, ...) 1159 : attr(r, "header") <- 1160 : sprintf('%d x %d diagonal Matrix of class "%s"', 1161 : d[1], d[2], class(object)) 1162 : ## use ole' S3 technology for such a simple case 1163 : class(r) <- c("diagSummary", class(r)) 1164 : r 1165 : }) 1166 : 1167 : print.diagSummary <- function (x, ...) { 1168 : cat(attr(x, "header"),"\n") 1169 : class(x) <- class(x)[-1] 1170 : print(x, ...) 1171 : invisible(x) 1172 : }

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