# SCM Repository

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

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

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

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