# SCM Repository

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

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

Revision 1635 - (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 : 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 : for(i in seq(along = mlist)) { 51 : 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 1109 setAs("diagonalMatrix", "triangularMatrix", 62 : function(from) { 63 : n <- from@Dim[1] 64 : i <- seq(length = n) 65 : x <- from@x 66 : maechler 1575 new(paste(.M.kind(from), "tTMatrix", sep=''), 67 : maechler 1109 diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames, 68 : x = x, i = i, j = i) 69 : }) 70 : 71 : setAs("diagonalMatrix", "matrix", 72 : function(from) { 73 : n <- from@Dim[1] 74 : diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE 75 : } else from@x, 76 : nrow = n, ncol = n) 77 : }) 78 : 79 : maechler 1174 setAs("diagonalMatrix", "generalMatrix", 80 : maechler 1575 function(from) as(from, paste(.M.kind(from), "geMatrix", sep=''))) 81 : maechler 1174 82 : maechler 1295 setAs("ddiMatrix", "dgTMatrix", 83 : function(from) { 84 : n <- from@Dim[1] 85 : i <- seq(length = n) - 1:1 86 : new("dgTMatrix", i = i, j = i, 87 : x = if(from@diag == "U") rep(1,n) else from@x, 88 : Dim = c(n,n), Dimnames = from@Dimnames) }) 89 : 90 : setAs("ddiMatrix", "dgCMatrix", 91 : function(from) as(as(from, "dgTMatrix"), "dgCMatrix")) 92 : 93 : setAs("ldiMatrix", "lgTMatrix", 94 : function(from) { 95 : n <- from@Dim[1] 96 : maechler 1575 if(from@diag == "U") { # unit-diagonal 97 : x <- rep.int(TRUE, n) 98 : i <- seq(length = n) 99 : } else { # "normal" 100 : nz <- nz.NA(from@x, na. = TRUE) 101 : x <- from@x[nz] 102 : i <- which(nz) - 1:1 103 : } 104 : new("lgTMatrix", i = i, j = i, x = x, 105 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) }) 106 : 107 : setAs("ldiMatrix", "lgCMatrix", 108 : function(from) as(as(from, "lgTMatrix"), "lgCMatrix")) 109 : 110 : setAs("diagonalMatrix", "sparseMatrix", 111 : function(from) 112 : as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix")) 113 : 114 : maechler 1447 if(FALSE) # now have faster "ddense" -> "dge" 115 : maechler 1174 setAs("ddiMatrix", "dgeMatrix", 116 : function(from) as(as(from, "matrix"), "dgeMatrix")) 117 : 118 : maechler 1109 setAs("matrix", "diagonalMatrix", 119 : function(from) { 120 : maechler 1295 d <- dim(from) 121 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix") 122 : if(any(from[row(from) != col(from)] != 0)) 123 : stop("has non-zero off-diagonal entries") 124 : maechler 1295 x <- diag(from) 125 : if(is.logical(x)) { 126 : cl <- "ldiMatrix" 127 : uni <- all(x) 128 : } else { 129 : cl <- "ddiMatrix" 130 : uni <- all(x == 1) 131 : storage.mode(x) <- "double" 132 : maechler 1575 } ## TODO: complex 133 : maechler 1295 new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 134 : x = if(uni) x[FALSE] else x) 135 : maechler 1109 }) 136 : 137 : ## ``generic'' coercion to diagonalMatrix : build on isDiagonal() and diag() 138 : setAs("Matrix", "diagonalMatrix", 139 : function(from) { 140 : d <- dim(from) 141 : if(d[1] != (n <- d[2])) stop("non-square matrix") 142 : if(!isDiagonal(from)) stop("matrix is not diagonal") 143 : ## else: 144 : x <- diag(from) 145 : if(is.logical(x)) { 146 : cl <- "ldiMatrix" 147 : uni <- all(x) 148 : } else { 149 : cl <- "ddiMatrix" 150 : uni <- all(x == 1) 151 : storage.mode(x) <- "double" 152 : } 153 : new(cl, Dim = c(n,n), diag = if(uni) "U" else "N", 154 : x = if(uni) x[FALSE] else x) 155 : }) 156 : 157 : maechler 1617 ## When you assign to a diagonalMatrix, the result should be 158 : ## diagonal or sparse 159 : setReplaceMethod("[", signature(x = "diagonalMatrix", 160 : i = "ANY", j = "ANY", value = "ANY"), 161 : function(x, i, j, value) { 162 : r <- callGeneric(x = as(x,"sparseMatrix"), 163 : i=i, j=j, value=value) 164 : if(isDiagonal(r)) as(r, "diagonalMatrix") else r 165 : }) 166 : 167 : 168 : maechler 1109 setMethod("t", signature(x = "diagonalMatrix"), 169 : function(x) { x@Dimnames <- x@Dimnames[2:1] ; x }) 170 : 171 : maechler 1331 setMethod("isDiagonal", signature(object = "diagonalMatrix"), 172 : function(object) TRUE) 173 : setMethod("isTriangular", signature(object = "diagonalMatrix"), 174 : function(object) TRUE) 175 : maechler 1109 setMethod("isSymmetric", signature(object = "diagonalMatrix"), 176 : function(object) TRUE) 177 : 178 : setMethod("diag", signature(x = "diagonalMatrix"), 179 : function(x = 1, nrow, ncol = n) { 180 : if(x@diag == "U") 181 : rep.int(if(is.logical(x@x)) TRUE else 1, x@Dim[1]) 182 : else x@x 183 : }) 184 : 185 : setMethod("!", "ldiMatrix", function(e1) { 186 : if(e1@diag == "N") 187 : e1@x <- !e1@x 188 : else { ## "U" 189 : e1@diag <- "N" 190 : e1@x <- rep.int(FALSE, e1@Dim[1]) 191 : } 192 : x 193 : }) 194 : 195 : ## Basic Matrix Multiplication {many more to add} 196 : 197 : ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix" 198 : diagdiagprod <- function(x, y) { 199 : if(any(dim(x) != dim(y))) stop("non-matching dimensions") 200 : if(x@diag != "U") { 201 : if(y@diag != "U") x@x <- x@x * y@x 202 : return(x) 203 : } else ## x is unit diagonal 204 : return(y) 205 : } 206 : 207 : setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"), 208 : diagdiagprod, valueClass = "ddiMatrix") 209 : 210 : formals(diagdiagprod) <- alist(x=, y=NULL) 211 : setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"), 212 : diagdiagprod, valueClass = "ddiMatrix") 213 : setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"), 214 : diagdiagprod, valueClass = "ddiMatrix") 215 : 216 : 217 : diagmatprod <- function(x, y) { 218 : dx <- dim(x) 219 : dy <- dim(y) 220 : if(dx[2] != dy[1]) stop("non-matching dimensions") 221 : n <- dx[1] 222 : as(if(x@diag == "U") y else x@x * y, "Matrix") 223 : } 224 : 225 : setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), 226 : diagmatprod) 227 : formals(diagmatprod) <- alist(x=, y=NULL) 228 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"), 229 : diagmatprod) 230 : 231 : diagdgeprod <- function(x, y) { 232 : dx <- dim(x) 233 : dy <- dim(y) 234 : if(dx[2] != dy[1]) stop("non-matching dimensions") 235 : if(x@diag != "U") 236 : y@x <- x@x * y@x 237 : y 238 : } 239 : setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"), 240 : diagdgeprod, valueClass = "dgeMatrix") 241 : formals(diagdgeprod) <- alist(x=, y=NULL) 242 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), 243 : diagdgeprod, valueClass = "dgeMatrix") 244 : 245 : setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), 246 : function(x, y) { 247 : maechler 1635 dx <- dim(x) 248 : dy <- dim(y) 249 : if(dx[2] != dy[1]) stop("non-matching dimensions") 250 : as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix") 251 : }) 252 : maechler 1109 253 : setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"), 254 : function(x, y) { 255 : maechler 1635 dx <- dim(x) 256 : dy <- dim(y) 257 : if(dx[2] != dy[1]) stop("non-matching dimensions") 258 : if(y@diag == "N") 259 : x@x <- x@x * rep(y@x, each = dx[1]) 260 : x 261 : }) 262 : maechler 1109 263 : maechler 1295 ## crossprod {more of these} 264 : maechler 1109 265 : maechler 1295 ## tcrossprod --- all are not yet there: do the dense ones here: 266 : maechler 1109 267 : maechler 1295 ## FIXME: 268 : ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"), 269 : ## function(x, y = NULL) { 270 : ## }) 271 : maechler 1109 272 : maechler 1295 ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"), 273 : ## function(x, y = NULL) { 274 : ## }) 275 : maechler 1109 276 : maechler 1295 277 : ### ---------------- diagonal o sparse ----------------------------- 278 : 279 : ## These are cheap implementations via coercion 280 : 281 : ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1() 282 : 283 : setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"), 284 : function(x, y) as(x, "sparseMatrix") %*% y) 285 : 286 : setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"), 287 : function(x, y) x %*% as(y, "sparseMatrix")) 288 : 289 : setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 290 : function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() }) 291 : 292 : setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 293 : function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() }) 294 : 295 : setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"), 296 : function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() }) 297 : 298 : setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"), 299 : function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() }) 300 : 301 : 302 : 303 : 304 : maechler 1109 ## similar to prTriang() in ./Auxiliaries.R : 305 : prDiag <- 306 : function(x, digits = getOption("digits"), justify = "none", right = TRUE) 307 : { 308 : cf <- array(".", dim = x@Dim, dimnames = x@Dimnames) 309 : cf[row(cf) == col(cf)] <- 310 : sapply(diag(x), format, digits = digits, justify = justify) 311 : print(cf, quote = FALSE, right = right) 312 : invisible(x) 313 : } 314 : 315 : setMethod("show", signature(object = "diagonalMatrix"), 316 : maechler 1592 function(object) { 317 : d <- dim(object) 318 : cl <- class(object) 319 : cat(sprintf('%d x %d diagonal matrix of class "%s"\n', 320 : d[1], d[2], cl)) 321 : prDiag(object) 322 : })

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