# SCM Repository

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

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

Revision 2508 - (view) (download)
Original Path: pkg/R/diagMatrix.R

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

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