SCM

SCM Repository

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

Annotation of /pkg/R/diagMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1575 - (view) (download)

1 : maechler 1109 ## Purpose: Constructor of diagonal matrices -- ~= diag() ,
2 :     ## but *not* diag() extractor!
3 :     Diagonal <- function(n, x = NULL)
4 :     {
5 : maechler 1575 ## Allow Diagonal(4) and Diagonal(x=1:5)
6 : maechler 1109 if(missing(n))
7 : maechler 1575 n <- length(x)
8 : maechler 1109 else {
9 : maechler 1575 stopifnot(length(n) == 1, n == as.integer(n), n >= 0)
10 :     n <- as.integer(n)
11 : maechler 1109 }
12 :    
13 :     if(missing(x)) # unit diagonal matrix
14 : maechler 1575 new("ddiMatrix", Dim = c(n,n), diag = "U")
15 : maechler 1109 else {
16 : maechler 1575 stopifnot(length(x) == n)
17 :     if(is.logical(x))
18 :     cl <- "ldiMatrix"
19 :     else if(is.numeric(x)) {
20 :     cl <- "ddiMatrix"
21 :     x <- as.numeric(x)
22 :     }
23 :     else if(is.complex(x)) {
24 :     cl <- "zdiMatrix" # will not yet work
25 :     } else stop("'x' has invalid data type")
26 :     new(cl, Dim = c(n,n), diag = "N", x = x)
27 : maechler 1109 }
28 :     }
29 :    
30 :     setAs("diagonalMatrix", "triangularMatrix",
31 :     function(from) {
32 :     n <- from@Dim[1]
33 :     i <- seq(length = n)
34 :     x <- from@x
35 : maechler 1575 new(paste(.M.kind(from), "tTMatrix", sep=''),
36 : maechler 1109 diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
37 :     x = x, i = i, j = i)
38 :     })
39 :    
40 :     setAs("diagonalMatrix", "matrix",
41 :     function(from) {
42 :     n <- from@Dim[1]
43 :     diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE
44 :     } else from@x,
45 :     nrow = n, ncol = n)
46 :     })
47 :    
48 : maechler 1174 setAs("diagonalMatrix", "generalMatrix",
49 : maechler 1575 function(from) as(from, paste(.M.kind(from), "geMatrix", sep='')))
50 : maechler 1174
51 : maechler 1295 setAs("ddiMatrix", "dgTMatrix",
52 :     function(from) {
53 :     n <- from@Dim[1]
54 :     i <- seq(length = n) - 1:1
55 :     new("dgTMatrix", i = i, j = i,
56 :     x = if(from@diag == "U") rep(1,n) else from@x,
57 :     Dim = c(n,n), Dimnames = from@Dimnames) })
58 :    
59 :     setAs("ddiMatrix", "dgCMatrix",
60 :     function(from) as(as(from, "dgTMatrix"), "dgCMatrix"))
61 :    
62 :     setAs("ldiMatrix", "lgTMatrix",
63 :     function(from) {
64 :     n <- from@Dim[1]
65 : maechler 1575 if(from@diag == "U") { # unit-diagonal
66 :     x <- rep.int(TRUE, n)
67 :     i <- seq(length = n)
68 :     } else { # "normal"
69 :     nz <- nz.NA(from@x, na. = TRUE)
70 :     x <- from@x[nz]
71 :     i <- which(nz) - 1:1
72 :     }
73 :     new("lgTMatrix", i = i, j = i, x = x,
74 : maechler 1295 Dim = c(n,n), Dimnames = from@Dimnames) })
75 :    
76 :     setAs("ldiMatrix", "lgCMatrix",
77 :     function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))
78 :    
79 :     setAs("diagonalMatrix", "sparseMatrix",
80 :     function(from)
81 :     as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))
82 :    
83 : maechler 1447 if(FALSE) # now have faster "ddense" -> "dge"
84 : maechler 1174 setAs("ddiMatrix", "dgeMatrix",
85 :     function(from) as(as(from, "matrix"), "dgeMatrix"))
86 :    
87 : maechler 1109 setAs("matrix", "diagonalMatrix",
88 :     function(from) {
89 : maechler 1295 d <- dim(from)
90 : maechler 1109 if(d[1] != (n <- d[2])) stop("non-square matrix")
91 :     if(any(from[row(from) != col(from)] != 0))
92 :     stop("has non-zero off-diagonal entries")
93 : maechler 1295 x <- diag(from)
94 :     if(is.logical(x)) {
95 :     cl <- "ldiMatrix"
96 :     uni <- all(x)
97 :     } else {
98 :     cl <- "ddiMatrix"
99 :     uni <- all(x == 1)
100 :     storage.mode(x) <- "double"
101 : maechler 1575 } ## TODO: complex
102 : maechler 1295 new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
103 :     x = if(uni) x[FALSE] else x)
104 : maechler 1109 })
105 :    
106 :     ## ``generic'' coercion to diagonalMatrix : build on isDiagonal() and diag()
107 :     setAs("Matrix", "diagonalMatrix",
108 :     function(from) {
109 :     d <- dim(from)
110 :     if(d[1] != (n <- d[2])) stop("non-square matrix")
111 :     if(!isDiagonal(from)) stop("matrix is not diagonal")
112 :     ## else:
113 :     x <- diag(from)
114 :     if(is.logical(x)) {
115 :     cl <- "ldiMatrix"
116 :     uni <- all(x)
117 :     } else {
118 :     cl <- "ddiMatrix"
119 :     uni <- all(x == 1)
120 :     storage.mode(x) <- "double"
121 :     }
122 :     new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
123 :     x = if(uni) x[FALSE] else x)
124 :     })
125 :    
126 :     setMethod("t", signature(x = "diagonalMatrix"),
127 :     function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
128 :    
129 : maechler 1331 setMethod("isDiagonal", signature(object = "diagonalMatrix"),
130 :     function(object) TRUE)
131 :     setMethod("isTriangular", signature(object = "diagonalMatrix"),
132 :     function(object) TRUE)
133 : maechler 1109 setMethod("isSymmetric", signature(object = "diagonalMatrix"),
134 :     function(object) TRUE)
135 :    
136 :     setMethod("diag", signature(x = "diagonalMatrix"),
137 :     function(x = 1, nrow, ncol = n) {
138 :     if(x@diag == "U")
139 :     rep.int(if(is.logical(x@x)) TRUE else 1, x@Dim[1])
140 :     else x@x
141 :     })
142 :    
143 :     setMethod("!", "ldiMatrix", function(e1) {
144 :     if(e1@diag == "N")
145 :     e1@x <- !e1@x
146 :     else { ## "U"
147 :     e1@diag <- "N"
148 :     e1@x <- rep.int(FALSE, e1@Dim[1])
149 :     }
150 :     x
151 :     })
152 :    
153 :     ## Basic Matrix Multiplication {many more to add}
154 :    
155 :     ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix"
156 :     diagdiagprod <- function(x, y) {
157 :     if(any(dim(x) != dim(y))) stop("non-matching dimensions")
158 :     if(x@diag != "U") {
159 :     if(y@diag != "U") x@x <- x@x * y@x
160 :     return(x)
161 :     } else ## x is unit diagonal
162 :     return(y)
163 :     }
164 :    
165 :     setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"),
166 :     diagdiagprod, valueClass = "ddiMatrix")
167 :    
168 :     formals(diagdiagprod) <- alist(x=, y=NULL)
169 :     setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),
170 :     diagdiagprod, valueClass = "ddiMatrix")
171 :     setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),
172 :     diagdiagprod, valueClass = "ddiMatrix")
173 :    
174 :    
175 :     diagmatprod <- function(x, y) {
176 :     dx <- dim(x)
177 :     dy <- dim(y)
178 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
179 :     n <- dx[1]
180 :     as(if(x@diag == "U") y else x@x * y, "Matrix")
181 :     }
182 :    
183 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
184 :     diagmatprod)
185 :     formals(diagmatprod) <- alist(x=, y=NULL)
186 :     setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
187 :     diagmatprod)
188 :    
189 :     diagdgeprod <- function(x, y) {
190 :     dx <- dim(x)
191 :     dy <- dim(y)
192 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
193 :     if(x@diag != "U")
194 :     y@x <- x@x * y@x
195 :     y
196 :     }
197 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),
198 :     diagdgeprod, valueClass = "dgeMatrix")
199 :     formals(diagdgeprod) <- alist(x=, y=NULL)
200 :     setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
201 :     diagdgeprod, valueClass = "dgeMatrix")
202 :    
203 :     setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
204 :     function(x, y) {
205 :     dx <- dim(x)
206 :     dy <- dim(y)
207 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
208 :     as(if(y@diag == "U") x else x * rep.int(y@x, dx[1]), "Matrix")
209 :     })
210 :    
211 :     setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),
212 :     function(x, y) {
213 :     dx <- dim(x)
214 :     dy <- dim(y)
215 :     if(dx[2] != dy[1]) stop("non-matching dimensions")
216 :     if(y@diag == "N")
217 :     x@x <- x@x * rep.int(y@x, dx[1])
218 :     x
219 :     })
220 :    
221 : maechler 1295 ## crossprod {more of these}
222 : maechler 1109
223 : maechler 1295 ## tcrossprod --- all are not yet there: do the dense ones here:
224 : maechler 1109
225 : maechler 1295 ## FIXME:
226 :     ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
227 :     ## function(x, y = NULL) {
228 :     ## })
229 : maechler 1109
230 : maechler 1295 ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),
231 :     ## function(x, y = NULL) {
232 :     ## })
233 : maechler 1109
234 : maechler 1295
235 :     ### ---------------- diagonal o sparse -----------------------------
236 :    
237 :     ## These are cheap implementations via coercion
238 :    
239 :     ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
240 :    
241 :     setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
242 :     function(x, y) as(x, "sparseMatrix") %*% y)
243 :    
244 :     setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
245 :     function(x, y) x %*% as(y, "sparseMatrix"))
246 :    
247 :     setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
248 :     function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })
249 :    
250 :     setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
251 :     function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })
252 :    
253 :     setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
254 :     function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })
255 :    
256 :     setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
257 :     function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })
258 :    
259 :    
260 :    
261 :    
262 : maechler 1109 ## similar to prTriang() in ./Auxiliaries.R :
263 :     prDiag <-
264 :     function(x, digits = getOption("digits"), justify = "none", right = TRUE)
265 :     {
266 :     cf <- array(".", dim = x@Dim, dimnames = x@Dimnames)
267 :     cf[row(cf) == col(cf)] <-
268 :     sapply(diag(x), format, digits = digits, justify = justify)
269 :     print(cf, quote = FALSE, right = right)
270 :     invisible(x)
271 :     }
272 :    
273 :     setMethod("show", signature(object = "diagonalMatrix"),
274 :     function(object) prDiag(object))
275 :    

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge