https://github.com/cran/Matrix
Tip revision: 891e875c01bd4f2885bed8274b0bec3e2926f839 authored by Martin Maechler on 26 March 2014, 00:00:00 UTC
version 1.1-3
version 1.1-3
Tip revision: 891e875
indexing.Rout.save
R Under development (unstable) (2014-01-26 r64897) -- "Unsuffered Consequences"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> #### For both 'Extract' ("[") and 'Replace' ("[<-") Method testing
> #### aka subsetting and subassignment
>
> #### suppressPackageStartupMessages(...) as we have an *.Rout.save to Rdiff against
> stopifnot(suppressPackageStartupMessages(require(Matrix)))
>
> source(system.file("test-tools.R", package = "Matrix"), keep.source = FALSE)
> ##-> identical3() etc
> cat("doExtras:",doExtras,"\n")
doExtras: FALSE
>
> if(interactive()) {
+ options(error = recover, warn = 1)
+ } else if(FALSE) { ## MM @ testing *manually* only
+ options(error = recover, Matrix.verbose = TRUE, warn = 1)
+ } else {
+ options(Matrix.verbose = TRUE, warn = 1)
+ }
>
>
> ### Dense Matrices
>
> m <- Matrix(1:28 +0, nrow = 7)
> validObject(m)
[1] TRUE
> stopifnot(identical(m, m[]),
+ identical(m[2, 3], 16), # simple number
+ identical(m[2, 3:4], c(16,23)), # simple numeric of length 2
+ identical(m[NA,NA], as(Matrix(NA, 7,4), "dMatrix")))
>
> m[2, 3:4, drop=FALSE] # sub matrix of class 'dgeMatrix'
1 x 2 Matrix of class "dgeMatrix"
[,1] [,2]
[1,] 16 23
> m[-(4:7), 3:4] # ditto; the upper right corner of 'm'
3 x 2 Matrix of class "dgeMatrix"
[,1] [,2]
[1,] 15 22
[2,] 16 23
[3,] 17 24
>
> ## rows or columns only:
> m[1,] # first row, as simple numeric vector
[1] 1 8 15 22
> m[,2] # 2nd column
[1] 8 9 10 11 12 13 14
> m[,1:2] # sub matrix of first two columns
7 x 2 Matrix of class "dgeMatrix"
[,1] [,2]
[1,] 1 8
[2,] 2 9
[3,] 3 10
[4,] 4 11
[5,] 5 12
[6,] 6 13
[7,] 7 14
> m[-(1:6),, drop=FALSE] # not the first 6 rows, i.e. only the 7th
1 x 4 Matrix of class "dgeMatrix"
[,1] [,2] [,3] [,4]
[1,] 7 14 21 28
> m[integer(0),] #-> 0 x 4 Matrix
0 x 4 Matrix of class "dgeMatrix"
[,1] [,2] [,3] [,4]
> m[2:4, numeric(0)] #-> 3 x 0 Matrix
3 x 0 Matrix of class "dgeMatrix"
[1,]
[2,]
[3,]
>
> ## logical indexing
> stopifnot(identical(m[2,3], m[(1:nrow(m)) == 2, (1:ncol(m)) == 3]),
+ identical(m[2,], m[(1:nrow(m)) == 2, ]),
+ identical(m[,3:4], m[, (1:4) >= 3]))
>
> ## dimnames indexing:
> mn <- m
> dimnames(mn) <- list(paste("r",letters[1:nrow(mn)],sep=""),
+ LETTERS[1:ncol(mn)])
> checkMatrix(mn)
norm(m [7 x 4]) : 1 I F M ok
Summary: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
> mn["rd", "D"]
[1] 25
> msr <- ms <- as(mn,"sparseMatrix")
> mnr <- mn
> v <- rev(as(ms, "vector"))
> mnr[] <- v
> msr[] <- v # [<- "sparse" -- not very sensical; did fail w/o a message
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> z <- msr; z[] <- 0
> zz <- as(array(0, dim(z)), "sparseMatrix")
> a.m <- as(mnr,"matrix")
> stopifnot(identical(mn["rc", "D"], mn[3,4]), mn[3,4] == 24,
+ identical(mn[, "A"], mn[,1]), mn[,1] == 1:7,
+ identical(mn[c("re", "rb"), "B"], mn[c(5,2), 2]),
+ identical(ms["rc", "D"], ms[3,4]), ms[3,4] == 24,
+ identical(ms[, "A"], ms[,1]), ms[,1] == 1:7,
+ identical(ms[ci <- c("re", "rb"), "B"], ms[c(5,2), 2]),
+ identical(rownames(mn[ci, ]), ci),
+ identical(rownames(ms[ci, ]), ci),
+ identical(colnames(mn[,cj <- c("B","D")]), cj),
+ identical(colnames(ms[,cj]), cj),
+ identical(a.m, as(msr,"matrix")),
+ identical(unname(z), zz),
+ identical(a.m, array(v, dim=dim(mn), dimnames=dimnames(mn)))
+ )
> showProc.time()
Time elapsed: 0.232 0.002 0.237
>
> ## Printing sparse colnames:
> ms[sample(28, 20)] <- 0
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> ms <- t(rbind2(ms, 3*ms))
> cnam1 <- capture.output(show(ms))[2] ; op <- options("sparse.colnames" = "abb3")
[[ suppressing 14 column names 'ra', 'rb', 'rc' ... ]]
> cnam2 <- capture.output(show(ms))[2] ; options(op) # revert
> stopifnot(## sparse printing
+ grep("^ +$", cnam1) == 1, # cnam1 is empty
+ identical(cnam2,
+ paste(" ", paste(rep(rownames(mn), 2), collapse=" "))))
>
> mo <- m
> m[2,3] <- 100
> m[1:2, 4] <- 200
> m[, 1] <- -1
> m[1:3,]
3 x 4 Matrix of class "dgeMatrix"
[,1] [,2] [,3] [,4]
[1,] -1 8 15 200
[2,] -1 9 100 200
[3,] -1 10 17 24
>
> m. <- as.matrix(m)
>
> ## m[ cbind(i,j) ] indexing:
> iN <- ij <- cbind(1:6, 2:3)
> iN[2:3,] <- iN[5,2] <- NA
> stopifnot(identical(m[ij], m.[ij]),
+ identical(m[iN], m.[iN]))
>
> ## testing operations on logical Matrices rather more than indexing:
> g10 <- m [ m > 10 ]
> stopifnot(18 == length(g10))
> stopifnot(10 == length(m[ m <= 10 ]))
> sel <- (20 < m) & (m < 150)
> sel.<- (20 < m.)& (m.< 150)
> nsel <-(20 >= m) | (m >= 150)
> (ssel <- as(sel, "sparseMatrix"))
7 x 4 sparse Matrix of class "lgCMatrix"
[1,] . . . .
[2,] . . | .
[3,] . . . |
[4,] . . . |
[5,] . . . |
[6,] . . . |
[7,] . . | |
> stopifnot(is(sel, "lMatrix"), is(ssel, "lsparseMatrix"),
+ identical3(as.mat(sel.), as.mat(sel), as.mat(ssel)),
+ identical3(!sel, !ssel, nsel), # !<sparse> is typically dense
+ identical3(m[ sel], m[ ssel], as.matrix(m)[as.matrix( ssel)]),
+ identical3(m[!sel], m[!ssel], as.matrix(m)[as.matrix(!ssel)])
+ )
> showProc.time()
Time elapsed: 0.122 0.001 0.123
>
> ## more sparse Matrices --------------------------------------
>
> ##' @title Check sparseMatrix sub-assignment m[i,j] <- v
> ##' @param ms sparse Matrix
> ##' @param mm its [traditional matrix]-equivalent
> ##' @param k (approximate) length of index vectors (i,j)
> ##' @param n.uniq (approximate) number of unique values in i,j
> ##' @param vRNG function(n) for random 'v' generation
> ##' @param show logical; if TRUE, it will not stop on error
> ##' @return
> ##' @author Martin Maechler
> chkAssign <- function(ms, mm = as(ms, "matrix"),
+ k = min(20,dim(mm)), n.uniq = k %/% 3,
+ vRNG = { if(is.numeric(mm) || is.complex(mm))
+ function(n) rpois(n,lambda= 0.75)# <- about 47% zeros
+ else ## logical
+ function(n) runif(n) > 0.8 }, ## 80% zeros
+ showOnly=FALSE)
+ {
+ stopifnot(is(ms,"sparseMatrix"))
+ d <- dim(ms)
+ s1 <- function(n) sample(n, pmin(n, pmax(1, rpois(1, n.uniq))))
+ i <- sample(s1(d[1]), k/2+ rpois(1, k/2), replace = TRUE)
+ j <- sample(s1(d[2]), k/2+ rpois(1, k/2), replace = TRUE)
+ assert.EQ.mat(ms[i,j], mm[i,j])
+ ms2 <- ms. <- ms; mm. <- mm # save
+ ## now sub*assign* to these repeated indices, and then compare -----
+ v <- vRNG(length(i) * length(j))
+ mm[i,j] <- v
+ ms[i,j] <- v
+ ## useful to see (ii,ij), but confusing R/ESS when additionally debugging:
+ ## if(!showOnly && interactive()) { op <- options(error = recover); on.exit(options(op)) }
+ assert.EQ.mat(ms, mm, show=showOnly)
+ ## vector indexing m[cbind(i,j)] == m[i + N(j-1)] , N = nrow(.)
+ ii <- seq_len(min(length(i), length(j)))
+ i <- i[ii]
+ j <- j[ii]
+ ij <- cbind(i, j)
+ ii <- i + nrow(ms)*(j - 1)
+ ord.i <- order(ii)
+ iio <- ii[ord.i]
+ ui <- unique(iio) # compare these with :
+ neg.ii <- - setdiff(seq_len(prod(d)), ii)
+ stopifnot(identical(mm[ii], mm[ij]),
+ identical(ms.[ui], ms.[neg.ii]),
+ ms.[ij] == mm.[ii], ## M[ cbind(i,j) ] was partly broken; now checking
+ ms.[ii] == mm.[ii])
+ v <- v[seq_len(length(i))]
+ if(is(ms,"nMatrix")) v <- as.logical(v) # !
+ mm.[ij] <- v
+ ms.[ii] <- v
+ nodup <- (length(ui) == length(ii)) ## <==> ! anyDuplicated(iio)
+ if(nodup) { cat("[-]") # rare, unfortunately
+ ms2[neg.ii] <- v[ord.i]
+ stopifnot(identical(ms2, ms.))
+ }
+ assert.EQ.mat(ms., mm., show=showOnly)
+ } ##{chkAssign}
>
> ## Get duplicated index {because these are "hard" (and rare)
> getDuplIndex <- function(n, k) {
+ repeat {
+ i <- sample(n, k, replace=TRUE) # 3 4 6 9 2 9 : 9 is twice
+ if(anyDuplicated(i)) break
+ }
+ i
+ }
>
> ## From package 'sfsmisc':
> repChar <- function (char, no) paste(rep.int(char, no), collapse = "")
>
> m <- 1:800
> set.seed(101) ; m[sample(800, 600)] <- 0
> m0 <- Matrix(m, nrow = 40)
> m1 <- add.simpleDimnames(m0)
> for(kind in c("n", "l", "d")) {
+ for(m in list(m0,m1)) { ## -- with and without dimnames -------------------------
+ kClass <-paste0(kind, "Matrix" )
+ Ckind <- paste0(kind, "gCMatrix")
+ Tkind <- paste0(kind, "gTMatrix")
+ str(mC <- as(m, Ckind))
+ str(mT <- as(as(as(m, kClass), "TsparseMatrix"), Tkind))
+ mm <- as(mC, "matrix") # also logical or double
+ IDENT <- if(kind == "n") function(x,y) Q.eq2(x,y, tol=0) else identical
+ stopifnot(identical(mT, as(as(mC, "TsparseMatrix"), Tkind)),
+ identical(mC, as(mT, Ckind)),
+ Qidentical(mC[0,0], new(Ckind)),
+ Qidentical(mT[0,0], new(Tkind)),
+ identical(unname(mT[0,]), new(Tkind, Dim = c(0L,ncol(m)))),
+ identical(unname(mT[,0]), new(Tkind, Dim = c(nrow(m),0L))),
+ IDENT(mC[0,], as(mT[FALSE,], Ckind)),
+ IDENT(mC[,0], as(mT[,FALSE], Ckind)),
+ sapply(pmin(min(dim(mC)), c(0:2, 5:10)),
+ function(k) {i <- seq_len(k); all(mC[i,i] == mT[i,i])}),
+ TRUE)
+ cat("ok\n")
+ show(mC[,1])
+ show(mC[1:2,])
+ show(mC[7, drop = FALSE])
+ assert.EQ.mat(mC[1:2,], mm[1:2,])
+ assert.EQ.mat(mC[0,], mm[0,])
+ assert.EQ.mat(mC[,FALSE], mm[,FALSE])
+ ##
+ ## *repeated* (aka 'duplicated') indices - did not work at all ...
+ i <- pmin(nrow(mC), rep(8:10,2))
+ j <- c(2:4, 4:3)
+ assert.EQ.mat(mC[i,], mm[i,])
+ assert.EQ.mat(mC[,j], mm[,j])
+ ## FIXME? assert.EQ.mat(mC[,NA], mm[,NA]) -- mC[,NA] is all 0 "instead" of all NA
+ ## MM currently thinks we should NOT allow <sparse>[ <NA> ]
+ assert.EQ.mat(mC[i, 2:1], mm[i, 2:1])
+ assert.EQ.mat(mC[c(4,1,2:1), j], mm[c(4,1,2:1), j])
+ assert.EQ.mat(mC[i,j], mm[i,j])
+ ##
+ ## set.seed(7)
+ op <- options(Matrix.verbose = FALSE)
+ cat(" for(): ")
+ for(n in 1:(if(doExtras) 50 else 5)) {
+ chkAssign(mC, mm)
+ chkAssign(mC[-3,-2], mm[-3,-2])
+ cat(".")
+ }
+ options(op)
+ cat(sprintf("\n[Ok]%s\n\n", repChar("-", 64)))
+ }
+ cat(sprintf("\nok( %s )\n== ###%s\n\n", kind, repChar("=", 70)))
+ }## end{for}---------------------------------------------------------------
Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ p : int [1:21] 0 8 22 28 37 41 50 63 71 81 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ factors : list()
Formal class 'ngTMatrix' [package "Matrix"] with 5 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ j : int [1:200] 0 0 0 0 0 0 0 0 1 1 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ factors : list()
Note: method with signature 'nsparseMatrix#sparseMatrix' chosen for function '==',
target signature 'ngCMatrix#ngTMatrix'.
"nMatrix#nMatrix", "sparseMatrix#nsparseMatrix" would also be valid
ok
[1] FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
[25] TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[37] FALSE TRUE TRUE FALSE
2 x 20 sparse Matrix of class "ngCMatrix"
[1,] . . . | . . | . . . . | . . | . | . . .
[2,] . | . . . | . . . . . . | | . . . . | .
[1] TRUE
for(): .....
[Ok]----------------------------------------------------------------
Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ p : int [1:21] 0 8 22 28 37 41 50 63 71 81 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
.. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
..@ factors : list()
Formal class 'ngTMatrix' [package "Matrix"] with 5 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ j : int [1:200] 0 0 0 0 0 0 0 0 1 1 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
.. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
..@ factors : list()
ok
r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 r13
FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE
r14 r15 r16 r17 r18 r19 r20 r21 r22 r23 r24 r25 r26
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE
r27 r28 r29 r30 r31 r32 r33 r34 r35 r36 r37 r38 r39
FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE
r40
FALSE
2 x 20 sparse Matrix of class "ngCMatrix"
[[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
r1 . . . | . . | . . . . | . . | . | . . .
r2 . | . . . | . . . . . . | | . . . . | .
[1] TRUE
for(): .....
[Ok]----------------------------------------------------------------
ok( n )
== ###======================================================================
Formal class 'lgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ p : int [1:21] 0 8 22 28 37 41 50 63 71 81 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
..@ factors : list()
Formal class 'lgTMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ j : int [1:200] 0 0 0 0 0 0 0 0 1 1 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
..@ factors : list()
ok
[1] FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
[25] TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[37] FALSE TRUE TRUE FALSE
2 x 20 sparse Matrix of class "lgCMatrix"
[1,] . . . | . . | . . . . | . . | . | . . .
[2,] . | . . . | . . . . . . | | . . . . | .
[1] TRUE
for(): .....
[Ok]----------------------------------------------------------------
Formal class 'lgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ p : int [1:21] 0 8 22 28 37 41 50 63 71 81 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
.. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
..@ x : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
..@ factors : list()
Formal class 'lgTMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ j : int [1:200] 0 0 0 0 0 0 0 0 1 1 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
.. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
..@ x : logi [1:200] TRUE TRUE TRUE TRUE TRUE TRUE ...
..@ factors : list()
ok
r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 r13
FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE
r14 r15 r16 r17 r18 r19 r20 r21 r22 r23 r24 r25 r26
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE
r27 r28 r29 r30 r31 r32 r33 r34 r35 r36 r37 r38 r39
FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE
r40
FALSE
2 x 20 sparse Matrix of class "lgCMatrix"
[[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
r1 . . . | . . | . . . . | . . | . | . . .
r2 . | . . . | . . . . . . | | . . . . | .
[1] TRUE
for(): .....
[Ok]----------------------------------------------------------------
ok( l )
== ###======================================================================
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ p : int [1:21] 0 8 22 28 37 41 50 63 71 81 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:200] 3 7 12 22 25 30 38 39 42 45 ...
..@ factors : list()
Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ j : int [1:200] 0 0 0 0 0 0 0 0 1 1 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:200] 3 7 12 22 25 30 38 39 42 45 ...
..@ factors : list()
ok
[1] 0 0 3 0 0 0 7 0 0 0 0 12 0 0 0 0 0 0 0 0 0 22 0 0 25
[26] 0 0 0 0 30 0 0 0 0 0 0 0 38 39 0
2 x 20 sparse Matrix of class "dgCMatrix"
[1,] . . . 121 . . 241 . . . . 441 . . 561 . 641 . . .
[2,] . 42 . . . 202 . . . . . . 482 522 . . . . 722 .
[1] 7
for(): .....
[Ok]----------------------------------------------------------------
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ p : int [1:21] 0 8 22 28 37 41 50 63 71 81 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
.. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
..@ x : num [1:200] 3 7 12 22 25 30 38 39 42 45 ...
..@ factors : list()
Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:200] 2 6 11 21 24 29 37 38 1 4 ...
..@ j : int [1:200] 0 0 0 0 0 0 0 0 1 1 ...
..@ Dim : int [1:2] 40 20
..@ Dimnames:List of 2
.. ..$ : chr [1:40] "r1" "r2" "r3" "r4" ...
.. ..$ : chr [1:20] "c1" "c2" "c3" "c4" ...
..@ x : num [1:200] 3 7 12 22 25 30 38 39 42 45 ...
..@ factors : list()
ok
r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 r13 r14 r15 r16 r17 r18 r19 r20
0 0 3 0 0 0 7 0 0 0 0 12 0 0 0 0 0 0 0 0
r21 r22 r23 r24 r25 r26 r27 r28 r29 r30 r31 r32 r33 r34 r35 r36 r37 r38 r39 r40
0 22 0 0 25 0 0 0 0 30 0 0 0 0 0 0 0 38 39 0
2 x 20 sparse Matrix of class "dgCMatrix"
[[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
r1 . . . 121 . . 241 . . . . 441 . . 561 . 641 . . .
r2 . 42 . . . 202 . . . . . . 482 522 . . . . 722 .
[1] 7
for(): .....
[Ok]----------------------------------------------------------------
ok( d )
== ###======================================================================
> showProc.time()
Time elapsed: 0.842 0.001 0.846
>
> if(doExtras) {### {was ./AAA_index.R, MM-only}
+ ## an nsparse-example
+ A <- Matrix(c(rep(c(1,0,0),2), rep(c(2,0),7), c(0,0,2), rep(0,4)), 3,9)
+ i <- c(3,1:2)
+ j <- c(3, 5, 9, 5, 9)
+ vv <- logical(length(i)*length(j)); vv[6:9] <- TRUE
+
+ print(An <- as(A,"nMatrix")); an <- as(An, "matrix")
+ assert.EQ.mat(An, an)
+ An[i, j] <- vv
+ an[i, j] <- vv
+ assert.EQ.mat(An, an)## error
+ if(!all(An == an)) show(drop0(An - an))
+ ## all are +1
+
+ options("Matrix.subassign.verbose" = TRUE)# output from C
+ An <- as(A,"nMatrix"); An[i, j] <- vv
+ ## and compare with this:
+ Al <- as(A,"lMatrix"); Al[i, j] <- vv
+ options("Matrix.subassign.verbose" = FALSE)
+
+ ##--- An interesting not small not large example for M[i,j] <- v ------------
+ ##
+ M <- Matrix(c(1, rep(0,7), 1:4), 3,4)
+ N0 <- kronecker(M,M)
+ mkN1 <- function(M) {
+ stopifnot(length(d <- dim(M)) == 2)
+ isC <- is(M,"CsparseMatrix")
+ M[,d[2]] <- c(0,2,0)
+ N <- kronecker(diag(x=1:2), M)## remains sparse if 'M' is
+ if(isC) N <- as(N, "CsparseMatrix")
+ diag(N[-1,]) <- -2
+ N[9,] <- 1:4 # is recycled
+ N[,12] <- -7:-9 # ditto
+ N
+ }
+
+ show(N1 <- t(N <- mkN1(N0))) # transpose {for display reasons}
+ C1 <- t(C <- mkN1(as(N0,"CsparseMatrix")))
+ stopifnot(all(C == N))
+ assert.EQ.mat(C, mkN1(as.matrix(N0)))
+
+ C. <- C1
+ show(N <- N1) ; n <- as.matrix(N); str(N)
+ sort(i <- c(6,8,19,11,21,20,10,7,12,9,5,18,17,22,13))## == c(5:13, 17:22))
+ sort(j <- c(3,8,6,15,10,4,14,13,16,2,11,17,7,5))## == c(2:8, 10:11, 13:17)
+ val <- v.l <- 5*c(0,6,0,7,0,0,8:9, 0,0)
+ show(spv <- as(val, "sparseVector")); str(spv)
+
+ n [i,j] <- v.l
+ N [i,j] <- val# is recycled, too
+ C.[i,j] <- val
+ assert.EQ.mat(N,n) ; stopifnot(all(C. == N))
+ ## and the same *again*:
+ n [i,j] <- v.l
+ N [i,j] <- val
+ C.[i,j] <- val
+ assert.EQ.mat(N,n)
+ stopifnot(all(C. == N))
+
+ print(load(system.file("external", "symA.rda", package="Matrix"))) # "As"
+ stopifnot(isValid(As, "dsCMatrix"), identical(As@factors, list()))
+ R. <- drop0(chol(As))
+ stopifnot(1:32 == sort(diag(R.)), ## !
+ R.@x == as.integer(R.@x),## so it is an integer-valued chol-decomp !
+ ## shows that (1) As is *not* singular (2) the matrix is not random
+ all.equal(crossprod(R.), As, tolerance =1e-15))
+ print(summary(evA <- eigen(As, only.values=TRUE)$val))
+ print(tail(evA)) ## largest three ~= 10^7, smallest two *negative*
+ print(rcond(As)) # 1.722 e-21 == very bad !
+ ##-> this *is* a border line case, i.e. very close to singular !
+ ## and also determinant(.) is rather random here!
+ cc0 <- Cholesky(As)# no problem
+ try({
+ cc <- Cholesky(As, super=TRUE)
+ ## gives --on 32-bit only--
+ ## Cholmod error 'matrix not positive definite' at file:../Supernodal/t_cholmod_super_numeric.c, line 613
+ ecc <- expand(cc)
+ L.P <- with(ecc, crossprod(L,P)) ## == L'P
+ ## crossprod(L.P) == (L'P)' L'P == P'LL'P
+ stopifnot( all.equal(crossprod(L.P), As) )
+ })
+ ##---- end{ eigen( As ) -----------
+
+ } ## only if(doExtras)
>
>
> ##---- Symmetric indexing of symmetric Matrix ----------
> m. <- mC
> m.[, c(2, 7:12)] <- 0
> isValid(S <- crossprod(add.simpleDimnames(m.) %% 100), "dsCMatrix")
[1] TRUE
> ss <- as(S, "matrix")
> ds <- as(S, "denseMatrix")
> ## NA-indexing of *dense* Matrices: should work as traditionally
> assert.EQ.mat(ds[NA,NA], ss[NA,NA])
> assert.EQ.mat(ds[NA, ], ss[NA,])
> assert.EQ.mat(ds[ ,NA], ss[,NA])
> T <- as(S, "TsparseMatrix")
> stopifnot(identical(ds[2 ,NA], ss[2,NA]),
+ identical(ds[NA, 1], ss[NA, 1]),
+ identical(S, as(T, "CsparseMatrix")) )
>
> ## non-repeated indices:
> i <- c(7:5, 2:4);assert.EQ.mat(T[i,i], ss[i,i])
> ## NA in indices -- check that we get a helpful error message:
> i[2] <- NA
> er <- tryCatch(T[i,i], error = function(e)e)
> stopifnot(as.logical(grep("indices.*sparse Matrices", er$message)))
>
> N <- nrow(T)
> set.seed(11)
> for(n in 1:(if(doExtras) 50 else 3)) {
+ i <- sample(N, max(2, sample(N,1)), replace = FALSE)
+ validObject(Tii <- T[i,i]) ; tTi <- t(T)[i,i]
+ stopifnot(is(Tii, "dsTMatrix"), # remained symmetric Tsparse
+ is(tTi, "dsTMatrix"), # may not be identical when *sorted* differently
+ identical(as(t(Tii),"CsparseMatrix"), as(tTi,"CsparseMatrix")))
+ assert.EQ.mat(Tii, ss[i,i])
+ }
>
> b <- diag(1:2)[,c(1,1,2,2)]
> cb <- crossprod(b)
> cB <- crossprod(Matrix(b, sparse=TRUE))
> a <- matrix(0, 6, 6)
> a[1:4, 1:4] <- cb
> A1 <- A2 <- Matrix(0, 6, 6)#-> sparse
> A1[1:4, 1:4] <- cb
> A2[1:4, 1:4] <- cB
> assert.EQ.mat(A1, a)# indeed
> stopifnot(identical(A1, A2), is(A1, "dsCMatrix"))
>
> ## repeated ones ``the challenge'' (to do smartly):
> j <- c(4, 4, 9, 12, 9, 4, 17, 3, 18, 4, 12, 18, 4, 9)
> assert.EQ.mat(T[j,j], ss[j,j])
> ## and another two sets (a, A) & (a., A.) :
> a <- matrix(0, 6,6)
> a[upper.tri(a)] <- (utr <- c(2, 0,-1, 0,0,5, 7,0,0,0, 0,0,-2,0,8))
> ta <- t(a); ta[upper.tri(a)] <- utr; a <- t(ta)
> diag(a) <- c(0,3,0,4,6,0)
> A <- as(Matrix(a), "TsparseMatrix")
> A. <- A
> diag(A.) <- 10 * (1:6)
> a. <- as(A., "matrix")
> ## More testing {this was not working for a long time..}
> set.seed(1)
> for(n in 1:(if(doExtras) 100 else 6)) {
+ i <- sample(1:nrow(A), 3+2*rpois(1, lam=3), replace=TRUE)
+ Aii <- A[i,i]
+ A.ii <- A.[i,i]
+ stopifnot(class(Aii) == class(A),
+ class(A.ii) == class(A.))
+ assert.EQ.mat(Aii , a [i,i])
+ assert.EQ.mat(A.ii, a.[i,i])
+ assert.EQ.mat(T[i,i], ss[i,i])
+ }
> showProc.time()
Time elapsed: 0.265 0 0.265
>
> stopifnot(all.equal(mC[,3], mm[,3]),
+ identical(mC[ij], mC[ij + 0.4]),
+ identical(mC[ij], mm[ij]),
+ identical(mC[iN], mm[iN]))
> ## out of bound indexing must be detected:
> assertError(mC[cbind(ij[,1] - 5, ij[,2])])
> assertError(mC[cbind(ij[,1], ij[,2] + ncol(mC))])
>
> assert.EQ.mat(mC[7, , drop=FALSE], mm[7, , drop=FALSE])
> identical (mC[7, drop=FALSE], mm[7, drop=FALSE]) # *vector* indexing
[1] TRUE
>
> stopifnot(dim(mC[numeric(0), ]) == c(0,20), # used to give warnings
+ dim(mC[, integer(0)]) == c(40,0),
+ identical(mC[, integer(0)], mC[, FALSE]))
> validObject(print(mT[,c(2,4)]))
40 x 2 sparse Matrix of class "dgTMatrix"
c2 c4
r1 . 121
r2 42 .
r3 . .
r4 . .
r5 45 .
r6 . .
r7 . .
r8 . 128
r9 . 129
r10 50 .
r11 . .
r12 52 132
r13 . 133
r14 . .
r15 55 .
r16 . .
r17 . .
r18 . 138
r19 . .
r20 . .
r21 . 141
r22 . 142
r23 63 .
r24 . .
r25 65 .
r26 . .
r27 67 .
r28 68 .
r29 . .
r30 . .
r31 71 .
r32 72 .
r33 . .
r34 74 .
r35 . .
r36 76 .
r37 . .
r38 . .
r39 . 159
r40 80 .
[1] TRUE
> stopifnot(all.equal(mT[2,], mm[2,]),
+ ## row or column indexing in combination with t() :
+ Q.C.identical(mT[2,], t(mT)[,2]),
+ Q.C.identical(mT[-2,], t(t(mT)[,-2])),
+ Q.C.identical(mT[c(2,5),], t(t(mT)[,c(2,5)])) )
> assert.EQ.mat(mT[4,, drop = FALSE], mm[4,, drop = FALSE])
> stopifnot(identical3(mm[,1], mC[,1], mT[,1]),
+ identical3(mm[3,], mC[3,], mT[3,]),
+ identical3(mT[2,3], mC[2,3], 0),
+ identical(mT[], mT),
+ identical4( mm[c(3,7), 2:4], as.mat( m[c(3,7), 2:4]),
+ as.mat(mT[c(3,7), 2:4]), as.mat(mC[c(3,7), 2:4]))
+ )
>
> x.x <- crossprod(mC)
> stopifnot(class(x.x) == "dsCMatrix",
+ class(x.x. <- round(x.x / 10000)) == "dsCMatrix",
+ identical(x.x[cbind(2:6, 2:6)],
+ diag(x.x [2:6, 2:6])))
> head(x.x.) # Note the *non*-structural 0's printed as "0"
6 x 20 sparse Matrix of class "dgCMatrix"
[[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
c1 1 0 . 1 . 1 1 3 . 3 2 1 6 1 . 2 4 6 5 1
c2 0 6 2 1 3 5 7 5 12 14 14 9 11 16 12 13 17 19 19 10
c3 . 2 6 . 4 2 5 3 8 12 5 16 9 11 23 . . 6 7 7
c4 1 1 . 17 . 8 10 13 8 6 18 18 29 35 14 8 25 10 19 21
c5 . 3 4 . 14 4 10 . . 29 8 9 19 11 11 . . 26 26 16
c6 1 5 2 8 4 42 5 19 14 9 8 10 42 56 50 27 29 32 64 16
> tail(x.x., -3) # all but the first three lines
17 x 20 sparse Matrix of class "dgCMatrix"
[[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
c4 1 1 . 17 . 8 10 13 8 6 18 18 29 35 14 8 25 10 19 21
c5 . 3 4 . 14 4 10 . . 29 8 9 19 11 11 . . 26 26 16
c6 1 5 2 8 4 42 5 19 14 9 8 10 42 56 50 27 29 32 64 16
c7 1 7 5 10 10 5 87 14 9 31 77 47 79 43 28 17 67 110 36 121
c8 3 5 3 13 . 19 14 70 10 24 37 13 59 62 34 19 58 21 64 44
c9 . 12 8 8 . 14 9 10 116 41 58 33 33 72 78 43 69 72 75 25
c10 3 14 12 6 29 9 31 24 41 167 69 56 99 44 70 24 105 82 85 32
c11 2 14 5 18 8 8 77 37 58 69 267 80 86 139 49 105 194 119 122 129
c12 1 9 16 18 9 10 47 13 33 56 80 194 70 77 81 . 90 32 . 106
c13 6 11 9 29 19 42 79 59 33 99 86 70 324 157 55 . 69 142 144 155
c14 1 16 11 35 11 56 43 62 72 44 139 77 157 375 123 102 145 39 196 81
c15 . 12 23 14 11 50 28 34 78 70 49 81 55 123 368 71 112 41 41 86
c16 2 13 . 8 . 27 17 19 43 24 105 . . 102 71 233 124 44 139 .
c17 4 17 . 25 . 29 67 58 69 105 194 90 69 145 112 124 523 141 245 100
c18 6 19 6 10 26 32 110 21 72 82 119 32 142 39 41 44 141 497 104 111
c19 5 19 7 19 26 64 36 64 75 85 122 . 144 196 41 139 245 104 542 55
c20 1 10 7 21 16 16 121 44 25 32 129 106 155 81 86 . 100 111 55 541
>
> lx.x <- as(x.x, "lsCMatrix") # FALSE only for "structural" 0
> (l10 <- lx.x[1:10, 1:10])# "lsC"
10 x 10 sparse Matrix of class "lsCMatrix"
[[ suppressing 10 column names 'c1', 'c2', 'c3' ... ]]
c1 | | . | . | | | . |
c2 | | | | | | | | | |
c3 . | | . | | | | | |
c4 | | . | . | | | | |
c5 . | | . | | | . . |
c6 | | | | | | | | | |
c7 | | | | | | | | | |
c8 | | | | . | | | | |
c9 . | | | . | | | | |
c10 | | | | | | | | | |
> (l3 <- lx.x[1:3, ])
3 x 20 sparse Matrix of class "lgCMatrix"
[[ suppressing 20 column names 'c1', 'c2', 'c3' ... ]]
c1 | | . | . | | | . | | | | | . | | | | |
c2 | | | | | | | | | | | | | | | | | | | |
c3 . | | . | | | | | | | | | | | . . | | |
> m.x <- as.mat(x.x) # as.mat() *drops* (NULL,NULL) dimnames
> stopifnot(class(l10) == "lsCMatrix", # symmetric indexing -> symmetric !
+ identical(as.mat(lx.x), m.x != 0),
+ identical(as.logical(lx.x), as.logical(m.x)),
+ identical(as.mat(l10), m.x[1:10, 1:10] != 0),
+ identical(as.mat(l3 ), m.x[1:3, ] != 0)
+ )
>
> ##-- Sub*assignment* with repeated / duplicated index:
> A <- Matrix(0,4,3) ; A[c(1,2,1), 2] <- 1 ; A
4 x 3 sparse Matrix of class "dgCMatrix"
[1,] . 1 .
[2,] . 1 .
[3,] . . .
[4,] . . .
> B <- A; B[c(1,2,1), 2] <- 1:3; B; B. <- B
4 x 3 sparse Matrix of class "dgCMatrix"
[1,] . 3 .
[2,] . 2 .
[3,] . . .
[4,] . . .
> B.[3,] <- rbind(4:2)
> ## change the diagonal and the upper and lower subdiagonal :
> diag(B.) <- 10 * diag(B.)
> diag(B.[,-1]) <- 5* diag(B.[,-1])
> diag(B.[-1,]) <- 4* diag(B.[-1,]) ; B.
4 x 3 sparse Matrix of class "dgCMatrix"
[1,] . 15 .
[2,] . 20 .
[3,] 4 12 20
[4,] . . .
> C <- B.; C[,2] <- C[,2]; C[1,] <- C[1,]; C[2:3,2:1] <- C[2:3,2:1]
> stopifnot(identical(unname(as.matrix(A)),
+ local({a <- matrix(0,4,3); a[c(1,2,1), 2] <- 1 ; a})),
+ identical(unname(as.matrix(B)),
+ local({a <- matrix(0,4,3); a[c(1,2,1), 2] <- 1:3; a})),
+ identical(C, drop0(B.)))
> ## <sparse>[<logicalSparse>] <- v failed in the past
> T <- as(C,"TsparseMatrix"); C. <- C
> T[T>0] <- 21
Note: method with signature 'TsparseMatrix#Matrix#missing#replValue' chosen for function '[<-',
target signature 'dgTMatrix#lgTMatrix#missing#numeric'.
"Matrix#lsparseMatrix#missing#replValue" would also be valid
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> C[C>0] <- 21
Note: method with signature 'CsparseMatrix#Matrix#missing#replValue' chosen for function '[<-',
target signature 'dgCMatrix#lgCMatrix#missing#numeric'.
"Matrix#lsparseMatrix#missing#replValue" would also be valid
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> a. <- local({a <- as.matrix(C.); a[a>0] <- 21; a})
> assert.EQ.mat(C, a.)
> stopifnot(identical(C, as(T, "CsparseMatrix")))
>
> ## used to fail
> n <- 5 ## or much larger
> sm <- new("dsTMatrix", i=1L, j=1L, Dim=as.integer(c(n,n)), x = 1)
> (cm <- as(sm, "CsparseMatrix"))
5 x 5 sparse Matrix of class "dsCMatrix"
[1,] . . . . .
[2,] . 1 . . .
[3,] . . . . .
[4,] . . . . .
[5,] . . . . .
> sm[2,]
[1] 0 1 0 0 0
> stopifnot(sm[2,] == c(0:1, rep.int(0,ncol(sm)-2)),
+ sm[2,] == cm[2,],
+ sm[,3] == sm[3,],
+ all(sm[,-(1:3)] == t(sm[-(1:3),])), # all(<lge.>)
+ all(sm[,-(1:3)] == 0)
+ )
> showProc.time()
Time elapsed: 0.149 0.001 0.151
>
> ##--- "nsparse*" sub-assignment :----------
> M <- Matrix(c(1, rep(0,7), 1:4), 3,4)
> N0 <- kronecker(M,M)
> Nn <- as(N0, "nMatrix"); nn <- as(Nn,"matrix")
> (Nn00 <- Nn0 <- Nn); nn00 <- nn0 <- nn
9 x 16 sparse Matrix of class "ngTMatrix"
[1,] | . . | . . . . . . . . | . . |
[2,] . . . | . . . . . . . . . . . |
[3,] . . | | . . . . . . . . . . | |
[4,] . . . . . . . . . . . . | . . |
[5,] . . . . . . . . . . . . . . . |
[6,] . . . . . . . . . . . . . . | |
[7,] . . . . . . . . | . . | | . . |
[8,] . . . . . . . . . . . | . . . |
[9,] . . . . . . . . . . | | . . | |
>
> set.seed(1)
> Nn0 <- Nn00; nn0 <- nn00
> for(i in 1:(if(doExtras) 200 else 25)) {
+ Nn <- Nn0
+ nn <- nn0
+ i. <- getDuplIndex(nrow(N0), 6)
+ j. <- getDuplIndex(ncol(N0), 4)
+ vv <- sample(c(FALSE,TRUE),
+ length(i.)*length(j.), replace=TRUE)
+ cat(",")
+ Nn[i., j.] <- vv
+ nn[i., j.] <- vv
+ assert.EQ.mat(Nn, nn)
+ if(!all(Nn == nn)) {
+ cat("i=",i,":\n i. <- "); dput(i.)
+ cat("j. <- "); dput(j.)
+ cat("which(vv): "); dput(which(vv))
+ cat("Difference matrix:\n")
+ show(drop0(Nn - nn))
+ }
+ cat("k")
+ ## sub-assign double precision to logical sparseMatrices: now *with* warning:
+ ## {earlier: gave *no* warning}:
+ assertWarning(Nn[1:2,] <- -pi)
+ assertWarning(Nn[, 5] <- -pi)
+ assertWarning(Nn[2:4, 5:8] <- -pi)
+ stopifnot(isValid(Nn,"nsparseMatrix"))
+ ##
+ cat(".")
+ if(i %% 10 == 0) cat("\n")
+ if(i == 100) {
+ Nn0 <- as(Nn0, "CsparseMatrix")
+ cat("Now: class", class(Nn0)," :\n~~~~~~~~~~~~~~~~~\n")
+ }
+ }
,k.,k.,k.,k.,k.,k.,k.,k.,k.,k.
,k.,k.,k.,k.,k.,k.,k.,k.,k.,k.
,k.,k.,k.,k.,k.> showProc.time()
Time elapsed: 0.198 0 0.199
> Nn <- Nn0
> ## Check that NA is interpreted as TRUE (with a warning), for "nsparseMatrix":
> assertWarning(Nn[ii <- 3 ] <- NA); stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii])
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> assertWarning(Nn[ii <- 22:24] <- NA); stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii])
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> assertWarning(Nn[ii <- -(1:99)] <- NA); stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii])
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> assertWarning(Nn[ii <- 3:4 ] <- c(0,NA))
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii] == 0:1)
> assertWarning(Nn[ii <- 25:27] <- c(0,1,NA))
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> stopifnot(isValid(Nn,"nsparseMatrix"), Nn[ii] == c(FALSE,TRUE,TRUE))
>
> m0 <- Diagonal(5)
> stopifnot(identical(m0[2,], m0[,2]),
+ identical(m0[,1], c(1,0,0,0,0)))
> ### Diagonal -- Sparse:
> (m1 <- as(m0, "TsparseMatrix")) # dtTMatrix unitriangular
5 x 5 sparse Matrix of class "dtTMatrix" (unitriangular)
[1,] 1 . . . .
[2,] . 1 . . .
[3,] . . 1 . .
[4,] . . . 1 .
[5,] . . . . 1
> (m2 <- as(m0, "CsparseMatrix")) # dtCMatrix unitriangular
5 x 5 sparse Matrix of class "dtCMatrix" (unitriangular)
[1,] 1 . . . .
[2,] . 1 . . .
[3,] . . 1 . .
[4,] . . . 1 .
[5,] . . . . 1
> m1g <- as(m1, "generalMatrix")
> stopifnot(is(m1g, "dgTMatrix"))
> assert.EQ.mat(m2[1:3,], diag(5)[1:3,])
> assert.EQ.mat(m2[,c(4,1)], diag(5)[,c(4,1)])
> stopifnot(identical(m2[1:3,], as(m1[1:3,], "CsparseMatrix")),
+ identical(Matrix:::uniqTsparse(m1[, c(4,2)]),
+ Matrix:::uniqTsparse(as(m2[, c(4,2)], "TsparseMatrix")))
+ )## failed in 0.9975-11
>
> (uTr <- new("dtTMatrix", Dim = c(3L,3L), diag="U"))
3 x 3 sparse Matrix of class "dtTMatrix" (unitriangular)
[1,] 1 . .
[2,] . 1 .
[3,] . . 1
> uTr[1,] <- 0
> assert.EQ.mat(uTr, cbind(0, rbind(0,diag(2))))
>
> M <- m0; M[1,] <- 0
> Z <- m0; Z[] <- 0; z <- array(0, dim(M))
> stopifnot(identical(M, Diagonal(x=c(0, rep(1,4)))),
+ all(Z == 0), Qidentical(as(Z, "matrix"), z))
> M <- m0; M[,3] <- 3 ; M ; stopifnot(is(M, "sparseMatrix"), M[,3] == 3)
5 x 5 sparse Matrix of class "dgCMatrix"
[1,] 1 . 3 . .
[2,] . 1 3 . .
[3,] . . 3 . .
[4,] . . 3 1 .
[5,] . . 3 . 1
> checkMatrix(M)
Note: method with signature 'sparseMatrix#ldiMatrix' chosen for function '==',
target signature 'nsCMatrix#ldiMatrix'.
"nsparseMatrix#sparseMatrix", "nMatrix#lMatrix" would also be valid
norm(m [5 x 5]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: ok; determinant(): ok
> M <- m0; M[1:3, 3] <- 0 ;M
5 x 5 diagonal matrix of class "ddiMatrix"
[,1] [,2] [,3] [,4] [,5]
[1,] 1 . . . .
[2,] . 1 . . .
[3,] . . 0 . .
[4,] . . . 1 .
[5,] . . . . 1
> T <- m0; T[1:3, 3] <- 10
> stopifnot(identical(M, Diagonal(x=c(1,1, 0, 1,1))),
+ isValid(T, "triangularMatrix"), identical(T[,3], c(10,10,10,0,0)))
>
> M <- m1; M[1,] <- 0 ; M ; assert.EQ.mat(M, diag(c(0,rep(1,4))), tol=0)
5 x 5 sparse Matrix of class "dtTMatrix"
[1,] . . . . .
[2,] . 1 . . .
[3,] . . 1 . .
[4,] . . . 1 .
[5,] . . . . 1
> M <- m1; M[,3] <- 3 ; stopifnot(is(M,"sparseMatrix"), M[,3] == 3)
M[i,j] <- v : coercing symmetric M[] into non-symmetric
> Z <- m1; Z[] <- 0
> checkMatrix(M)
norm(m [5 x 5]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: Note: method with signature 'sparseMatrix#ldiMatrix' chosen for function '&',
target signature 'nsCMatrix#ldiMatrix'.
"nsparseMatrix#sparseMatrix", "nMatrix#lMatrix" would also be valid
ok
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: ok; determinant(): ok
> M <- m1; M[1:3, 3] <- 0 ;M
5 x 5 sparse Matrix of class "dtTMatrix"
[1,] 1 . . . .
[2,] . 1 . . .
[3,] . . . . .
[4,] . . . 1 .
[5,] . . . . 1
> assert.EQ.mat(M, diag(c(1,1, 0, 1,1)), tol=0)
> T <- m1; T[1:3, 3] <- 10; checkMatrix(T)
norm(m [5 x 5]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: ok; determinant(): ok
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
as(<triangular (ge)matrix>, dtTMatrix): valid: TRUE
> stopifnot(is(T, "triangularMatrix"), identical(T[,3], c(10,10,10,0,0)),
+ Qidentical(as(Z, "matrix"), z))
>
> M <- m2; M[1,] <- 0 ; M ; assert.EQ.mat(M, diag(c(0,rep(1,4))), tol=0)
5 x 5 sparse Matrix of class "dtCMatrix"
[1,] . . . . .
[2,] . 1 . . .
[3,] . . 1 . .
[4,] . . . 1 .
[5,] . . . . 1
> M <- m2; M[,3] <- 3 ; stopifnot(is(M,"sparseMatrix"), M[,3] == 3)
> checkMatrix(M)
norm(m [5 x 5]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: ok; determinant(): ok
> Z <- m2; Z[] <- 0
> M <- m2; M[1:3, 3] <- 0 ;M
5 x 5 sparse Matrix of class "dtCMatrix"
[1,] 1 . . . .
[2,] . 1 . . .
[3,] . . . . .
[4,] . . . 1 .
[5,] . . . . 1
> assert.EQ.mat(M, diag(c(1,1, 0, 1,1)), tol=0)
> T <- m2; T[1:3, 3] <- 10; checkMatrix(T)
norm(m [5 x 5]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: ok; determinant(): ok
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
as(<triangular (ge)matrix>, dtCMatrix): valid: TRUE
> stopifnot(is(T, "dtCMatrix"), identical(T[,3], c(10,10,10,0,0)),
+ Qidentical(as(Z, "matrix"), z))
> showProc.time()
Time elapsed: 0.75 0.002 0.756
>
>
> ## "Vector indices" -------------------
> asLogical <- function(x) {
+ stopifnot(is.atomic(x))
+ storage.mode(x) <- "logical"
+ x
+ }
> .iniDiag.example <- expression({
+ D <- Diagonal(6)
+ M <- as(D,"dgeMatrix")
+ m <- as(D,"matrix")
+ s <- as(D,"TsparseMatrix"); N <- as(s,"nMatrix")
+ S <- as(s,"CsparseMatrix"); C <- as(S,"nMatrix")
+ })
> eval(.iniDiag.example)
> i <- c(3,1,6); v <- c(10,15,20)
> ## (logical,value) which both are recycled:
> L <- c(TRUE, rep(FALSE,8)) ; z <- c(50,99)
>
> ## vector subassignment, both with integer & logical
> ## these now work correctly {though not very efficiently; hence warnings}
> m[i] <- v # the role model: only first column is affected
> M[i] <- v; assert.EQ.mat(M,m) # dge
> D[i] <- v; assert.EQ.mat(D,m) # ddi -> dtT -> dgT
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
> s[i] <- v; assert.EQ.mat(s,m) # dtT -> dgT
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
> S[i] <- v; assert.EQ.mat(S,m); S # dtC -> dtT -> dgT -> dgC
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
6 x 6 sparse Matrix of class "dgCMatrix"
[1,] 15 . . . . .
[2,] . 1 . . . .
[3,] 10 . 1 . . .
[4,] . . . 1 . .
[5,] . . . . 1 .
[6,] 20 . . . . 1
> m.L <- asLogical(m)
> C[i] <- v; assert.EQ.mat(C,m.L); validObject(C)
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
Warning in `[<-`(`*tmp*`, i, value = c(10, 15, 20)) :
x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
[1] TRUE
> N[i] <- v; assert.EQ.mat(N,m.L); validObject(N)
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
Warning in `[<-`(`*tmp*`, i, value = c(10, 15, 20)) :
x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
[1] TRUE
> stopifnot(Q.C.identical(D,s, checkClass=FALSE))
> ## logical *vector* indexing
> eval(.iniDiag.example)
> m[L] <- z; m.L <- asLogical(m)
> M[L] <- z; assert.EQ.mat(M,m)
> D[L] <- z; assert.EQ.mat(D,m)
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
> s[L] <- z; assert.EQ.mat(s,m)
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
> S[L] <- z; assert.EQ.mat(S,m) ; S
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
6 x 6 sparse Matrix of class "dgCMatrix"
[1,] 50 . . 50 . .
[2,] . 1 . . . .
[3,] . . 1 . . .
[4,] . 99 . 1 99 .
[5,] . . . . 1 .
[6,] . . . . . 1
> C[L] <- z; assert.EQ.mat(C,m.L)
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
Warning in `[<-`(`*tmp*`, i, value = c(50, 99)) :
x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
> N[L] <- z; assert.EQ.mat(N,m.L)
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
Warning in `[<-`(`*tmp*`, L, value = c(50, 99)) :
x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
>
>
> ## indexing [i] vs [i,] --- now ok
> eval(.iniDiag.example)
> stopifnot(identical5(m[i], M[i], D[i], s[i], S[i]), identical3(as.logical(m[i]), C[i], N[i]),
+ identical5(m[L], M[L], D[L], s[L], S[L]), identical3(as.logical(m[L]), C[L], N[L]))
<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
> ## bordercase ' drop = .' *vector* indexing {failed till 2009-04-..)
> stopifnot(identical5(m[i,drop=FALSE], M[i,drop=FALSE], D[i,drop=FALSE],
+ s[i,drop=FALSE], S[i,drop=FALSE]),
+ identical3(as.logical(m[i,drop=FALSE]),
+ C[i,drop=FALSE], N[i,drop=FALSE]))
> stopifnot(identical5(m[L,drop=FALSE], M[L,drop=FALSE], D[L,drop=FALSE],
+ s[L,drop=FALSE], S[L,drop=FALSE]),
+ identical3(as.logical(m[L,drop=FALSE]),
+ C[L,drop=FALSE], N[L,drop=FALSE]))
> ## using L for row-indexing should give an error
> assertError(m[L,]); assertError(m[L,, drop=FALSE])
> ## these did not signal an error, upto (including) 0.999375-30:
> assertError(s[L,]); assertError(s[L,, drop=FALSE])
> assertError(S[L,]); assertError(S[L,, drop=FALSE])
> assertError(N[L,]); assertError(N[L,, drop=FALSE])
>
> ## row indexing:
> assert.EQ.mat(D[i,], m[i,])
> assert.EQ.mat(M[i,], m[i,])
> assert.EQ.mat(s[i,], m[i,])
> assert.EQ.mat(S[i,], m[i,])
> assert.EQ.mat(C[i,], asLogical(m[i,]))
> assert.EQ.mat(N[i,], asLogical(m[i,]))
> ## column indexing:
> assert.EQ.mat(D[,i], m[,i])
> assert.EQ.mat(M[,i], m[,i])
> assert.EQ.mat(s[,i], m[,i])
> assert.EQ.mat(S[,i], m[,i])
> assert.EQ.mat(C[,i], asLogical(m[,i]))
> assert.EQ.mat(N[,i], asLogical(m[,i]))
>
>
> ### --- negative indices ----------
>
> ## 1) negative *vector* indexing
> eval(.iniDiag.example)
> i <- -(2:30)
> stopifnot(identical5(m[i], M[i], D[i], s[i], S[i]),
+ identical3(as.logical(m[i]), C[i], N[i]))
> ## negative vector subassignment :
> v <- seq_along(m[i])
> m[i] <- v; m.L <- asLogical(m)
> M[i] <- v; assert.EQ.mat(M,m) # dge
> D[i] <- v; assert.EQ.mat(D,m) # ddi -> dtT -> dgT
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
> s[i] <- v; assert.EQ.mat(s,m) # dtT -> dgT
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
> S[i] <- v; assert.EQ.mat(S,m); S # dtC -> dtT -> dgT -> dgC
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
6 x 6 sparse Matrix of class "dgCMatrix"
[1,] 1 . . . . 2
[2,] . 1 . . . 3
[3,] . . 1 . . 4
[4,] . . . 1 . 5
[5,] . . . . 1 6
[6,] . . . . . 7
> N[i] <- v; assert.EQ.mat(N,m.L); N #
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
Warning in `[<-`(`*tmp*`, i, value = 1:7) :
x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
6 x 6 sparse Matrix of class "ngTMatrix"
[1,] | . . . . |
[2,] . | . . . |
[3,] . . | . . |
[4,] . . . | . |
[5,] . . . . | |
[6,] . . . . . |
> C[i] <- v; assert.EQ.mat(C,m.L); C #
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ntTMatrix to ngTMatrix
Warning in `[<-`(`*tmp*`, i, value = 1:7) :
x[.] <- val: x is "ngTMatrix", val not in {TRUE, FALSE} is coerced.
6 x 6 sparse Matrix of class "ngCMatrix"
[1,] | . . . . |
[2,] . | . . . |
[3,] . . | . . |
[4,] . . . | . |
[5,] . . . . | |
[6,] . . . . . |
>
>
> ## 2) negative [i,j] indices
> mc <- mC[1:5, 1:7]
> mt <- mT[1:5, 1:7]
> ## sub matrix
> assert.EQ.mat(mC[1:2, 0:3], mm[1:2, 0:3]) # test 0-index
> stopifnot(identical(mc[-(3:5), 0:2], mC[1:2, 0:2]),
+ identical(mt[-(3:5), 0:2], mT[1:2, 0:2]),
+ identical(mC[2:3, 4], mm[2:3, 4]))
> assert.EQ.mat(mC[1:2,], mm[1:2,])
> ## sub vector
> stopifnot(identical4(mc[-(1:4), ], mC[5, 1:7],
+ mt[-(1:4), ], mT[5, 1:7]))
> stopifnot(identical4(mc[-(1:4), -(2:4)], mC[5, c(1,5:7)],
+ mt[-(1:4), -(2:4)], mT[5, c(1,5:7)]))
>
> ## mixing of negative and positive must give error
> assertError(mT[-1:1,])
> showProc.time()
Time elapsed: 0.183 0 0.184
>
> ## Sub *Assignment* ---- now works (partially):
> mt0 <- mt
> nt <- as(mt, "nMatrix")
> mt[1, 4] <- -99
> mt[2:3, 1:6] <- 0
> mt
5 x 7 sparse Matrix of class "dgTMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . . -99 . . 241
r2 . . . . . . .
r3 . . . . . . 243
r4 . . . . . . .
r5 . 45 . . . . .
> m2 <- mt+mt
> m2[1,4] <- -200
> m2[c(1,3), c(5:6,2)] <- 1:6
> stopifnot(m2[1,4] == -200,
+ as.vector(m2[c(1,3), c(5:6,2)]) == 1:6)
> mt[,3] <- 30
> mt[2:3,] <- 250
> mt[1:5 %% 2 == 1, 3] <- 0
> mt[3:1, 1:7 > 5] <- 0
> mt
5 x 7 sparse Matrix of class "dgTMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . . -99 . . .
r2 250 250 250 250 250 . .
r3 250 250 . 250 250 . .
r4 . . 30 . . . .
r5 . 45 . . . . .
>
> tt <- as(mt,"matrix")
> ii <- c(0,2,5)
> jj <- c(2:3,5)
> tt[ii, jj] <- 1:6 # 0 is just "dropped"
> mt[ii, jj] <- 1:6
> assert.EQ.mat(mt, tt)
>
> mt[1:5, 2:6]
5 x 5 sparse Matrix of class "dgTMatrix"
c2 c3 c4 c5 c6
r1 . . -99 . .
r2 1 3 250 5 .
r3 250 . 250 250 .
r4 . 30 . . .
r5 2 4 . 6 .
> as((mt0 - mt)[1:5,], "dsparseMatrix")# [1,5] and lines 2:3
5 x 7 sparse Matrix of class "dgCMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . . 220 . . 241
r2 -250 41 -3 -250 -5 202 .
r3 -247 -250 . -250 -250 . 243
r4 . . -30 . . . .
r5 . 43 -4 . -6 . .
>
> mt[c(2,4), ] <- 0; stopifnot(as(mt[c(2,4), ],"matrix") == 0)
> mt[2:3, 4:7] <- 33
> checkMatrix(mt)
norm(m [5 x 7]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: ok
m >= m for all: ok
m < m for none: ok
> mt
5 x 7 sparse Matrix of class "dgTMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . . -99 . . .
r2 . . . 33 33 33 33
r3 250 250 . 33 33 33 33
r4 . . . . . . .
r5 . 2 4 . 6 . .
>
> mc[1,4] <- -99 ; stopifnot(mc[1,4] == -99)
> mc[1,4] <- 00 ; stopifnot(mc[1,4] == 00)
> mc[1,4] <- -99 ; stopifnot(mc[1,4] == -99)
> mc[1:2,4:3] <- 4:1; stopifnot(as.matrix(mc[1:2,4:3]) == 4:1)
>
> mc[-1, 3] <- -2:1 # 0 should not be entered; 'value' recycled
> mt[-1, 3] <- -2:1
> stopifnot(mc@x != 0, mt@x != 0,
+ mc[-1,3] == -2:1, mt[-1,3] == -2:1) ## failed earlier
>
> mc0 <- mc
> mt0 <- as(mc0, "TsparseMatrix")
> m0 <- as(mc0, "matrix")
> set.seed(1); options(Matrix.verbose = FALSE)
> for(i in 1:(if(doExtras) 50 else 4)) {
+ mc <- mc0; mt <- mt0 ; m <- m0
+ ev <- 1:5 %% 2 == round(runif(1))# 0 or 1
+ j <- sample(ncol(mc), 1 + round(runif(1)))
+ nv <- rpois(sum(ev) * length(j), lambda = 1)
+ mc[ev, j] <- nv
+ m[ev, j] <- nv
+ mt[ev, j] <- nv
+ if(i %% 10 == 1) print(mc[ev,j, drop = FALSE])
+ stopifnot(as.vector(mc[ev, j]) == nv, ## failed earlier...
+ as.vector(mt[ev, j]) == nv)
+ validObject(mc) ; assert.EQ.mat(mc, m)
+ validObject(mt) ; assert.EQ.mat(mt, m)
+ }
2 x 1 sparse Matrix of class "dgCMatrix"
c5
r2 2
r4 .
> showProc.time()
Time elapsed: 0.157 0.001 0.158
> options(Matrix.verbose = TRUE)
>
> mc # no longer has non-structural zeros
5 x 7 sparse Matrix of class "dgCMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . 2 4 . . 241
r2 . 42 . 3 . 202 .
r3 3 . -1 . . . 243
r4 . . 1 . . . .
r5 . 45 1 . . . .
> mc[ii, jj] <- 1:6
> mc[c(2,5), c(3,5)] <- 3.2
> checkMatrix(mc)
norm(m [5 x 7]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
> m. <- mc
> mc[4,] <- 0
> mc
5 x 7 sparse Matrix of class "dgCMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . 2.0 4 . . 241
r2 . 1 3.2 3 3.2 202 .
r3 3 . -1.0 . . . 243
r4 . . . . . . .
r5 . 2 3.2 . 3.2 . .
>
> S <- as(Diagonal(5),"TsparseMatrix")
> H <- Hilbert(9)
> Hc <- as(round(H, 3), "dsCMatrix")# a sparse matrix with no 0 ...
> (trH <- tril(Hc[1:5, 1:5]))
5 x 5 sparse Matrix of class "dtCMatrix"
[1,] 1.000 . . . .
[2,] 0.500 0.333 . . .
[3,] 0.333 0.250 0.200 . .
[4,] 0.250 0.200 0.167 0.143 .
[5,] 0.200 0.167 0.143 0.125 0.111
> stopifnot(is(trH, "triangularMatrix"), trH@uplo == "L",
+ is(S, "triangularMatrix"))
>
> ## triangular assignment
> ## the slick (but inefficient in case of sparse!) way to assign sub-diagonals:
> ## equivalent to tmp <- `diag<-`(S[,-1], -2:1); S[,-1] <- tmp
> ## which dispatches to (x="TsparseMatrix", i="missing",j="index", value="replValue")
> diag(S[,-1]) <- -2:1 # used to give a wrong warning
M[i,j] <- v : coercing symmetric M[] into non-symmetric
> S <- as(S,"triangularMatrix")
> assert.EQ.mat(S, local({s <- diag(5); diag(s[,-1]) <- -2:1; s}))
>
> trH[c(1:2,4), c(2:3,5)] <- 0 # gave an *error* upto Jan.2008
> trH[ lower.tri(trH) ] <- 0 # ditto, because of callNextMethod()
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class dtTMatrix to dgTMatrix
>
> m <- Matrix(0+1:28, nrow = 4)
> m[-3,c(2,4:5,7)] <- m[ 3, 1:4] <- m[1:3, 6] <- 0
> mT <- as(m, "dgTMatrix")
> stopifnot(identical(mT[lower.tri(mT)],
+ m [lower.tri(m) ]))
<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient
> lM <- upper.tri(mT, diag=TRUE)
> mT[lM] <- 0
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> m[lM] <- 0
> assert.EQ.mat(mT, as(m,"matrix"))
> mT[lM] <- -1:0
diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1)
> m[lM] <- -1:0
> assert.EQ.mat(mT, as(m,"matrix"))
> (mT <- drop0(mT))
4 x 7 sparse Matrix of class "dgCMatrix"
[1,] -1 . . -1 -1 -1 -1
[2,] 2 -1 -1 . . . .
[3,] . . . -1 -1 -1 -1
[4,] 4 . 12 . . . .
>
> i <- c(1:2, 4, 6:7); j <- c(2:4,6)
> H[i,j] <- 0
> (H. <- round(as(H, "sparseMatrix"), 3)[ , 2:7])
9 x 6 sparse Matrix of class "dgCMatrix"
[1,] . . . 0.200 . 0.143
[2,] . . . 0.167 . 0.125
[3,] 0.250 0.200 0.167 0.143 0.125 0.111
[4,] . . . 0.125 . 0.100
[5,] 0.167 0.143 0.125 0.111 0.100 0.091
[6,] . . . 0.100 . 0.083
[7,] . . . 0.091 . 0.077
[8,] 0.111 0.100 0.091 0.083 0.077 0.071
[9,] 0.100 0.091 0.083 0.077 0.071 0.067
> Hc. <- Hc
> Hc.[i,j] <- 0 ## now "works", but setting "non-structural" 0s
> stopifnot(as.matrix(Hc.[i,j]) == 0)
> Hc.[, 1:6]
9 x 6 sparse Matrix of class "dgCMatrix"
[1,] 1.000 . . . 0.200 .
[2,] 0.500 . . . 0.167 .
[3,] 0.333 0.250 0.200 0.167 0.143 0.125
[4,] 0.250 . . . 0.125 .
[5,] 0.200 0.167 0.143 0.125 0.111 0.100
[6,] 0.167 . . . 0.100 .
[7,] 0.143 . . . 0.091 .
[8,] 0.125 0.111 0.100 0.091 0.083 0.077
[9,] 0.111 0.100 0.091 0.083 0.077 0.071
>
> ## an example that failed for a long time
> sy3 <- new("dsyMatrix", Dim = as.integer(c(2, 2)), x = c(14, -1, 2, -7))
> checkMatrix(dm <- kronecker(Diagonal(2), sy3))# now sparse with new kronecker
Note: method with signature 'sparseMatrix#ANY' chosen for function 'kronecker',
target signature 'dtTMatrix#dsyMatrix'.
"ANY#Matrix" would also be valid
norm(m [4 x 4]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: ok
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: ok; determinant(): ok
> dm <- Matrix(as.matrix(dm))# -> "dsyMatrix"
> (s2 <- as(dm, "sparseMatrix"))
4 x 4 sparse Matrix of class "dsCMatrix"
[1,] 14 2 . .
[2,] 2 -7 . .
[3,] . . 14 2
[4,] . . 2 -7
> checkMatrix(st <- as(s2, "TsparseMatrix"))
norm(m [4 x 4]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: suboptimal 'Arith' implementation of 'dsC* o dsC*'
identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: suboptimal 'Arith' implementation of 'dsC* o dsC*'
ok; determinant(): ok
> stopifnot(is(s2, "symmetricMatrix"),
+ is(st, "symmetricMatrix"))
> checkMatrix(s.32 <- st[1:3,1:2]) ## 3 x 2 - and *not* dsTMatrix
norm(m [3 x 2]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: ok
m >= m for all: ok
m < m for none: ok
> checkMatrix(s2.32 <- s2[1:3,1:2])
norm(m [3 x 2]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: identical
m >= m for all: ok
m < m for none: ok
> I <- c(1,4:3)
> stopifnot(is(s2.32, "generalMatrix"),
+ is(s.32, "generalMatrix"),
+ identical(as.mat(s.32), as.mat(s2.32)),
+ identical3(dm[1:3,-1], asD(s2[1:3,-1]), asD(st[1:3,-1])),
+ identical4(2, dm[4,3], s2[4,3], st[4,3]),
+ identical3(diag(dm), diag(s2), diag(st)),
+ is((cI <- s2[I,I]), "dsCMatrix"),
+ is((tI <- st[I,I]), "dsTMatrix"),
+ identical4(as.mat(dm)[I,I], as.mat(dm[I,I]), as.mat(tI), as.mat(cI))
+ )
>
> ## now sub-assign and check for consistency
> ## symmetric subassign should keep symmetry
> st[I,I] <- 0; checkMatrix(st); stopifnot(is(st,"symmetricMatrix"))
norm(m [4 x 4]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: suboptimal 'Arith' implementation of 'dsC* o dsC*'
identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: suboptimal 'Arith' implementation of 'dsC* o dsC*'
ok; determinant(): ok
> s2[I,I] <- 0; checkMatrix(s2); stopifnot(is(s2,"symmetricMatrix"))
norm(m [4 x 4]) : 1 I F M ok
Summary: ok
as(., "nMatrix") giving full nonzero-pattern: ok
2*m =?= m+m: suboptimal 'Arith' implementation of 'dsC* o dsC*'
identical
m >= m for all: ok
m < m for none: ok
symmpart(m) + skewpart(m) == m: suboptimal 'Arith' implementation of 'dsC* o dsC*'
ok; determinant(): ok
> ##
> m <- as.mat(st)
> m[2:1,2:1] <- 4:1
> st[2:1,2:1] <- 4:1
M[i,j] <- v : coercing symmetric M[] into non-symmetric
> s2[2:1,2:1] <- 4:1
> stopifnot(identical(m, as.mat(st)),
+ 1:4 == as.vector(s2[1:2,1:2]),
+ identical(m, as.mat(s2)))
>
> ## now a slightly different situation for 's2' (had bug)
> s2 <- as(dm, "sparseMatrix")
> s2[I,I] <- 0; diag(s2)[2:3] <- -(1:2)
> stopifnot(is(s2,"symmetricMatrix"), diag(s2) == c(0:-2,0))
> t2 <- as(s2, "TsparseMatrix")
> m <- as.mat(s2)
> s2[2:1,2:1] <- 4:1
> t2[2:1,2:1] <- 4:1
M[i,j] <- v : coercing symmetric M[] into non-symmetric
> m[2:1,2:1] <- 4:1
> assert.EQ.mat(t2, m)
> assert.EQ.mat(s2, m)
> ## and the same (for a different s2 !)
> s2[2:1,2:1] <- 4:1
> t2[2:1,2:1] <- 4:1
> assert.EQ.mat(t2, m)# ok
> assert.EQ.mat(s2, m)# failed in 0.9975-8
> showProc.time()
Time elapsed: 0.709 0.001 0.713
>
>
> ## m[cbind(i,j)] <- value: (2-column matrix subassignment):
> m.[ cbind(3:5, 1:3) ] <- 1:3
> stopifnot(m.[3,1] == 1, m.[4,2] == 2)
> nt. <- nt ; nt[rbind(2:3, 3:4, c(3,3))] <- FALSE
> s. <- m. ; m.[cbind(3,4:6)] <- 0 ## assigning 0 where there *is* 0 ..
> stopifnot(identical(nt.,nt), ## should not have changed
+ identical(s., m.))
> x.x[ cbind(2:6, 2:6)] <- 12:16
> stopifnot(isValid(x.x, "dsCMatrix"),
+ 12:16 == as.mat(x.x)[cbind(2:6, 2:6)])
> (ne1 <- (mc - m.) != 0)
5 x 7 sparse Matrix of class "lgCMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 . . : : . . :
r2 . : : : : : .
r3 | . : . . . :
r4 . | | . . . .
r5 . : | . : . .
> stopifnot(identical(ne1, 0 != abs(mc - m.)))
> (ge <- m. >= mc) # contains "=" -> result is dense
5 x 7 Matrix of class "lgeMatrix"
c1 c2 c3 c4 c5 c6 c7
r1 TRUE TRUE TRUE TRUE TRUE TRUE TRUE
r2 TRUE TRUE TRUE TRUE TRUE TRUE TRUE
r3 FALSE TRUE TRUE TRUE TRUE TRUE TRUE
r4 TRUE TRUE TRUE TRUE TRUE TRUE TRUE
r5 TRUE TRUE FALSE TRUE TRUE TRUE TRUE
> ne. <- mc != m. # was wrong (+ warning)
> stopifnot(identical(!(m. < mc), m. >= mc),
+ identical(m. < mc, as(!ge, "sparseMatrix")),
+ identical(ne., drop0(ne1)))
>
> d6 <- Diagonal(6)
> ii <- c(1:2, 4:5)
> d6[cbind(ii,ii)] <- 7*ii
> stopifnot(is(d6, "ddiMatrix"), identical(d6, Diagonal(x=c(7*1:2,1,7*4:5,1))))
>
> for(j in 3:6) { ## even and odd j used to behave differently
+ M <- Matrix(0, j,j); m <- matrix(0, j,j)
+ T <- as(M, "TsparseMatrix")
+ TG <- as(T, "generalMatrix")
+ G <- as(M, "generalMatrix")
+ id <- cbind(1:j,1:j)
+ i2 <- cbind(1:j,j:1)
+ m[id] <- 1:j
+ M[id] <- 1:j ; stopifnot(is(M,"symmetricMatrix"))
+ T[id] <- 1:j ; stopifnot(is(T,"symmetricMatrix"))
+ G[id] <- 1:j
+ TG[id]<- 1:j
+ m[i2] <- 10
+ M[i2] <- 10 ; stopifnot(is(M,"symmetricMatrix"))
+ T[i2] <- 10 ; stopifnot(is(T,"symmetricMatrix"))
+ G[i2] <- 10
+ TG[i2]<- 10
+ ##
+ assert.EQ.mat(M, m)
+ assert.EQ.mat(T, m)
+ assert.EQ.mat(G, m)
+ assert.EQ.mat(TG,m)
+ }
>
>
> ## drop, triangular, ...
> (M3 <- Matrix(upper.tri(matrix(, 3, 3)))) # ltC; indexing used to fail
3 x 3 sparse Matrix of class "ltCMatrix"
[1,] . | |
[2,] . . |
[3,] . . .
> T3 <- as(M3, "TsparseMatrix")
> stopifnot(identical(drop(M3), M3),
+ identical4(drop(M3[,2, drop = FALSE]), M3[,2, drop = TRUE],
+ drop(T3[,2, drop = FALSE]), T3[,2, drop = TRUE]),
+ is(T3, "triangularMatrix"),
+ !is(T3[,2, drop=FALSE], "triangularMatrix")
+ )
>
> (T6 <- as(as(kronecker(Matrix(c(0,0,1,0),2,2), t(T3)), "lMatrix"),
+ "triangularMatrix"))
6 x 6 sparse Matrix of class "ltTMatrix"
[1,] . . . . . .
[2,] . . . | . .
[3,] . . . | | .
[4,] . . . . . .
[5,] . . . . . .
[6,] . . . . . .
> T6[1:4, -(1:3)] # failed (trying to coerce back to ltTMatrix)
4 x 3 sparse Matrix of class "lgTMatrix"
[1,] . . .
[2,] | . .
[3,] | | .
[4,] . . .
> stopifnot(identical(T6[1:4, -(1:3)][2:3, -3],
+ spMatrix(2,2, i=c(1,2,2), j=c(1,1,2), x=rep(TRUE,3))))
>
> M <- Diagonal(4); M[1,2] <- 2
> M. <- as(M, "CsparseMatrix")
> (R <- as(M., "RsparseMatrix"))
4 x 4 sparse Matrix of class "dtRMatrix" (unitriangular)
[1,] 1 2 . .
[2,] . 1 . .
[3,] . . 1 .
[4,] . . . 1
> (Ms <- symmpart(M.))
4 x 4 sparse Matrix of class "dsCMatrix"
[1,] 1 1 . .
[2,] 1 1 . .
[3,] . . 1 .
[4,] . . . 1
> Rs <- as(Ms, "RsparseMatrix")
> stopifnot(isValid(M, "triangularMatrix"),
+ isValid(M.,"triangularMatrix"),
+ isValid(Ms, "dsCMatrix"),
+ isValid(R, "dtRMatrix"),
+ isValid(Rs, "dsRMatrix") )
> stopifnot(dim(M[2:3, FALSE]) == c(2,0),
+ dim(R[2:3, FALSE]) == c(2,0),
+ identical(M [2:3,TRUE], M [2:3,]),
+ identical(M.[2:3,TRUE], M.[2:3,]),
+ identical(R [2:3,TRUE], R [2:3,]),
+ dim(R[FALSE, FALSE]) == c(0,0))
>
> n <- 50000L
> Lrg <- new("dgTMatrix", Dim = c(n,n))
> diag(Lrg) <- 1:n
> dLrg <- as(Lrg, "diagonalMatrix")
> stopifnot(identical(Diagonal(x = 1:n), dLrg))
> diag(dLrg) <- 1 + diag(dLrg)
> Clrg <- as(Lrg,"CsparseMatrix")
> Ctrg <- as(Clrg, "triangularMatrix")
> diag(Ctrg) <- 1 + diag(Ctrg)
> stopifnot(identical(Diagonal(x = 1+ 1:n), dLrg),
+ identical(Ctrg, as(dLrg,"CsparseMatrix")))
>
> cc <- capture.output(show(dLrg))# show(<diag>) used to error for large n
> showProc.time()
Time elapsed: 0.687 0.003 0.692
>
> ## Large Matrix indexing / subassignment
> ## ------------------------------------- (from ex. by Imran Rashid)
> n <- 7000000
> m <- 100000
> nnz <- 20000
>
> set.seed(12)
> f <- sparseMatrix(i = sample(n, size=nnz, replace=TRUE),
+ j = sample(m, size=nnz, replace=TRUE))
> str(f)
Formal class 'ngCMatrix' [package "Matrix"] with 5 slots
..@ i : int [1:20000] 6692226 4657233 4490801 3688935 344371 6380246 2797160 3584813 6553304 2327896 ...
..@ p : int [1:99993] 0 1 1 1 1 1 1 1 1 1 ...
..@ Dim : int [1:2] 6999863 99992
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ factors : list()
> dim(f) # 6999863 x 99992
[1] 6999863 99992
> prod(dim(f)) # 699930301096 == 699'930'301'096 (~ 700'000 millions)
[1] 699930301096
> str(thisCol <- f[,5000])# logi [~ 7 mio....]
logi [1:6999863] FALSE FALSE FALSE FALSE FALSE FALSE ...
> sv <- as(thisCol, "sparseVector")
> str(sv) ## "empty" !
Formal class 'lsparseVector' [package "Matrix"] with 3 slots
..@ x : logi(0)
..@ length: int 6999863
..@ i : int(0)
> validObject(spCol <- f[,5000, drop=FALSE])
[1] TRUE
> ##
> ## *not* identical(): as(spCol, "sparseVector")@length is "double"prec:
> stopifnot(all.equal(as(spCol, "sparseVector"),
+ as(sv, "nsparseVector"), tolerance=0))
> if(doExtras) {#-----------------------------------------------------------------
+ f[,5762] <- thisCol # now "fine" <<<<<<<<<< FIXME uses LARGE objects -- slow --
+ ## is using replCmat() in ../R/Csparse.R, then
+ ## replTmat() in ../R/Tsparse.R
+
+ fx <- sparseMatrix(i = sample(n, size=nnz, replace=TRUE),
+ j = sample(m, size=nnz, replace=TRUE),
+ x = round(10*rnorm(nnz)))
+ class(fx)## dgCMatrix
+ fx[,6000] <- (tC <- rep(thisCol, length=nrow(fx)))# slow (as above)
+ thCol <- fx[,2000]
+ fx[,5762] <- thCol# slow
+ stopifnot(is(f, "ngCMatrix"), is(fx, "dgCMatrix"),
+ identical(thisCol, f[,5762]),# perfect
+ identical(as.logical(fx[,6000]), tC),
+ identical(thCol, fx[,5762]))
+
+ showProc.time()
+ ##
+ cat("checkMatrix() of all: \n---------\n")
+ Sys.setlocale("LC_COLLATE", "C")# to keep ls() reproducible
+ for(nm in ls()) if(is(.m <- get(nm), "Matrix")) {
+ cat(nm, "\n")
+ checkMatrix(.m, verbose = FALSE
+ , doDet = nm != "As" ## <- "As" almost singular <=> det() "ill posed"
+ )
+ }
+ showProc.time()
+ }#--------------end if(doExtras) -----------------------------------------------
>
> if(!interactive()) warnings()
NULL
>
> proc.time()
user system elapsed
6.942 0.118 7.167
