1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
|
#### Index Matrices -- Coercion and Methods (--> ../man/indMatrix-class.Rd )
## The typical 'constructor' : coerce from 'perm'
setAs("integer", "indMatrix",
function(from)
new("indMatrix", Dim = c(length(from), max(from)),
Dimnames = list(names(from), NULL),
perm = from))
setAs("numeric", "indMatrix",
function(from)
if(all(from == (i <- as.integer(from)))) as(i, "indMatrix")
else stop("coercion to \"indMatrix\" only works from integer numeric"))
## A constructor from a list giving the index ('perm') and the number of columns
## (need this for cases in which the value(s) represented by the last
## column(s) has no observations):
.list2indMatrix <- function(from) {
if(length(from) == 2 &&
all(from[[1]] == (i <- as.integer(from[[1]]))) &&
from[[2]] == (d <- as.integer(from[[2]])) &&
length(d) == 1 && d >= max(i)) {
new("indMatrix", perm = i, Dim = c(length(i), d))
} else
stop("coercion from list(i1,...,ik, d) to \"indMatrix\" failed.
All entries must be integer valued and the number of columns, d, not smaller
than the maximal index i*.")
}
setAs("list", "indMatrix", .list2indMatrix)
setAs("indMatrix", "matrix",
function(from) {
fp <- from@perm
r <- ldiag(n = from@Dim[2])[fp,]
if(.has.DN(from)) dimnames(r) <- from@Dimnames
r
})
## coerce to 0/1 sparse matrix, i.e. sparse pattern
.ind2ngT <- function(from) {
d <- from@Dim
new("ngTMatrix", i = seq_len(d[1]) - 1L, j = from@perm - 1L,
Dim = d, Dimnames = from@Dimnames)
}
setAs("indMatrix", "ngTMatrix", .ind2ngT)
setAs("indMatrix", "TsparseMatrix", .ind2ngT)
setAs("indMatrix", "nMatrix", .ind2ngT)
setAs("indMatrix", "lMatrix", function(from) as(.ind2ngT(from), "lMatrix"))
setAs("indMatrix", "dMatrix", function(from) as(.ind2ngT(from), "dMatrix"))
setAs("indMatrix", "dsparseMatrix", function(from) as(from, "dMatrix"))
setAs("indMatrix", "lsparseMatrix", function(from) as(from, "lMatrix"))
setAs("indMatrix", "nsparseMatrix", .ind2ngT)
setAs("indMatrix", "CsparseMatrix",
function(from) as(.ind2ngT(from), "CsparseMatrix"))
setAs("indMatrix", "ngeMatrix", function(from) as(.ind2ngT(from),"ngeMatrix"))
setAs("nMatrix", "indMatrix",
function(from) {
from <- as(as(from, "TsparseMatrix"), "ngTMatrix")
n <- (d <- from@Dim)[1]
if(length(i <- from@i) != n)
stop("the number of non-zero entries differs from nrow(.)")
if((need.sort <- is.unsorted(i))) {
ii <- sort.list(i)
i <- i[ii]
}
if(n >= 1 && !identical(i, 0:(n - 1)))
stop("must have exactly one non-zero entry per row")
new("indMatrix", ## validity checking checks the 'perm' slot:
perm = 1L + if(need.sort) from@j[ii] else from@j,
Dim = d, Dimnames = from@Dimnames)
})
setAs("matrix", "indMatrix", function(from) as(as(from, "nMatrix"), "indMatrix"))
setAs("indMatrix", "matrix", function(from) as(.ind2ngT(from), "matrix"))
setAs("sparseMatrix", "indMatrix", function(from)
as(as(from, "nsparseMatrix"), "indMatrix"))
setMethod("is.na", signature(x = "indMatrix"), is.na_nsp)
setMethod("is.infinite", signature(x = "indMatrix"), is.na_nsp)
setMethod("is.finite", signature(x = "indMatrix"), allTrueMatrix)
setMethod("t", signature(x = "indMatrix"), function(x) t(.ind2ngT(x)))
setMethod("isSymmetric", signature(object = "indMatrix"),
function(object, ...) {
d <- dim(object)
if(d[1L] != d[2L])
FALSE
else ## using "!=" (instead of "==") as the former is typically sparse
!any(object != t(object))
})
setMethod("%*%", signature(x = "matrix", y = "indMatrix"),
function(x, y) x %*% as(y, "lMatrix"))
setMethod("%*%", signature(x = "Matrix", y = "indMatrix"),
function(x, y) x %*% as(y, "lMatrix"))
setMethod("%*%", signature(x = "indMatrix", y = "matrix"),
function(x, y) { mmultCheck(x,y); y[x@perm ,] })
setMethod("%*%", signature(x = "indMatrix", y = "Matrix"),
function(x, y) { mmultCheck(x,y); y[x@perm ,] })
setMethod("crossprod", signature(x = "indMatrix", y = "matrix"),
function(x, y) as(t(x), "lMatrix") %*% y)
setMethod("crossprod", signature(x = "indMatrix", y = "Matrix"),
function(x, y) as(t(x), "lMatrix") %*% y)
setMethod("crossprod", signature(x = "indMatrix", y = "indMatrix"),
function(x, y) {
mmultCheck(x,y, 2L)
## xy <- interaction(x@perm, y@perm)
## this is wrong if any of the columns in X or Y are empty because interaction()
## drops non-occuring levels from a non-factor. Explicitly defining a factor with
## levels 1:ncol(<indMatrix>) avoids that.
nx <- x@Dim[2L]
ny <- y@Dim[2L]
## xy <- interaction(factor(x@perm, levels=seq_len(nx)),
## factor(y@perm, levels=seq_len(ny)))
## much faster (notably for large x,y):
xy <- x@perm + nx*as.double(y@perm-1L)
Matrix(tabulate(xy, nbins = nx*ny), nrow = nx, ncol = ny,
dimnames = list(x@Dimnames[[2L]], y@Dimnames[[2L]]))
})
setMethod("tcrossprod", signature(x = "matrix", y = "indMatrix"),
function(x, y) { mmultCheck(x,y, 3L); x[, y@perm] })
setMethod("tcrossprod", signature(x = "Matrix", y = "indMatrix"),
function(x, y) { mmultCheck(x,y, 3L); x[, y@perm] })
setMethod("tcrossprod", signature(x = "indMatrix", y = "indMatrix"),
function(x, y) { mmultCheck(x,y, 3L); x[, y@perm] })
setMethod("crossprod", signature(x = "indMatrix", y = "missing"),
function(x, y=NULL) Diagonal(x = tabulate(x@perm, nbins=x@Dim[2L])))
setMethod("tcrossprod", signature(x = "indMatrix", y = "missing"),
function(x, y=NULL) x[,x@perm])
setMethod("kronecker", signature(X = "indMatrix", Y = "indMatrix"),
function (X, Y, FUN = "*", make.dimnames = FALSE, ...) {
if (FUN != "*") stop("kronecker method must use default 'FUN'")
if(any(as.double(X@Dim)*Y@Dim >= .Machine$integer.max))
stop("resulting matrix dimension would be too large")
## Explicitly defining a factor with levels 1:ncol(.) avoids that
## interaction() drops non-occuring levels when any of the
## columns in X or Y are empty:
## perm <- as.integer(interaction(factor(rep(X@perm, each =Y@Dim[1]),
## levels=seq_len(X@Dim[2])),
## factor(rep.int(Y@perm, times=X@Dim[1]),
## levels=seq_len(Y@Dim[2])),
## lex.order=TRUE))
## much faster (notably for large X, Y):
fX <- rep (X@perm-1L, each = Y@Dim[1])
fY <- rep.int(Y@perm-1L, times = X@Dim[1])
new("indMatrix", perm = 1L + fY + Y@Dim[2] * fX,
Dim = X@Dim*Y@Dim)
})
setMethod("[", signature(x = "indMatrix", i = "index", j = "missing",
drop = "logical"),
function (x, i, j, ..., drop)
{
n <- length(newperm <- x@perm[i])
if(drop && n == 1) { ## -> logical unit vector
newperm == seq_len(x@Dim[2])
} else { ## stay matrix
if(!is.null((DN <- x@Dimnames)[[1]])) DN[[1]] <- DN[[1]][i]
new("indMatrix", perm = newperm,
Dim = c(n, x@Dim[2]), Dimnames = DN)
}
})
.indMat.nosense <- function (x, i, j, ..., value)
stop('replacing "indMatrix" entries is not allowed, as rarely sensible')
setReplaceMethod("[", signature(x = "indMatrix", i = "index"), .indMat.nosense)
setReplaceMethod("[", signature(x = "indMatrix", i = "missing", j = "index"),
.indMat.nosense) ## explicit ^^^^^^^^^^^^ for disambiguation
setReplaceMethod("[", signature(x = "indMatrix", i = "missing", j = "missing"),
.indMat.nosense)
### rbind2() method: --> bind2.R
|
