SCM

SCM Repository

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

Annotation of /pkg/Matrix/R/Matrix.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3020 - (view) (download)

1 : maechler 607 #### Toplevel ``virtual'' class "Matrix"
2 :    
3 : maechler 1331
4 :     ### Virtual coercions -- via smart "helpers" (-> ./Auxiliaries.R)
5 :    
6 : maechler 1799 setAs("Matrix", "sparseMatrix", function(from) as(from, "CsparseMatrix"))
7 :     setAs("Matrix", "CsparseMatrix", function(from) as_Csparse(from))
8 : maechler 1331 setAs("Matrix", "denseMatrix", function(from) as_dense(from))
9 :    
10 : maechler 1799 ## Maybe TODO:
11 :     ## setAs("Matrix", "nMatrix", function(from) ....)
12 :    
13 : maechler 1618 ## Most of these work; this is a last resort:
14 : maechler 1548 setAs(from = "Matrix", to = "matrix", # do *not* call base::as.matrix() here:
15 :     function(from) .bail.out.2("coerce", class(from), class(to)))
16 : maechler 1618 setAs(from = "matrix", to = "Matrix", function(from) Matrix(from))
17 : maechler 1548
18 : maechler 956 ## ## probably not needed eventually:
19 :     ## setAs(from = "ddenseMatrix", to = "matrix",
20 :     ## function(from) {
21 :     ## if(length(d <- dim(from)) != 2) stop("dim(.) has not length 2")
22 :     ## array(from@x, dim = d, dimnames = dimnames(from))
23 :     ## })
24 : maechler 607
25 : maechler 676 ## should propagate to all subclasses:
26 : maechler 658 setMethod("as.matrix", signature(x = "Matrix"), function(x) as(x, "matrix"))
27 : maechler 956 ## for 'Matrix' objects, as.array() should be equivalent:
28 :     setMethod("as.array", signature(x = "Matrix"), function(x) as(x, "matrix"))
29 : mmaechler 2676 ## Such that also base functions dispatch properly on our classes:
30 :     as.array.Matrix <- as.matrix.Matrix <- function(x, ...) as(x, "matrix")
31 : maechler 658
32 : maechler 1347 ## head and tail apply to all Matrix objects for which subscripting is allowed:
33 : maechler 1513 setMethod("head", signature(x = "Matrix"), utils::head.matrix)
34 :     setMethod("tail", signature(x = "Matrix"), utils::tail.matrix)
35 : maechler 1390
36 : maechler 1751 setMethod("drop", signature(x = "Matrix"),
37 :     function(x) if(all(dim(x) != 1)) x else drop(as(x, "matrix")))
38 :    
39 : maechler 956 ## slow "fall back" method {subclasses should have faster ones}:
40 :     setMethod("as.vector", signature(x = "Matrix", mode = "missing"),
41 : mmaechler 2175 function(x, mode) as.vector(as(x, "matrix"), mode))
42 : mmaechler 2676 ## so base functions calling as.vector() work too:
43 :     as.vector.Matrix <- function(x, mode) as.vector(as(x, "matrix"), mode)
44 : maechler 956
45 : mmaechler 2363 setAs("Matrix", "vector", function(from) as.vector (as(from, "matrix")))
46 :     setAs("Matrix", "numeric", function(from) as.numeric(as(from, "matrix")))
47 :     setAs("Matrix", "logical", function(from) as.logical(as(from, "matrix")))
48 :     setAs("Matrix", "integer", function(from) as.integer(as(from, "matrix")))
49 :     setAs("Matrix", "complex", function(from) as.complex(as(from, "matrix")))
50 :    
51 : maechler 1389 ## mainly need these for "dMatrix" or "lMatrix" respectively, but why not general:
52 :     setMethod("as.numeric", signature(x = "Matrix"),
53 :     function(x, ...) as.numeric(as.vector(x)))
54 :     setMethod("as.logical", signature(x = "Matrix"),
55 :     function(x, ...) as.logical(as.vector(x)))
56 : maechler 956
57 : mmaechler 2553 setMethod("mean", signature(x = "sparseMatrix"),
58 :     function(x, ...) mean(as(x,"sparseVector"), ...))
59 :     setMethod("mean", signature(x = "sparseVector"),
60 :     function(x, trim = 0, na.rm = FALSE, ...)
61 :     {
62 :     if (na.rm) # remove NAs such that new length() is ok
63 :     x <- x[!is.na(x)] # remains sparse!
64 :     if(is0(trim)) sum(x) / length(x)
65 :     else {
66 :     ## fast trimmed mean for sparseVector:
67 :     ## ---> we'd need fast & sparse sort(<sparseV>).
68 :     ## Normally this means to define a xtfrm() method;
69 :     ## however, that plus x[order(x, ..)] will NOT be sparse
70 :     ## TODO: sortSparseVector(.)
71 :     warning("trimmed mean of 'sparseVector' -- suboptimally using as.numeric(.)")
72 :     mean(as.numeric(x), trim=trim)
73 :     }
74 :     })
75 :     ## for the non-"sparseMatrix" ones:
76 : mmaechler 2175 setMethod("mean", signature(x = "Matrix"),
77 : mmaechler 2553 function(x, trim = 0, na.rm = FALSE, ...)
78 :     {
79 :     if (na.rm)
80 :     x <- x[!is.na(x)]
81 :     if(is0(trim)) sum(x) / length(x)
82 :     else mean(as.numeric(x), trim=trim)
83 :     })
84 : mmaechler 2175
85 : mmaechler 2553
86 : mmaechler 2906 ## for non-"sparseMatrix" :
87 : maechler 2096 setMethod("cov2cor", signature(V = "Matrix"),
88 : maechler 2115 function(V) { ## was as(cov2cor(as(V, "matrix")), "dpoMatrix"))
89 :     r <- V
90 :     p <- (d <- dim(V))[1]
91 :     if(p != d[2]) stop("'V' is not a square matrix")
92 :     Is <- sqrt(1/diag(V)) # diag( 1/sigma_i )
93 :     if(any(!is.finite(Is)))
94 :     warning("diag(.) had 0 or NA entries; non-finite result is doubtful")
95 :     Is <- Diagonal(x = Is)
96 :     r <- Is %*% V %*% Is
97 : mmaechler 2411 r[cbind(1:p,1:p)] <- 1 # exact in diagonal
98 : maechler 2115 as(forceSymmetric(r), "dpoMatrix")
99 :     })
100 : maechler 1389
101 : maechler 1467 ## "base" has an isSymmetric() S3-generic since R 2.3.0
102 : maechler 1108 setMethod("isSymmetric", signature(object = "symmetricMatrix"),
103 : maechler 2112 function(object, ...) TRUE)
104 : maechler 1108 setMethod("isSymmetric", signature(object = "triangularMatrix"),
105 : maechler 2112 ## TRUE iff diagonal:
106 :     function(object, ...) isDiagonal(object))
107 : maechler 1189
108 : maechler 1238 setMethod("isTriangular", signature(object = "matrix"), isTriMat)
109 :    
110 : maechler 1174 setMethod("isDiagonal", signature(object = "matrix"), .is.diagonal)
111 : maechler 886
112 : maechler 2112 ## The "catch all" methods -- far from optimal:
113 : mmaechler 3020 setMethod("symmpart", signature(x = "Matrix"), function(x)
114 :     as(symmetrizeDimnames(x + t(x))/2, "symmetricMatrix"))
115 :     setMethod("skewpart", signature(x = "Matrix"), function(x) symmetrizeDimnames(x - t(x))/2)
116 : maechler 1174
117 : maechler 2112 ## FIXME: do this (similarly as for "ddense.." in C
118 : mmaechler 3020 setMethod("symmpart", signature(x = "matrix"), function(x) symmetrizeDimnames(x + t(x))/2)
119 :     setMethod("skewpart", signature(x = "matrix"), function(x) symmetrizeDimnames(x - t(x))/2)
120 : maechler 1174
121 : maechler 2112
122 : mmaechler 2980 if(getRversion() >= "3.1.0")
123 : mmaechler 2963 ## NB: ./nsparseMatrix.R and ./sparseVector.R have extra methods
124 : mmaechler 2980 setMethod("anyNA", signature(x = "xMatrix"),
125 :     function(x) anyNA(x@x))
126 : maechler 2112
127 : mmaechler 2980
128 : maechler 579 setMethod("dim", signature(x = "Matrix"),
129 : maechler 698 function(x) x@Dim, valueClass = "integer")
130 : maechler 1705
131 : maechler 1738 setMethod("length", "Matrix", function(x) prod(dim(x)))
132 :    
133 : maechler 579 setMethod("dimnames", signature(x = "Matrix"), function(x) x@Dimnames)
134 : maechler 1705
135 :    
136 : maechler 607 ## not exported but used more than once for "dimnames<-" method :
137 :     ## -- or do only once for all "Matrix" classes ??
138 :     dimnamesGets <- function (x, value) {
139 :     d <- dim(x)
140 :     if (!is.list(value) || length(value) != 2 ||
141 : maechler 698 !(is.null(v1 <- value[[1]]) || length(v1) == d[1]) ||
142 :     !(is.null(v2 <- value[[2]]) || length(v2) == d[2]))
143 : mmaechler 2863 stop(gettextf("invalid dimnames given for %s object", dQuote(class(x))),
144 :     domain=NA)
145 : mmaechler 3020 x@Dimnames <- # preserve names(value)!
146 :     lapply(value, function(v) if(!is.null(v)) as.character(v))
147 : maechler 607 x
148 :     }
149 :     setMethod("dimnames<-", signature(x = "Matrix", value = "list"),
150 : maechler 698 dimnamesGets)
151 : maechler 579
152 : mmaechler 2193 setMethod("dimnames<-", signature(x = "Matrix", value = "NULL"),
153 :     function(x, value) {
154 : mmaechler 2863 message("dimnames(.) <- NULL: translated to \ndimnames(.) <- list(NULL,NULL) <==> unname(.)")
155 : mmaechler 2193 x@Dimnames <- list(NULL,NULL)
156 :     x
157 :     })
158 :    
159 : maechler 658 setMethod("unname", signature("Matrix", force="missing"),
160 : maechler 698 function(obj) { obj@Dimnames <- list(NULL,NULL); obj})
161 : maechler 607
162 : maechler 2015
163 : maechler 2115 Matrix <- function (data = NA, nrow = 1, ncol = 1, byrow = FALSE,
164 : mmaechler 2725 dimnames = NULL, sparse = NULL,
165 :     doDiag = TRUE, forceCheck = FALSE)
166 : bates 10 {
167 : maechler 1665 sparseDefault <- function(m) prod(dim(m)) > 2*sum(isN0(as(m, "matrix")))
168 : maechler 954
169 :     i.M <- is(data, "Matrix")
170 : mmaechler 2375 sM <- FALSE
171 : mmaechler 2341 if(i.M) {
172 :     if(is(data, "diagonalMatrix")) return(data) # in all cases
173 : mmaechler 2375 sV <- FALSE
174 : mmaechler 2341 } else if(inherits(data, "table")) # special treatment
175 : maechler 2043 class(data) <- "matrix" # "matrix" first for S4 dispatch
176 : mmaechler 2375 else if(is(data, "sparseVector")) {
177 :     data <- spV2M(data, nrow, ncol, byrow=byrow)
178 :     i.M <- sparse <- forceCheck <- sM <- sV <- TRUE
179 :     }
180 : maechler 1618 if(is.null(sparse1 <- sparse) && (i.M || is(data, "matrix")))
181 : maechler 1174 sparse <- sparseDefault(data)
182 : maechler 1455 doDN <- TRUE
183 : maechler 1551 if (i.M) {
184 : mmaechler 2375 if (!sV) {
185 :     if(!missing(nrow) || !missing(ncol)|| !missing(byrow))
186 :     warning("'nrow', 'ncol', etc, are disregarded when 'data' is \"Matrix\" already")
187 :     sM <- is(data,"sparseMatrix")
188 :     if(!forceCheck && ((sparse && sM) || (!sparse && !sM)))
189 :     return(data)
190 :     ## else : convert dense <-> sparse -> at end
191 :     }
192 : maechler 954 }
193 :     else if (!is.matrix(data)) { ## cut & paste from "base::matrix" :
194 : mmaechler 2627 ## avoid copying to strip attributes in simple cases
195 :     if (is.object(data) || !is.atomic(data)) data <- as.vector(data)
196 : maechler 1618 if(length(data) == 1 && is0(data) && !identical(sparse, FALSE)) {
197 : maechler 2115 ## Matrix(0, ...) : always sparse unless "sparse = FALSE":
198 : maechler 1618 if(is.null(sparse)) sparse1 <- sparse <- TRUE
199 : maechler 2115 i.M <- sM <- TRUE
200 : mmaechler 2627 if (missing(nrow)) nrow <- ceiling(1/ncol) else
201 :     if (missing(ncol)) ncol <- ceiling(1/nrow)
202 : maechler 2115 isSym <- nrow == ncol
203 : maechler 1315 ## will be sparse: do NOT construct full matrix!
204 : mmaechler 2837 data <- new(paste0(if(is.numeric(data)) "d" else
205 :     if(is.logical(data)) "l" else
206 :     stop("invalid 'data'"),
207 :     if(isSym) "s" else "g", "CMatrix"),
208 : maechler 2115 p = rep.int(0L, ncol+1L),
209 : maechler 1315 Dim = as.integer(c(nrow,ncol)),
210 : mmaechler 2865 Dimnames = if(is.null.DN(dimnames)) list(NULL,NULL)
211 : maechler 1332 else dimnames)
212 : mmaechler 2837 } else { ## normal case
213 : mmaechler 2845 ## Now 'forbidden' :
214 : mmaechler 2843 ## data <- .Internal(matrix(data, nrow, ncol, byrow, dimnames,
215 : mmaechler 2845 ## missing(nrow), missing(ncol)))
216 : mmaechler 2843 data <- .External(Mmatrix,
217 : mmaechler 2845 data, nrow, ncol, byrow, dimnames,
218 :     missing(nrow), missing(ncol))
219 : maechler 1315 if(is.null(sparse))
220 :     sparse <- sparseDefault(data)
221 :     }
222 : maechler 1455 doDN <- FALSE
223 : maechler 1751 } else if(!missing(nrow) || !missing(ncol)|| !missing(byrow))
224 : maechler 1747 warning("'nrow', 'ncol', etc, are disregarded for matrix 'data'")
225 :    
226 : maechler 1455 ## 'data' is now a "matrix" or "Matrix"
227 :     if (doDN && !is.null(dimnames))
228 : maechler 1189 dimnames(data) <- dimnames
229 : maechler 954
230 : maechler 1174 ## check for symmetric / triangular / diagonal :
231 :     isSym <- isSymmetric(data)
232 :     if((isTri <- !isSym))
233 :     isTri <- isTriangular(data)
234 :     isDiag <- isSym # cannot be diagonal if it isn't symmetric
235 : mmaechler 2341 if(isDiag) # do not *build* 1 x 1 diagonalMatrix
236 : mmaechler 2725 isDiag <- doDiag && !isTRUE(sparse1) && nrow(data) > 1 && isDiagonal(data)
237 : maechler 1174
238 : maechler 2115 ## try to coerce ``via'' virtual classes
239 :     if(isDiag) { ## diagonal is preferred to sparse !
240 :     data <- as(data, "diagonalMatrix")
241 :     isSym <- FALSE
242 :     } else if(sparse && !sM)
243 :     data <- as(data, "sparseMatrix")
244 :     else if(!sparse) {
245 :     if(i.M) { ## data is 'Matrix'
246 :     if(!is(data, "denseMatrix"))
247 :     data <- as(data, "denseMatrix")
248 :     } else { ## data is "matrix" (and result "dense" -> go via "general"
249 :     ctype <- typeof(data)
250 :     if (ctype == "complex")
251 : maechler 1174 stop("complex matrices not yet implemented in Matrix package")
252 : maechler 2115 if (ctype == "integer") ## integer Matrices not yet implemented
253 :     storage.mode(data) <- "double"
254 : mmaechler 2778 data <- new(paste0(.M.kind(data), "geMatrix"),
255 : maechler 2115 Dim = dim(data),
256 :     Dimnames = .M.DN(data),
257 :     x = c(data))
258 : maechler 1174 }
259 : maechler 2115 }
260 : maechler 1174
261 : maechler 2115 if(isTri && !is(data, "triangularMatrix")) {
262 :     data <- if(attr(isTri,"kind") == "L") tril(data) else triu(data)
263 :     #was as(data, "triangularMatrix")
264 :     } else if(isSym && !is(data, "symmetricMatrix"))
265 :     data <- forceSymmetric(data) #was as(data, "symmetricMatrix")
266 :    
267 :     data
268 : bates 10 }
269 :    
270 : bates 683 ## Methods for operations where one argument is numeric
271 : maechler 512
272 : maechler 2043 ## maybe not 100% optimal, but elegant:
273 : maechler 1665 setMethod("solve", signature(a = "Matrix", b = "missing"),
274 :     function(a, b, ...) solve(a, Diagonal(nrow(a))))
275 :    
276 : bates 683 setMethod("solve", signature(a = "Matrix", b = "numeric"),
277 : maechler 2036 function(a, b, ...) callGeneric(a, Matrix(b)))
278 :     setMethod("solve", signature(a = "Matrix", b = "matrix"),
279 :     function(a, b, ...) callGeneric(a, Matrix(b)))
280 :     setMethod("solve", signature(a = "matrix", b = "Matrix"),
281 :     function(a, b, ...) callGeneric(Matrix(a), b))
282 :    
283 : mmaechler 2236 setMethod("solve", signature(a = "Matrix", b = "diagonalMatrix"),
284 :     function(a, b, ...) callGeneric(a, as(b,"CsparseMatrix")))
285 :    
286 : maechler 1659 ## when no sub-class method is found, bail out
287 : mmaechler 2256 setMethod("solve", signature(a = "Matrix", b = "ANY"),
288 : maechler 1747 function(a, b, ...) .bail.out.2("solve", class(a), class(b)))
289 : mmaechler 2256 setMethod("solve", signature(a = "ANY", b = "Matrix"),
290 :     function(a, b, ...) .bail.out.2("solve", class(a), class(b)))
291 : maechler 848
292 : mmaechler 2472 setMethod("chol2inv", signature(x = "denseMatrix"),
293 : mmaechler 2485 function (x, ...) chol2inv(as(as(x, "dMatrix"), "dtrMatrix"), ...))
294 : mmaechler 2770 setMethod("chol2inv", signature(x = "diagonalMatrix"),
295 :     function (x, ...) {
296 : mmaechler 2992 chk.s(..., which.call=-2)
297 : mmaechler 2770 tcrossprod(solve(x))
298 :     })
299 : mmaechler 2485 setMethod("chol2inv", signature(x = "sparseMatrix"),
300 :     function (x, ...) {
301 : mmaechler 2992 chk.s(..., which.call=-2)
302 : mmaechler 2485 ## for now:
303 :     tcrossprod(solve(as(x,"triangularMatrix")))
304 :     })
305 :    
306 : mmaechler 2872 ## There are special sparse methods in ./kronecker.R ; this is a "fall back":
307 : maechler 956 setMethod("kronecker", signature(X = "Matrix", Y = "ANY",
308 : maechler 1805 FUN = "ANY", make.dimnames = "ANY"),
309 :     function(X, Y, FUN, make.dimnames, ...) {
310 : maechler 1825 if(is(X, "sparseMatrix"))
311 :     warning("using slow kronecker() method")
312 : maechler 1805 X <- as(X, "matrix") ; Matrix(callGeneric()) })
313 :    
314 : maechler 956 setMethod("kronecker", signature(X = "ANY", Y = "Matrix",
315 : maechler 1805 FUN = "ANY", make.dimnames = "ANY"),
316 :     function(X, Y, FUN, make.dimnames, ...) {
317 : maechler 1825 if(is(Y, "sparseMatrix"))
318 :     warning("using slow kronecker() method")
319 : maechler 1805 Y <- as(Y, "matrix") ; Matrix(callGeneric()) })
320 : maechler 956
321 :    
322 : mmaechler 2175 setMethod("determinant", signature(x = "Matrix", logarithm = "missing"),
323 :     function(x, logarithm, ...)
324 :     determinant(x, logarithm = TRUE, ...))
325 :    
326 : mmaechler 2872 ## The ``Right Thing'' to do :
327 : mmaechler 2175 ## base::det() calls [base::]determinant();
328 :     ## our det() should call our determinant() :
329 :     det <- base::det
330 : mmaechler 2520 environment(det) <- environment()## == asNamespace("Matrix")
331 : mmaechler 2175
332 : mmaechler 2442 setMethod("Cholesky", signature(A = "Matrix"),
333 :     function(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...)
334 :     stop(gettextf("Cholesky(A) called for 'A' of class \"%s\";\n\t it is currently defined for sparseMatrix only; consider using chol() instead",
335 : mmaechler 2857 class(A)), call. = FALSE, domain=NA))
336 : mmaechler 2442
337 : maechler 1654 ## FIXME: All of these should never be called
338 :     setMethod("chol", signature(x = "Matrix"),
339 : mmaechler 2440 function(x, pivot, ...) .bail.out.1("chol", class(x)))
340 : mmaechler 2175 setMethod("determinant", signature(x = "Matrix", logarithm = "logical"),
341 : mmaechler 2185 function(x, logarithm, ...)
342 :     determinant(as(x,"dMatrix"), logarithm=logarithm, ...))
343 : maechler 1654
344 : maechler 1315 setMethod("diag", signature(x = "Matrix"),
345 : mmaechler 2440 function(x, nrow, ncol) .bail.out.1("diag", class(x)))
346 : mmaechler 2811 if(FALSE)## TODO: activate later
347 :     setMethod("diag<-", signature(x = "Matrix"),
348 :     function(x, value) .bail.out.1("diag", class(x)))
349 : maechler 946 setMethod("t", signature(x = "Matrix"),
350 :     function(x) .bail.out.1(.Generic, class(x)))
351 : maechler 925
352 : maechler 2005 setMethod("norm", signature(x = "Matrix", type = "character"),
353 : mmaechler 2440 function(x, type, ...) .bail.out.1("norm", class(x)))
354 : maechler 2158 setMethod("rcond", signature(x = "Matrix", norm = "character"),
355 : mmaechler 2440 function(x, norm, ...) .bail.out.1("rcond", class(x)))
356 : maechler 2005
357 :    
358 :     ## for all :
359 :     setMethod("norm", signature(x = "ANY", type = "missing"),
360 :     function(x, type, ...) norm(x, type = "O", ...))
361 : maechler 2158 setMethod("rcond", signature(x = "ANY", norm = "missing"),
362 :     function(x, norm, ...) rcond(x, norm = "O", ...))
363 : maechler 2005
364 : mmaechler 2534 setMethod("lu", "matrix", function(x, warnSing = TRUE, ...)
365 : mmaechler 2960 lu(..2dge(x), warnSing=warnSing, ...))
366 : maechler 2005
367 :    
368 :    
369 : mmaechler 2245 ## We want to use all.equal.numeric() *and* make sure that uses
370 :     ## not just base::as.vector but the generic with our methods:
371 :     all.equal_num <- base::all.equal.numeric ## from <R>/src/library/base/R/all.equal.R
372 :     environment(all.equal_num) <- environment()## == as.environment("Matrix")
373 : mmaechler 2255
374 : mmaechler 2445 all.equal_Mat <- function(target, current, check.attributes = TRUE,
375 :     factorsCheck = FALSE, ...)
376 : mmaechler 2255 {
377 : mmaechler 2445 msg <- attr.all_Mat(target, current, check.attributes=check.attributes,
378 :     factorsCheck=factorsCheck, ...)
379 : mmaechler 2906 if(is.list(msg)) msg[[1]]
380 :     else .a.e.comb(msg,
381 :     all.equal_num(as.vector(target), as.vector(current),
382 :     check.attributes=check.attributes, ...))
383 : mmaechler 2255 }
384 : mmaechler 2245 ## The all.equal() methods for dense matrices (and fallback):
385 :     setMethod("all.equal", c(target = "Matrix", current = "Matrix"),
386 : mmaechler 2256 all.equal_Mat)
387 : mmaechler 2245 setMethod("all.equal", c(target = "Matrix", current = "ANY"),
388 : mmaechler 2256 all.equal_Mat)
389 : mmaechler 2245 setMethod("all.equal", c(target = "ANY", current = "Matrix"),
390 : mmaechler 2256 all.equal_Mat)
391 : mmaechler 2245 ## -> ./sparseMatrix.R, ./sparseVector.R have specific methods
392 : maechler 2005
393 : mmaechler 2245
394 :    
395 : maechler 1799 ## MM: More or less "Cut & paste" from
396 :     ## --- diff.default() from R/src/library/base/R/diff.R :
397 :     setMethod("diff", signature(x = "Matrix"),
398 :     function(x, lag = 1, differences = 1, ...) {
399 :     if (length(lag) > 1 || length(differences) > 1 ||
400 :     lag < 1 || differences < 1)
401 :     stop("'lag' and 'differences' must be integers >= 1")
402 :     xlen <- nrow(x)
403 :     if (lag * differences >= xlen)
404 :     return(x[,FALSE][0]) # empty of proper mode
405 :    
406 :     i1 <- -1:-lag
407 :     for (i in 1:differences)
408 :     x <- x[i1, , drop = FALSE] -
409 :     x[-nrow(x):-(nrow(x)-lag+1), , drop = FALSE]
410 :     x
411 :     })
412 :    
413 : maechler 1829 setMethod("image", "Matrix",
414 :     function(x, ...) { # coercing to sparse is not inefficient,
415 :     ## since we need 'i' and 'j' for levelplot()
416 : mmaechler 2482 x <- as(as(x, "sparseMatrix"), "dsparseMatrix")
417 :     ## note that "ddiMatrix" is "sparse*" and "d*", but *not* dsparse
418 : maechler 1829 callGeneric()
419 :     })
420 : maechler 1799
421 : maechler 1829
422 : maechler 956 ## Group Methods
423 :    
424 : mmaechler 2464 ## NOTE: "&" and "|" are now in group "Logic" c "Ops" --> ./Ops.R
425 :     ## "!" is in ./not.R
426 : mmaechler 2345
427 : mmaechler 2175 ## Further, see ./Ops.R
428 :     ## ~~~~~
429 :    
430 :    
431 : maechler 871 ### --------------------------------------------------------------------------
432 :     ###
433 :     ### Subsetting "[" and
434 :     ### SubAssign "[<-" : The "missing" cases can be dealt with here, "at the top":
435 : maechler 868
436 : maechler 886 ## Using "index" for indices should allow
437 : maechler 875 ## integer (numeric), logical, or character (names!) indices :
438 :    
439 : maechler 868 ## "x[]":
440 :     setMethod("[", signature(x = "Matrix",
441 :     i = "missing", j = "missing", drop = "ANY"),
442 : maechler 2056 function (x, i, j, ..., drop) x)
443 : maechler 1225
444 : maechler 868 ## missing 'drop' --> 'drop = TRUE'
445 :     ## -----------
446 : maechler 2111 ## select rows __ or __ vector indexing:
447 : maechler 886 setMethod("[", signature(x = "Matrix", i = "index", j = "missing",
448 : maechler 868 drop = "missing"),
449 : maechler 2056 function(x,i,j, ..., drop) {
450 : mmaechler 2363 Matrix.msg("M[i,m,m] : nargs()=",nargs(), .M.level = 2)
451 : maechler 2098 if(nargs() == 2) { ## e.g. M[0] , M[TRUE], M[1:2]
452 : mmaechler 2363 .M.vectorSub(x,i)
453 : maechler 2098 } else {
454 :     callGeneric(x, i=i, , drop=TRUE)
455 :     ## ^^
456 :     }
457 :     })
458 : maechler 1799
459 : maechler 868 ## select columns
460 : maechler 886 setMethod("[", signature(x = "Matrix", i = "missing", j = "index",
461 : maechler 868 drop = "missing"),
462 : mmaechler 2363 function(x,i,j, ..., drop) {
463 :     Matrix.msg("M[m,i,m] : nargs()=",nargs(), .M.level = 2)
464 :     callGeneric(x, j=j, drop= TRUE)
465 :     })
466 : maechler 886 setMethod("[", signature(x = "Matrix", i = "index", j = "index",
467 : mmaechler 2363 drop = "missing"),
468 :     function(x,i,j, ..., drop) {
469 :     Matrix.msg("M[i,i,m] : nargs()=",nargs(), .M.level = 2)
470 :     callGeneric(x, i=i, j=j, drop= TRUE)
471 :     })
472 : maechler 868
473 : maechler 875 ## bail out if any of (i,j,drop) is "non-sense"
474 :     setMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", drop = "ANY"),
475 : maechler 2056 function(x,i,j, ..., drop)
476 : mmaechler 2363 stop("invalid or not-yet-implemented 'Matrix' subsetting"))
477 : maechler 875
478 : maechler 1472 ## logical indexing, such as M[ M >= 7 ] *BUT* also M[ M[,1] >= 3,],
479 :     ## The following is *both* for M [ <logical> ]
480 :     ## and also for M [ <logical> , ]
481 : maechler 2056 .M.sub.i.logical <- function (x, i, j, ..., drop)
482 : maechler 1472 {
483 : mmaechler 2331 nA <- nargs() # counts 'M[i]' as 2 arguments, 'M[i,]' as 3
484 : maechler 1472 if(nA == 2) { ## M [ M >= 7 ]
485 : maechler 2112 ## FIXME: when both 'x' and 'i' are sparse, this can be very inefficient
486 :     if(is(x, "sparseMatrix"))
487 :     message("<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient")
488 :     toC <- geClass(x)
489 :     if(canCoerce(x, toC)) as(x, toC)@x[as.vector(i)]
490 :     else as(as(as(x, "generalMatrix"), "denseMatrix"), toC)@x[as.vector(i)]
491 : maechler 1472 ## -> error when lengths don't match
492 : mmaechler 2331 }
493 :     else if(nA == 3) { ## M[i, ] e.g., M [ M[,1, drop=FALSE] >= 7, ]
494 : maechler 1472
495 : mmaechler 2331 ## Note: current method dispatch seems not to call this ever
496 :    
497 :     if(!any(is.na(i)) && all(i)) ## select everything
498 :     x
499 :     else ## not selecting all -> result is *NOT* diagonal/triangular/symmetric/..
500 :     ## keep j missing, but drop = "logical"
501 :     callGeneric(as(x,"generalMatrix"), i = i, , drop = TRUE)
502 :    
503 : mmaechler 2387 } else stop(gettextf(
504 :     "nargs() = %d. Extraneous illegal arguments inside '[ .. ]' (i.logical)?",
505 : mmaechler 2857 nA), domain=NA)
506 : maechler 1472 }
507 : maechler 1225
508 : mmaechler 2472 ## instead of using 'drop = "ANY"' {against ambiguity notices}:
509 :     for(ii in c("lMatrix", "logical"))
510 :     setMethod("[", signature(x = "Matrix", i = ii, j = "missing", drop = "missing"),
511 :     .M.sub.i.logical)
512 : mmaechler 2898 rm(ii)
513 : maechler 1225
514 : mmaechler 2472
515 : mmaechler 2950 ##' x[ ij ] where ij is (i,j) 2-column matrix
516 :     ##' @note only called from .M.sub.i.2col(x, i) below
517 : maechler 2112 subset.ij <- function(x, ij) {
518 :     m <- nrow(ij)
519 :     if(m > 3) {
520 :     cld <- getClassDef(class(x))
521 :     sym.x <- extends(cld, "symmetricMatrix")
522 :     if(sym.x) {
523 :     W <- if(x@uplo == "U") # stored only [i,j] with i <= j
524 :     ij[,1] > ij[,2] else ij[,1] < ij[,2]
525 :     if(any(W))
526 :     ij[W,] <- ij[W, 2:1]
527 :     }
528 :     if(extends(cld, "sparseMatrix")) {
529 :     ## do something smarter:
530 : mmaechler 2203 di <- dim(x)
531 : maechler 2112 if(!extends(cld, "CsparseMatrix")) {
532 :     x <- as(x, "CsparseMatrix") # simpler; our standard
533 :     cld <- getClassDef(class(x))
534 :     }
535 :     tri.x <- extends(cld, "triangularMatrix")
536 :     if(tri.x) {
537 :     ## need these for the 'x' slot in any case
538 :     if (x@diag == "U") x <- .Call(Csparse_diagU2N, x)
539 :     ## slightly more efficient than non0.i() or non0ind():
540 :     ij.x <- .Call(compressed_non_0_ij, x, isC=TRUE)
541 : mmaechler 2323 } else { ## symmetric / general : for symmetric, only "existing" part
542 : maechler 2112 ij.x <- non0.i(x, cld)
543 :     }
544 :    
545 : mmaechler 2987 m1 <- .Call(m_encodeInd, ij.x, di, orig1=FALSE, checkBounds=FALSE)
546 :     m2 <- .Call(m_encodeInd, ij, di, orig1= TRUE, checkBounds= TRUE)
547 : mmaechler 2950 mi <- match(m2, m1, nomatch=0)
548 :     mmi <- mi != 0L ## == (m2 %in% m1)
549 :     ## Result: all FALSE or 0 apart from where we match non-zero entries
550 : maechler 2112 ans <- vector(mode = .type.kind[.M.kindC(cld)], length = m)
551 :     ## those that are *not* zero:
552 : mmaechler 2950 ans[mmi] <- if(extends(cld, "nsparseMatrix")) TRUE else x@x[mi[mmi]]
553 : mmaechler 2323 if(any(ina <- is.na(m2))) # has one or two NA in that (i,j) row
554 :     is.na(ans) <- ina
555 : maechler 2112 ans
556 :    
557 :     } else { ## non-sparse : dense
558 : maechler 2115 ##---- NEVER happens: 'denseMatrix' has its own setMethod(.) !
559 : maechler 2112 message("m[ <ij-matrix> ]: inefficiently indexing single elements")
560 :     i1 <- ij[,1]
561 :     i2 <- ij[,2]
562 : maechler 2115 ## very inefficient for large m
563 : maechler 2112 unlist(lapply(seq_len(m), function(j) x[i1[j], i2[j]]))
564 :     }
565 :     } else { # 1 <= m <= 3
566 :     i1 <- ij[,1]
567 :     i2 <- ij[,2]
568 :     unlist(lapply(seq_len(m), function(j) x[i1[j], i2[j]]))
569 :     }
570 :     }
571 :    
572 : maechler 2110 ## A[ ij ] where ij is (i,j) 2-column matrix -- but also when that is logical mat!
573 : maechler 2056 .M.sub.i.2col <- function (x, i, j, ..., drop)
574 : maechler 1673 {
575 :     nA <- nargs()
576 : mmaechler 2950 if(nA != 2)
577 :     stop(domain=NA, gettextf(
578 :     "nargs() = %d. Extraneous illegal arguments inside '[ .. ]' (i.2col)?", nA))
579 :     ## else: (nA == 2): M [ cbind(ii,jj) ] or M [ <logical matrix> ]
580 :     if(!is.integer(nc <- ncol(i)))
581 :     stop(".M.sub.i.2col(): 'i' has no integer column number;\n should never happen; please report")
582 :     if(is.logical(i))
583 :     return(.M.sub.i.logical(x, i=i)) # call with 2 args!
584 :     else if(!is.numeric(i) || nc != 2)
585 :     stop("such indexing must be by logical or 2-column numeric matrix")
586 :     if(!nrow(i)) return(vector(mode = .type.kind[.M.kind(x)]))
587 :     ## else
588 :     subset.ij(x, i)
589 : maechler 868
590 : maechler 1673 }
591 :     setMethod("[", signature(x = "Matrix", i = "matrix", j = "missing"),# drop="ANY"
592 :     .M.sub.i.2col)
593 : mmaechler 2472 ## just against ambiguity notices :
594 :     setMethod("[", signature(x = "Matrix", i = "matrix", j = "missing", drop="missing"),
595 :     .M.sub.i.2col)
596 : maechler 868
597 : maechler 1472
598 : maechler 871 ### "[<-" : -----------------
599 : maechler 868
600 : maechler 1673 ## A[ ij ] <- value, where ij is (i,j) 2-column matrix :
601 : maechler 2056 ## ----------------
602 : mmaechler 2950 ## The cheap general method, now only used for "pMatrix","indMatrix"
603 :     ## sparse all use .TM.repl.i.mat()
604 : maechler 2056 ## NOTE: need '...' below such that setMethod() does
605 :     ## not use .local() such that nargs() will work correctly:
606 :     .M.repl.i.2col <- function (x, i, j, ..., value)
607 : maechler 1673 {
608 :     nA <- nargs()
609 : maechler 2110 if(nA == 3) { ## M [ cbind(ii,jj) ] <- value or M [ Lmat ] <- value
610 : maechler 1673 if(!is.integer(nc <- ncol(i)))
611 : mmaechler 2863 stop(".M.repl.i.2col(): 'i' has no integer column number;\n should never happen; please report")
612 : maechler 1724 else if(!is.numeric(i) || nc != 2)
613 :     stop("such indexing must be by logical or 2-column numeric matrix")
614 : maechler 1673 if(is.logical(i)) {
615 : maechler 1738 message(".M.repl.i.2col(): drop 'matrix' case ...")
616 : maechler 2110 ## c(i) : drop "matrix" to logical vector
617 :     return( callGeneric(x, i=c(i), value=value) )
618 :     }
619 : maechler 1724 if(!is.integer(i)) storage.mode(i) <- "integer"
620 :     if(any(i < 0))
621 :     stop("negative values are not allowed in a matrix subscript")
622 :     if(any(is.na(i)))
623 :     stop("NAs are not allowed in subscripted assignments")
624 :     if(any(i0 <- (i == 0))) # remove them
625 :     i <- i[ - which(i0, arr.ind = TRUE)[,"row"], ]
626 :     ## now have integer i >= 1
627 : maechler 1673 m <- nrow(i)
628 : maechler 1724 ## mod.x <- .type.kind[.M.kind(x)]
629 : maechler 1673 if(length(value) > 0 && m %% length(value) != 0)
630 :     warning("number of items to replace is not a multiple of replacement length")
631 :     ## recycle:
632 : mmaechler 2945 value <- rep_len(value, m)
633 : maechler 1673 i1 <- i[,1]
634 :     i2 <- i[,2]
635 : maechler 2112 if(m > 2)
636 :     message("m[ <ij-matrix> ] <- v: inefficiently treating single elements")
637 : maechler 1673 ## inefficient -- FIXME -- (also loses "symmetry" unnecessarily)
638 :     for(k in seq_len(m))
639 :     x[i1[k], i2[k]] <- value[k]
640 : maechler 1724
641 : maechler 1673 x
642 : mmaechler 2387 } else stop(gettextf(
643 :     "nargs() = %d. Extraneous illegal arguments inside '[ .. ]' ?",
644 : mmaechler 2857 nA), domain=NA)
645 : maechler 1673 }
646 : maechler 1724
647 : maechler 1673 setReplaceMethod("[", signature(x = "Matrix", i = "matrix", j = "missing",
648 :     value = "replValue"),
649 : mmaechler 2268 .M.repl.i.2col)
650 : maechler 1673
651 : mmaechler 2207 ## Three catch-all methods ... would be very inefficient for sparse*
652 :     ## --> extra methods in ./sparseMatrix.R
653 : maechler 2062 setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "ANY",
654 : maechler 1615 value = "Matrix"),
655 : maechler 2062 function (x, i, j, ..., value)
656 : maechler 2098 callGeneric(x=x, , j=j, value = as.vector(value)))
657 : maechler 1705
658 : maechler 2062 setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "missing",
659 :     value = "Matrix"),
660 :     function (x, i, j, ..., value)
661 : maechler 2098 callGeneric(x=x, i=i, , value = as.vector(value)))
662 : maechler 2056
663 : maechler 1705 setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY",
664 :     value = "Matrix"),
665 : maechler 2062 function (x, i, j, ..., value)
666 : maechler 1615 callGeneric(x=x, i=i, j=j, value = as.vector(value)))
667 : maechler 1705
668 : mmaechler 2207
669 : maechler 2062 setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "ANY",
670 :     value = "matrix"),
671 :     function (x, i, j, ..., value)
672 : maechler 2098 callGeneric(x=x, , j=j, value = c(value)))
673 : maechler 2062
674 :     setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "missing",
675 :     value = "matrix"),
676 :     function (x, i, j, ..., value)
677 : maechler 2098 callGeneric(x=x, i=i, , value = c(value)))
678 : maechler 2062
679 : maechler 1615 setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY",
680 :     value = "matrix"),
681 :     function (x, i, j, value)
682 :     callGeneric(x=x, i=i, j=j, value = c(value)))
683 :    
684 : mmaechler 2549 ## M [ <lMatrix> ] <- value; used notably for x = "CsparseMatrix" -------------------
685 :     .repl.i.lDMat <- function (x, i, j, ..., value)
686 :     {
687 :     ## nA <- nargs()
688 : mmaechler 2857 ## if(nA != 3) stop(gettextf("nargs() = %d should never happen; please report.", nA), domain=NA)
689 : mmaechler 2549 ## else: nA == 3 i.e., M [ Lmat ] <- value
690 :     ## x[i] <- value ; return(x)
691 :     `[<-`(x, i=which(as.vector(i)), value=value)
692 :     }
693 :     setReplaceMethod("[", signature(x = "Matrix", i = "ldenseMatrix", j = "missing",
694 :     value = "replValue"), .repl.i.lDMat)
695 :     setReplaceMethod("[", signature(x = "Matrix", i = "ndenseMatrix", j = "missing",
696 :     value = "replValue"), .repl.i.lDMat)
697 :     .repl.i.lSMat <- function (x, i, j, ..., value)
698 :     {
699 :     ## nA <- nargs()
700 : mmaechler 2857 ## if(nA != 3) stop(gettextf("nargs() = %d should never happen; please report.", nA), domain=NA)
701 : mmaechler 2549 ## else: nA == 3 i.e., M [ Lmat ] <- value
702 :     ## x[i] <- value ; return(x)
703 :     `[<-`(x, i=which(as(i, "sparseVector")), value=value)
704 :     }
705 :     setReplaceMethod("[", signature(x = "Matrix", i = "lsparseMatrix", j = "missing",
706 :     value = "replValue"), .repl.i.lSMat)
707 :     setReplaceMethod("[", signature(x = "Matrix", i = "nsparseMatrix", j = "missing",
708 :     value = "replValue"), .repl.i.lSMat)
709 :    
710 : maechler 1189 ## (ANY,ANY,ANY) is used when no `real method' is implemented :
711 : maechler 871 setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY",
712 :     value = "ANY"),
713 : maechler 1189 function (x, i, j, value) {
714 :     if(!is.atomic(value))
715 : mmaechler 2387 stop(gettextf(
716 :     "RHS 'value' (class %s) matches 'ANY', but must match matrix class %s",
717 : mmaechler 2857 class(value), class(x)), domain=NA)
718 : maechler 1189 else stop("not-yet-implemented 'Matrix[<-' method")
719 :     })

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge