SCM

SCM Repository

[matrix] Annotation of /pkg/R/sparseMatrix.R
ViewVC logotype

Annotation of /pkg/R/sparseMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2207 - (view) (download)

1 : bates 684 ### Define Methods that can be inherited for all subclasses
2 :    
3 : maechler 925 ### Idea: Coercion between *VIRTUAL* classes -- as() chooses "closest" classes
4 :     ### ---- should also work e.g. for dense-triangular --> sparse-triangular !
5 : maechler 868
6 : maechler 1472 ##-> see als ./dMatrix.R, ./ddenseMatrix.R and ./lMatrix.R
7 : maechler 868
8 : maechler 1472 setAs("ANY", "sparseMatrix", function(from) as(from, "CsparseMatrix"))
9 : maechler 868
10 : maechler 1845 setAs("sparseMatrix", "generalMatrix", as_gSparse)
11 : maechler 1472
12 : maechler 1845 setAs("sparseMatrix", "symmetricMatrix", as_sSparse)
13 :    
14 :     setAs("sparseMatrix", "triangularMatrix", as_tSparse)
15 :    
16 : maechler 1932 spMatrix <- function(nrow, ncol,
17 :     i = integer(), j = integer(), x = numeric())
18 :     {
19 : maechler 1852 dim <- c(as.integer(nrow), as.integer(ncol))
20 :     ## The conformability of (i,j,x) with itself and with 'dim'
21 :     ## is checked automatically by internal "validObject()" inside new(.):
22 :     kind <- .M.kind(x)
23 :     new(paste(kind, "gTMatrix", sep=''), Dim = dim,
24 :     x = if(kind == "d") as.double(x) else x,
25 :     ## our "Tsparse" Matrices use 0-based indices :
26 :     i = as.integer(i - 1L),
27 :     j = as.integer(j - 1L))
28 :     }
29 :    
30 :    
31 : maechler 871 ## "graph" coercions -- this needs the graph package which is currently
32 :     ## ----- *not* required on purpose
33 :     ## Note: 'undirected' graph <==> 'symmetric' matrix
34 :    
35 : maechler 1271 ## Add some utils that may no longer be needed in future versions of the 'graph' package
36 :     graph.has.weights <- function(g) "weight" %in% names(edgeDataDefaults(g))
37 :    
38 :     graph.wgtMatrix <- function(g)
39 :     {
40 :     ## Purpose: work around "graph" package's as(g, "matrix") bug
41 :     ## ----------------------------------------------------------------------
42 :     ## Arguments: g: an object inheriting from (S4) class "graph"
43 :     ## ----------------------------------------------------------------------
44 :     ## Author: Martin Maechler, based on Seth Falcon's code; Date: 12 May 2006
45 :    
46 :     ## MM: another buglet for the case of "no edges":
47 :     if(numEdges(g) == 0) {
48 :     p <- length(nd <- nodes(g))
49 :     return( matrix(0, p,p, dimnames = list(nd, nd)) )
50 :     }
51 :     ## Usual case, when there are edges:
52 :     has.w <- "weight" %in% names(edgeDataDefaults(g))
53 :     if(has.w) {
54 :     w <- unlist(edgeData(g, attr = "weight"))
55 :     has.w <- any(w != 1)
56 :     } ## now 'has.w' is TRUE iff there are weights != 1
57 :     m <- as(g, "matrix")
58 :     ## now is a 0/1 - matrix (instead of 0/wgts) with the 'graph' bug
59 :     if(has.w) { ## fix it if needed
60 :     tm <- t(m)
61 :     tm[tm != 0] <- w
62 :     t(tm)
63 :     }
64 :     else m
65 :     }
66 :    
67 :    
68 :     setAs("graphAM", "sparseMatrix",
69 : bates 862 function(from) {
70 : maechler 1271 symm <- edgemode(from) == "undirected" && isSymmetric(from@adjMat)
71 :     ## This is only ok if there are no weights...
72 :     if(graph.has.weights(from)) {
73 :     as(graph.wgtMatrix(from),
74 :     if(symm) "dsTMatrix" else "dgTMatrix")
75 :     }
76 :     else { ## no weights: 0/1 matrix -> logical
77 :     as(as(from, "matrix"),
78 : maechler 1548 if(symm) "nsTMatrix" else "ngTMatrix")
79 : maechler 1271 }
80 : bates 862 })
81 : maechler 1271
82 : bates 1476 setAs("graph", "CsparseMatrix",
83 : bates 1479 function(from) as(as(from, "graphNEL"), "CsparseMatrix"))
84 : maechler 687
85 : bates 1474 setAs("graphNEL", "CsparseMatrix",
86 : maechler 1565 function(from) as(as(from, "TsparseMatrix"), "CsparseMatrix"))
87 :    
88 :     setAs("graphNEL", "TsparseMatrix",
89 : maechler 1271 function(from) {
90 : maechler 1565 nd <- nodes(from)
91 : bates 1474 dm <- rep.int(length(nd), 2)
92 : maechler 1271 symm <- edgemode(from) == "undirected"
93 : bates 1474
94 : maechler 1565 if(graph.has.weights(from)) {
95 :     eWts <- edgeWeights(from)
96 :     lens <- unlist(lapply(eWts, length))
97 :     i <- rep.int(0:(dm[1]-1), lens) # column indices (0-based)
98 :     To <- unlist(lapply(eWts, names))
99 : maechler 1832 j <- as.integer(match(To,nd) - 1L) # row indices (0-based)
100 : maechler 1565 ## symm <- symm && <weights must also be symmetric>: improbable
101 :     ## if(symm) new("dsTMatrix", .....) else
102 :     new("dgTMatrix", i = i, j = j, x = unlist(eWts),
103 :     Dim = dm, Dimnames = list(nd, nd))
104 :     }
105 :     else { ## no weights: 0/1 matrix -> logical
106 :     edges <- lapply(from@edgeL[nd], "[[", "edges")
107 :     lens <- unlist(lapply(edges, length))
108 :     ## nnz <- sum(unlist(lens)) # number of non-zeros
109 :     i <- rep.int(0:(dm[1]-1), lens) # column indices (0-based)
110 :     j <- as.integer(unlist(edges) - 1) # row indices (0-based)
111 :     if(symm) { # symmetric: ensure upper triangle
112 :     tmp <- i
113 :     flip <- i > j
114 :     i[flip] <- j[flip]
115 :     j[flip] <- tmp[flip]
116 :     new("nsTMatrix", i = i, j = j, Dim = dm,
117 :     Dimnames = list(nd, nd), uplo = "U")
118 :     } else {
119 :     new("ngTMatrix", i = i, j = j, Dim = dm,
120 :     Dimnames = list(nd, nd))
121 :     }
122 : bates 1474 }
123 : maechler 1271 })
124 : maechler 687
125 : maechler 871 setAs("sparseMatrix", "graph", function(from) as(from, "graphNEL"))
126 :     setAs("sparseMatrix", "graphNEL",
127 : mmaechler 2207 ## since have specific method for Tsparse below, are Csparse or Rsparse,
128 :     ## i.e. do not need to "uniquify" the T* matrix:
129 :     function(from) Tsp2grNEL(as(from, "TsparseMatrix"), need.uniq=FALSE))
130 : maechler 908
131 : mmaechler 2207 Tsp2grNEL <- function(from, need.uniq = is_not_uniqT(from)) {
132 : maechler 1348 d <- dim(from)
133 :     if(d[1] != d[2])
134 : maechler 1513 stop("only square matrices can be used as incidence matrices for graphs")
135 : maechler 1348 n <- d[1]
136 :     if(n == 0) return(new("graphNEL"))
137 :     if(is.null(rn <- dimnames(from)[[1]]))
138 :     rn <- as.character(1:n)
139 : mmaechler 2207 if(need.uniq) ## Need to 'uniquify' the triplets!
140 :     from <- uniq(from)
141 : maechler 908
142 : maechler 1348 if(isSymmetric(from)) { # either "symmetricMatrix" or otherwise
143 :     ##-> undirected graph: every edge only once!
144 :     if(!is(from, "symmetricMatrix")) {
145 :     ## a general matrix which happens to be symmetric
146 :     ## ==> remove the double indices
147 :     from <- tril(from)
148 :     }
149 : maechler 1565 eMode <- "undirected"
150 :     } else {
151 :     eMode <- "directed"
152 : maechler 1348 }
153 : maechler 1565 ## every edge is there only once, either upper or lower triangle
154 : maechler 1832 ft1 <- cbind(rn[from@i + 1L], rn[from@j + 1L])
155 : maechler 1565 ## not yet: graph::ftM2graphNEL(.........)
156 :     ftM2graphNEL(ft1, W = from@x, V= rn, edgemode= eMode)
157 : maechler 871
158 : maechler 1348 }
159 : mmaechler 2207 setAs("TsparseMatrix", "graphNEL", function(from) Tsp2grNEL(from))
160 : maechler 871
161 : maechler 1348
162 : maechler 868 ### Subsetting -- basic things (drop = "missing") are done in ./Matrix.R
163 : maechler 687
164 : maechler 925 ### FIXME : we defer to the "*gT" -- conveniently, but not efficient for gC !
165 : maechler 687
166 : maechler 925 ## [dl]sparse -> [dl]gT -- treat both in one via superclass
167 :     ## -- more useful when have "z" (complex) and even more
168 : maechler 687
169 : maechler 925 setMethod("[", signature(x = "sparseMatrix", i = "index", j = "missing",
170 : maechler 868 drop = "logical"),
171 : maechler 2098 function (x, i,j, ..., drop) {
172 : maechler 1725 cld <- getClassDef(class(x))
173 : maechler 2113 ##> why should this be needed; can still happen in <Tsparse>[..]:
174 :     ##> if(!extends(cld, "generalMatrix")) x <- as(x, "generalMatrix")
175 :     ## viaCl <- paste(.M.kind(x, cld), "gTMatrix", sep='')
176 :     x <- as(x, "TsparseMatrix")[i, , drop=drop]
177 :     ##simpler than x <- callGeneric(x = as(x, "TsparseMatrix"), i=i, drop=drop)
178 : maechler 1751 ## try_as(x, c(cl, sub("T","C", viaCl)))
179 :     if(is(x, "Matrix") && extends(cld, "CsparseMatrix"))
180 :     as(x, "CsparseMatrix") else x
181 :     })
182 : maechler 687
183 : maechler 925 setMethod("[", signature(x = "sparseMatrix", i = "missing", j = "index",
184 : maechler 868 drop = "logical"),
185 : maechler 2098 function (x,i,j, ..., drop) {
186 : maechler 1725 cld <- getClassDef(class(x))
187 : maechler 2113 ##> why should this be needed; can still happen in <Tsparse>[..]:
188 :     ##> if(!extends(cld, "generalMatrix")) x <- as(x, "generalMatrix")
189 :     ## viaCl <- paste(.M.kind(x, cld), "gTMatrix", sep='')
190 :    
191 :     x <- as(x, "TsparseMatrix")[, j, drop=drop]
192 :     ##simpler than x <- callGeneric(x = as(x, "TsparseMatrix"), j=j, drop=drop)
193 : maechler 1751 if(is(x, "Matrix") && extends(cld, "CsparseMatrix"))
194 :     as(x, "CsparseMatrix") else x
195 :     })
196 : maechler 868
197 : maechler 925 setMethod("[", signature(x = "sparseMatrix",
198 : maechler 886 i = "index", j = "index", drop = "logical"),
199 : maechler 2056 function (x, i, j, ..., drop) {
200 : maechler 1725 cld <- getClassDef(class(x))
201 : maechler 1665 ## be smart to keep symmetric indexing of <symm.Mat.> symmetric:
202 : maechler 2113 ##> doSym <- (extends(cld, "symmetricMatrix") &&
203 :     ##> length(i) == length(j) && all(i == j))
204 :     ##> why should this be needed; can still happen in <Tsparse>[..]:
205 :     ##> if(!doSym && !extends(cld, "generalMatrix"))
206 :     ##> x <- as(x, "generalMatrix")
207 :     ## viaCl <- paste(.M.kind(x, cld),
208 :     ## if(doSym) "sTMatrix" else "gTMatrix", sep='')
209 :     x <- as(x, "TsparseMatrix")[i, j, drop=drop]
210 : maechler 1725 if(is(x, "Matrix") && extends(cld, "CsparseMatrix"))
211 : maechler 1751 as(x, "CsparseMatrix") else x
212 : maechler 1665 })
213 : maechler 868
214 :    
215 : maechler 1673 ## setReplaceMethod("[", .........)
216 :     ## -> ./Tsparse.R
217 : mmaechler 2207 ## & ./Csparse.R & ./Rsparse.R {those go via Tsparse}
218 :     ##
219 :     ## Do not use as.vector() (see ./Matrix.R ) for "scarce" matrices :
220 :     setReplaceMethod("[", signature(x = "sparseMatrix", i = "missing", j = "ANY",
221 :     value = "scarceMatrix"),
222 :     function (x, i, j, ..., value)
223 :     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
224 : maechler 1226
225 : mmaechler 2207 setReplaceMethod("[", signature(x = "sparseMatrix", i = "ANY", j = "missing",
226 :     value = "scarceMatrix"),
227 :     function (x, i, j, ..., value)
228 :     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
229 : maechler 1226
230 : mmaechler 2207 setReplaceMethod("[", signature(x = "sparseMatrix", i = "ANY", j = "ANY",
231 :     value = "scarceMatrix"),
232 :     function (x, i, j, ..., value)
233 :     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
234 :    
235 :    
236 :    
237 : maechler 1714 ## Group Methods
238 : maechler 1226
239 : maechler 1737 setMethod("Math",
240 :     signature(x = "sparseMatrix"),
241 :     function(x) callGeneric(as(x, "CsparseMatrix")))
242 : maechler 868
243 : mmaechler 2175 ## further group methods -> see ./Ops.R {"Summary": ./dMatrix.R }
244 : maechler 1472
245 :    
246 : maechler 1737
247 : maechler 687 ### --- show() method ---
248 :    
249 : maechler 1389 ## FIXME(?) -- ``merge this'' (at least ``synchronize'') with
250 :     ## - - - prMatrix() from ./Auxiliaries.R
251 : maechler 1705 ## FIXME: prTriang() in ./Auxiliaries.R should also get align = "fancy"
252 : mmaechler 2203 ##
253 : maechler 1947 printSpMatrix <- function(x, digits = getOption("digits"),
254 :     maxp = getOption("max.print"), zero.print = ".",
255 :     col.names, note.dropping.colnames = TRUE,
256 :     col.trailer = '', align = c("fancy", "right"))
257 : maechler 687 {
258 : maechler 1903 cl <- getClassDef(class(x))
259 : maechler 1737 stopifnot(extends(cl, "sparseMatrix"))
260 : mmaechler 2175 validObject(x) # have seen seg.faults for invalid objects
261 : maechler 1903 d <- dim(x)
262 : maechler 1389 if(prod(d) > maxp) { # "Large" => will be "cut"
263 :     ## only coerce to dense that part which won't be cut :
264 :     nr <- maxp %/% d[2]
265 : maechler 1903 m <- as(x[1:max(1, nr), ,drop=FALSE], "Matrix")
266 : maechler 1389 } else {
267 : maechler 1903 m <- as(x, "matrix")
268 : maechler 1389 }
269 : maechler 1947 dn <- dimnames(m) ## will be === dimnames(cx)
270 : maechler 1737 logi <- extends(cl,"lsparseMatrix") || extends(cl,"nsparseMatrix")
271 : maechler 1315 if(logi)
272 : maechler 1947 cx <- array("N", dim(m), dimnames=dn)
273 : mmaechler 2203 else { ## numeric (or --not yet implemented-- complex):
274 : maechler 1903 cx <- apply(m, 2, format)
275 :     if(is.null(dim(cx))) {# e.g. in 1 x 1 case
276 :     dim(cx) <- dim(m)
277 : maechler 1947 dimnames(cx) <- dn
278 : maechler 1315 }
279 : maechler 687 }
280 : maechler 1947 if (missing(col.names))
281 :     col.names <- {
282 :     if(!is.null(cc <- getOption("sparse.colnames")))
283 :     cc
284 :     else if(is.null(dn[[2]]))
285 :     FALSE
286 :     else { # has column names == dn[[2]]
287 :     ncol(x) < 10
288 :     }
289 :     }
290 :     if(identical(col.names, FALSE))
291 : maechler 1903 cx <- emptyColnames(cx, msg.if.not.empty = note.dropping.colnames)
292 : maechler 1947 else if(is.character(col.names)) {
293 :     stopifnot(length(col.names) == 1)
294 :     cn <- col.names
295 :     switch(substr(cn, 1,3),
296 :     "abb" = {
297 :     iarg <- as.integer(sub("^[^0-9]*", '', cn))
298 :     colnames(cx) <- abbreviate(colnames(cx), minlength = iarg)
299 :     },
300 :     "sub" = {
301 :     iarg <- as.integer(sub("^[^0-9]*", '', cn))
302 :     colnames(cx) <- substr(colnames(cx), 1, iarg)
303 :     },
304 :     stop("invalid 'col.names' string: ", cn))
305 :     }
306 :     ## else: nothing to do for col.names == TRUE
307 : maechler 687 if(is.logical(zero.print))
308 :     zero.print <- if(zero.print) "0" else " "
309 : maechler 1315 if(logi) {
310 : maechler 1903 cx[!m] <- zero.print
311 :     cx[m] <- "|"
312 : maechler 1315 } else { # non logical
313 :     ## show only "structural" zeros as 'zero.print', not all of them..
314 :     ## -> cannot use 'm'
315 : maechler 1903 d <- dim(cx)
316 : mmaechler 2203 ne <- length(iN0 <- 1L + .Call(m_encodeInd, non0ind(x, cl), di = d))
317 : maechler 1705 if(0 < ne && ne < prod(d)) {
318 :     align <- match.arg(align)
319 : maechler 2098 if(align == "fancy" && !is.integer(m)) {
320 : maechler 1705 fi <- apply(m, 2, format.info) ## fi[3,] == 0 <==> not expo.
321 :     ## now 'format' the zero.print by padding it with ' ' on the right:
322 :     ## case 1: non-exponent: fi[2,] + as.logical(fi[2,] > 0)
323 :     ## the column numbers of all 'zero' entries -- (*large*)
324 : maechler 1832 cols <- 1L + (0:(prod(d)-1L))[-iN0] %/% d[1]
325 : maechler 1705 pad <-
326 :     ifelse(fi[3,] == 0,
327 :     fi[2,] + as.logical(fi[2,] > 0),
328 :     ## exponential:
329 :     fi[2,] + fi[3,] + 4)
330 : maechler 1737 ## now be efficient ; sprintf() is relatively slow
331 :     ## and pad is much smaller than 'cols'; instead of "simply"
332 :     ## zero.print <- sprintf("%-*s", pad[cols] + 1, zero.print)
333 :     if(any(doP <- pad > 0)) {#
334 :     ## only pad those that need padding - *before* expanding
335 :     z.p.pad <- rep.int(zero.print, length(pad))
336 :     z.p.pad[doP] <- sprintf("%-*s", pad[doP] + 1, zero.print)
337 :     zero.print <- z.p.pad[cols]
338 :     }
339 :     else
340 :     zero.print <- rep.int(zero.print, length(cols))
341 : maechler 1705 } ## else "right" : nothing to do
342 :    
343 : maechler 1903 cx[-iN0] <- zero.print
344 : maechler 1705 } else if (ne == 0)# all zeroes
345 : maechler 1903 cx[] <- zero.print
346 : maechler 1315 }
347 : maechler 1870 if(col.trailer != '')
348 : maechler 1903 cx <- cbind(cx, col.trailer, deparse.level = 0)
349 : maechler 1705 ## right = TRUE : cheap attempt to get better "." alignment
350 : maechler 1903 print(cx, quote = FALSE, right = TRUE, max = maxp)
351 :     invisible(x)
352 : mmaechler 2203 } ## printSpMatrix()
353 : maechler 687
354 : mmaechler 2203 printSpMatrix2 <- function(x, digits = getOption("digits"),
355 :     maxp = getOption("max.print"), zero.print = ".",
356 :     col.names, note.dropping.colnames = TRUE,
357 :     suppRows = NULL, suppCols = NULL,
358 :     col.trailer = if(suppCols) "......" else "",
359 :     align = c("fancy", "right"))
360 :     {
361 :     d <- dim(x)
362 :     if((identical(suppRows,FALSE) && identical(suppCols, FALSE)) ||
363 :     (!isTRUE(suppRows) && !isTRUE(suppCols) && prod(d) <= maxp))
364 :     {
365 :     if(missing(col.trailer) && is.null(suppCols))
366 :     suppCols <- FALSE # for 'col.trailer'
367 :     printSpMatrix(x, digits=digits, maxp=maxp,
368 :     zero.print=zero.print, col.names=col.names,
369 :     note.dropping.colnames=note.dropping.colnames,
370 :     col.trailer=col.trailer, align=align)
371 :     } else { ## d[1] > maxp / d[2] >= nr : -- this needs [,] working:
372 :     nR <- d[1] ## nrow
373 :     useW <- getOption("width") - (format.info(nR)[1] + 3+1)
374 :     ## space for "[<last>,] "
375 : maechler 1903
376 : mmaechler 2203 ## --> suppress rows and/or columns in printing ...
377 : maechler 1737
378 : mmaechler 2203 if(is.null(suppCols)) suppCols <- (d[2] * 2 > useW)
379 :     nc <- if(suppCols) (useW - (1 + nchar(col.trailer))) %/% 2 else d[2]
380 :     nr <- maxp %/% nc
381 :     if(is.null(suppRows)) suppRows <- (nr < nR)
382 : maechler 1845
383 : mmaechler 2203 sTxt <- c("in show(); maybe adjust 'options(max.print= *)'",
384 :     "\n ..............................\n")
385 :     if(suppRows) {
386 :     if(suppCols)
387 :     x <- x[ , 1:nc, drop = FALSE]
388 :     n2 <- ceiling(nr / 2)
389 :     printSpMatrix(x[seq_len(min(nR, max(1, n2))), , drop=FALSE],
390 :     digits=digits, maxp=maxp,
391 :     zero.print=zero.print, col.names=col.names,
392 :     note.dropping.colnames=note.dropping.colnames,
393 :     col.trailer = col.trailer, align=align)
394 :     cat("\n ..............................",
395 :     "\n ........suppressing rows ", sTxt, "\n", sep='')
396 :     ## tail() automagically uses "[..,]" rownames:
397 :     printSpMatrix(tail(x, max(1, nr-n2)),
398 :     digits=digits, maxp=maxp,
399 :     zero.print=zero.print, col.names=col.names,
400 :     note.dropping.colnames=note.dropping.colnames,
401 :     col.trailer = col.trailer, align=align)
402 :     }
403 :     else if(suppCols) {
404 :     printSpMatrix(x[ , 1:nc , drop = FALSE],
405 :     digits=digits, maxp=maxp,
406 :     zero.print=zero.print, col.names=col.names,
407 :     note.dropping.colnames=note.dropping.colnames,
408 :     col.trailer = col.trailer, align=align)
409 :     cat("\n .....suppressing columns ", sTxt, sep='')
410 :     }
411 :     else stop("logic programming error in printSpMatrix2(), please report")
412 : maechler 1845
413 : mmaechler 2203 invisible(x)
414 :     }
415 :     } ## printSpMatrix2 ()
416 : maechler 1737
417 : mmaechler 2203 setMethod("print", signature(x = "sparseMatrix"), printSpMatrix)
418 : maechler 1737
419 : mmaechler 2203 setMethod("show", signature(object = "sparseMatrix"),
420 :     function(object) {
421 :     d <- dim(object)
422 :     cl <- class(object)
423 :     cat(sprintf('%d x %d sparse Matrix of class "%s"\n',
424 :     d[1], d[2], cl))
425 :     printSpMatrix2(object)
426 :     })
427 : maechler 886
428 :    
429 : maechler 1852 ## For very large and very sparse matrices, the above show()
430 :     ## is not really helpful; Use summary() as an alternative:
431 :    
432 :     setMethod("summary", signature(object = "sparseMatrix"),
433 :     function(object, ...) {
434 :     d <- dim(object)
435 :     T <- as(object, "TsparseMatrix")
436 :     ## return a data frame (int, int, {double|logical|...}) :
437 : mmaechler 2175 r <- if(is(object,"nsparseMatrix"))
438 :     data.frame(i = T@i + 1L, j = T@j + 1L)
439 :     else data.frame(i = T@i + 1L, j = T@j + 1L, x = T@x)
440 : maechler 1852 attr(r, "header") <-
441 : maechler 1855 sprintf('%d x %d sparse Matrix of class "%s", with %d entries',
442 : mmaechler 2175 d[1], d[2], class(object), length(T@i))
443 : maechler 1852 ## use ole' S3 technology for such a simple case
444 :     class(r) <- c("sparseSummary", class(r))
445 :     r
446 :     })
447 :    
448 :     print.sparseSummary <- function (x, ...) {
449 :     cat(attr(x, "header"),"\n")
450 :     print.data.frame(x, ...)
451 :     invisible(x)
452 :     }
453 :    
454 : maechler 1108 setMethod("isSymmetric", signature(object = "sparseMatrix"),
455 : maechler 2113 function(object, tol = 100*.Machine$double.eps, ...) {
456 : maechler 886 ## pretest: is it square?
457 :     d <- dim(object)
458 :     if(d[1] != d[2]) return(FALSE)
459 : maechler 2115
460 :     ## else slower test using t() --
461 :    
462 :     ## FIXME (for tol = 0): use cholmod_symmetry(A, 1, ...)
463 :     ## for tol > 0 should modify cholmod_symmetry(..) to work with tol
464 :    
465 :     ## or slightly simpler, rename and export is_sym() in ../src/cs_utils.c
466 :    
467 :    
468 : maechler 973 if (is(object, "dMatrix"))
469 : maechler 886 ## use gC; "T" (triplet) is *not* unique!
470 : maechler 1659 isTRUE(all.equal(.as.dgC.0.factors( object),
471 : maechler 2113 .as.dgC.0.factors(t(object)),
472 :     tol = tol, ...))
473 : maechler 973 else if (is(object, "lMatrix"))
474 : maechler 886 ## test for exact equality; FIXME(?): identical() too strict?
475 : maechler 2115 identical(as( object, "lgCMatrix"),
476 : maechler 886 as(t(object), "lgCMatrix"))
477 : maechler 1548 else if (is(object, "nMatrix"))
478 :     ## test for exact equality; FIXME(?): identical() too strict?
479 : maechler 2115 identical(as( object, "ngCMatrix"),
480 : maechler 1548 as(t(object), "ngCMatrix"))
481 : maechler 886 else stop("not yet implemented")
482 :     })
483 : maechler 1174
484 : maechler 1238
485 : maechler 2110 setMethod("isTriangular", signature(object = "CsparseMatrix"), isTriC)
486 :     setMethod("isTriangular", signature(object = "TsparseMatrix"), isTriT)
487 : maechler 1174
488 :     setMethod("isDiagonal", signature(object = "sparseMatrix"),
489 :     function(object) {
490 : maechler 1799 d <- dim(object)
491 :     if(d[1] != d[2]) return(FALSE)
492 :     ## else
493 : maechler 1174 gT <- as(object, "TsparseMatrix")
494 :     all(gT@i == gT@j)
495 :     })
496 :    
497 : maechler 1290
498 : mmaechler 2175 setMethod("determinant", signature(x = "sparseMatrix", logarithm = "missing"),
499 :     function(x, logarithm, ...)
500 :     determinant(x, logarithm = TRUE, ...))
501 :     setMethod("determinant", signature(x = "sparseMatrix", logarithm = "logical"),
502 :     function(x, logarithm = TRUE, ...)
503 :     determinant(as(x,"dsparseMatrix"), logarithm, ...))
504 :    
505 :    
506 : maechler 1472 setMethod("diag", signature(x = "sparseMatrix"),
507 :     function(x, nrow, ncol = n) diag(as(x, "CsparseMatrix")))
508 :    
509 : maechler 1845 setMethod("dim<-", signature(x = "sparseMatrix", value = "ANY"),
510 :     function(x, value) {
511 :     if(!is.numeric(value) || length(value) != 2)
512 :     stop("dim(.) value must be numeric of length 2")
513 :     if(prod(dim(x)) != prod(value <- as.integer(value)))
514 :     stop("dimensions don't match the number of cells")
515 :     ## be careful to keep things sparse
516 :     as(spV2M(as(x, "sparseVector"), nrow=value[1], ncol=value[2]),
517 :     class(x))
518 :     })
519 :    
520 :    
521 : maechler 2005 setMethod("norm", signature(x = "sparseMatrix", type = "character"),
522 :     function(x, type, ...) {
523 :     type <- toupper(substr(type[1], 1, 1))
524 :     switch(type, ## max(<empty>, 0) |--> 0
525 :     "O" = ,
526 :     "1" = max(colSums(abs(x)), 0), ## One-norm (L_1)
527 :     "I" = max(rowSums(abs(x)), 0), ## L_Infinity
528 :     "F" = sqrt(sum(x^2)), ## Frobenius
529 :     "M" = max(abs(x), 0), ## Maximum modulus of all
530 :     ## otherwise:
531 :     stop("invalid 'type'"))
532 :     })
533 :    
534 : maechler 2158 setMethod("rcond", signature(x = "sparseMatrix", norm = "character"),
535 :     function(x, norm, ...) {
536 : maechler 2072 d <- dim(x)
537 : maechler 2102 ## FIXME: qr.R(qr(.)) warns about differing R (permutation!)
538 :     ## really fix qr.R() *or* go via dense in any cases
539 : maechler 2072 rcond(if(d[1] == d[2]) {
540 : maechler 2102 warning("rcond(.) via sparse -> dense coercion")
541 : maechler 2072 as(x, "denseMatrix")
542 :     } else if(d[1] > d[2]) qr.R(qr(x)) else qr.R(qr(t(x))),
543 : maechler 2158 norm = norm, ...)
544 : maechler 2072 })
545 : maechler 2048
546 : maechler 2096 setMethod("cov2cor", signature(V = "sparseMatrix"),
547 :     function(V) {
548 :     ## like stats::cov2cor() but making sure all matrices stay sparse
549 :     p <- (d <- dim(V))[1]
550 :     if (p != d[2])
551 :     stop("'V' is not a *square* matrix")
552 :     if(!is(V, "dMatrix"))
553 :     V <- as(V, "dMatrix")# actually "dsparseMatrix"
554 :     Is <- sqrt(1/diag(V))
555 :     if (any(!is.finite(Is))) ## original had 0 or NA
556 :     warning("diag(.) had 0 or NA entries; non-finite result is doubtful")
557 :     ## TODO: if <diagonal> %*% <sparse> was implemented more efficiently
558 :     ## we'd rather use that!
559 :     Is <- as(Diagonal(x = Is), "sparseMatrix")
560 :     r <- Is %*% V %*% Is
561 :     r[cbind(1L:p,1L:p)] <- 1 # exact in diagonal
562 :     as(r, "symmetricMatrix")
563 :     })
564 : maechler 2072
565 : mmaechler 2175 setMethod("is.na", signature(x = "sparseMatrix"),## NB: nsparse* have own method!
566 : maechler 2144 function(x) {
567 :     if(any((inax <- is.na(x@x)))) {
568 : mmaechler 2175 cld <- getClassDef(class(x))
569 :     if(extends(cld, "triangularMatrix") && x@diag == "U")
570 :     inax <- is.na((x <- .diagU2N(x, cld))@x)
571 :     r <- as(x, "lMatrix") # will be "lsparseMatrix" - *has* x slot
572 : maechler 2144 r@x <- inax
573 : mmaechler 2185 if(!extends(cld, "CsparseMatrix"))
574 :     r <- as(r, "CsparseMatrix")
575 :     as(.Call(Csparse_drop, r, 0), "nMatrix") # a 'pattern matrix
576 : maechler 2144 }
577 : mmaechler 2175 else is.na_nsp(x)
578 : maechler 2144 })
579 : maechler 2096
580 :    
581 : mmaechler 2199 ### Keep this namespace-hidden: Would need to return a classed object
582 : bates 1895 lm.fit.sparse <-
583 :     function(x, y, offset = NULL, method = c("qr", "cholesky"),
584 :     tol = 1e-7, singular.ok = TRUE, transpose = FALSE, ...)
585 : maechler 2043 ### Fit a linear model, __ given __ a sparse model matrix 'x'
586 :     ### using a sparse QR or a sparse Cholesky factorization
587 : bates 1895 {
588 :     stopifnot(is(x, "dsparseMatrix"))
589 : maechler 2005 ## if(!is(x, "dsparseMatrix"))
590 :     ## x <- as(x, "dsparseMatrix")
591 : bates 1895 yy <- as.numeric(y)
592 :     if (!is.null(offset)) {
593 : maechler 2005 stopifnot(length(offset) == length(y))
594 :     yy <- yy - as.numeric(offset)
595 : bates 1895 }
596 :     ans <- switch(as.character(method)[1],
597 : maechler 2005 cholesky =
598 :     .Call(dgCMatrix_cholsol,
599 :     as(if (transpose) x else t(x), "dgCMatrix"), yy),
600 :     qr =
601 :     .Call(dgCMatrix_qrsol,
602 :     as(if (transpose) t(x) else x, "dgCMatrix"), yy),
603 :     ## otherwise:
604 :     stop("unknown method ", dQuote(method))
605 :     )
606 : bates 1895 ans
607 :     }
608 :    
609 : maechler 1911 fac2sparse <- function(from, to = c("d","i","l","n","z"))
610 :     {
611 :     ## factor(-like) --> sparseMatrix {also works for integer, character}
612 :     levs <- levels(fact <- factor(from)) # drop unused levels
613 :     n <- length(fact)
614 :     to <- match.arg(to)
615 :     res <- new(paste(to, "gCMatrix", sep=''))
616 :     res@i <- as.integer(fact) - 1L # 0-based
617 :     res@p <- 0:n
618 :     res@Dim <- c(length(levs), n)
619 :     res@Dimnames <- list(levs, NULL)
620 :     if(to != "n")
621 :     res@x <- rep.int(switch(to,
622 :     "d" = 1., "i" = 1L, "l" = TRUE, "z" = 1+0i),
623 :     n)
624 :     res
625 :     }
626 : bates 1895
627 : maechler 1911 setAs("factor", "sparseMatrix", function(from) fac2sparse(from, to = "d"))
628 : bates 2094
629 : maechler 2095 ## xtabs returning a sparse matrix. This is cut'n'paste
630 :     ## of xtabs() in <Rsrc>/src/library/stats/R/xtabs.R ;
631 :     ## with the new argument 'sparse'
632 :     xtabs <- function(formula = ~., data = parent.frame(), subset, sparse = FALSE,
633 :     na.action, exclude = c(NA, NaN), drop.unused.levels = FALSE)
634 : bates 2094 {
635 : maechler 2095 if (missing(formula) && missing(data))
636 :     stop("must supply either 'formula' or 'data'")
637 :     if(!missing(formula)){
638 :     ## We need to coerce the formula argument now, but model.frame
639 :     ## will coerce the original version later.
640 :     formula <- as.formula(formula)
641 :     if (!inherits(formula, "formula"))
642 :     stop("'formula' missing or incorrect")
643 : bates 2094 }
644 : maechler 2095 if (any(attr(terms(formula, data = data), "order") > 1))
645 :     stop("interactions are not allowed")
646 : bates 2094 m <- match.call(expand.dots = FALSE)
647 : maechler 2095 if (is.matrix(eval(m$data, parent.frame())))
648 :     m$data <- as.data.frame(data)
649 :     m$... <- m$exclude <- m$drop.unused.levels <- m$sparse <- NULL
650 : bates 2094 m[[1]] <- as.name("model.frame")
651 :     mf <- eval(m, parent.frame())
652 :     if(length(formula) == 2) {
653 :     by <- mf
654 :     y <- NULL
655 :     }
656 :     else {
657 :     i <- attr(attr(mf, "terms"), "response")
658 :     by <- mf[-i]
659 :     y <- mf[[i]]
660 :     }
661 :     by <- lapply(by, function(u) {
662 :     if(!is.factor(u)) u <- factor(u, exclude = exclude)
663 :     u[ , drop = drop.unused.levels]
664 :     })
665 : maechler 2095 if(!sparse) { ## identical to stats::xtabs
666 :     x <-
667 :     if(is.null(y))
668 :     do.call("table", by)
669 :     else if(NCOL(y) == 1)
670 :     tapply(y, by, sum)
671 :     else {
672 :     z <- lapply(as.data.frame(y), tapply, by, sum)
673 :     array(unlist(z),
674 :     dim = c(dim(z[[1]]), length(z)),
675 :     dimnames = c(dimnames(z[[1]]), list(names(z))))
676 :     }
677 :     x[is.na(x)] <- 0
678 :     class(x) <- c("xtabs", "table")
679 :     attr(x, "call") <- match.call()
680 :     x
681 :    
682 :     } else { ## sparse
683 :     if (length(by) != 2)
684 :     stop("xtabs(*, sparse=TRUE) applies only to two-way tables")
685 :     rows <- by[[1]]
686 :     cols <- by[[2]]
687 :     rl <- levels(rows)
688 :     cl <- levels(cols)
689 :     if (is.null(y))
690 :     y <- rep.int(1, length(rows))
691 :     as(new("dgTMatrix",
692 :     i = as.integer(rows) - 1L,
693 :     j = as.integer(cols) - 1L,
694 : bates 2097 x = as.double(y),
695 : maechler 2095 Dim = c(length(rl), length(cl)),
696 :     Dimnames = list(rl, cl)), "CsparseMatrix")
697 :     }
698 : bates 2094 }

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