SCM

SCM Repository

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

Annotation of /pkg/R/Auxiliaries.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1227 - (view) (download)

1 : maechler 632 #### "Namespace private" Auxiliaries such as method functions
2 :     #### (called from more than one place --> need to be defined early)
3 :    
4 : maechler 656 ## For %*% (M = Matrix; v = vector (double or integer {complex maybe?}):
5 :     .M.v <- function(x, y) callGeneric(x, as.matrix(y))
6 :     .v.M <- function(x, y) callGeneric(rbind(x), y)
7 : maechler 632
8 : maechler 656 .has.DN <- ## has non-trivial Dimnames slot?
9 :     function(x) !identical(list(NULL,NULL), x@Dimnames)
10 :    
11 : maechler 949 .bail.out.1 <- function(fun, cl) {
12 :     stop(gettextf('not-yet-implemented method for %s(<%s>)', fun, cl),
13 :     call. = FALSE)
14 :     }
15 :     .bail.out.2 <- function(fun, cl1, cl2) {
16 :     stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)',
17 :     fun, cl1, cl2), call. = FALSE)
18 :     }
19 :    
20 : maechler 632 ## chol() via "dpoMatrix"
21 :     cholMat <- function(x, pivot, LINPACK) {
22 :     px <- as(x, "dpoMatrix")
23 : bates 703 if (isTRUE(validObject(px, test=TRUE))) chol(px)
24 : maechler 632 else stop("'x' is not positive definite -- chol() undefined.")
25 :     }
26 : maechler 908
27 : maechler 954 dimCheck <- function(a, b) {
28 :     da <- dim(a)
29 :     db <- dim(b)
30 :     if(any(da != db))
31 :     stop(gettextf("Matrices must have same dimensions in %s",
32 :     deparse(sys.call(sys.parent()))),
33 :     call. = FALSE)
34 :     da
35 :     }
36 :    
37 : maechler 956 dimNamesCheck <- function(a, b) {
38 :     ## assume dimCheck() has happened before
39 :     nullDN <- list(NULL,NULL)
40 :     h.a <- !identical(nullDN, dna <- dimnames(a))
41 :     h.b <- !identical(nullDN, dnb <- dimnames(b))
42 :     if(h.a || h.b) {
43 : maechler 1084 if (!h.b) dna
44 :     else if(!h.a) dnb
45 : maechler 956 else { ## both have non-trivial dimnames
46 :     r <- dna # "default" result
47 :     for(j in 1:2) {
48 :     dn <- dnb[[j]]
49 :     if(is.null(r[[j]]))
50 :     r[[j]] <- dn
51 :     else if (!is.null(dn) && any(r[[j]] != dn))
52 :     warning(gettextf("dimnames [%d] mismatch in %s", j,
53 :     deparse(sys.call(sys.parent()))),
54 :     call. = FALSE)
55 :     }
56 :     r
57 :     }
58 :     }
59 :     else
60 :     nullDN
61 :     }
62 :    
63 : maechler 908 rowCheck <- function(a, b) {
64 :     da <- dim(a)
65 :     db <- dim(b)
66 :     if(da[1] != db[1])
67 :     stop(gettextf("Matrices must have same number of rows in %s",
68 :     deparse(sys.call(sys.parent()))),
69 :     call. = FALSE)
70 :     ## return the common nrow()
71 :     da[1]
72 :     }
73 :    
74 :     colCheck <- function(a, b) {
75 :     da <- dim(a)
76 :     db <- dim(b)
77 :     if(da[2] != db[2])
78 :     stop(gettextf("Matrices must have same number of columns in %s",
79 :     deparse(sys.call(sys.parent()))),
80 :     call. = FALSE)
81 :     ## return the common ncol()
82 :     da[2]
83 :     }
84 :    
85 : maechler 1227 isPacked <- function(x)
86 :     {
87 :     ## Is 'x' a packed (dense) matrix ?
88 :     is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
89 :     }
90 :    
91 : maechler 956 emptyColnames <- function(x)
92 :     {
93 :     ## Useful for compact printing of (parts) of sparse matrices
94 :     ## possibly dimnames(x) "==" NULL :
95 :     dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2]))
96 :     x
97 :     }
98 : maechler 908
99 : maechler 919 prTriang <- function(x, digits = getOption("digits"),
100 :     justify = "none", right = TRUE)
101 :     {
102 :     ## modeled along stats:::print.dist
103 :     diag <- TRUE
104 :     upper <- x@uplo == "U"
105 :    
106 :     m <- as(x, "matrix")
107 :     cf <- format(m, digits = digits, justify = justify)
108 :     if(upper)
109 :     cf[row(cf) > col(cf)] <- "."
110 :     else
111 :     cf[row(cf) < col(cf)] <- "."
112 :     print(cf, quote = FALSE, right = right)
113 :     invisible(x)
114 :     }
115 :    
116 :     prMatrix <- function(x, digits = getOption("digits")) {
117 :     d <- dim(x)
118 :     cl <- class(x)
119 :     cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
120 :     maxp <- getOption("max.print")
121 :     if(prod(d) <= maxp) {
122 :     if(is(x, "triangularMatrix"))
123 :     prTriang(x, digits = digits)
124 :     else
125 :     print(as(x, "matrix"), digits = digits)
126 :     }
127 :     else { ## d[1] > maxp / d[2] >= nr :
128 :     m <- as(x, "matrix")
129 :     nr <- maxp %/% d[2]
130 :     n2 <- ceiling(nr / 2)
131 :     print(head(m, max(1, n2)))
132 :     cat("\n ..........\n\n")
133 :     print(tail(m, max(1, nr - n2)))
134 :     }
135 :     ## DEBUG: cat("str(.):\n") ; str(x)
136 :     invisible(x)# as print() S3 methods do
137 :     }
138 :    
139 :     ## For sparseness handling
140 :     non0ind <- function(x) {
141 :     if(is.numeric(x))
142 : maechler 1226 return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))
143 :     ## else
144 : maechler 919 stopifnot(is(x, "sparseMatrix"))
145 : maechler 1226 ## return a 2-column (i,j) matrix of
146 :     ## 0-based indices of non-zero entries :
147 : maechler 925 if(is(x, "TsparseMatrix"))
148 : maechler 954 return(unique(cbind(x@i,x@j)))
149 : maechler 1226 ## else:
150 :     isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
151 :     .Call("compressed_non_0_ij", x, isC, PACKAGE = "Matrix")
152 :     }
153 : maechler 954
154 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings:
155 :     ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
156 :     encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
157 :     decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
158 :    
159 :     complementInd <- function(ij, dim)
160 :     {
161 :     ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
162 :     ## but as 1-based indices
163 :     n <- prod(dim)
164 :     if(n == 0) return(integer(0))
165 :     ii <- 1:n
166 :     ii[-(1 + encodeInd(ij, nr = dim[1]))]
167 : maechler 919 }
168 :    
169 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
170 :    
171 :     intersectInd <- function(ij1, ij2, nrow) {
172 :     ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
173 :     ## return only the *common* entries
174 :     decodeInd(intersect(encodeInd(ij1, nrow),
175 :     encodeInd(ij2, nrow)), nrow)
176 :     }
177 :    
178 :     WhichintersectInd <- function(ij1, ij2, nrow) {
179 :     ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
180 :     ## find *where* common entries are in ij1 & ij2
181 :     m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
182 :     ni <- !is.na(m1)
183 :     list(which(ni), m1[ni])
184 :     }
185 :    
186 :    
187 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R !
188 : maechler 919 uniq <- function(x) {
189 : maechler 925 if(is(x, "TsparseMatrix")) {
190 : maechler 919 ## Purpose: produce a *unique* triplet representation:
191 :     ## by having (i,j) sorted and unique
192 :     ## -----------------------------------------------------------
193 : maechler 1226 ## The following is not quite efficient {but easy to program,
194 :     ## and both as() are based on C code
195 : maechler 919 if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
196 :     else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
197 :     else stop("not implemented for class", class(x))
198 :    
199 : maechler 1226 } else x ## not 'gT' ; i.e. "uniquely" represented in any case
200 : maechler 919 }
201 :    
202 :     if(FALSE) ## try an "efficient" version
203 :     uniq_gT <- function(x)
204 :     {
205 :     ## Purpose: produce a *unique* triplet representation:
206 :     ## by having (i,j) sorted and unique
207 : maechler 1226 ## ------------------------------------------------------------------
208 : maechler 919 ## Arguments: a "gT" Matrix
209 :     stopifnot(is(x, "gTMatrix"))
210 :     if((n <- length(x@i)) == 0) return(x)
211 :     ii <- order(x@i, x@j)
212 :     if(any(ii != 1:n)) {
213 :     x@i <- x@i[ii]
214 :     x@j <- x@j[ii]
215 :     x@x <- x@x[ii]
216 :     }
217 :     ij <- x@i + nrow(x) * x@j
218 :     if(any(dup <- duplicated(ij))) {
219 :    
220 :     }
221 :     ### We should use a .Call() based utility for this!
222 :    
223 :     }
224 :    
225 : maechler 946 t_geMatrix <- function(x) {
226 :     x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
227 :     x@Dim <- x@Dim[2:1]
228 :     x@Dimnames <- x@Dimnames[2:1]
229 :     ## FIXME: how to set factors?
230 :     x
231 :     }
232 :    
233 :     ## t( [dl]trMatrix ) and t( [dl]syMatrix ) :
234 :     t_trMatrix <- function(x) {
235 :     x@x <- as.vector(t(as(x, "matrix")))
236 :     x@Dim <- x@Dim[2:1]
237 :     x@Dimnames <- x@Dimnames[2:1]
238 :     x@uplo <- if (x@uplo == "U") "L" else "U"
239 :     # and keep x@diag
240 :     x
241 :     }
242 : maechler 956
243 :     fixupDense <- function(m, from) {
244 :     if(is(m, "triangularMatrix")) {
245 :     m@uplo <- from@uplo
246 :     m@diag <- from@diag
247 :     } else if(is(m, "symmetricMatrix")) {
248 :     m@uplo <- from@uplo
249 :     }
250 :     m
251 :     }
252 :    
253 :     ## -> ./ldenseMatrix.R :
254 :     l2d_Matrix <- function(from) {
255 :     stopifnot(is(from, "lMatrix"))
256 :     fixupDense(new(sub("^l", "d", class(from)),
257 :     x = as.double(from@x),
258 : maechler 1198 Dim = from@Dim, Dimnames = from@Dimnames),
259 : maechler 956 from)
260 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
261 : maechler 956 }
262 :    
263 :     if(FALSE)# unused
264 :     l2d_meth <- function(x) {
265 :     cl <- class(x)
266 :     as(callGeneric(as(x, sub("^l", "d", cl))), cl)
267 :     }
268 :    
269 : maechler 1226 dClass2 <- function(dClass, kind = "l") {
270 :     ## Find "corresponding" class for a dMatrix;
271 :     # since pos.def. matrices have no pendant:
272 :     if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
273 :     else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
274 :     else sub("^d", kind, dClass)
275 :     }
276 :    
277 :     geClass <- function(x) {
278 :     if(is(x, "dMatrix")) "dgeMatrix"
279 :     else if(is(x, "lMatrix")) "lgeMatrix"
280 :     else stop("general Matrix class not yet implemented for",
281 :     class(x))
282 :     }
283 :    
284 : maechler 956 ## -> ./ddenseMatrix.R :
285 :     d2l_Matrix <- function(from) {
286 :     stopifnot(is(from, "dMatrix"))
287 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
288 : maechler 1198 Dim = from@Dim, Dimnames = from@Dimnames),
289 : maechler 956 from)
290 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
291 : maechler 956 }
292 : maechler 973
293 :    
294 :     try_as <- function(x, classes, tryAnyway = FALSE) {
295 :     if(!tryAnyway && !is(x, "Matrix"))
296 :     return(x)
297 :     ## else
298 :     ok <- canCoerce(x, classes[1])
299 :     while(!ok && length(classes <- classes[-1])) {
300 :     ok <- canCoerce(x, classes[1])
301 :     }
302 :     if(ok) as(x, classes[1]) else x
303 :     }
304 :    
305 : maechler 1174 if(paste(R.version$major, R.version$minor, sep=".") < "2.3")
306 :     ## This will be in R 2.3.0
307 : maechler 973 canCoerce <- function(object, Class) {
308 :     ## Purpose: test if 'object' is coercable to 'Class', i.e.,
309 :     ## as(object, Class) will {typically} work
310 :     ## ----------------------------------------------------------------------
311 :     ## Author: John Chambers, Date: 6 Oct 2005
312 :     is(object, Class) ||
313 :     !is.null(selectMethod("coerce", c(class(object), Class),
314 :     optional = TRUE,
315 :     useInherited = c(from = TRUE, to = FALSE)))
316 :     }
317 :    
318 : maechler 1226 .is.triangular <- function(object, upper = NA) {
319 : maechler 1174 ## pretest: is it square?
320 :     d <- dim(object)
321 :     if(d[1] != d[2]) return(FALSE)
322 :     ## else slower test
323 :     if(!is.matrix(object))
324 : maechler 1226 object <- as(object,"matrix")
325 : maechler 1174 ## == 0 even works for logical & complex:
326 : maechler 1226 if(is.na(upper)) {
327 :     if(all(object[lower.tri(object)] == 0))
328 :     structure(TRUE, kind = "U")
329 :     else if(all(object[upper.tri(object)] == 0))
330 :     structure(TRUE, kind = "L")
331 :     else FALSE
332 :     } else if(upper)
333 :     all(object[lower.tri(object)] == 0)
334 :     else ## upper is FALSE
335 :     all(object[upper.tri(object)] == 0)
336 : maechler 1174 }
337 :    
338 :     .is.diagonal <- function(object) {
339 :     d <- dim(object)
340 :     if(d[1] != (n <- d[2])) FALSE
341 :     else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)
342 :     }

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