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 1280 - (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 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 : maechler 1238 ## possibly dimnames(x) "==" NULL :
95 : maechler 956 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 : maechler 1238 justify = "none", right = TRUE)
101 : maechler 919 {
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 : maechler 1238 cf[row(cf) > col(cf)] <- "."
110 : maechler 919 else
111 : maechler 1238 cf[row(cf) < col(cf)] <- "."
112 : maechler 919 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 : maechler 1238 if(is(x, "triangularMatrix"))
123 :     prTriang(x, digits = digits)
124 :     else
125 :     print(as(x, "matrix"), digits = digits)
126 : maechler 919 }
127 :     else { ## d[1] > maxp / d[2] >= nr :
128 : maechler 1238 m <- as(x, "matrix")
129 : maechler 919 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 : maechler 1280 .Call(compressed_non_0_ij, x, isC)
152 : maechler 1226 }
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 : maechler 1238 ## but as 1-based indices
163 : maechler 1226 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 1270
190 :     ### Note: maybe, using
191 : maechler 1280 ### ---- xj <- .Call(Matrix_expand_pointers, x@p)
192 : maechler 1270 ### would be slightly more efficient than as( <dgC> , "dgTMatrix")
193 :     ### but really efficient would be to use only one .Call(.) for uniq(.) !
194 :     ### Try to do it particularly fast for the case where 'x' is already a 'uniq' <dgT>
195 :    
196 : maechler 925 if(is(x, "TsparseMatrix")) {
197 : maechler 1238 ## Purpose: produce a *unique* triplet representation:
198 :     ## by having (i,j) sorted and unique
199 :     ## -----------------------------------------------------------
200 :     ## The following is not quite efficient {but easy to program,
201 :     ## and both as() are based on C code
202 :     if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
203 :     else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
204 :     else stop("not implemented for class", class(x))
205 : maechler 919
206 : maechler 1226 } else x ## not 'gT' ; i.e. "uniquely" represented in any case
207 : maechler 919 }
208 :    
209 :     if(FALSE) ## try an "efficient" version
210 :     uniq_gT <- function(x)
211 :     {
212 :     ## Purpose: produce a *unique* triplet representation:
213 :     ## by having (i,j) sorted and unique
214 : maechler 1226 ## ------------------------------------------------------------------
215 : maechler 919 ## Arguments: a "gT" Matrix
216 :     stopifnot(is(x, "gTMatrix"))
217 :     if((n <- length(x@i)) == 0) return(x)
218 :     ii <- order(x@i, x@j)
219 :     if(any(ii != 1:n)) {
220 : maechler 1238 x@i <- x@i[ii]
221 :     x@j <- x@j[ii]
222 :     x@x <- x@x[ii]
223 : maechler 919 }
224 :     ij <- x@i + nrow(x) * x@j
225 :     if(any(dup <- duplicated(ij))) {
226 :    
227 :     }
228 :     ### We should use a .Call() based utility for this!
229 :    
230 :     }
231 :    
232 : maechler 946 t_geMatrix <- function(x) {
233 :     x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
234 :     x@Dim <- x@Dim[2:1]
235 :     x@Dimnames <- x@Dimnames[2:1]
236 :     ## FIXME: how to set factors?
237 :     x
238 :     }
239 :    
240 :     ## t( [dl]trMatrix ) and t( [dl]syMatrix ) :
241 :     t_trMatrix <- function(x) {
242 :     x@x <- as.vector(t(as(x, "matrix")))
243 :     x@Dim <- x@Dim[2:1]
244 :     x@Dimnames <- x@Dimnames[2:1]
245 :     x@uplo <- if (x@uplo == "U") "L" else "U"
246 :     # and keep x@diag
247 :     x
248 :     }
249 : maechler 956
250 :     fixupDense <- function(m, from) {
251 :     if(is(m, "triangularMatrix")) {
252 : maechler 1238 m@uplo <- from@uplo
253 :     m@diag <- from@diag
254 : maechler 956 } else if(is(m, "symmetricMatrix")) {
255 : maechler 1238 m@uplo <- from@uplo
256 : maechler 956 }
257 :     m
258 :     }
259 :    
260 :     ## -> ./ldenseMatrix.R :
261 :     l2d_Matrix <- function(from) {
262 :     stopifnot(is(from, "lMatrix"))
263 :     fixupDense(new(sub("^l", "d", class(from)),
264 : maechler 1238 x = as.double(from@x),
265 :     Dim = from@Dim, Dimnames = from@Dimnames),
266 :     from)
267 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
268 : maechler 956 }
269 :    
270 :     if(FALSE)# unused
271 :     l2d_meth <- function(x) {
272 :     cl <- class(x)
273 :     as(callGeneric(as(x, sub("^l", "d", cl))), cl)
274 :     }
275 :    
276 : maechler 1226 dClass2 <- function(dClass, kind = "l") {
277 :     ## Find "corresponding" class for a dMatrix;
278 :     # since pos.def. matrices have no pendant:
279 :     if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
280 :     else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
281 :     else sub("^d", kind, dClass)
282 :     }
283 :    
284 :     geClass <- function(x) {
285 :     if(is(x, "dMatrix")) "dgeMatrix"
286 :     else if(is(x, "lMatrix")) "lgeMatrix"
287 :     else stop("general Matrix class not yet implemented for",
288 :     class(x))
289 :     }
290 :    
291 : maechler 956 ## -> ./ddenseMatrix.R :
292 :     d2l_Matrix <- function(from) {
293 :     stopifnot(is(from, "dMatrix"))
294 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
295 : maechler 1238 Dim = from@Dim, Dimnames = from@Dimnames),
296 :     from)
297 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
298 : maechler 956 }
299 : maechler 973
300 :    
301 :     try_as <- function(x, classes, tryAnyway = FALSE) {
302 :     if(!tryAnyway && !is(x, "Matrix"))
303 : maechler 1238 return(x)
304 : maechler 973 ## else
305 :     ok <- canCoerce(x, classes[1])
306 :     while(!ok && length(classes <- classes[-1])) {
307 : maechler 1238 ok <- canCoerce(x, classes[1])
308 : maechler 973 }
309 :     if(ok) as(x, classes[1]) else x
310 :     }
311 :    
312 : maechler 1174 if(paste(R.version$major, R.version$minor, sep=".") < "2.3")
313 :     ## This will be in R 2.3.0
314 : maechler 973 canCoerce <- function(object, Class) {
315 :     ## Purpose: test if 'object' is coercable to 'Class', i.e.,
316 :     ## as(object, Class) will {typically} work
317 :     ## ----------------------------------------------------------------------
318 :     ## Author: John Chambers, Date: 6 Oct 2005
319 :     is(object, Class) ||
320 :     !is.null(selectMethod("coerce", c(class(object), Class),
321 :     optional = TRUE,
322 :     useInherited = c(from = TRUE, to = FALSE)))
323 :     }
324 :    
325 : maechler 1238 ## For *dense* matrices
326 :     isTriMat <- function(object, upper = NA) {
327 : maechler 1174 ## pretest: is it square?
328 :     d <- dim(object)
329 :     if(d[1] != d[2]) return(FALSE)
330 :     ## else slower test
331 :     if(!is.matrix(object))
332 : maechler 1226 object <- as(object,"matrix")
333 : maechler 1174 ## == 0 even works for logical & complex:
334 : maechler 1226 if(is.na(upper)) {
335 :     if(all(object[lower.tri(object)] == 0))
336 :     structure(TRUE, kind = "U")
337 :     else if(all(object[upper.tri(object)] == 0))
338 :     structure(TRUE, kind = "L")
339 :     else FALSE
340 :     } else if(upper)
341 :     all(object[lower.tri(object)] == 0)
342 :     else ## upper is FALSE
343 :     all(object[upper.tri(object)] == 0)
344 : maechler 1174 }
345 :    
346 : maechler 1238 ## For Csparse matrices
347 :     isTriC <- function(x, upper = NA) {
348 :     ## pretest: is it square?
349 :     d <- dim(x)
350 :     if(d[1] != d[2]) return(FALSE)
351 :     ## else
352 :     if(d[1] == 0) return(TRUE)
353 :     ni <- 1:d[1]
354 :     ## the row indices split according to column:
355 :     ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
356 :     lil <- unlist(lapply(ilist, length), use.names = FALSE)
357 :     if(any(lil == 0)) {
358 :     pos <- lil > 0
359 :     if(!any(pos)) ## matrix of all 0's
360 :     return(TRUE)
361 :     ilist <- ilist[pos]
362 :     ni <- ni[pos]
363 :     }
364 :     if(is.na(upper)) {
365 :     if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni))
366 :     structure(TRUE, kind = "U")
367 :     else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni))
368 :     structure(TRUE, kind = "L")
369 :     else FALSE
370 :     } else if(upper) {
371 :     all(sapply(ilist, max, USE.NAMES = FALSE) <= ni)
372 :     } else { ## 'lower'
373 :     all(sapply(ilist, min, USE.NAMES = FALSE) >= ni)
374 :     }
375 :     }
376 :    
377 : maechler 1174 .is.diagonal <- function(object) {
378 :     d <- dim(object)
379 :     if(d[1] != (n <- d[2])) FALSE
380 :     else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)
381 :     }
382 : maechler 1253
383 :     diagU2N <- function(x)
384 :     {
385 :     ## Purpose: Transform a *unit diagonal* triangular matrix
386 :     ## into one with explicit diagonal entries '1'
387 :     xT <- as(x, "dgTMatrix")
388 :     ## leave it as T* - the caller can always coerce to C* if needed:
389 :     new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
390 :     Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
391 :     }

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