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 1315 - (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 : maechler 1238 call. = FALSE)
14 : maechler 949 }
15 :     .bail.out.2 <- function(fun, cl1, cl2) {
16 :     stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>)',
17 : maechler 1238 fun, cl1, cl2), call. = FALSE)
18 : maechler 949 }
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 1285 ## Note: !isPacked(.) i.e. `full' still contains
86 :     ## ---- "*sy" and "*tr" which have "undefined" lower or upper part
87 : maechler 1227 isPacked <- function(x)
88 :     {
89 :     ## Is 'x' a packed (dense) matrix ?
90 :     is(x,"Matrix") && !is.null(x@x) && length(x@x) < prod(dim(x))
91 :     }
92 :    
93 : maechler 956 emptyColnames <- function(x)
94 :     {
95 :     ## Useful for compact printing of (parts) of sparse matrices
96 : maechler 1238 ## possibly dimnames(x) "==" NULL :
97 : maechler 956 dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2]))
98 :     x
99 :     }
100 : maechler 908
101 : maechler 919 prTriang <- function(x, digits = getOption("digits"),
102 : maechler 1238 justify = "none", right = TRUE)
103 : maechler 919 {
104 :     ## modeled along stats:::print.dist
105 :     diag <- TRUE
106 :     upper <- x@uplo == "U"
107 :    
108 :     m <- as(x, "matrix")
109 :     cf <- format(m, digits = digits, justify = justify)
110 :     if(upper)
111 : maechler 1238 cf[row(cf) > col(cf)] <- "."
112 : maechler 919 else
113 : maechler 1238 cf[row(cf) < col(cf)] <- "."
114 : maechler 919 print(cf, quote = FALSE, right = right)
115 :     invisible(x)
116 :     }
117 :    
118 :     prMatrix <- function(x, digits = getOption("digits")) {
119 :     d <- dim(x)
120 :     cl <- class(x)
121 :     cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
122 :     maxp <- getOption("max.print")
123 :     if(prod(d) <= maxp) {
124 : maechler 1238 if(is(x, "triangularMatrix"))
125 :     prTriang(x, digits = digits)
126 :     else
127 :     print(as(x, "matrix"), digits = digits)
128 : maechler 919 }
129 :     else { ## d[1] > maxp / d[2] >= nr :
130 : maechler 1238 m <- as(x, "matrix")
131 : maechler 919 nr <- maxp %/% d[2]
132 :     n2 <- ceiling(nr / 2)
133 :     print(head(m, max(1, n2)))
134 :     cat("\n ..........\n\n")
135 :     print(tail(m, max(1, nr - n2)))
136 :     }
137 :     ## DEBUG: cat("str(.):\n") ; str(x)
138 :     invisible(x)# as print() S3 methods do
139 :     }
140 :    
141 :     ## For sparseness handling
142 :     non0ind <- function(x) {
143 :     if(is.numeric(x))
144 : maechler 1226 return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))
145 :     ## else
146 : maechler 919 stopifnot(is(x, "sparseMatrix"))
147 : maechler 1226 ## return a 2-column (i,j) matrix of
148 :     ## 0-based indices of non-zero entries :
149 : maechler 925 if(is(x, "TsparseMatrix"))
150 : maechler 954 return(unique(cbind(x@i,x@j)))
151 : maechler 1315 if(is(x, "pMatrix"))
152 :     return(cbind(seq(length=nrow(x)), x@perm) - 1:1)
153 : maechler 1226 ## else:
154 :     isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
155 : maechler 1280 .Call(compressed_non_0_ij, x, isC)
156 : maechler 1226 }
157 : maechler 954
158 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings:
159 :     ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
160 :     encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
161 :     decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
162 :    
163 :     complementInd <- function(ij, dim)
164 :     {
165 :     ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
166 : maechler 1238 ## but as 1-based indices
167 : maechler 1226 n <- prod(dim)
168 :     if(n == 0) return(integer(0))
169 :     ii <- 1:n
170 :     ii[-(1 + encodeInd(ij, nr = dim[1]))]
171 : maechler 919 }
172 :    
173 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
174 :    
175 :     intersectInd <- function(ij1, ij2, nrow) {
176 :     ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
177 :     ## return only the *common* entries
178 :     decodeInd(intersect(encodeInd(ij1, nrow),
179 :     encodeInd(ij2, nrow)), nrow)
180 :     }
181 :    
182 :     WhichintersectInd <- function(ij1, ij2, nrow) {
183 :     ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
184 :     ## find *where* common entries are in ij1 & ij2
185 :     m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
186 :     ni <- !is.na(m1)
187 :     list(which(ni), m1[ni])
188 :     }
189 :    
190 :    
191 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R !
192 : maechler 919 uniq <- function(x) {
193 : maechler 1270
194 :     ### Note: maybe, using
195 : maechler 1280 ### ---- xj <- .Call(Matrix_expand_pointers, x@p)
196 : maechler 1270 ### would be slightly more efficient than as( <dgC> , "dgTMatrix")
197 :     ### but really efficient would be to use only one .Call(.) for uniq(.) !
198 :     ### Try to do it particularly fast for the case where 'x' is already a 'uniq' <dgT>
199 :    
200 : maechler 925 if(is(x, "TsparseMatrix")) {
201 : maechler 1238 ## Purpose: produce a *unique* triplet representation:
202 :     ## by having (i,j) sorted and unique
203 :     ## -----------------------------------------------------------
204 :     ## The following is not quite efficient {but easy to program,
205 :     ## and both as() are based on C code
206 :     if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
207 :     else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
208 :     else stop("not implemented for class", class(x))
209 : maechler 919
210 : maechler 1226 } else x ## not 'gT' ; i.e. "uniquely" represented in any case
211 : maechler 919 }
212 :    
213 :     if(FALSE) ## try an "efficient" version
214 :     uniq_gT <- function(x)
215 :     {
216 :     ## Purpose: produce a *unique* triplet representation:
217 :     ## by having (i,j) sorted and unique
218 : maechler 1226 ## ------------------------------------------------------------------
219 : maechler 919 ## Arguments: a "gT" Matrix
220 :     stopifnot(is(x, "gTMatrix"))
221 :     if((n <- length(x@i)) == 0) return(x)
222 :     ii <- order(x@i, x@j)
223 :     if(any(ii != 1:n)) {
224 : maechler 1238 x@i <- x@i[ii]
225 :     x@j <- x@j[ii]
226 :     x@x <- x@x[ii]
227 : maechler 919 }
228 :     ij <- x@i + nrow(x) * x@j
229 :     if(any(dup <- duplicated(ij))) {
230 :    
231 :     }
232 :     ### We should use a .Call() based utility for this!
233 :    
234 :     }
235 :    
236 : maechler 946 t_geMatrix <- function(x) {
237 :     x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
238 :     x@Dim <- x@Dim[2:1]
239 :     x@Dimnames <- x@Dimnames[2:1]
240 :     ## FIXME: how to set factors?
241 :     x
242 :     }
243 :    
244 :     ## t( [dl]trMatrix ) and t( [dl]syMatrix ) :
245 :     t_trMatrix <- function(x) {
246 :     x@x <- as.vector(t(as(x, "matrix")))
247 :     x@Dim <- x@Dim[2:1]
248 :     x@Dimnames <- x@Dimnames[2:1]
249 :     x@uplo <- if (x@uplo == "U") "L" else "U"
250 :     # and keep x@diag
251 :     x
252 :     }
253 : maechler 956
254 :     fixupDense <- function(m, from) {
255 :     if(is(m, "triangularMatrix")) {
256 : maechler 1238 m@uplo <- from@uplo
257 :     m@diag <- from@diag
258 : maechler 956 } else if(is(m, "symmetricMatrix")) {
259 : maechler 1238 m@uplo <- from@uplo
260 : maechler 956 }
261 :     m
262 :     }
263 :    
264 :     ## -> ./ldenseMatrix.R :
265 :     l2d_Matrix <- function(from) {
266 :     stopifnot(is(from, "lMatrix"))
267 :     fixupDense(new(sub("^l", "d", class(from)),
268 : maechler 1238 x = as.double(from@x),
269 :     Dim = from@Dim, Dimnames = from@Dimnames),
270 :     from)
271 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
272 : maechler 956 }
273 :    
274 :     if(FALSE)# unused
275 :     l2d_meth <- function(x) {
276 :     cl <- class(x)
277 :     as(callGeneric(as(x, sub("^l", "d", cl))), cl)
278 :     }
279 :    
280 : maechler 1226 dClass2 <- function(dClass, kind = "l") {
281 :     ## Find "corresponding" class for a dMatrix;
282 :     # since pos.def. matrices have no pendant:
283 :     if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
284 :     else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
285 :     else sub("^d", kind, dClass)
286 :     }
287 :    
288 :     geClass <- function(x) {
289 :     if(is(x, "dMatrix")) "dgeMatrix"
290 :     else if(is(x, "lMatrix")) "lgeMatrix"
291 :     else stop("general Matrix class not yet implemented for",
292 :     class(x))
293 :     }
294 :    
295 : maechler 956 ## -> ./ddenseMatrix.R :
296 :     d2l_Matrix <- function(from) {
297 :     stopifnot(is(from, "dMatrix"))
298 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
299 : maechler 1238 Dim = from@Dim, Dimnames = from@Dimnames),
300 :     from)
301 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
302 : maechler 956 }
303 : maechler 973
304 :    
305 :     try_as <- function(x, classes, tryAnyway = FALSE) {
306 :     if(!tryAnyway && !is(x, "Matrix"))
307 : maechler 1238 return(x)
308 : maechler 973 ## else
309 :     ok <- canCoerce(x, classes[1])
310 :     while(!ok && length(classes <- classes[-1])) {
311 : maechler 1238 ok <- canCoerce(x, classes[1])
312 : maechler 973 }
313 :     if(ok) as(x, classes[1]) else x
314 :     }
315 :    
316 :    
317 : maechler 1238 ## For *dense* matrices
318 :     isTriMat <- 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 : maechler 1238 ## For Csparse matrices
339 :     isTriC <- function(x, upper = NA) {
340 :     ## pretest: is it square?
341 :     d <- dim(x)
342 :     if(d[1] != d[2]) return(FALSE)
343 :     ## else
344 :     if(d[1] == 0) return(TRUE)
345 : maechler 1315 ni <- 1:d[2]
346 : maechler 1238 ## the row indices split according to column:
347 :     ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
348 :     lil <- unlist(lapply(ilist, length), use.names = FALSE)
349 :     if(any(lil == 0)) {
350 :     pos <- lil > 0
351 :     if(!any(pos)) ## matrix of all 0's
352 :     return(TRUE)
353 :     ilist <- ilist[pos]
354 :     ni <- ni[pos]
355 :     }
356 : maechler 1315 ni0 <- ni - 1:1 # '0-based ni'
357 : maechler 1238 if(is.na(upper)) {
358 : maechler 1315 if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
359 : maechler 1238 structure(TRUE, kind = "U")
360 : maechler 1315 else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
361 : maechler 1238 structure(TRUE, kind = "L")
362 :     else FALSE
363 :     } else if(upper) {
364 : maechler 1315 all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
365 : maechler 1238 } else { ## 'lower'
366 : maechler 1315 all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
367 : maechler 1238 }
368 :     }
369 :    
370 : maechler 1174 .is.diagonal <- function(object) {
371 :     d <- dim(object)
372 :     if(d[1] != (n <- d[2])) FALSE
373 :     else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)
374 :     }
375 : maechler 1253
376 :     diagU2N <- function(x)
377 :     {
378 :     ## Purpose: Transform a *unit diagonal* triangular matrix
379 :     ## into one with explicit diagonal entries '1'
380 :     xT <- as(x, "dgTMatrix")
381 :     ## leave it as T* - the caller can always coerce to C* if needed:
382 :     new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
383 :     Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
384 :     }
385 : maechler 1290
386 : maechler 1295 .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
387 :     x <- as(x, "dgCMatrix")
388 :     callGeneric()
389 :     }
390 :    
391 :     .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
392 :     x <- as(x, "dgTMatrix")
393 :     callGeneric()
394 :     }
395 :    
396 :    
397 : maechler 1290 ### Fast much simplified version of tapply()
398 :     tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
399 :     sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
400 :     }
401 :    
402 :     ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
403 :     ## tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
404 :     ## }
405 :    

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