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 1349 - (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 1332 .M.DN <- function(x) if(!is.null(dn <- dimnames(x))) dn else list(NULL,NULL)
9 :    
10 : maechler 656 .has.DN <- ## has non-trivial Dimnames slot?
11 :     function(x) !identical(list(NULL,NULL), x@Dimnames)
12 :    
13 : maechler 949 .bail.out.1 <- function(fun, cl) {
14 :     stop(gettextf('not-yet-implemented method for %s(<%s>)', fun, cl),
15 : maechler 1238 call. = FALSE)
16 : maechler 949 }
17 :     .bail.out.2 <- function(fun, cl1, cl2) {
18 :     stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)',
19 : maechler 1238 fun, cl1, cl2), call. = FALSE)
20 : maechler 949 }
21 :    
22 : maechler 632 ## chol() via "dpoMatrix"
23 :     cholMat <- function(x, pivot, LINPACK) {
24 :     px <- as(x, "dpoMatrix")
25 : bates 703 if (isTRUE(validObject(px, test=TRUE))) chol(px)
26 : maechler 632 else stop("'x' is not positive definite -- chol() undefined.")
27 :     }
28 : maechler 908
29 : maechler 954 dimCheck <- function(a, b) {
30 :     da <- dim(a)
31 :     db <- dim(b)
32 :     if(any(da != db))
33 :     stop(gettextf("Matrices must have same dimensions in %s",
34 :     deparse(sys.call(sys.parent()))),
35 :     call. = FALSE)
36 :     da
37 :     }
38 :    
39 : maechler 956 dimNamesCheck <- function(a, b) {
40 :     ## assume dimCheck() has happened before
41 :     nullDN <- list(NULL,NULL)
42 :     h.a <- !identical(nullDN, dna <- dimnames(a))
43 :     h.b <- !identical(nullDN, dnb <- dimnames(b))
44 :     if(h.a || h.b) {
45 : maechler 1084 if (!h.b) dna
46 :     else if(!h.a) dnb
47 : maechler 956 else { ## both have non-trivial dimnames
48 :     r <- dna # "default" result
49 :     for(j in 1:2) {
50 :     dn <- dnb[[j]]
51 :     if(is.null(r[[j]]))
52 :     r[[j]] <- dn
53 :     else if (!is.null(dn) && any(r[[j]] != dn))
54 :     warning(gettextf("dimnames [%d] mismatch in %s", j,
55 :     deparse(sys.call(sys.parent()))),
56 :     call. = FALSE)
57 :     }
58 :     r
59 :     }
60 :     }
61 :     else
62 :     nullDN
63 :     }
64 :    
65 : maechler 908 rowCheck <- function(a, b) {
66 :     da <- dim(a)
67 :     db <- dim(b)
68 :     if(da[1] != db[1])
69 :     stop(gettextf("Matrices must have same number of rows in %s",
70 :     deparse(sys.call(sys.parent()))),
71 :     call. = FALSE)
72 :     ## return the common nrow()
73 :     da[1]
74 :     }
75 :    
76 :     colCheck <- function(a, b) {
77 :     da <- dim(a)
78 :     db <- dim(b)
79 :     if(da[2] != db[2])
80 :     stop(gettextf("Matrices must have same number of columns in %s",
81 :     deparse(sys.call(sys.parent()))),
82 :     call. = FALSE)
83 :     ## return the common ncol()
84 :     da[2]
85 :     }
86 :    
87 : maechler 1285 ## Note: !isPacked(.) i.e. `full' still contains
88 :     ## ---- "*sy" and "*tr" which have "undefined" lower or upper part
89 : maechler 1227 isPacked <- function(x)
90 :     {
91 :     ## Is 'x' a packed (dense) matrix ?
92 :     is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
93 :     }
94 :    
95 : maechler 956 emptyColnames <- function(x)
96 :     {
97 :     ## Useful for compact printing of (parts) of sparse matrices
98 : maechler 1238 ## possibly dimnames(x) "==" NULL :
99 : maechler 956 dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2]))
100 :     x
101 :     }
102 : maechler 908
103 : maechler 919 prTriang <- function(x, digits = getOption("digits"),
104 : maechler 1238 justify = "none", right = TRUE)
105 : maechler 919 {
106 :     ## modeled along stats:::print.dist
107 :     diag <- TRUE
108 :     upper <- x@uplo == "U"
109 :    
110 :     m <- as(x, "matrix")
111 :     cf <- format(m, digits = digits, justify = justify)
112 :     if(upper)
113 : maechler 1238 cf[row(cf) > col(cf)] <- "."
114 : maechler 919 else
115 : maechler 1238 cf[row(cf) < col(cf)] <- "."
116 : maechler 919 print(cf, quote = FALSE, right = right)
117 :     invisible(x)
118 :     }
119 :    
120 :     prMatrix <- function(x, digits = getOption("digits")) {
121 :     d <- dim(x)
122 :     cl <- class(x)
123 :     cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
124 :     maxp <- getOption("max.print")
125 :     if(prod(d) <= maxp) {
126 : maechler 1238 if(is(x, "triangularMatrix"))
127 :     prTriang(x, digits = digits)
128 :     else
129 :     print(as(x, "matrix"), digits = digits)
130 : maechler 919 }
131 :     else { ## d[1] > maxp / d[2] >= nr :
132 : maechler 1238 m <- as(x, "matrix")
133 : maechler 919 nr <- maxp %/% d[2]
134 :     n2 <- ceiling(nr / 2)
135 :     print(head(m, max(1, n2)))
136 :     cat("\n ..........\n\n")
137 :     print(tail(m, max(1, nr - n2)))
138 :     }
139 :     ## DEBUG: cat("str(.):\n") ; str(x)
140 :     invisible(x)# as print() S3 methods do
141 :     }
142 :    
143 :     ## For sparseness handling
144 :     non0ind <- function(x) {
145 :     if(is.numeric(x))
146 : maechler 1226 return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))
147 :     ## else
148 : maechler 919 stopifnot(is(x, "sparseMatrix"))
149 : maechler 1226 ## return a 2-column (i,j) matrix of
150 :     ## 0-based indices of non-zero entries :
151 : maechler 1319 non0.i <- function(M) {
152 :     if(is(M, "TsparseMatrix"))
153 :     return(unique(cbind(M@i,M@j)))
154 :     if(is(M, "pMatrix"))
155 :     return(cbind(seq(length=nrow(M)), M@perm) - 1:1)
156 :     ## else:
157 :     isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
158 :     .Call(compressed_non_0_ij, M, isC)
159 :     }
160 :    
161 :     if(is(x, "symmetricMatrix")) { # also get "other" triangle
162 :     ij <- non0.i(x)
163 :     notdiag <- ij[,1] != ij[,2]# but not the diagonals again
164 :     rbind(ij, ij[notdiag, 2:1])
165 :     }
166 :     else
167 :     non0.i(x)
168 : maechler 1226 }
169 : maechler 954
170 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings:
171 :     ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
172 :     encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
173 :     decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
174 :    
175 :     complementInd <- function(ij, dim)
176 :     {
177 :     ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
178 : maechler 1238 ## but as 1-based indices
179 : maechler 1226 n <- prod(dim)
180 :     if(n == 0) return(integer(0))
181 :     ii <- 1:n
182 :     ii[-(1 + encodeInd(ij, nr = dim[1]))]
183 : maechler 919 }
184 :    
185 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
186 :    
187 :     intersectInd <- function(ij1, ij2, nrow) {
188 :     ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
189 :     ## return only the *common* entries
190 :     decodeInd(intersect(encodeInd(ij1, nrow),
191 :     encodeInd(ij2, nrow)), nrow)
192 :     }
193 :    
194 :     WhichintersectInd <- function(ij1, ij2, nrow) {
195 :     ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
196 :     ## find *where* common entries are in ij1 & ij2
197 :     m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
198 :     ni <- !is.na(m1)
199 :     list(which(ni), m1[ni])
200 :     }
201 :    
202 :    
203 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R !
204 : maechler 1270
205 : maechler 1331 uniqTsparse <- function(x, class.x = c(class(x))) {
206 :     ## Purpose: produce a *unique* triplet representation:
207 :     ## by having (i,j) sorted and unique
208 :     ## -----------------------------------------------------------
209 :     ## The following is not quite efficient {but easy to program,
210 :     ## and as() are based on C code (all of them?)
211 :     ##
212 :     ## FIXME: Do it fast for the case where 'x' is already 'uniq'
213 : maechler 1270
214 : maechler 1331 switch(class.x,
215 :     "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
216 :     "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
217 :     "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
218 :     ## do we need this for "logical" ones, there's no sum() there!
219 :     "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
220 :     "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
221 :     "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
222 :     ## otherwise:
223 :     stop("not yet implemented for class ", clx))
224 :     }
225 : maechler 919
226 : maechler 1331 ## Note: maybe, using
227 :     ## ---- xj <- .Call(Matrix_expand_pointers, x@p)
228 :     ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
229 :     ## but really efficient would be to use only one .Call(.) for uniq(.) !
230 :    
231 :     uniq <- function(x) {
232 :     if(is(x, "TsparseMatrix")) uniqTsparse(x) else x
233 :     ## else: not 'Tsparse', i.e. "uniquely" represented in any case
234 : maechler 919 }
235 :    
236 :     if(FALSE) ## try an "efficient" version
237 :     uniq_gT <- function(x)
238 :     {
239 :     ## Purpose: produce a *unique* triplet representation:
240 :     ## by having (i,j) sorted and unique
241 : maechler 1226 ## ------------------------------------------------------------------
242 : maechler 919 ## Arguments: a "gT" Matrix
243 :     stopifnot(is(x, "gTMatrix"))
244 :     if((n <- length(x@i)) == 0) return(x)
245 :     ii <- order(x@i, x@j)
246 :     if(any(ii != 1:n)) {
247 : maechler 1238 x@i <- x@i[ii]
248 :     x@j <- x@j[ii]
249 :     x@x <- x@x[ii]
250 : maechler 919 }
251 :     ij <- x@i + nrow(x) * x@j
252 :     if(any(dup <- duplicated(ij))) {
253 :    
254 :     }
255 :     ### We should use a .Call() based utility for this!
256 :    
257 :     }
258 :    
259 : maechler 946 t_geMatrix <- function(x) {
260 :     x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
261 :     x@Dim <- x@Dim[2:1]
262 :     x@Dimnames <- x@Dimnames[2:1]
263 :     ## FIXME: how to set factors?
264 :     x
265 :     }
266 :    
267 :     ## t( [dl]trMatrix ) and t( [dl]syMatrix ) :
268 :     t_trMatrix <- function(x) {
269 :     x@x <- as.vector(t(as(x, "matrix")))
270 :     x@Dim <- x@Dim[2:1]
271 :     x@Dimnames <- x@Dimnames[2:1]
272 :     x@uplo <- if (x@uplo == "U") "L" else "U"
273 :     # and keep x@diag
274 :     x
275 :     }
276 : maechler 956
277 :     fixupDense <- function(m, from) {
278 :     if(is(m, "triangularMatrix")) {
279 : maechler 1238 m@uplo <- from@uplo
280 :     m@diag <- from@diag
281 : maechler 956 } else if(is(m, "symmetricMatrix")) {
282 : maechler 1238 m@uplo <- from@uplo
283 : maechler 956 }
284 :     m
285 :     }
286 :    
287 :     ## -> ./ldenseMatrix.R :
288 :     l2d_Matrix <- function(from) {
289 :     stopifnot(is(from, "lMatrix"))
290 :     fixupDense(new(sub("^l", "d", class(from)),
291 : maechler 1238 x = as.double(from@x),
292 :     Dim = from@Dim, Dimnames = from@Dimnames),
293 :     from)
294 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
295 : maechler 956 }
296 :    
297 :     if(FALSE)# unused
298 :     l2d_meth <- function(x) {
299 :     cl <- class(x)
300 :     as(callGeneric(as(x, sub("^l", "d", cl))), cl)
301 :     }
302 :    
303 : maechler 1331 ## return "d" or "l" or "z"
304 :     .M.kind <- function(x, clx = class(x)) {
305 :     if(is.matrix(x)) { ## 'old style matrix'
306 :     if (is.numeric(x)) "d"
307 :     else if(is.logical(x)) "l"
308 :     else if(is.complex(x)) "z"
309 :     else stop("not yet implemented for matrix w/ typeof ", typeof(x))
310 :     }
311 :     else if(extends(clx, "dMatrix")) "d"
312 :     else if(extends(clx, "lMatrix")) "l"
313 :     else if(extends(clx, "zMatrix")) "z"
314 :     else stop(" not yet be implemented for ", clx)
315 :     }
316 :    
317 :     .M.shape <- function(x, clx = class(x)) {
318 :     if(is.matrix(x)) { ## 'old style matrix'
319 :     if (isDiagonal (x)) "d"
320 :     else if(isTriangular(x)) "t"
321 :     else if(isSymmetric (x)) "s"
322 :     else "g" # general
323 :     }
324 :     else if(extends(clx, "diagonalMatrix")) "d"
325 :     else if(extends(clx, "triangularMatrix"))"t"
326 :     else if(extends(clx, "symmetricMatrix")) "s"
327 :     else "g"
328 :     }
329 :    
330 :    
331 : maechler 1329 class2 <- function(cl, kind = "l", do.sub = TRUE) {
332 :     ## Find "corresponding" class; since pos.def. matrices have no pendant:
333 :     if (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
334 :     else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
335 :     else if(do.sub) sub("^d", kind, cl)
336 :     else cl
337 : maechler 1226 }
338 :    
339 :     geClass <- function(x) {
340 : maechler 1331 if (is(x, "dMatrix")) "dgeMatrix"
341 : maechler 1226 else if(is(x, "lMatrix")) "lgeMatrix"
342 : maechler 1331 else if(is(x, "zMatrix")) "zgeMatrix"
343 : maechler 1329 else stop("general Matrix class not yet implemented for ",
344 : maechler 1226 class(x))
345 :     }
346 : maechler 1331
347 :     .dense.prefixes <- c("d" = "di",
348 :     "t" = "tr",
349 :     "s" = "sy",
350 :     "g" = "ge")
351 :    
352 : maechler 1349 .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
353 :     "t" = "t",
354 :     "s" = "s",
355 :     "g" = "g")
356 : maechler 1331
357 : maechler 1329 ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
358 : maechler 1331 as_dense <- function(x) {
359 :     as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
360 :     }
361 :    
362 :     as_Csparse <- function(x) {
363 : maechler 1349 as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
364 : maechler 1331 }
365 : maechler 1349 as_Rsparse <- function(x) {
366 :     as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
367 :     }
368 :     as_Tsparse <- function(x) {
369 :     as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
370 :     }
371 : maechler 1331
372 : maechler 1329 as_geClass <- function(x, cl) {
373 : maechler 1331 if (extends(cl, "diagonalMatrix") && isDiagonal(x))
374 :     as(x, cl)
375 :     else if(extends(cl, "symmetricMatrix") && isSymmetric(x))
376 :     as(x, class2(cl, kind, do.sub= kind != "d"))
377 :     else if(extends(cl, "triangularMatrix") && isTriangular(x))
378 :     as(x, cl)
379 :     else
380 :     as(x, paste(.M.kind(x), "geMatrix", sep=''))
381 :     }
382 : maechler 1226
383 : maechler 1331 as_CspClass <- function(x, cl) {
384 :     if ((extends(cl, "diagonalMatrix") && isDiagonal(x)) ||
385 :     (extends(cl, "symmetricMatrix") && isSymmetric(x)) ||
386 :     (extends(cl, "triangularMatrix")&& isTriangular(x)))
387 :     as(x, cl)
388 :     else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
389 : maechler 1329 }
390 :    
391 : maechler 1331
392 : maechler 956 ## -> ./ddenseMatrix.R :
393 :     d2l_Matrix <- function(from) {
394 :     stopifnot(is(from, "dMatrix"))
395 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
396 : maechler 1238 Dim = from@Dim, Dimnames = from@Dimnames),
397 :     from)
398 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
399 : maechler 956 }
400 : maechler 973
401 :    
402 :     try_as <- function(x, classes, tryAnyway = FALSE) {
403 :     if(!tryAnyway && !is(x, "Matrix"))
404 : maechler 1238 return(x)
405 : maechler 973 ## else
406 :     ok <- canCoerce(x, classes[1])
407 :     while(!ok && length(classes <- classes[-1])) {
408 : maechler 1238 ok <- canCoerce(x, classes[1])
409 : maechler 973 }
410 :     if(ok) as(x, classes[1]) else x
411 :     }
412 :    
413 :    
414 : maechler 1238 ## For *dense* matrices
415 :     isTriMat <- function(object, upper = NA) {
416 : maechler 1174 ## pretest: is it square?
417 :     d <- dim(object)
418 :     if(d[1] != d[2]) return(FALSE)
419 :     ## else slower test
420 :     if(!is.matrix(object))
421 : maechler 1226 object <- as(object,"matrix")
422 : maechler 1174 ## == 0 even works for logical & complex:
423 : maechler 1226 if(is.na(upper)) {
424 :     if(all(object[lower.tri(object)] == 0))
425 :     structure(TRUE, kind = "U")
426 :     else if(all(object[upper.tri(object)] == 0))
427 :     structure(TRUE, kind = "L")
428 :     else FALSE
429 :     } else if(upper)
430 :     all(object[lower.tri(object)] == 0)
431 :     else ## upper is FALSE
432 :     all(object[upper.tri(object)] == 0)
433 : maechler 1174 }
434 :    
435 : maechler 1238 ## For Csparse matrices
436 :     isTriC <- function(x, upper = NA) {
437 :     ## pretest: is it square?
438 :     d <- dim(x)
439 :     if(d[1] != d[2]) return(FALSE)
440 :     ## else
441 :     if(d[1] == 0) return(TRUE)
442 : maechler 1315 ni <- 1:d[2]
443 : maechler 1238 ## the row indices split according to column:
444 :     ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
445 :     lil <- unlist(lapply(ilist, length), use.names = FALSE)
446 :     if(any(lil == 0)) {
447 :     pos <- lil > 0
448 :     if(!any(pos)) ## matrix of all 0's
449 :     return(TRUE)
450 :     ilist <- ilist[pos]
451 :     ni <- ni[pos]
452 :     }
453 : maechler 1315 ni0 <- ni - 1:1 # '0-based ni'
454 : maechler 1238 if(is.na(upper)) {
455 : maechler 1315 if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
456 : maechler 1238 structure(TRUE, kind = "U")
457 : maechler 1315 else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
458 : maechler 1238 structure(TRUE, kind = "L")
459 :     else FALSE
460 :     } else if(upper) {
461 : maechler 1315 all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
462 : maechler 1238 } else { ## 'lower'
463 : maechler 1315 all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
464 : maechler 1238 }
465 :     }
466 :    
467 : maechler 1174 .is.diagonal <- function(object) {
468 :     d <- dim(object)
469 :     if(d[1] != (n <- d[2])) FALSE
470 :     else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)
471 :     }
472 : maechler 1253
473 :     diagU2N <- function(x)
474 :     {
475 :     ## Purpose: Transform a *unit diagonal* triangular matrix
476 :     ## into one with explicit diagonal entries '1'
477 :     xT <- as(x, "dgTMatrix")
478 :     ## leave it as T* - the caller can always coerce to C* if needed:
479 :     new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
480 :     Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
481 :     }
482 : maechler 1290
483 : maechler 1349 ## FIXME: this should probably be dropped / replaced by as_Csparse
484 : maechler 1295 .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
485 :     x <- as(x, "dgCMatrix")
486 :     callGeneric()
487 :     }
488 :    
489 :     .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
490 : maechler 1349 ## used e.g. inside colSums() etc methods
491 : maechler 1295 x <- as(x, "dgTMatrix")
492 :     callGeneric()
493 :     }
494 :    
495 :    
496 : maechler 1290 ### Fast much simplified version of tapply()
497 :     tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
498 :     sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
499 :     }
500 :    
501 :     ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
502 :     ## tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
503 :     ## }
504 :    

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