1 |
#### Toplevel ``virtual'' class "Matrix" |
#### Toplevel ``virtual'' class "Matrix" |
2 |
|
|
3 |
|
|
4 |
|
### Virtual coercions -- via smart "helpers" (-> ./Auxiliaries.R) |
5 |
|
|
6 |
|
setAs("Matrix", "sparseMatrix", function(from) as_Csparse(from)) |
7 |
|
setAs("Matrix", "denseMatrix", function(from) as_dense(from)) |
8 |
|
|
9 |
## ## probably not needed eventually: |
## ## probably not needed eventually: |
10 |
## setAs(from = "ddenseMatrix", to = "matrix", |
## setAs(from = "ddenseMatrix", to = "matrix", |
11 |
## function(from) { |
## function(from) { |
18 |
## for 'Matrix' objects, as.array() should be equivalent: |
## for 'Matrix' objects, as.array() should be equivalent: |
19 |
setMethod("as.array", signature(x = "Matrix"), function(x) as(x, "matrix")) |
setMethod("as.array", signature(x = "Matrix"), function(x) as(x, "matrix")) |
20 |
|
|
21 |
|
## head and tail apply to all Matrix objects for which subscripting is allowed: |
22 |
|
## if(paste(R.version$major, R.version$minor, sep=".") < "2.4") { |
23 |
|
setMethod("head", signature(x = "Matrix"), utils:::head.matrix) |
24 |
|
setMethod("tail", signature(x = "Matrix"), utils:::tail.matrix) |
25 |
|
## } else { # R 2.4.0 and newer |
26 |
|
## setMethod("head", signature(x = "Matrix"), utils::head.matrix) |
27 |
|
## setMethod("tail", signature(x = "Matrix"), utils::tail.matrix) |
28 |
|
## } |
29 |
|
|
30 |
## slow "fall back" method {subclasses should have faster ones}: |
## slow "fall back" method {subclasses should have faster ones}: |
31 |
setMethod("as.vector", signature(x = "Matrix", mode = "missing"), |
setMethod("as.vector", signature(x = "Matrix", mode = "missing"), |
32 |
function(x) as.vector(as(x, "matrix"))) |
function(x) as.vector(as(x, "matrix"))) |
33 |
|
|
34 |
|
## mainly need these for "dMatrix" or "lMatrix" respectively, but why not general: |
35 |
|
setMethod("as.numeric", signature(x = "Matrix"), |
36 |
|
function(x, ...) as.numeric(as.vector(x))) |
37 |
|
setMethod("as.logical", signature(x = "Matrix"), |
38 |
|
function(x, ...) as.logical(as.vector(x))) |
39 |
|
|
40 |
|
|
41 |
## Note that isSymmetric is *not* exported --- |
## Note that isSymmetric is *not* exported |
42 |
### but also note that "base" eigen may get an isSymmetric() that *would* be exported! |
## but that "base" has an isSymmetric() S3-generic since R 2.3.0 |
43 |
setMethod("isSymmetric", signature(object = "symmetricMatrix", tol="ANY"), |
setMethod("isSymmetric", signature(object = "symmetricMatrix"), |
44 |
function(object,tol) TRUE) |
function(object,tol) TRUE) |
45 |
setMethod("isSymmetric", signature(object = "triangularMatrix", tol="ANY"), |
setMethod("isSymmetric", signature(object = "triangularMatrix"), |
46 |
## FIXME: 'TRUE' if *diagonal*, i.e. return(isDiagonal(object)) |
## TRUE iff diagonal: |
47 |
function(object,tol) FALSE) |
function(object,tol) isDiagonal(object)) |
48 |
|
|
49 |
setMethod("isDiagonal", signature(object = "sparseMatrix"), |
if(paste(R.version$major, R.version$minor, sep=".") < "2.3") |
50 |
function(object) { |
## need a "matrix" method as in R 2.3 and later |
51 |
gT <- as(object, "TsparseMatrix") |
setMethod("isSymmetric", signature(object = "matrix"), |
52 |
all(gT@i == gT@j) |
function(object, tol = 100*.Machine$double.eps, ...) |
53 |
|
{ |
54 |
|
## pretest: is it square? |
55 |
|
d <- dim(object) |
56 |
|
if(d[1] != d[2]) return(FALSE) |
57 |
|
## for `broken' all.equal in R <= 2.2.x: |
58 |
|
dn <- dimnames(object) |
59 |
|
if(!identical(dn[1], dn[2])) return(FALSE) |
60 |
|
test <- |
61 |
|
if(is.complex(object)) |
62 |
|
all.equal.numeric(object, Conj(t(object)), tol = tol, ...) |
63 |
|
else # numeric, character, .. |
64 |
|
all.equal(object, t(object), tol = tol, ...) |
65 |
|
isTRUE(test) |
66 |
}) |
}) |
67 |
|
|
68 |
|
|
69 |
|
setMethod("isTriangular", signature(object = "triangularMatrix"), |
70 |
|
function(object, ...) TRUE) |
71 |
|
|
72 |
|
setMethod("isTriangular", signature(object = "matrix"), isTriMat) |
73 |
|
|
74 |
|
setMethod("isDiagonal", signature(object = "matrix"), .is.diagonal) |
75 |
|
|
76 |
|
|
77 |
|
|
78 |
setMethod("dim", signature(x = "Matrix"), |
setMethod("dim", signature(x = "Matrix"), |
79 |
function(x) x@Dim, valueClass = "integer") |
function(x) x@Dim, valueClass = "integer") |
80 |
setMethod("dimnames", signature(x = "Matrix"), function(x) x@Dimnames) |
setMethod("dimnames", signature(x = "Matrix"), function(x) x@Dimnames) |
98 |
|
|
99 |
Matrix <- |
Matrix <- |
100 |
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL, |
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL, |
101 |
sparse = NULL) |
sparse = NULL, forceCheck = FALSE) |
102 |
{ |
{ |
103 |
sparseDefault <- function(m) |
sparseDefault <- function(m) |
104 |
prod(dim(m)) > 2*sum(as(m, "matrix") != 0) |
prod(dim(m)) > 2*sum(as(m, "matrix") != 0) |
107 |
if(is.null(sparse) && (i.M || is(data, "matrix"))) |
if(is.null(sparse) && (i.M || is(data, "matrix"))) |
108 |
sparse <- sparseDefault(data) |
sparse <- sparseDefault(data) |
109 |
|
|
110 |
if (i.M) { |
doDN <- TRUE |
111 |
|
if (i.M && !forceCheck) { |
112 |
sM <- is(data,"sparseMatrix") |
sM <- is(data,"sparseMatrix") |
113 |
if((sparse && sM) || (!sparse && !sM)) |
if((sparse && sM) || (!sparse && !sM)) |
114 |
return(data) |
return(data) |
119 |
nrow <- ceiling(length(data)/ncol) |
nrow <- ceiling(length(data)/ncol) |
120 |
else if (missing(ncol)) |
else if (missing(ncol)) |
121 |
ncol <- ceiling(length(data)/nrow) |
ncol <- ceiling(length(data)/nrow) |
122 |
|
if(length(data) == 1 && data == 0 && !identical(sparse,FALSE)) { |
123 |
|
if(is.null(sparse)) sparse <- TRUE |
124 |
|
## will be sparse: do NOT construct full matrix! |
125 |
|
data <- new(if(is.numeric(data)) "dgTMatrix" else |
126 |
|
if(is.logical(data)) "lgTMatrix" else |
127 |
|
stop("invalid 'data'"), |
128 |
|
Dim = as.integer(c(nrow,ncol)), |
129 |
|
Dimnames = if(is.null(dimnames)) list(NULL,NULL) |
130 |
|
else dimnames) |
131 |
|
} else { ## normal case |
132 |
data <- .Internal(matrix(data, nrow, ncol, byrow)) |
data <- .Internal(matrix(data, nrow, ncol, byrow)) |
133 |
if(is.null(sparse)) |
if(is.null(sparse)) |
134 |
sparse <- sparseDefault(data) |
sparse <- sparseDefault(data) |
135 |
dimnames(data) <- dimnames |
dimnames(data) <- dimnames |
136 |
} |
} |
137 |
|
doDN <- FALSE |
138 |
|
} |
139 |
## 'data' is now a "matrix" or "Matrix" |
## 'data' is now a "matrix" or "Matrix" |
140 |
## FIXME: consider it's type (logical,....) |
if (doDN && !is.null(dimnames)) |
141 |
## ctype <- substr(class(data), 1,1) # "d", "l", ... |
dimnames(data) <- dimnames |
142 |
## FIXME(2): check for symmetric / triangular / ... |
|
143 |
|
## check for symmetric / triangular / diagonal : |
144 |
|
isSym <- isSymmetric(data) |
145 |
|
if((isTri <- !isSym)) |
146 |
|
isTri <- isTriangular(data) |
147 |
|
isDiag <- isSym # cannot be diagonal if it isn't symmetric |
148 |
|
if(isDiag) |
149 |
|
isDiag <- isDiagonal(data) |
150 |
|
|
151 |
### TODO: Compare with as.Matrix() and its tests in ./dgeMatrix.R |
### TODO: Compare with as.Matrix() and its tests in ./dgeMatrix.R |
152 |
if(sparse) |
|
153 |
as(data, "dgCMatrix") |
## Find proper matrix class 'cl' |
154 |
else |
cl <- |
155 |
as(data, "dgeMatrix") |
if(isDiag) |
156 |
|
"diagonalMatrix" # -> will automatically check for type |
157 |
|
else { |
158 |
|
## consider it's type |
159 |
|
ctype <- |
160 |
|
if(is(data,"Matrix")) class(data) |
161 |
|
else { |
162 |
|
if("complex" == (ctype <- typeof(data))) |
163 |
|
"z" else ctype |
164 |
|
} |
165 |
|
ctype <- substr(ctype, 1,1) # "d", "l", "i" or "z" |
166 |
|
if(ctype == "z") |
167 |
|
stop("complex matrices not yet implemented in Matrix package") |
168 |
|
if(ctype == "i") { |
169 |
|
warning("integer matrices not yet implemented in 'Matrix'; ", |
170 |
|
"using 'double' ones'") |
171 |
|
ctype <- "d" |
172 |
|
} |
173 |
|
paste(ctype, |
174 |
|
if(sparse) { |
175 |
|
if(isSym) "sCMatrix" else |
176 |
|
if(isTri) "tCMatrix" else "gCMatrix" |
177 |
|
} else { ## dense |
178 |
|
if(isSym) "syMatrix" else |
179 |
|
if(isTri) "trMatrix" else "geMatrix" |
180 |
|
}, sep="") |
181 |
|
} |
182 |
|
|
183 |
|
## Now coerce and return |
184 |
|
as(data, cl) |
185 |
} |
} |
186 |
|
|
187 |
## Methods for operations where one argument is numeric |
## Methods for operations where one argument is numeric |
197 |
|
|
198 |
setMethod("crossprod", signature(x = "Matrix", y = "numeric"), |
setMethod("crossprod", signature(x = "Matrix", y = "numeric"), |
199 |
function(x, y = NULL) callGeneric(x, as.matrix(y))) |
function(x, y = NULL) callGeneric(x, as.matrix(y))) |
|
|
|
200 |
setMethod("crossprod", signature(x = "numeric", y = "Matrix"), |
setMethod("crossprod", signature(x = "numeric", y = "Matrix"), |
201 |
function(x, y = NULL) callGeneric(rbind(x), y)) |
function(x, y = NULL) callGeneric(as.matrix(x), y)) |
202 |
|
|
203 |
|
## The as.matrix() promotion seems illogical to MM, |
204 |
|
## but is according to help(tcrossprod, package = "base") : |
205 |
|
setMethod("tcrossprod", signature(x = "Matrix", y = "numeric"), |
206 |
|
function(x, y = NULL) callGeneric(x, as.matrix(y))) |
207 |
|
setMethod("tcrossprod", signature(x = "numeric", y = "Matrix"), |
208 |
|
function(x, y = NULL) callGeneric(as.matrix(x), y)) |
209 |
|
|
210 |
setMethod("solve", signature(a = "Matrix", b = "numeric"), |
setMethod("solve", signature(a = "Matrix", b = "numeric"), |
211 |
function(a, b, ...) callGeneric(a, as.matrix(b))) |
function(a, b, ...) callGeneric(a, as.matrix(b))) |
220 |
function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y))) |
function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y))) |
221 |
setMethod("crossprod", signature(x = "ANY", y = "Matrix"), |
setMethod("crossprod", signature(x = "ANY", y = "Matrix"), |
222 |
function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y))) |
function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y))) |
223 |
|
setMethod("tcrossprod", signature(x = "Matrix", y = "ANY"), |
224 |
|
function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y))) |
225 |
|
setMethod("tcrossprod", signature(x = "ANY", y = "Matrix"), |
226 |
|
function (x, y = NULL) .bail.out.2(.Generic, class(x), class(y))) |
227 |
|
|
228 |
|
## cheap fallbacks |
229 |
|
setMethod("crossprod", signature(x = "Matrix", y = "Matrix"), |
230 |
|
function(x, y = NULL) t(x) %*% y) |
231 |
|
setMethod("tcrossprod", signature(x = "Matrix", y = "Matrix"), |
232 |
|
function(x, y = NULL) x %*% t(y)) |
233 |
|
|
234 |
## There are special sparse methods; this is a "fall back": |
## There are special sparse methods; this is a "fall back": |
235 |
setMethod("kronecker", signature(X = "Matrix", Y = "ANY", |
setMethod("kronecker", signature(X = "Matrix", Y = "ANY", |
242 |
Y <- as(Y, "matrix") ; Matrix(callGeneric()) }) |
Y <- as(Y, "matrix") ; Matrix(callGeneric()) }) |
243 |
|
|
244 |
|
|
245 |
|
setMethod("diag", signature(x = "Matrix"), |
246 |
|
function(x, nrow, ncol) .bail.out.1(.Generic, class(x))) |
247 |
setMethod("t", signature(x = "Matrix"), |
setMethod("t", signature(x = "Matrix"), |
248 |
function(x) .bail.out.1(.Generic, class(x))) |
function(x) .bail.out.1(.Generic, class(x))) |
249 |
|
|
281 |
setMethod("[", signature(x = "Matrix", |
setMethod("[", signature(x = "Matrix", |
282 |
i = "missing", j = "missing", drop = "ANY"), |
i = "missing", j = "missing", drop = "ANY"), |
283 |
function (x, i, j, drop) x) |
function (x, i, j, drop) x) |
284 |
|
|
285 |
## missing 'drop' --> 'drop = TRUE' |
## missing 'drop' --> 'drop = TRUE' |
286 |
## ----------- |
## ----------- |
287 |
## select rows |
## select rows |
301 |
function(x,i,j, drop) |
function(x,i,j, drop) |
302 |
stop("invalid or not-yet-implemented 'Matrix' subsetting")) |
stop("invalid or not-yet-implemented 'Matrix' subsetting")) |
303 |
|
|
304 |
|
## "logical *vector* indexing, such as M [ M >= 10 ] : |
305 |
|
setMethod("[", signature(x = "Matrix", i = "lMatrix", j = "missing", |
306 |
|
drop = "ANY"), |
307 |
|
function (x, i, j, drop) { |
308 |
|
as(x, geClass(x))@x[as.vector(i)] |
309 |
|
## -> error when lengths don't match |
310 |
|
}) |
311 |
|
|
312 |
|
## FIXME: The following is good for M [ <logical> ] |
313 |
|
## *BUT* it also triggers for M [ <logical> , ] where it is *WRONG* |
314 |
|
## using nargs() does not help: it gives '3' for both cases |
315 |
|
if(FALSE) |
316 |
|
setMethod("[", signature(x = "Matrix", i = "logical", j = "missing", |
317 |
|
drop = "ANY"), |
318 |
|
function (x, i, j, drop) { |
319 |
|
## DEBUG |
320 |
|
cat("[(Matrix,i,..): nargs=", nargs(),"\n") |
321 |
|
as(x, geClass(x))@x[i] }) |
322 |
|
|
323 |
|
|
324 |
## "FIXME:" |
## "FIXME:" |
325 |
## How can we get at A[ ij ] where ij is (i,j) 2-column matrix? |
## How can we get at A[ ij ] where ij is (i,j) 2-column matrix? |
326 |
## and A[ LL ] where LL is a logical *vector* |
## and A[ LL ] where LL is a logical *vector* |
331 |
|
|
332 |
## x[] <- value : |
## x[] <- value : |
333 |
setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "missing", |
setReplaceMethod("[", signature(x = "Matrix", i = "missing", j = "missing", |
334 |
value = "index"),## double/logical/... |
value = "ANY"),## double/logical/... |
335 |
function (x, value) { x@x <- value ; validObject(x); x }) |
function (x, value) { |
336 |
|
x@x <- value |
337 |
|
validObject(x)# check if type and lengths above match |
338 |
|
x |
339 |
|
}) |
340 |
|
|
341 |
## Otherwise (value is not "index"): bail out |
## Method for all 'Matrix' kinds (rather than incomprehensible error messages); |
342 |
|
## (ANY,ANY,ANY) is used when no `real method' is implemented : |
343 |
setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", |
setReplaceMethod("[", signature(x = "Matrix", i = "ANY", j = "ANY", |
344 |
value = "ANY"), |
value = "ANY"), |
345 |
function (x, i, j, value) |
function (x, i, j, value) { |
346 |
if(!is(value,"index")) |
if(!is.atomic(value)) |
347 |
stop("RHS 'value' must be of class \"index\"") |
stop("RHS 'value' must match matrix class ", class(x)) |
348 |
else stop("not-yet-implemented 'Matrix[<-' method")) |
else stop("not-yet-implemented 'Matrix[<-' method") |
349 |
|
}) |
|
|
|
|
|
|
|
## NOTE: the following only works for R 2.2.x (and later) --- |
|
|
## ---- *and* 'Matrix' must have been *installed* by R >= 2.2.x |
|
350 |
|
|
|
if(paste(R.version$major, R.version$minor, sep=".") >= "2.2") { |
|
351 |
|
|
352 |
## The trivial methods : |
## The trivial methods : |
353 |
setMethod("cbind2", signature(x = "Matrix", y = "NULL"), |
setMethod("cbind2", signature(x = "Matrix", y = "NULL"), |
379 |
colCheck(x,y) |
colCheck(x,y) |
380 |
t(cbind2(t(x), t(y))) |
t(cbind2(t(x), t(y))) |
381 |
}) |
}) |
|
|
|
|
}## R-2.2.x and newer |
|