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 1226 - (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 956 emptyColnames <- function(x)
86 :     {
87 :     ## Useful for compact printing of (parts) of sparse matrices
88 :     ## possibly dimnames(x) "==" NULL :
89 :     dimnames(x) <- list(dimnames(x)[[1]], rep("", dim(x)[2]))
90 :     x
91 :     }
92 : maechler 908
93 : maechler 919 prTriang <- function(x, digits = getOption("digits"),
94 :     justify = "none", right = TRUE)
95 :     {
96 :     ## modeled along stats:::print.dist
97 :     diag <- TRUE
98 :     upper <- x@uplo == "U"
99 :    
100 :     m <- as(x, "matrix")
101 :     cf <- format(m, digits = digits, justify = justify)
102 :     if(upper)
103 :     cf[row(cf) > col(cf)] <- "."
104 :     else
105 :     cf[row(cf) < col(cf)] <- "."
106 :     print(cf, quote = FALSE, right = right)
107 :     invisible(x)
108 :     }
109 :    
110 :     prMatrix <- function(x, digits = getOption("digits")) {
111 :     d <- dim(x)
112 :     cl <- class(x)
113 :     cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
114 :     maxp <- getOption("max.print")
115 :     if(prod(d) <= maxp) {
116 :     if(is(x, "triangularMatrix"))
117 :     prTriang(x, digits = digits)
118 :     else
119 :     print(as(x, "matrix"), digits = digits)
120 :     }
121 :     else { ## d[1] > maxp / d[2] >= nr :
122 :     m <- as(x, "matrix")
123 :     nr <- maxp %/% d[2]
124 :     n2 <- ceiling(nr / 2)
125 :     print(head(m, max(1, n2)))
126 :     cat("\n ..........\n\n")
127 :     print(tail(m, max(1, nr - n2)))
128 :     }
129 :     ## DEBUG: cat("str(.):\n") ; str(x)
130 :     invisible(x)# as print() S3 methods do
131 :     }
132 :    
133 :     ## For sparseness handling
134 :     non0ind <- function(x) {
135 :     if(is.numeric(x))
136 : maechler 1226 return(if((n <- length(x))) (0:(n-1))[x != 0] else integer(0))
137 :     ## else
138 : maechler 919 stopifnot(is(x, "sparseMatrix"))
139 : maechler 1226 ## return a 2-column (i,j) matrix of
140 :     ## 0-based indices of non-zero entries :
141 : maechler 925 if(is(x, "TsparseMatrix"))
142 : maechler 954 return(unique(cbind(x@i,x@j)))
143 : maechler 1226 ## else:
144 :     isC <- any("i" == slotNames(x))# is Csparse (not Rsparse)
145 :     .Call("compressed_non_0_ij", x, isC, PACKAGE = "Matrix")
146 :     }
147 : maechler 954
148 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings:
149 :     ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
150 :     encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
151 :     decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
152 :    
153 :     complementInd <- function(ij, dim)
154 :     {
155 :     ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
156 :     ## but as 1-based indices
157 :     n <- prod(dim)
158 :     if(n == 0) return(integer(0))
159 :     ii <- 1:n
160 :     ii[-(1 + encodeInd(ij, nr = dim[1]))]
161 : maechler 919 }
162 :    
163 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
164 :    
165 :     intersectInd <- function(ij1, ij2, nrow) {
166 :     ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
167 :     ## return only the *common* entries
168 :     decodeInd(intersect(encodeInd(ij1, nrow),
169 :     encodeInd(ij2, nrow)), nrow)
170 :     }
171 :    
172 :     WhichintersectInd <- function(ij1, ij2, nrow) {
173 :     ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
174 :     ## find *where* common entries are in ij1 & ij2
175 :     m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
176 :     ni <- !is.na(m1)
177 :     list(which(ni), m1[ni])
178 :     }
179 :    
180 :    
181 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R !
182 : maechler 919 uniq <- function(x) {
183 : maechler 925 if(is(x, "TsparseMatrix")) {
184 : maechler 919 ## Purpose: produce a *unique* triplet representation:
185 :     ## by having (i,j) sorted and unique
186 :     ## -----------------------------------------------------------
187 : maechler 1226 ## The following is not quite efficient {but easy to program,
188 :     ## and both as() are based on C code
189 : maechler 919 if(is(x, "dgTMatrix")) as(as(x, "dgCMatrix"), "dgTMatrix")
190 :     else if(is(x, "lgTMatrix")) as(as(x, "lgCMatrix"), "lgTMatrix")
191 :     else stop("not implemented for class", class(x))
192 :    
193 : maechler 1226 } else x ## not 'gT' ; i.e. "uniquely" represented in any case
194 : maechler 919 }
195 :    
196 :     if(FALSE) ## try an "efficient" version
197 :     uniq_gT <- function(x)
198 :     {
199 :     ## Purpose: produce a *unique* triplet representation:
200 :     ## by having (i,j) sorted and unique
201 : maechler 1226 ## ------------------------------------------------------------------
202 : maechler 919 ## Arguments: a "gT" Matrix
203 :     stopifnot(is(x, "gTMatrix"))
204 :     if((n <- length(x@i)) == 0) return(x)
205 :     ii <- order(x@i, x@j)
206 :     if(any(ii != 1:n)) {
207 :     x@i <- x@i[ii]
208 :     x@j <- x@j[ii]
209 :     x@x <- x@x[ii]
210 :     }
211 :     ij <- x@i + nrow(x) * x@j
212 :     if(any(dup <- duplicated(ij))) {
213 :    
214 :     }
215 :     ### We should use a .Call() based utility for this!
216 :    
217 :     }
218 :    
219 : maechler 946 t_geMatrix <- function(x) {
220 :     x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
221 :     x@Dim <- x@Dim[2:1]
222 :     x@Dimnames <- x@Dimnames[2:1]
223 :     ## FIXME: how to set factors?
224 :     x
225 :     }
226 :    
227 :     ## t( [dl]trMatrix ) and t( [dl]syMatrix ) :
228 :     t_trMatrix <- function(x) {
229 :     x@x <- as.vector(t(as(x, "matrix")))
230 :     x@Dim <- x@Dim[2:1]
231 :     x@Dimnames <- x@Dimnames[2:1]
232 :     x@uplo <- if (x@uplo == "U") "L" else "U"
233 :     # and keep x@diag
234 :     x
235 :     }
236 : maechler 956
237 :     fixupDense <- function(m, from) {
238 :     if(is(m, "triangularMatrix")) {
239 :     m@uplo <- from@uplo
240 :     m@diag <- from@diag
241 :     } else if(is(m, "symmetricMatrix")) {
242 :     m@uplo <- from@uplo
243 :     }
244 :     m
245 :     }
246 :    
247 :     ## -> ./ldenseMatrix.R :
248 :     l2d_Matrix <- function(from) {
249 :     stopifnot(is(from, "lMatrix"))
250 :     fixupDense(new(sub("^l", "d", class(from)),
251 :     x = as.double(from@x),
252 : maechler 1198 Dim = from@Dim, Dimnames = from@Dimnames),
253 : maechler 956 from)
254 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
255 : maechler 956 }
256 :    
257 :     if(FALSE)# unused
258 :     l2d_meth <- function(x) {
259 :     cl <- class(x)
260 :     as(callGeneric(as(x, sub("^l", "d", cl))), cl)
261 :     }
262 :    
263 : maechler 1226 dClass2 <- function(dClass, kind = "l") {
264 :     ## Find "corresponding" class for a dMatrix;
265 :     # since pos.def. matrices have no pendant:
266 :     if(dClass == "dpoMatrix") paste(kind,"syMatrix", sep='')
267 :     else if(dClass == "dppMatrix") paste(kind,"spMatrix", sep='')
268 :     else sub("^d", kind, dClass)
269 :     }
270 :    
271 :     geClass <- function(x) {
272 :     if(is(x, "dMatrix")) "dgeMatrix"
273 :     else if(is(x, "lMatrix")) "lgeMatrix"
274 :     else stop("general Matrix class not yet implemented for",
275 :     class(x))
276 :     }
277 :    
278 : maechler 956 ## -> ./ddenseMatrix.R :
279 :     d2l_Matrix <- function(from) {
280 :     stopifnot(is(from, "dMatrix"))
281 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
282 : maechler 1198 Dim = from@Dim, Dimnames = from@Dimnames),
283 : maechler 956 from)
284 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
285 : maechler 956 }
286 : maechler 973
287 :    
288 :     try_as <- function(x, classes, tryAnyway = FALSE) {
289 :     if(!tryAnyway && !is(x, "Matrix"))
290 :     return(x)
291 :     ## else
292 :     ok <- canCoerce(x, classes[1])
293 :     while(!ok && length(classes <- classes[-1])) {
294 :     ok <- canCoerce(x, classes[1])
295 :     }
296 :     if(ok) as(x, classes[1]) else x
297 :     }
298 :    
299 : maechler 1174 if(paste(R.version$major, R.version$minor, sep=".") < "2.3")
300 :     ## This will be in R 2.3.0
301 : maechler 973 canCoerce <- function(object, Class) {
302 :     ## Purpose: test if 'object' is coercable to 'Class', i.e.,
303 :     ## as(object, Class) will {typically} work
304 :     ## ----------------------------------------------------------------------
305 :     ## Author: John Chambers, Date: 6 Oct 2005
306 :     is(object, Class) ||
307 :     !is.null(selectMethod("coerce", c(class(object), Class),
308 :     optional = TRUE,
309 :     useInherited = c(from = TRUE, to = FALSE)))
310 :     }
311 :    
312 : maechler 1226 .is.triangular <- function(object, upper = NA) {
313 : maechler 1174 ## pretest: is it square?
314 :     d <- dim(object)
315 :     if(d[1] != d[2]) return(FALSE)
316 :     ## else slower test
317 :     if(!is.matrix(object))
318 : maechler 1226 object <- as(object,"matrix")
319 : maechler 1174 ## == 0 even works for logical & complex:
320 : maechler 1226 if(is.na(upper)) {
321 :     if(all(object[lower.tri(object)] == 0))
322 :     structure(TRUE, kind = "U")
323 :     else if(all(object[upper.tri(object)] == 0))
324 :     structure(TRUE, kind = "L")
325 :     else FALSE
326 :     } else if(upper)
327 :     all(object[lower.tri(object)] == 0)
328 :     else ## upper is FALSE
329 :     all(object[upper.tri(object)] == 0)
330 : maechler 1174 }
331 :    
332 :     .is.diagonal <- function(object) {
333 :     d <- dim(object)
334 :     if(d[1] != (n <- d[2])) FALSE
335 :     else all(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)] == 0)
336 :     }

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