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 1659 - (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 : maechler 1654 stop(gettextf('not-yet-implemented method for %s(<%s>).\n ->> Ask the package authors to implement the missing feature.', fun, cl),
26 : maechler 1238 call. = FALSE)
27 : maechler 949 }
28 :     .bail.out.2 <- function(fun, cl1, cl2) {
29 : maechler 1654 stop(gettextf('not-yet-implemented method for %s(<%s>, <%s>).\n ->> Ask the package authors to implement the missing feature.',
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 1592 ### TODO: write in C and port to base (or 'utils') R
125 :     indTri <- function(n, upper = TRUE) {
126 :     ## == which(upper.tri(diag(n)) or
127 :     ## which(lower.tri(diag(n)) -- but much more efficiently for largish 'n'
128 :     stopifnot(length(n) == 1, n == (n. <- as.integer(n)), (n <- n.) >= 0)
129 :     if(n <= 2)
130 :     return(if(n == 2) as.integer(if(upper) n+1 else n) else integer(0))
131 :     ## First, compute the 'diff(.)' fast. Use integers
132 :     one <- 1:1 ; two <- 2:2
133 :     n1 <- n - one
134 :     n2 <- n1 - one
135 :     r <- rep.int(one, n*n1/two - one)
136 :     r[cumsum(if(upper) 1:n2 else c(n1, if(n >= 4) n2:two))] <- if(upper) n:3 else 3:n
137 :     ## now have "dliu" difference; revert to "liu":
138 :     cumsum(c(if(upper) n+one else two, r))
139 :     }
140 :    
141 :    
142 : maechler 919 prTriang <- function(x, digits = getOption("digits"),
143 : maechler 1389 maxp = getOption("max.print"),
144 : maechler 1238 justify = "none", right = TRUE)
145 : maechler 919 {
146 :     ## modeled along stats:::print.dist
147 :     upper <- x@uplo == "U"
148 :    
149 :     m <- as(x, "matrix")
150 :     cf <- format(m, digits = digits, justify = justify)
151 :     if(upper)
152 : maechler 1238 cf[row(cf) > col(cf)] <- "."
153 : maechler 919 else
154 : maechler 1238 cf[row(cf) < col(cf)] <- "."
155 : maechler 1472 if(.isR_24)
156 :     print(cf, quote = FALSE, right = right, max = maxp)
157 :     else print(cf, quote = FALSE, right = right)
158 : maechler 919 invisible(x)
159 :     }
160 :    
161 : maechler 1389 prMatrix <- function(x, digits = getOption("digits"),
162 :     maxp = getOption("max.print")) {
163 : maechler 919 d <- dim(x)
164 :     cl <- class(x)
165 :     cat(sprintf('%d x %d Matrix of class "%s"\n', d[1], d[2], cl))
166 :     if(prod(d) <= maxp) {
167 : maechler 1238 if(is(x, "triangularMatrix"))
168 : maechler 1472 prTriang(x, digits = digits, maxp = maxp)
169 :     else {
170 :     if(.isR_24)
171 :     print(as(x, "matrix"), digits = digits, max = maxp)
172 :     else print(as(x, "matrix"), digits = digits)
173 :     }
174 : maechler 919 }
175 :     else { ## d[1] > maxp / d[2] >= nr :
176 : maechler 1238 m <- as(x, "matrix")
177 : maechler 919 nr <- maxp %/% d[2]
178 :     n2 <- ceiling(nr / 2)
179 :     print(head(m, max(1, n2)))
180 :     cat("\n ..........\n\n")
181 :     print(tail(m, max(1, nr - n2)))
182 :     }
183 :     ## DEBUG: cat("str(.):\n") ; str(x)
184 :     invisible(x)# as print() S3 methods do
185 :     }
186 :    
187 : maechler 1575 nonFALSE <- function(x) {
188 :     ## typically used for lMatrices: (TRUE,NA,FALSE) |-> (TRUE,FALSE)
189 :     if(any(ix <- is.na(x))) x[ix] <- TRUE
190 :     x
191 :     }
192 :    
193 :     nz.NA <- function(x, na.value) {
194 :     ## Non-Zeros of x
195 :     ## na.value: TRUE: NA's give TRUE, they are not 0
196 :     ## NA: NA's are not known ==> result := NA
197 :     ## FALSE: NA's give FALSE, could be 0
198 :     stopifnot(is.logical(na.value) && length(na.value) == 1)
199 :     if(is.na(na.value)) x != 0
200 :     else if(na.value) isN0(x)
201 :     else x != 0 & !is.na(x)
202 :     }
203 :    
204 : maechler 1592 ## Number of non-zeros :
205 :     ## FIXME? -- make this into a generic function (?)
206 : maechler 1575 nnzero <- function(x, na.counted = NA) {
207 :     ## na.counted: TRUE: NA's are counted, they are not 0
208 :     ## NA: NA's are not known (0 or not) ==> result := NA
209 :     ## FALSE: NA's are omitted before counting
210 : maechler 1548 cl <- class(x)
211 :     if(!extends(cl, "Matrix"))
212 : maechler 1575 sum(nz.NA(x, na.counted))
213 : maechler 1548 else if(extends(cl, "sparseMatrix"))
214 :     ## NOTA BENE: The number of *structural* non-zeros {could have other '0'}!
215 :     switch(.sp.class(cl),
216 :     "CsparseMatrix" = length(x@i),
217 :     "TsparseMatrix" = length(x@i),
218 :     "RsparseMatrix" = length(x@j))
219 :     else ## denseMatrix
220 : maechler 1575 sum(nz.NA(as_geClass(x, cl)@x, na.counted))
221 : maechler 1548 }
222 :    
223 : maechler 919 ## For sparseness handling
224 : maechler 1467 ## return a 2-column (i,j) matrix of
225 :     ## 0-based indices of non-zero entries :
226 : maechler 919 non0ind <- function(x) {
227 : maechler 1467
228 : maechler 919 if(is.numeric(x))
229 : maechler 1575 return(if((n <- length(x))) (0:(n-1))[isN0(x)] else integer(0))
230 : maechler 1226 ## else
231 : maechler 919 stopifnot(is(x, "sparseMatrix"))
232 : maechler 1319 non0.i <- function(M) {
233 :     if(is(M, "TsparseMatrix"))
234 :     return(unique(cbind(M@i,M@j)))
235 :     if(is(M, "pMatrix"))
236 : maechler 1654 return(cbind(seq_len(nrow(M)), M@perm) - 1:1)
237 : maechler 1319 ## else:
238 :     isC <- any("i" == slotNames(M)) # is Csparse (not Rsparse)
239 :     .Call(compressed_non_0_ij, M, isC)
240 :     }
241 :    
242 :     if(is(x, "symmetricMatrix")) { # also get "other" triangle
243 :     ij <- non0.i(x)
244 :     notdiag <- ij[,1] != ij[,2]# but not the diagonals again
245 :     rbind(ij, ij[notdiag, 2:1])
246 :     }
247 : maechler 1389 else if(is(x, "triangularMatrix")) { # check for "U" diag
248 :     if(x@diag == "U") {
249 : maechler 1654 i <- seq_len(dim(x)[1]) - 1:1
250 : maechler 1389 rbind(non0.i(x), cbind(i,i))
251 :     } else non0.i(x)
252 :     }
253 : maechler 1319 else
254 :     non0.i(x)
255 : maechler 1226 }
256 : maechler 954
257 : maechler 1226 ## nr= nrow: since i in {0,1,.., nrow-1} these are 1:1 "decimal" encodings:
258 :     ## Further, these map to and from the usual "Fortran-indexing" (but 0-based)
259 :     encodeInd <- function(ij, nr) ij[,1] + ij[,2] * nr
260 :     decodeInd <- function(code, nr) cbind(code %% nr, code %/% nr)
261 :    
262 :     complementInd <- function(ij, dim)
263 :     {
264 :     ## Purpose: Compute the complement of the 2-column 0-based ij-matrix
265 : maechler 1238 ## but as 1-based indices
266 : maechler 1226 n <- prod(dim)
267 :     if(n == 0) return(integer(0))
268 :     ii <- 1:n
269 :     ii[-(1 + encodeInd(ij, nr = dim[1]))]
270 : maechler 919 }
271 :    
272 : maechler 1226 unionInd <- function(ij1, ij2) unique(rbind(ij1, ij2))
273 :    
274 :     intersectInd <- function(ij1, ij2, nrow) {
275 :     ## from 2-column (i,j) matrices where i in {0,.., nrow-1},
276 :     ## return only the *common* entries
277 :     decodeInd(intersect(encodeInd(ij1, nrow),
278 :     encodeInd(ij2, nrow)), nrow)
279 :     }
280 :    
281 :     WhichintersectInd <- function(ij1, ij2, nrow) {
282 :     ## from 2-column (i,j) matrices where i \in {0,.., nrow-1},
283 :     ## find *where* common entries are in ij1 & ij2
284 :     m1 <- match(encodeInd(ij1, nrow), encodeInd(ij2, nrow))
285 :     ni <- !is.na(m1)
286 :     list(which(ni), m1[ni])
287 :     }
288 :    
289 :    
290 : maechler 925 ### There is a test on this in ../tests/dgTMatrix.R !
291 : maechler 1270
292 : maechler 1331 uniqTsparse <- function(x, class.x = c(class(x))) {
293 :     ## Purpose: produce a *unique* triplet representation:
294 :     ## by having (i,j) sorted and unique
295 :     ## -----------------------------------------------------------
296 :     ## The following is not quite efficient {but easy to program,
297 :     ## and as() are based on C code (all of them?)
298 :     ##
299 :     ## FIXME: Do it fast for the case where 'x' is already 'uniq'
300 : maechler 1270
301 : maechler 1331 switch(class.x,
302 :     "dgTMatrix" = as(as(x, "dgCMatrix"), "dgTMatrix"),
303 :     "dsTMatrix" = as(as(x, "dsCMatrix"), "dsTMatrix"),
304 :     "dtTMatrix" = as(as(x, "dtCMatrix"), "dtTMatrix"),
305 :     ## do we need this for "logical" ones, there's no sum() there!
306 :     "lgTMatrix" = as(as(x, "lgCMatrix"), "lgTMatrix"),
307 :     "lsTMatrix" = as(as(x, "lsCMatrix"), "lsTMatrix"),
308 :     "ltTMatrix" = as(as(x, "ltCMatrix"), "ltTMatrix"),
309 : maechler 1548 ## do we need this for "logical" ones, there's no sum() there!
310 :     "ngTMatrix" = as(as(x, "ngCMatrix"), "ngTMatrix"),
311 :     "nsTMatrix" = as(as(x, "nsCMatrix"), "nsTMatrix"),
312 :     "ntTMatrix" = as(as(x, "ntCMatrix"), "ntTMatrix"),
313 : maechler 1331 ## otherwise:
314 : maechler 1472 stop("not yet implemented for class ", class.x))
315 : maechler 1331 }
316 : maechler 919
317 : maechler 1331 ## Note: maybe, using
318 :     ## ---- xj <- .Call(Matrix_expand_pointers, x@p)
319 :     ## would be slightly more efficient than as( <dgC> , "dgTMatrix")
320 :     ## but really efficient would be to use only one .Call(.) for uniq(.) !
321 :    
322 : maechler 1654 drop0 <- function(x, clx = c(class(x))) {
323 :     ## FIXME: Csparse_drop should do this (not losing symm./triang.):
324 :     ## Careful: 'Csparse_drop' also drops triangularity,...
325 :     ## .Call(Csparse_drop, as_CspClass(x, clx), 0)
326 :     as_CspClass(.Call(Csparse_drop, as_CspClass(x, clx), 0.),
327 :     clx)
328 :     }
329 :    
330 : maechler 1331 uniq <- function(x) {
331 : maechler 1654 if(is(x, "TsparseMatrix")) uniqTsparse(x) else
332 :     if(is(x, "sparseMatrix")) drop0(x) else x
333 : maechler 919 }
334 :    
335 : maechler 1472 asTuniq <- function(x) {
336 :     if(is(x, "TsparseMatrix")) uniqTsparse(x) else as(x,"TsparseMatrix")
337 :     }
338 :    
339 : maechler 919 if(FALSE) ## try an "efficient" version
340 :     uniq_gT <- function(x)
341 :     {
342 :     ## Purpose: produce a *unique* triplet representation:
343 :     ## by having (i,j) sorted and unique
344 : maechler 1226 ## ------------------------------------------------------------------
345 : maechler 919 ## Arguments: a "gT" Matrix
346 :     stopifnot(is(x, "gTMatrix"))
347 :     if((n <- length(x@i)) == 0) return(x)
348 :     ii <- order(x@i, x@j)
349 :     if(any(ii != 1:n)) {
350 : maechler 1238 x@i <- x@i[ii]
351 :     x@j <- x@j[ii]
352 :     x@x <- x@x[ii]
353 : maechler 919 }
354 :     ij <- x@i + nrow(x) * x@j
355 :     if(any(dup <- duplicated(ij))) {
356 :    
357 :     }
358 :     ### We should use a .Call() based utility for this!
359 :    
360 :     }
361 :    
362 : maechler 946 t_geMatrix <- function(x) {
363 :     x@x <- as.vector(t(array(x@x, dim = x@Dim))) # no dimnames here
364 :     x@Dim <- x@Dim[2:1]
365 :     x@Dimnames <- x@Dimnames[2:1]
366 :     ## FIXME: how to set factors?
367 :     x
368 :     }
369 :    
370 :     ## t( [dl]trMatrix ) and t( [dl]syMatrix ) :
371 :     t_trMatrix <- function(x) {
372 :     x@x <- as.vector(t(as(x, "matrix")))
373 :     x@Dim <- x@Dim[2:1]
374 :     x@Dimnames <- x@Dimnames[2:1]
375 :     x@uplo <- if (x@uplo == "U") "L" else "U"
376 :     # and keep x@diag
377 :     x
378 :     }
379 : maechler 956
380 :     fixupDense <- function(m, from) {
381 :     if(is(m, "triangularMatrix")) {
382 : maechler 1238 m@uplo <- from@uplo
383 :     m@diag <- from@diag
384 : maechler 956 } else if(is(m, "symmetricMatrix")) {
385 : maechler 1238 m@uplo <- from@uplo
386 : maechler 956 }
387 :     m
388 :     }
389 :    
390 :     ## -> ./ldenseMatrix.R :
391 :     l2d_Matrix <- function(from) {
392 :     stopifnot(is(from, "lMatrix"))
393 :     fixupDense(new(sub("^l", "d", class(from)),
394 : maechler 1238 x = as.double(from@x),
395 :     Dim = from@Dim, Dimnames = from@Dimnames),
396 :     from)
397 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
398 : maechler 956 }
399 :    
400 : maechler 1548 ## -> ./ndenseMatrix.R :
401 :     n2d_Matrix <- function(from) {
402 :     stopifnot(is(from, "nMatrix"))
403 :     fixupDense(new(sub("^n", "d", class(from)),
404 :     x = as.double(from@x),
405 :     Dim = from@Dim, Dimnames = from@Dimnames),
406 :     from)
407 :     ## FIXME: treat 'factors' smartly {not for triangular!}
408 :     }
409 :     n2l_spMatrix <- function(from) {
410 :     stopifnot(is(from, "nMatrix"))
411 :     new(sub("^n", "l", class(from)),
412 :     ##x = as.double(from@x),
413 :     Dim = from@Dim, Dimnames = from@Dimnames)
414 :     }
415 :    
416 : maechler 956 if(FALSE)# unused
417 :     l2d_meth <- function(x) {
418 :     cl <- class(x)
419 :     as(callGeneric(as(x, sub("^l", "d", cl))), cl)
420 :     }
421 :    
422 : maechler 1548 ## return "d" or "l" or "n" or "z"
423 : maechler 1331 .M.kind <- function(x, clx = class(x)) {
424 :     if(is.matrix(x)) { ## 'old style matrix'
425 :     if (is.numeric(x)) "d"
426 : maechler 1548 else if(is.logical(x)) "l" ## FIXME ? "n" if no NA ??
427 : maechler 1331 else if(is.complex(x)) "z"
428 :     else stop("not yet implemented for matrix w/ typeof ", typeof(x))
429 :     }
430 :     else if(extends(clx, "dMatrix")) "d"
431 : maechler 1548 else if(extends(clx, "nMatrix")) "n"
432 : maechler 1331 else if(extends(clx, "lMatrix")) "l"
433 :     else if(extends(clx, "zMatrix")) "z"
434 : maechler 1548 else if(extends(clx, "pMatrix")) "n" # permutation -> pattern
435 : maechler 1331 else stop(" not yet be implemented for ", clx)
436 :     }
437 :    
438 :     .M.shape <- function(x, clx = class(x)) {
439 :     if(is.matrix(x)) { ## 'old style matrix'
440 :     if (isDiagonal (x)) "d"
441 :     else if(isTriangular(x)) "t"
442 :     else if(isSymmetric (x)) "s"
443 :     else "g" # general
444 :     }
445 :     else if(extends(clx, "diagonalMatrix")) "d"
446 :     else if(extends(clx, "triangularMatrix"))"t"
447 :     else if(extends(clx, "symmetricMatrix")) "s"
448 :     else "g"
449 :     }
450 :    
451 :    
452 : maechler 1329 class2 <- function(cl, kind = "l", do.sub = TRUE) {
453 :     ## Find "corresponding" class; since pos.def. matrices have no pendant:
454 :     if (cl == "dpoMatrix") paste(kind, "syMatrix", sep='')
455 :     else if(cl == "dppMatrix") paste(kind, "spMatrix", sep='')
456 :     else if(do.sub) sub("^d", kind, cl)
457 :     else cl
458 : maechler 1226 }
459 :    
460 :     geClass <- function(x) {
461 : maechler 1331 if (is(x, "dMatrix")) "dgeMatrix"
462 : maechler 1226 else if(is(x, "lMatrix")) "lgeMatrix"
463 : maechler 1548 else if(is(x, "nMatrix")) "ngeMatrix"
464 : maechler 1331 else if(is(x, "zMatrix")) "zgeMatrix"
465 : maechler 1329 else stop("general Matrix class not yet implemented for ",
466 : maechler 1226 class(x))
467 :     }
468 : maechler 1331
469 :     .dense.prefixes <- c("d" = "di",
470 :     "t" = "tr",
471 :     "s" = "sy",
472 :     "g" = "ge")
473 :    
474 : maechler 1349 .sparse.prefixes <- c("d" = "t", ## map diagonal to triangular
475 :     "t" = "t",
476 :     "s" = "s",
477 :     "g" = "g")
478 : maechler 1331
479 : maechler 1329 ## Used, e.g. after subsetting: Try to use specific class -- if feasible :
480 : maechler 1331 as_dense <- function(x) {
481 :     as(x, paste(.M.kind(x), .dense.prefixes[.M.shape(x)], "Matrix", sep=''))
482 :     }
483 :    
484 : maechler 1548 .sp.class <- function(x) { ## find and return the "sparseness class"
485 :     if(!is.character(x)) x <- class(x)
486 :     for(cl in paste(c("C","T","R"), "sparseMatrix", sep=''))
487 :     if(extends(x, cl))
488 :     return(cl)
489 :     ## else (should rarely happen)
490 :     as.character(NA)
491 :     }
492 :    
493 : maechler 1331 as_Csparse <- function(x) {
494 : maechler 1349 as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "CMatrix", sep=''))
495 : maechler 1331 }
496 : maechler 1349 as_Rsparse <- function(x) {
497 :     as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "RMatrix", sep=''))
498 :     }
499 :     as_Tsparse <- function(x) {
500 :     as(x, paste(.M.kind(x), .sparse.prefixes[.M.shape(x)], "TMatrix", sep=''))
501 :     }
502 : maechler 1331
503 : maechler 1329 as_geClass <- function(x, cl) {
504 : maechler 1331 if (extends(cl, "diagonalMatrix") && isDiagonal(x))
505 :     as(x, cl)
506 : maechler 1357 else if(extends(cl, "symmetricMatrix") && isSymmetric(x)) {
507 :     kind <- .M.kind(x)
508 : maechler 1331 as(x, class2(cl, kind, do.sub= kind != "d"))
509 : maechler 1357 } else if(extends(cl, "triangularMatrix") && isTriangular(x))
510 : maechler 1331 as(x, cl)
511 :     else
512 :     as(x, paste(.M.kind(x), "geMatrix", sep=''))
513 :     }
514 : maechler 1226
515 : maechler 1331 as_CspClass <- function(x, cl) {
516 : maechler 1654 if (## diagonal is *not* sparse:
517 :     ##(extends(cl, "diagonalMatrix") && isDiagonal(x)) ||
518 :     (extends(cl, "symmetricMatrix") && isSymmetric(x)) ||
519 : maechler 1331 (extends(cl, "triangularMatrix")&& isTriangular(x)))
520 :     as(x, cl)
521 : maechler 1654 else if(is(x, "CsparseMatrix")) x
522 : maechler 1331 else as(x, paste(.M.kind(x), "gCMatrix", sep=''))
523 : maechler 1329 }
524 :    
525 : maechler 1331
526 : maechler 956 ## -> ./ddenseMatrix.R :
527 :     d2l_Matrix <- function(from) {
528 :     stopifnot(is(from, "dMatrix"))
529 : maechler 1226 fixupDense(new(sub("^d", "l", class(from)), # no need for dClass2 here
530 : maechler 1238 Dim = from@Dim, Dimnames = from@Dimnames),
531 :     from)
532 : maechler 1198 ## FIXME: treat 'factors' smartly {not for triangular!}
533 : maechler 956 }
534 : maechler 973
535 :    
536 :     try_as <- function(x, classes, tryAnyway = FALSE) {
537 :     if(!tryAnyway && !is(x, "Matrix"))
538 : maechler 1238 return(x)
539 : maechler 973 ## else
540 :     ok <- canCoerce(x, classes[1])
541 :     while(!ok && length(classes <- classes[-1])) {
542 : maechler 1238 ok <- canCoerce(x, classes[1])
543 : maechler 973 }
544 :     if(ok) as(x, classes[1]) else x
545 :     }
546 :    
547 :    
548 : maechler 1238 ## For *dense* matrices
549 :     isTriMat <- function(object, upper = NA) {
550 : maechler 1174 ## pretest: is it square?
551 :     d <- dim(object)
552 :     if(d[1] != d[2]) return(FALSE)
553 :     ## else slower test
554 :     if(!is.matrix(object))
555 : maechler 1226 object <- as(object,"matrix")
556 :     if(is.na(upper)) {
557 : maechler 1467 if(all0(object[lower.tri(object)]))
558 : maechler 1226 structure(TRUE, kind = "U")
559 : maechler 1467 else if(all0(object[upper.tri(object)]))
560 : maechler 1226 structure(TRUE, kind = "L")
561 :     else FALSE
562 :     } else if(upper)
563 : maechler 1467 all0(object[lower.tri(object)])
564 : maechler 1226 else ## upper is FALSE
565 : maechler 1467 all0(object[upper.tri(object)])
566 : maechler 1174 }
567 :    
568 : maechler 1238 ## For Csparse matrices
569 :     isTriC <- function(x, upper = NA) {
570 :     ## pretest: is it square?
571 :     d <- dim(x)
572 :     if(d[1] != d[2]) return(FALSE)
573 :     ## else
574 :     if(d[1] == 0) return(TRUE)
575 : maechler 1315 ni <- 1:d[2]
576 : maechler 1238 ## the row indices split according to column:
577 :     ilist <- split(x@i, factor(rep.int(ni, diff(x@p)), levels= ni))
578 :     lil <- unlist(lapply(ilist, length), use.names = FALSE)
579 :     if(any(lil == 0)) {
580 :     pos <- lil > 0
581 :     if(!any(pos)) ## matrix of all 0's
582 :     return(TRUE)
583 :     ilist <- ilist[pos]
584 :     ni <- ni[pos]
585 :     }
586 : maechler 1315 ni0 <- ni - 1:1 # '0-based ni'
587 : maechler 1238 if(is.na(upper)) {
588 : maechler 1315 if(all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0))
589 : maechler 1238 structure(TRUE, kind = "U")
590 : maechler 1315 else if(all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0))
591 : maechler 1238 structure(TRUE, kind = "L")
592 :     else FALSE
593 :     } else if(upper) {
594 : maechler 1315 all(sapply(ilist, max, USE.NAMES = FALSE) <= ni0)
595 : maechler 1238 } else { ## 'lower'
596 : maechler 1315 all(sapply(ilist, min, USE.NAMES = FALSE) >= ni0)
597 : maechler 1238 }
598 :     }
599 :    
600 : maechler 1174 .is.diagonal <- function(object) {
601 : maechler 1357 ## "matrix" or "denseMatrix" (but not "diagonalMatrix")
602 : maechler 1174 d <- dim(object)
603 :     if(d[1] != (n <- d[2])) FALSE
604 : maechler 1357 else if(is.matrix(object))
605 :     ## requires that "vector-indexing" works for 'object' :
606 : maechler 1467 all0(object[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
607 : maechler 1357 else ## "denseMatrix" -- packed or unpacked
608 :     if(is(object, "generalMatrix")) # "dge", "lge", ...
609 : maechler 1467 all0(object@x[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
610 : maechler 1357 else { ## "dense" but not {diag, general}, i.e. triangular or symmetric:
611 :     ## -> has 'uplo' differentiate between packed and unpacked
612 :    
613 :     ### .......... FIXME ...............
614 :    
615 :     packed <- isPacked(object)
616 :     if(object@uplo == "U") {
617 :     } else { ## uplo == "L"
618 :     }
619 :    
620 :     ### very cheap workaround
621 : maechler 1467 all0(as.matrix(object)[rep(c(FALSE, rep.int(TRUE,n)), length = n^2)])
622 : maechler 1357 }
623 : maechler 1174 }
624 : maechler 1253
625 : maechler 1357
626 : maechler 1592 ## FIXME? -- this should also work for "ltT", "ntT", ... :
627 : maechler 1253 diagU2N <- function(x)
628 :     {
629 : maechler 1592 ## Purpose: Transform a *unit diagonal* sparse triangular matrix
630 : maechler 1253 ## into one with explicit diagonal entries '1'
631 :     xT <- as(x, "dgTMatrix")
632 :     ## leave it as T* - the caller can always coerce to C* if needed:
633 :     new("dtTMatrix", x = xT@x, i = xT@i, j = xT@j, Dim = x@Dim,
634 :     Dimnames = x@Dimnames, uplo = x@uplo, diag = "N")
635 :     }
636 : maechler 1290
637 : maechler 1659 ## Needed, e.g., in ./Csparse.R for colSums() etc:
638 : maechler 1295 .as.dgC.Fun <- function(x, na.rm = FALSE, dims = 1) {
639 :     x <- as(x, "dgCMatrix")
640 :     callGeneric()
641 :     }
642 :    
643 : maechler 1659 .as.dgC.0.factors <- function(x) {
644 :     if(!is(x, "dgCMatrix"))
645 :     as(x, "dgCMatrix") # will not have 'factors'
646 :     else ## dgCMatrix
647 :     if(!length(x@factors)) x else { x@factors <- list() ; x }
648 :     }
649 :    
650 : maechler 1295 .as.dgT.Fun <- function(x, na.rm = FALSE, dims = 1) {
651 : maechler 1349 ## used e.g. inside colSums() etc methods
652 : maechler 1295 x <- as(x, "dgTMatrix")
653 :     callGeneric()
654 :     }
655 :    
656 :    
657 : maechler 1290 ### Fast much simplified version of tapply()
658 :     tapply1 <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
659 :     sapply(split(X, INDEX), FUN, ..., simplify = simplify, USE.NAMES = FALSE)
660 :     }
661 :    
662 :     ## tapply.x <- function (X, n, INDEX, FUN = NULL, ..., simplify = TRUE) {
663 :     ## tapply1(X, factor(INDEX, 0:(n-1)), FUN = FUN, ..., simplify = simplify)
664 :     ## }
665 :    

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