# SCM Repository

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

# Annotation of /pkg/R/diagMatrix.R

Revision 1832 - (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 : 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 1575 stopifnot(length(x) == n) 20 : if(is.logical(x)) 21 : cl <- "ldiMatrix" 22 : else if(is.numeric(x)) { 23 : cl <- "ddiMatrix" 24 : x <- as.numeric(x) 25 : } 26 : else if(is.complex(x)) { 27 : cl <- "zdiMatrix" # will not yet work 28 : } else stop("'x' has invalid data type") 29 : new(cl, Dim = c(n,n), diag = "N", x = x) 30 : maechler 1109 } 31 : } 32 : 33 : maechler 1617 ### This is modified from a post of Bert Gunter to R-help on 1 Sep 2005. 34 : ### Bert's code built on a post by Andy Liaw who most probably was influenced 35 : ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002 36 : ### who posted his bdiag() function written in December 1995. 37 : 38 : bdiag <- function(...) { 39 : if(nargs() == 0) return(new("dgCMatrix")) 40 : ## else : 41 : mlist <- if (nargs() == 1) as.list(...) else list(...) 42 : dims <- sapply(mlist, dim) 43 : ## make sure we had all matrices: 44 : if(!(is.matrix(dims) && nrow(dims) == 2)) 45 : stop("some arguments are not matrices") 46 : csdim <- rbind(rep.int(0:0, 2), 47 : apply(sapply(mlist, dim), 1, cumsum)) 48 : ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),])) 49 : add1 <- matrix(1:0, 2,2) 50 : maechler 1654 for(i in seq_along(mlist)) { 51 : maechler 1617 indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2]) 52 : if(is.null(dim(indx))) ## non-square matrix 53 : ret[indx[[1]],indx[[2]]] <- mlist[[i]] 54 : else ## square matrix 55 : ret[indx[,1],indx[,2]] <- mlist[[i]] 56 : } 57 : ## slightly debatable if we really should return Csparse.. : 58 : as(ret, "CsparseMatrix") 59 : } 60 : 61 : maechler 1654 diag2T <- function(from) { 62 : i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1:1 63 : new(paste(.M.kind(from), "tTMatrix", sep=''), 64 : diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames, 65 : x = from@x, # <- ok for diag = "U" and "N" (!) 66 : i = i, j = i) 67 : } 68 : maechler 1109 69 : maechler 1654 setAs("diagonalMatrix", "triangularMatrix", diag2T) 70 : setAs("diagonalMatrix", "sparseMatrix", diag2T) 71 : maechler 1805 ## needed too (otherwise -> Tsparse is taken): 72 : setAs("diagonalMatrix", "TsparseMatrix", diag2T) 73 : maechler 1654 ## is better than this: 74 : ## setAs("diagonalMatrix", "sparseMatrix", 75 : ## function(from) 76 : ## as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix")) 77 : setAs("diagonalMatrix", "CsparseMatrix", 78 : function(from) as(diag2T(from), "CsparseMatrix")) 79 : 80 : maechler 1109 setAs("diagonalMatrix", "matrix", 81 : function(from) { 82 : n <- from@Dim[1] 83 : diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE 84 : } else from@x, 85 : nrow = n, ncol = n) 86 : }) 87 : 88 : maechler 1654 setAs("diagonalMatrix", "generalMatrix", # prefer sparse: 89 : function(from) as(from, paste(.M.kind(from), "gCMatrix", sep=''))) 90 : maechler 1174 91 : maechler 1655 .diag.x <- function(m) { 92 : if(m@diag == "U") 93 : rep.int(if(is.numeric(m@x)) 1. else TRUE, 94 : m@Dim[1]) 95 : else m@x 96 : } 97 : 98 : .diag.2N <- function(m) { 99 : if(m@diag == "U") m@diag <- "N" 100 : m 101 : } 102 : 103 : maechler 1654 ## given the above, the following 4 coercions should be all unneeded; 104 : ## we prefer triangular to general: 105 : maechler 1295 setAs("ddiMatrix", "dgTMatrix", 106 : function(from) { 107 : maechler 1654 .Deprecated("as(, \"sparseMatrix\")") 108 : maechler 1295 n <- from@Dim[1] 109 : maechler 1654 i <- seq_len(n) - 1:1 110 : maechler 1655 new("dgTMatrix", i = i, j = i, x = .diag.x(from), 111 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) }) 112 : 113 : setAs("ddiMatrix", "dgCMatrix", 114 : maechler 1655 function(from) as(as(from, "sparseMatrix"), "dgCMatrix")) 115 : maechler 1295 116 : setAs("ldiMatrix", "lgTMatrix", 117 : function(from) { 118 : maechler 1654 .Deprecated("as(, \"sparseMatrix\")") 119 : maechler 1295 n <- from@Dim[1] 120 : maechler 1575 if(from@diag == "U") { # unit-diagonal 121 : x <- rep.int(TRUE, n) 122 : maechler 1654 i <- seq_len(n) - 1:1 123 : maechler 1575 } else { # "normal" 124 : nz <- nz.NA(from@x, na. = TRUE) 125 : x <- from@x[nz] 126 : i <- which(nz) - 1:1 127 : } 128 : new("lgTMatrix", i = i, j = i, x = x, 129 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) }) 130 : 131 : setAs("ldiMatrix", "lgCMatrix", 132 : function(from) as(as(from, "lgTMatrix"), "lgCMatrix")) 133 : 134 : 135 : maechler 1447 if(FALSE) # now have faster "ddense" -> "dge" 136 : maechler 1174 setAs("ddiMatrix", "dgeMatrix", 137 : function(from) as(as(from, "matrix"), "dgeMatrix")) 138 : 139 : maechler 1109 setAs("matrix", "diagonalMatrix", 140 : function(from) { 141 : maechler 1295 d <- dim(from) 142 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix") 143 : if(any(from[row(from) != col(from)] != 0)) 144 : stop("has non-zero off-diagonal entries") 145 : maechler 1295 x <- diag(from) 146 : if(is.logical(x)) { 147 : cl <- "ldiMatrix" 148 : uni <- all(x) 149 : } else { 150 : cl <- "ddiMatrix" 151 : uni <- all(x == 1) 152 : storage.mode(x) <- "double" 153 : maechler 1575 } ## TODO: complex 154 : maechler 1295 new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 155 : x = if(uni) x[FALSE] else x) 156 : maechler 1109 }) 157 : 158 : ## ``generic'' coercion to diagonalMatrix : build on isDiagonal() and diag() 159 : setAs("Matrix", "diagonalMatrix", 160 : function(from) { 161 : d <- dim(from) 162 : if(d[1] != (n <- d[2])) stop("non-square matrix") 163 : if(!isDiagonal(from)) stop("matrix is not diagonal") 164 : ## else: 165 : x <- diag(from) 166 : if(is.logical(x)) { 167 : cl <- "ldiMatrix" 168 : uni <- all(x) 169 : } else { 170 : cl <- "ddiMatrix" 171 : uni <- all(x == 1) 172 : storage.mode(x) <- "double" 173 : } 174 : new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 175 : x = if(uni) x[FALSE] else x) 176 : }) 177 : 178 : maechler 1708 179 : setMethod("diag", signature(x = "diagonalMatrix"), 180 : function(x = 1, nrow, ncol = n) .diag.x(x)) 181 : 182 : maechler 1799 183 : subDiag <- function(x, i, j, drop) { 184 : x <- as(x, "sparseMatrix") 185 : x <- if(missing(i)) 186 : x[, j, drop=drop] 187 : else if(missing(j)) 188 : x[i, , drop=drop] 189 : else 190 : x[i,j, drop=drop] 191 : if(isDiagonal(x)) as(x, "diagonalMatrix") else x 192 : } 193 : 194 : setMethod("[", signature(x = "diagonalMatrix", i = "index", 195 : j = "index", drop = "logical"), subDiag) 196 : setMethod("[", signature(x = "diagonalMatrix", i = "index", 197 : j = "missing", drop = "logical"), 198 : function(x, i, drop) subDiag(x, i=i, drop=drop)) 199 : setMethod("[", signature(x = "diagonalMatrix", i = "missing", 200 : j = "index", drop = "logical"), 201 : function(x, j, drop) subDiag(x, j=j, drop=drop)) 202 : 203 : maechler 1617 ## When you assign to a diagonalMatrix, the result should be 204 : maechler 1708 ## diagonal or sparse --- 205 : ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch 206 : maechler 1710 replDiag <- function(x, i, j, value) { 207 : x <- as(x, "sparseMatrix") 208 : if(missing(i)) 209 : x[, j] <- value 210 : else if(missing(j)) 211 : x[i, ] <- value 212 : else 213 : x[i,j] <- value 214 : if(isDiagonal(x)) as(x, "diagonalMatrix") else x 215 : } 216 : maechler 1617 217 : maechler 1710 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", 218 : j = "index", value = "replValue"), replDiag) 219 : setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", 220 : j = "missing", value = "replValue"), 221 : function(x, i, value) replDiag(x, i=i, value=value)) 222 : setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", 223 : j = "index", value = "replValue"), 224 : function(x, j, value) replDiag(x, j=j, value=value)) 225 : 226 : 227 : maechler 1109 setMethod("t", signature(x = "diagonalMatrix"), 228 : function(x) { x@Dimnames <- x@Dimnames[2:1] ; x }) 229 : 230 : maechler 1331 setMethod("isDiagonal", signature(object = "diagonalMatrix"), 231 : function(object) TRUE) 232 : setMethod("isTriangular", signature(object = "diagonalMatrix"), 233 : function(object) TRUE) 234 : maechler 1109 setMethod("isSymmetric", signature(object = "diagonalMatrix"), 235 : function(object) TRUE) 236 : 237 : maechler 1654 setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY" 238 : function(x, pivot) { 239 : if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries") 240 : x@x <- sqrt(x@x) 241 : x 242 : }) 243 : ## chol(L) is L for logical diagonal: 244 : setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x) 245 : 246 : maechler 1109 ## Basic Matrix Multiplication {many more to add} 247 : maechler 1654 ## --------------------- 248 : ## Note that "ldi" logical are treated as numeric 249 : maechler 1109 diagdiagprod <- function(x, y) { 250 : if(any(dim(x) != dim(y))) stop("non-matching dimensions") 251 : if(x@diag != "U") { 252 : maechler 1654 if(y@diag != "U") { 253 : nx <- x@x * y@x 254 : if(is.numeric(nx) && !is.numeric(x@x)) 255 : x <- as(x, "dMatrix") 256 : x@x <- as.numeric(nx) 257 : } 258 : return(x) 259 : maechler 1109 } else ## x is unit diagonal 260 : return(y) 261 : } 262 : 263 : maechler 1654 setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 264 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 265 : 266 : maechler 1654 formals(diagdiagprod) <- alist(x=, y=x) 267 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 268 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 269 : maechler 1654 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 270 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 271 : maechler 1654 setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"), 272 : diagdiagprod, valueClass = "ddiMatrix") 273 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), 274 : diagdiagprod, valueClass = "ddiMatrix") 275 : maechler 1109 276 : 277 : diagmatprod <- function(x, y) { 278 : dx <- dim(x) 279 : dy <- dim(y) 280 : if(dx[2] != dy[1]) stop("non-matching dimensions") 281 : n <- dx[1] 282 : as(if(x@diag == "U") y else x@x * y, "Matrix") 283 : } 284 : 285 : setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), 286 : maechler 1654 diagmatprod) 287 : maechler 1109 formals(diagmatprod) <- alist(x=, y=NULL) 288 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"), 289 : maechler 1654 diagmatprod) 290 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"), 291 : diagmatprod) 292 : maechler 1109 293 : diagdgeprod <- function(x, y) { 294 : dx <- dim(x) 295 : dy <- dim(y) 296 : if(dx[2] != dy[1]) stop("non-matching dimensions") 297 : if(x@diag != "U") 298 : y@x <- x@x * y@x 299 : y 300 : } 301 : setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"), 302 : diagdgeprod, valueClass = "dgeMatrix") 303 : formals(diagdgeprod) <- alist(x=, y=NULL) 304 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), 305 : diagdgeprod, valueClass = "dgeMatrix") 306 : 307 : setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), 308 : function(x, y) { 309 : maechler 1635 dx <- dim(x) 310 : dy <- dim(y) 311 : if(dx[2] != dy[1]) stop("non-matching dimensions") 312 : as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix") 313 : }) 314 : maechler 1109 315 : setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"), 316 : function(x, y) { 317 : maechler 1635 dx <- dim(x) 318 : dy <- dim(y) 319 : if(dx[2] != dy[1]) stop("non-matching dimensions") 320 : if(y@diag == "N") 321 : x@x <- x@x * rep(y@x, each = dx[1]) 322 : x 323 : }) 324 : maechler 1109 325 : maechler 1295 ## crossprod {more of these} 326 : maechler 1109 327 : maechler 1295 ## tcrossprod --- all are not yet there: do the dense ones here: 328 : maechler 1109 329 : maechler 1295 ## FIXME: 330 : ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"), 331 : ## function(x, y = NULL) { 332 : ## }) 333 : maechler 1109 334 : maechler 1295 ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"), 335 : ## function(x, y = NULL) { 336 : ## }) 337 : maechler 1109 338 : maechler 1799 setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 339 : function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() }) 340 : maechler 1295 341 : maechler 1799 setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 342 : function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() }) 343 : 344 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 345 : function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() }) 346 : 347 : setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 348 : function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() }) 349 : 350 : 351 : ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1() 352 : setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"), 353 : function(x, y) as(x, "sparseMatrix") %*% y) 354 : ## NB: The previous is *not* triggering for "ddi" o "dgC" (= distance 3) 355 : ## since there's a "ddense" o "Csparse" at dist. 2 => triggers first. 356 : ## ==> do this: 357 : setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"), 358 : function(x, y) as(x, "CsparseMatrix") %*% y) 359 : setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"), 360 : function(x, y) x %*% as(y, "CsparseMatrix")) 361 : ## NB: this is *not* needed for Tsparse & Rsparse 362 : ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal* 363 : ## do indeed work by going through sparse (and *not* ddense)! 364 : 365 : setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"), 366 : function(x, y) x %*% as(y, "sparseMatrix")) 367 : 368 : 369 : setMethod("solve", signature(a = "diagonalMatrix", b = "missing"), 370 : function(a, b, ...) { 371 : a@x <- 1/ a@x 372 : a@Dimnames <- a@Dimnames[2:1] 373 : a 374 : }) 375 : 376 : solveDiag <- function(a, b, ...) { 377 : if((n <- a@Dim[1]) != nrow(b)) 378 : stop("incompatible matrix dimensions") 379 : ## trivially invert a 'in place' and multiply: 380 : a@x <- 1/ a@x 381 : a@Dimnames <- a@Dimnames[2:1] 382 : a %*% b 383 : } 384 : setMethod("solve", signature(a = "diagonalMatrix", b = "matrix"), 385 : solveDiag) 386 : setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"), 387 : solveDiag) 388 : 389 : 390 : 391 : 392 : maechler 1654 ### ---------------- diagonal o sparse ----------------------------- 393 : maechler 1295 394 : maechler 1655 395 : ## Use function for several signatures, in order to evade 396 : ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.) 397 : diagOdiag <- function(e1,e2) { # result should also be diagonal 398 : r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible" 399 : if(is.numeric(r)) { 400 : if(is.numeric(e2@x)) { 401 : e2@x <- r; return(.diag.2N(e2)) } 402 : if(!is.numeric(e1@x)) 403 : ## e.g. e1, e2 are logical; 404 : e1 <- as(e1, "dMatrix") 405 : } 406 : else if(is.logical(r)) 407 : e1 <- as(e1, "lMatrix") 408 : else stop("intermediate 'r' is of type", typeof(r)) 409 : e1@x <- r 410 : .diag.2N(e1) 411 : } 412 : 413 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"), 414 : diagOdiag) 415 : ## These two are just for method disambiguation: 416 : setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"), 417 : diagOdiag) 418 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"), 419 : diagOdiag) 420 : 421 : ## For almost everything else, diag* shall be treated "as sparse" : 422 : maechler 1295 ## These are cheap implementations via coercion 423 : 424 : maechler 1655 ## for disambiguation 425 : maechler 1654 setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"), 426 : function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2)) 427 : setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"), 428 : function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix"))) 429 : maechler 1655 ## in general: 430 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"), 431 : function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2)) 432 : setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"), 433 : function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix"))) 434 : maechler 1654 435 : maechler 1655 436 : 437 : maechler 1109 ## similar to prTriang() in ./Auxiliaries.R : 438 : prDiag <- 439 : function(x, digits = getOption("digits"), justify = "none", right = TRUE) 440 : { 441 : cf <- array(".", dim = x@Dim, dimnames = x@Dimnames) 442 : cf[row(cf) == col(cf)] <- 443 : sapply(diag(x), format, digits = digits, justify = justify) 444 : print(cf, quote = FALSE, right = right) 445 : invisible(x) 446 : } 447 : 448 : setMethod("show", signature(object = "diagonalMatrix"), 449 : maechler 1592 function(object) { 450 : d <- dim(object) 451 : cl <- class(object) 452 : cat(sprintf('%d x %d diagonal matrix of class "%s"\n', 453 : d[1], d[2], cl)) 454 : prDiag(object) 455 : })

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