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

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