# SCM Repository

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

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

Revision 1707 - (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 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 : ## is better than this: 72 : ## setAs("diagonalMatrix", "sparseMatrix", 73 : ## function(from) 74 : ## as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix")) 75 : setAs("diagonalMatrix", "CsparseMatrix", 76 : function(from) as(diag2T(from), "CsparseMatrix")) 77 : 78 : maechler 1109 setAs("diagonalMatrix", "matrix", 79 : function(from) { 80 : n <- from@Dim[1] 81 : diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE 82 : } else from@x, 83 : nrow = n, ncol = n) 84 : }) 85 : 86 : maechler 1654 setAs("diagonalMatrix", "generalMatrix", # prefer sparse: 87 : function(from) as(from, paste(.M.kind(from), "gCMatrix", sep=''))) 88 : maechler 1174 89 : maechler 1655 .diag.x <- function(m) { 90 : if(m@diag == "U") 91 : rep.int(if(is.numeric(m@x)) 1. else TRUE, 92 : m@Dim[1]) 93 : else m@x 94 : } 95 : 96 : .diag.2N <- function(m) { 97 : if(m@diag == "U") m@diag <- "N" 98 : m 99 : } 100 : 101 : maechler 1654 ## given the above, the following 4 coercions should be all unneeded; 102 : ## we prefer triangular to general: 103 : maechler 1295 setAs("ddiMatrix", "dgTMatrix", 104 : function(from) { 105 : maechler 1654 .Deprecated("as(, \"sparseMatrix\")") 106 : maechler 1295 n <- from@Dim[1] 107 : maechler 1654 i <- seq_len(n) - 1:1 108 : maechler 1655 new("dgTMatrix", i = i, j = i, x = .diag.x(from), 109 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) }) 110 : 111 : setAs("ddiMatrix", "dgCMatrix", 112 : maechler 1655 function(from) as(as(from, "sparseMatrix"), "dgCMatrix")) 113 : maechler 1295 114 : setAs("ldiMatrix", "lgTMatrix", 115 : function(from) { 116 : maechler 1654 .Deprecated("as(, \"sparseMatrix\")") 117 : maechler 1295 n <- from@Dim[1] 118 : maechler 1575 if(from@diag == "U") { # unit-diagonal 119 : x <- rep.int(TRUE, n) 120 : maechler 1654 i <- seq_len(n) - 1:1 121 : maechler 1575 } else { # "normal" 122 : nz <- nz.NA(from@x, na. = TRUE) 123 : x <- from@x[nz] 124 : i <- which(nz) - 1:1 125 : } 126 : new("lgTMatrix", i = i, j = i, x = x, 127 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) }) 128 : 129 : setAs("ldiMatrix", "lgCMatrix", 130 : function(from) as(as(from, "lgTMatrix"), "lgCMatrix")) 131 : 132 : 133 : maechler 1447 if(FALSE) # now have faster "ddense" -> "dge" 134 : maechler 1174 setAs("ddiMatrix", "dgeMatrix", 135 : function(from) as(as(from, "matrix"), "dgeMatrix")) 136 : 137 : maechler 1109 setAs("matrix", "diagonalMatrix", 138 : function(from) { 139 : maechler 1295 d <- dim(from) 140 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix") 141 : if(any(from[row(from) != col(from)] != 0)) 142 : stop("has non-zero off-diagonal entries") 143 : maechler 1295 x <- diag(from) 144 : if(is.logical(x)) { 145 : cl <- "ldiMatrix" 146 : uni <- all(x) 147 : } else { 148 : cl <- "ddiMatrix" 149 : uni <- all(x == 1) 150 : storage.mode(x) <- "double" 151 : maechler 1575 } ## TODO: complex 152 : maechler 1295 new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 153 : x = if(uni) x[FALSE] else x) 154 : maechler 1109 }) 155 : 156 : ## ``generic'' coercion to diagonalMatrix : build on isDiagonal() and diag() 157 : setAs("Matrix", "diagonalMatrix", 158 : function(from) { 159 : d <- dim(from) 160 : if(d[1] != (n <- d[2])) stop("non-square matrix") 161 : if(!isDiagonal(from)) stop("matrix is not diagonal") 162 : ## else: 163 : x <- diag(from) 164 : if(is.logical(x)) { 165 : cl <- "ldiMatrix" 166 : uni <- all(x) 167 : } else { 168 : cl <- "ddiMatrix" 169 : uni <- all(x == 1) 170 : storage.mode(x) <- "double" 171 : } 172 : new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 173 : x = if(uni) x[FALSE] else x) 174 : }) 175 : 176 : maechler 1617 ## When you assign to a diagonalMatrix, the result should be 177 : ## diagonal or sparse 178 : setReplaceMethod("[", signature(x = "diagonalMatrix", 179 : i = "ANY", j = "ANY", value = "ANY"), 180 : function(x, i, j, value) { 181 : r <- callGeneric(x = as(x,"sparseMatrix"), 182 : i=i, j=j, value=value) 183 : if(isDiagonal(r)) as(r, "diagonalMatrix") else r 184 : }) 185 : 186 : 187 : maechler 1109 setMethod("t", signature(x = "diagonalMatrix"), 188 : function(x) { x@Dimnames <- x@Dimnames[2:1] ; x }) 189 : 190 : maechler 1331 setMethod("isDiagonal", signature(object = "diagonalMatrix"), 191 : function(object) TRUE) 192 : setMethod("isTriangular", signature(object = "diagonalMatrix"), 193 : function(object) TRUE) 194 : maechler 1109 setMethod("isSymmetric", signature(object = "diagonalMatrix"), 195 : function(object) TRUE) 196 : 197 : maechler 1654 setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY" 198 : function(x, pivot) { 199 : if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries") 200 : x@x <- sqrt(x@x) 201 : x 202 : }) 203 : ## chol(L) is L for logical diagonal: 204 : setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x) 205 : 206 : 207 : maechler 1109 setMethod("diag", signature(x = "diagonalMatrix"), 208 : function(x = 1, nrow, ncol = n) { 209 : if(x@diag == "U") 210 : rep.int(if(is.logical(x@x)) TRUE else 1, x@Dim[1]) 211 : else x@x 212 : }) 213 : 214 : setMethod("!", "ldiMatrix", function(e1) { 215 : if(e1@diag == "N") 216 : e1@x <- !e1@x 217 : else { ## "U" 218 : e1@diag <- "N" 219 : e1@x <- rep.int(FALSE, e1@Dim[1]) 220 : } 221 : x 222 : }) 223 : 224 : ## Basic Matrix Multiplication {many more to add} 225 : maechler 1654 ## --------------------- 226 : ## Note that "ldi" logical are treated as numeric 227 : maechler 1109 diagdiagprod <- function(x, y) { 228 : if(any(dim(x) != dim(y))) stop("non-matching dimensions") 229 : if(x@diag != "U") { 230 : maechler 1654 if(y@diag != "U") { 231 : nx <- x@x * y@x 232 : if(is.numeric(nx) && !is.numeric(x@x)) 233 : x <- as(x, "dMatrix") 234 : x@x <- as.numeric(nx) 235 : } 236 : return(x) 237 : maechler 1109 } else ## x is unit diagonal 238 : return(y) 239 : } 240 : 241 : maechler 1654 setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 242 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 243 : 244 : maechler 1654 formals(diagdiagprod) <- alist(x=, y=x) 245 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 246 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 247 : maechler 1654 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"), 248 : maechler 1109 diagdiagprod, valueClass = "ddiMatrix") 249 : maechler 1654 setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"), 250 : diagdiagprod, valueClass = "ddiMatrix") 251 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), 252 : diagdiagprod, valueClass = "ddiMatrix") 253 : maechler 1109 254 : 255 : diagmatprod <- function(x, y) { 256 : dx <- dim(x) 257 : dy <- dim(y) 258 : if(dx[2] != dy[1]) stop("non-matching dimensions") 259 : n <- dx[1] 260 : as(if(x@diag == "U") y else x@x * y, "Matrix") 261 : } 262 : 263 : setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), 264 : maechler 1654 diagmatprod) 265 : maechler 1109 formals(diagmatprod) <- alist(x=, y=NULL) 266 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"), 267 : maechler 1654 diagmatprod) 268 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"), 269 : diagmatprod) 270 : maechler 1109 271 : diagdgeprod <- function(x, y) { 272 : dx <- dim(x) 273 : dy <- dim(y) 274 : if(dx[2] != dy[1]) stop("non-matching dimensions") 275 : if(x@diag != "U") 276 : y@x <- x@x * y@x 277 : y 278 : } 279 : setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"), 280 : diagdgeprod, valueClass = "dgeMatrix") 281 : formals(diagdgeprod) <- alist(x=, y=NULL) 282 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), 283 : diagdgeprod, valueClass = "dgeMatrix") 284 : 285 : setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), 286 : function(x, y) { 287 : maechler 1635 dx <- dim(x) 288 : dy <- dim(y) 289 : if(dx[2] != dy[1]) stop("non-matching dimensions") 290 : as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix") 291 : }) 292 : maechler 1109 293 : setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"), 294 : function(x, y) { 295 : maechler 1635 dx <- dim(x) 296 : dy <- dim(y) 297 : if(dx[2] != dy[1]) stop("non-matching dimensions") 298 : if(y@diag == "N") 299 : x@x <- x@x * rep(y@x, each = dx[1]) 300 : x 301 : }) 302 : maechler 1109 303 : maechler 1295 ## crossprod {more of these} 304 : maechler 1109 305 : maechler 1295 ## tcrossprod --- all are not yet there: do the dense ones here: 306 : maechler 1109 307 : maechler 1295 ## FIXME: 308 : ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"), 309 : ## function(x, y = NULL) { 310 : ## }) 311 : maechler 1109 312 : maechler 1295 ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"), 313 : ## function(x, y = NULL) { 314 : ## }) 315 : maechler 1109 316 : maechler 1295 317 : maechler 1654 ### ---------------- diagonal o sparse ----------------------------- 318 : maechler 1295 319 : maechler 1655 320 : ## Use function for several signatures, in order to evade 321 : ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.) 322 : diagOdiag <- function(e1,e2) { # result should also be diagonal 323 : r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible" 324 : if(is.numeric(r)) { 325 : if(is.numeric(e2@x)) { 326 : e2@x <- r; return(.diag.2N(e2)) } 327 : if(!is.numeric(e1@x)) 328 : ## e.g. e1, e2 are logical; 329 : e1 <- as(e1, "dMatrix") 330 : } 331 : else if(is.logical(r)) 332 : e1 <- as(e1, "lMatrix") 333 : else stop("intermediate 'r' is of type", typeof(r)) 334 : e1@x <- r 335 : .diag.2N(e1) 336 : } 337 : 338 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"), 339 : diagOdiag) 340 : ## These two are just for method disambiguation: 341 : setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"), 342 : diagOdiag) 343 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"), 344 : diagOdiag) 345 : 346 : ## For almost everything else, diag* shall be treated "as sparse" : 347 : maechler 1295 ## These are cheap implementations via coercion 348 : 349 : maechler 1655 ## for disambiguation 350 : maechler 1654 setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"), 351 : function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2)) 352 : setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"), 353 : function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix"))) 354 : maechler 1655 ## in general: 355 : setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"), 356 : function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2)) 357 : setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"), 358 : function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix"))) 359 : maechler 1654 360 : maechler 1655 361 : 362 : maechler 1295 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1() 363 : setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"), 364 : function(x, y) as(x, "sparseMatrix") %*% y) 365 : maechler 1682 ## NB: The previous is *not* triggering for "ddi" o "dgC" (= distance 3) 366 : ## since there's a "ddense" o "Csparse" at dist. 2 => triggers first. 367 : ## ==> do this: 368 : setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"), 369 : function(x, y) as(x, "CsparseMatrix") %*% y) 370 : ## NB: this is *not* needed for Tsparse & Rsparse 371 : ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal* 372 : ## do indeed work by going throug sparse (and *not* ddense)! 373 : maechler 1295 374 : setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"), 375 : function(x, y) x %*% as(y, "sparseMatrix")) 376 : 377 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 378 : function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() }) 379 : 380 : setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 381 : function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() }) 382 : 383 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 384 : function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() }) 385 : 386 : setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 387 : function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() }) 388 : 389 : 390 : maechler 1707 ## cbind/rbind: preserve sparseness {not always optimal} 391 : maechler 1295 392 : maechler 1707 ## hack to suppress the obnoxious dispatch ambiguity warnings: 393 : diag2Sp <- function(x) suppressWarnings(as(x, "CsparseMatrix")) 394 : maechler 1295 395 : maechler 1707 ## in order to evade method dispatch ambiguity, but still remain "general" 396 : ## we use this hack instead of signature x = "diagonalMatrix" 397 : for(cls in names(getClass("diagonalMatrix")@subclasses)) { 398 : setMethod("cbind2", signature(x = cls, y = "matrix"), 399 : function(x,y) cbind2(diag2Sp(x), as_Csparse(y))) 400 : setMethod("cbind2", signature(x = "matrix", y = cls), 401 : function(x,y) cbind2(as_Csparse(x), diag2Sp(y))) 402 : setMethod("rbind2", signature(x = cls, y = "matrix"), 403 : function(x,y) rbind2(diag2Sp(x), as_Csparse(y))) 404 : setMethod("rbind2", signature(x = "matrix", y = cls), 405 : function(x,y) rbind2(as_Csparse(x), diag2Sp(y))) 406 : } 407 : 408 : 409 : maechler 1109 ## similar to prTriang() in ./Auxiliaries.R : 410 : prDiag <- 411 : function(x, digits = getOption("digits"), justify = "none", right = TRUE) 412 : { 413 : cf <- array(".", dim = x@Dim, dimnames = x@Dimnames) 414 : cf[row(cf) == col(cf)] <- 415 : sapply(diag(x), format, digits = digits, justify = justify) 416 : print(cf, quote = FALSE, right = right) 417 : invisible(x) 418 : } 419 : 420 : setMethod("show", signature(object = "diagonalMatrix"), 421 : maechler 1592 function(object) { 422 : d <- dim(object) 423 : cl <- class(object) 424 : cat(sprintf('%d x %d diagonal matrix of class "%s"\n', 425 : d[1], d[2], cl)) 426 : prDiag(object) 427 : })

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