https://github.com/cran/Matrix
Tip revision: d2fe595da814057e6d64fbb8c6fb4452e6cf522e authored by Doug and Martin on 28 July 2009, 00:00:00 UTC
version 0.999375-30
version 0.999375-30
Tip revision: d2fe595
indexing.Rout.save
R version 2.9.1 Patched (2009-07-18 r48959)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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
>
> library(Matrix)
Loading required package: lattice
Attaching package: 'Matrix'
The following object(s) are masked from package:stats :
contr.helmert,
contr.poly,
contr.SAS,
contr.sum,
contr.treatment,
xtabs
The following object(s) are masked from package:base :
rcond
>
> source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
>
> if(interactive()) {
+ options(error = recover, 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]))
Note: Method with signature "Matrix#index#missing#missing" chosen for function "[",
target signature "dgeMatrix#logical#missing#missing".
"Matrix#logical#missing#ANY" would also be valid
Note: Method with signature "denseMatrix#index#missing#logical" chosen for function "[",
target signature "dgeMatrix#logical#missing#logical".
"Matrix#logical#missing#ANY" would also be valid
>
> ## 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
> ## Printing sparse colnames:
> ms <- as(mn,"sparseMatrix")
> 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(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]),
+ ## 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)])
+ )
>
> ## more sparse Matrices --------------------------------------
>
> m <- 1:800
> set.seed(101) ; m[sample(800, 600)] <- 0
> m <- Matrix(m, nrow = 40)
> mm <- as(m, "matrix")
> dimnames(mm) <- NULL ## << workaround: as(<sparse>, "matrix") has NULL dimnames
> str(mC <- as(m, "dgCMatrix"))
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()
> str(mT <- as(m, "dgTMatrix"))
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()
> stopifnot(identical(mT, as(mC, "dgTMatrix")),
+ identical(mC, as(mT, "dgCMatrix")))
>
> mC[,1]
[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
> mC[1:2,]
2 x 20 sparse Matrix of class "dgCMatrix"
[1,] . . . 121 . . 241 . . . . 441 . . 561 . 641 . . .
[2,] . 42 . . . 202 . . . . . . 482 522 . . . . 722 .
> mC[7, drop = FALSE]
[1] 7
> assert.EQ.mat(mC[1:2,], mm[1:2,])
>
> ## *repeated* (aka 'duplicated') indices - did not work at all ...
> i <- 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)
> for(n in 1:50) {
+ i <- sample(sample(nrow(mC), 7), 20, replace = TRUE)
+ j <- sample(sample(ncol(mC), 6), 17, replace = TRUE)
+ assert.EQ.mat(mC[i,j], mm[i,j])
+ }
>
> ##---- Symmetric indexing of symmetric Matrix ----------
> m. <- mC; m.[, c(2, 7:12)] <- 0
> validObject(S <- crossprod(add.simpleDimnames(m.) %% 100))
[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,])
Note: Method with signature "Matrix#index#missing#missing" chosen for function "[",
target signature "dsyMatrix#logical#missing#missing".
"Matrix#logical#missing#ANY" would also be valid
Note: Method with signature "denseMatrix#index#missing#logical" chosen for function "[",
target signature "dsyMatrix#logical#missing#logical".
"Matrix#logical#missing#ANY" would also be valid
> assert.EQ.mat(ds[ ,NA], ss[,NA])
> stopifnot(identical(ds[2 ,NA], ss[2,NA]),
+ identical(ds[NA, 1], ss[NA, 1]))
> T <- as(S, "TsparseMatrix")
> ## 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:50) {
+ i <- sample(N, max(2, sample(N,1)), replace = FALSE)
+ validObject(Tii <- T[i,i])
+ stopifnot(is(Tii, "dsTMatrix"), # remained symmetric Tsparse
+ identical(t(Tii), t(T)[i,i]))
+ assert.EQ.mat(Tii, ss[i,i])
+ }
>
> ## 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:100) {
+ 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])
+ }
>
>
> stopifnot(all.equal(mC[,3], mm[,3]),
+ identical(mC[ij], mm[ij]),
+ identical(mC[iN], mm[iN]))
>
> 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"
[1,] . 121
[2,] 42 .
[3,] . .
[4,] . .
[5,] 45 .
[6,] . .
[7,] . .
[8,] . 128
[9,] . 129
[10,] 50 .
[11,] . .
[12,] 52 132
[13,] . 133
[14,] . .
[15,] 55 .
[16,] . .
[17,] . .
[18,] . 138
[19,] . .
[20,] . .
[21,] . 141
[22,] . 142
[23,] 63 .
[24,] . .
[25,] 65 .
[26,] . .
[27,] 67 .
[28,] 68 .
[29,] . .
[30,] . .
[31,] 71 .
[32,] 72 .
[33,] . .
[34,] 74 .
[35,] . .
[36,] 76 .
[37,] . .
[38,] . .
[39,] . 159
[40,] 80 .
[1] TRUE
> stopifnot(all.equal(mT[2,], mm[2,]),
+ ## row or column indexing in combination with t() :
+ identical(mT[2,], t(mT)[,2]),
+ identical(mT[-2,], t(t(mT)[,-2])),
+ 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"
[1,] 1 0 . 1 . 1 1 3 . 3 2 1 6 1 . 2 4 6 5 1
[2,] 0 6 2 1 3 5 7 5 12 14 14 9 11 16 12 13 17 19 19 10
[3,] . 2 6 . 4 2 5 3 8 12 5 16 9 11 23 . . 6 7 7
[4,] 1 1 . 17 . 8 10 13 8 6 18 18 29 35 14 8 25 10 19 21
[5,] . 3 4 . 14 4 10 . . 29 8 9 19 11 11 . . 26 26 16
[6,] 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"
[4,] 1 1 . 17 . 8 10 13 8 6 18 18 29 35 14 8 25 10 19 21
[5,] . 3 4 . 14 4 10 . . 29 8 9 19 11 11 . . 26 26 16
[6,] 1 5 2 8 4 42 5 19 14 9 8 10 42 56 50 27 29 32 64 16
[7,] 1 7 5 10 10 5 87 14 9 31 77 47 79 43 28 17 67 110 36 121
[8,] 3 5 3 13 . 19 14 70 10 24 37 13 59 62 34 19 58 21 64 44
[9,] . 12 8 8 . 14 9 10 116 41 58 33 33 72 78 43 69 72 75 25
[10,] 3 14 12 6 29 9 31 24 41 167 69 56 99 44 70 24 105 82 85 32
[11,] 2 14 5 18 8 8 77 37 58 69 267 80 86 139 49 105 194 119 122 129
[12,] 1 9 16 18 9 10 47 13 33 56 80 194 70 77 81 . 90 32 . 106
[13,] 6 11 9 29 19 42 79 59 33 99 86 70 324 157 55 . 69 142 144 155
[14,] 1 16 11 35 11 56 43 62 72 44 139 77 157 375 123 102 145 39 196 81
[15,] . 12 23 14 11 50 28 34 78 70 49 81 55 123 368 71 112 41 41 86
[16,] 2 13 . 8 . 27 17 19 43 24 105 . . 102 71 233 124 44 139 .
[17,] 4 17 . 25 . 29 67 58 69 105 194 90 69 145 112 124 523 141 245 100
[18,] 6 19 6 10 26 32 110 21 72 82 119 32 142 39 41 44 141 497 104 111
[19,] 5 19 7 19 26 64 36 64 75 85 122 . 144 196 41 139 245 104 542 55
[20,] 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"
[1,] | | . | . | | | . |
[2,] | | | | | | | | | |
[3,] . | | . | | | | | |
[4,] | | . | . | | | | |
[5,] . | | . | | | . . |
[6,] | | | | | | | | | |
[7,] | | | | | | | | | |
[8,] | | | | . | | | | |
[9,] . | | | . | | | | |
[10,] | | | | | | | | | |
> (l3 <- lx.x[1:3, ])
3 x 20 sparse Matrix of class "lgCMatrix"
[1,] | | . | . | | | . | | | | | . | | | | |
[2,] | | | | | | | | | | | | | | | | | | | |
[3,] . | | . | | | | | | | | | | | . . | | |
> 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)
> diag(B.) <- 10 * diag(B.)
> 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.)))
>
>
> ## used to fail
> n <- 5 ## or much larger
> sm <- new("dsTMatrix", i=as.integer(1),j=as.integer(1),
+ 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)
+ )
>
> 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
5 x 5 sparse Matrix of class "dtTMatrix"
[1,] 1 . . . .
[2,] . 1 . . .
[3,] . . 1 . .
[4,] . . . 1 .
[5,] . . . . 1
> (m2 <- as(m0, "CsparseMatrix")) # dtCMatrix
5 x 5 sparse Matrix of class "dtCMatrix"
[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"
[1,] 1 . .
[2,] . 1 .
[3,] . . 1
> uTr[1,] <- 0
> assert.EQ.mat(uTr, cbind(0, rbind(0,diag(2))))
>
> M <- m0; M[1,] <- 0
> stopifnot(identical(M, Diagonal(x=c(0, rep(1,4)))))
> M <- m0; M[,3] <- 3 ; M ; stopifnot(is(M, "sparseMatrix"), M[,3] == 3)
M[i,j] <- v : coercing symmetric M[] into non-symmetric
5 x 5 sparse Matrix of class "dgTMatrix"
[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: 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
det...(): 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))),
+ is(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
> 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: ok
m >= m for all: ok
m < m for none: ok
det...(): 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
det...(): ok
.TM.repl.i.2col(): drop 'matrix' case ...
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, "dtTMatrix"), identical(T[,3], c(10,10,10,0,0)))
>
> 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
det...(): ok
> 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
det...(): ok
.TM.repl.i.2col(): drop 'matrix' case ...
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)))
>
>
> ## "Vector indices" -------------------
> .iniDiag.example <- expression({
+ D <- Diagonal(6)
+ M <- as(D,"dgeMatrix")
+ m <- as(D,"matrix")
+ s <- as(D,"TsparseMatrix")
+ S <- as(s,"CsparseMatrix")
+ })
> 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
> stopifnot(identical(s,D))
> ## logical
> eval(.iniDiag.example)
> m[L] <- z
> 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
>
> ## indexing [i] vs [i,] --- now ok
> eval(.iniDiag.example)
> stopifnot(identical5(m[i], M[i], D[i], s[i], S[i]))
> stopifnot(identical5(m[L], M[L], D[L], s[L], S[L]))
Note: Method with signature "Matrix#index#missing#missing" chosen for function "[",
target signature "ddiMatrix#logical#missing#missing".
"Matrix#logical#missing#ANY" would also be valid
Note: Method with signature "Matrix#index#missing#missing" chosen for function "[",
target signature "dtTMatrix#logical#missing#missing".
"Matrix#logical#missing#ANY" would also be valid
Note: Method with signature "Matrix#index#missing#missing" chosen for function "[",
target signature "dtCMatrix#logical#missing#missing".
"Matrix#logical#missing#ANY" would also be valid
> ## 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]))
> stopifnot(identical5(m[L,drop=FALSE], M[L,drop=FALSE], D[L,drop=FALSE],
+ s[L,drop=FALSE], S[L,drop=FALSE]))
Note: Method with signature "diagonalMatrix#index#missing#logical" chosen for function "[",
target signature "ddiMatrix#logical#missing#logical".
"Matrix#logical#missing#ANY" would also be valid
Note: Method with signature "TsparseMatrix#index#missing#logical" chosen for function "[",
target signature "dtTMatrix#logical#missing#logical".
"Matrix#logical#missing#ANY" would also be valid
Note: Method with signature "sparseMatrix#index#missing#logical" chosen for function "[",
target signature "dtCMatrix#logical#missing#logical".
"Matrix#logical#missing#ANY" would also be valid
> ##
> 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(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])
>
>
> ## --- negative 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,])
>
> ## Sub *Assignment* ---- now works (partially):
> mt0 <- mt
> mt[1, 4] <- -99
> mt[2:3, 1:6] <- 0
> mt
5 x 7 sparse Matrix of class "dgTMatrix"
[1,] . . . -99 . . 241
[2,] . . . . . . .
[3,] . . . . . . 243
[4,] . . . . . . .
[5,] . 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"
[1,] . . . -99 . . .
[2,] 250 250 250 250 250 . .
[3,] 250 250 . 250 250 . .
[4,] . . 30 . . . .
[5,] . 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"
[1,] . . -99 . .
[2,] 1 3 250 5 .
[3,] 250 . 250 250 .
[4,] . 30 . . .
[5,] 2 4 . 6 .
> as((mt0 - mt)[1:5,], "dsparseMatrix")# [1,5] and lines 2:3
5 x 7 sparse Matrix of class "dgCMatrix"
[1,] . . . 220 . . 241
[2,] -250 41 -3 -250 -5 202 .
[3,] -247 -250 . -250 -250 . 243
[4,] . . -30 . . . .
[5,] . 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"
[1,] . . . -99 . . .
[2,] . . . 33 33 33 33
[3,] 250 250 . 33 33 33 33
[4,] . . . . . . .
[5,] . 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)
> for(i in 1:50) {
+ 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"
[1,] 2
[2,] .
2 x 1 sparse Matrix of class "dgCMatrix"
[1,] 2
[2,] 1
3 x 2 sparse Matrix of class "dgCMatrix"
[1,] 1 .
[2,] . .
[3,] 1 .
3 x 1 sparse Matrix of class "dgCMatrix"
[1,] 1
[2,] 1
[3,] 1
3 x 1 sparse Matrix of class "dgCMatrix"
[1,] .
[2,] 3
[3,] 1
>
> mc # no longer has non-structural zeros
5 x 7 sparse Matrix of class "dgCMatrix"
[1,] . . 2 4 . . 241
[2,] 1 42 -2 3 . 1 .
[3,] 3 . -1 . . . 243
[4,] 1 . . . . 1 .
[5,] . 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"
[1,] . . 2.0 4 . . 241
[2,] 1 1 3.2 3 3.2 1 .
[3,] 3 . -1.0 . . . 243
[4,] . . . . . . .
[5,] . 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()
.TM.repl.i.2col(): drop 'matrix' case ...
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
.TM.repl.i.2col(): drop 'matrix' case ...
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
.TM.repl.i.2col(): drop 'matrix' case ...
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 "diagonalMatrix#ANY" chosen for function "kronecker",
target signature "ddiMatrix#dsyMatrix".
"ANY#Matrix" would also be valid
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
det...(): 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 implementation of sparse 'symm. o symm.'
identical
m >= m for all: ok
m < m for none: ok
det...(): 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 implementation of sparse 'symm. o symm.'
identical
m >= m for all: ok
m < m for none: ok
det...(): 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 implementation of sparse 'symm. o symm.'
identical
m >= m for all: ok
m < m for none: ok
det...(): 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
>
>
> ## m[cbind(i,j)] <- value:
> m.[ cbind(3:5, 1:3) ] <- 1:3
> stopifnot(m.[3,1] == 1, m.[4,2] == 2)
> 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"
[1,] . . . . . . .
[2,] . . . . . . .
[3,] | . . . . . .
[4,] | | . . . | .
[5,] . . | . . . .
> stopifnot(identical(ne1, 0 != abs(mc - m.)))
> (ge <- m. >= mc) # contains "=" -> result is dense
5 x 7 Matrix of class "lgeMatrix"
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[2,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[3,] FALSE TRUE TRUE TRUE TRUE TRUE TRUE
[4,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[5,] 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"
[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
>
> ## 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"), tol=0))
> f[,5762] <- thisCol # now "fine" <<<<<<<<<< FIXME uses LARGE objects
> ## 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
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"
> fx[,6000] <- (tC <- rep(thisCol, length=nrow(fx)))
> thCol <- fx[,2000]
> fx[,5762] <- thCol
> stopifnot(is(f, "ngCMatrix"), is(fx, "dgCMatrix"),
+ identical(thisCol, f[,5762]),# perfect
+ identical(as.logical(fx[,6000]), tC),
+ identical(thCol, fx[,5762]))
>
> cat('Time elapsed: ', (.pt <- proc.time()),'\n') # "stats"
Time elapsed: 21.691 1.553 25.207 0 0
> ##
> cat("checkMatrix() of all: \n---------\n")
checkMatrix() of all:
---------
> Sys.setlocale("LC_COLLATE", "C")# to keep ls() reproducible
[1] "C"
> for(nm in ls()) if(is(.m <- get(nm), "Matrix")) {
+ cat(nm, "\n")
+ checkMatrix(.m, verbose = FALSE)
+ }
A
A.
suboptimal implementation of sparse 'symm. o symm.'
A.ii
suboptimal implementation of sparse 'symm. o symm.'
Aii
suboptimal implementation of sparse 'symm. o symm.'
B
B.
C
Clrg
Ctrg
D
exploding <diag> o <diag> into dense matrix
Note: Method with signature "nMatrix#Matrix" chosen for function "|",
target signature "nsyMatrix#nsCMatrix".
"Matrix#nMatrix" would also be valid
exploding <diag> o <diag> into dense matrix
G
H
Note: Method with signature "sparseMatrix#ldiMatrix" chosen for function "==",
target signature "ngCMatrix#ldiMatrix".
"nsparseMatrix#sparseMatrix", "nMatrix#lMatrix" would also be valid
H.
Hc
suboptimal implementation of sparse 'symm. o symm.'
Hc.
Lrg
M
.TM.repl.i.2col(): drop 'matrix' case ...
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.
.TM.repl.i.2col(): drop 'matrix' case ...
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
M3
.TM.repl.i.2col(): drop 'matrix' case ...
diagnosing replTmat(x,i,j,v): nargs()= 3; missing(i,j)= (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ltTMatrix to lgTMatrix
Ms
suboptimal implementation of sparse 'symm. o symm.'
R
Rs
suboptimal implementation of sparse 'symm. o symm.'
S
.TM.repl.i.2col(): drop 'matrix' case ...
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
T
suboptimal implementation of sparse 'symm. o symm.'
T3
.TM.repl.i.2col(): drop 'matrix' case ...
diagnosing replTmat(x,i,j,v): nargs()= 3; missing(i,j)= (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ltTMatrix to lgTMatrix
T6
.TM.repl.i.2col(): drop 'matrix' case ...
diagnosing replTmat(x,i,j,v): nargs()= 3; missing(i,j)= (0,1)
'sub-optimal sparse 'x[i] <- v' assignment: Coercing class ltTMatrix to lgTMatrix
TG
Tii
suboptimal implementation of sparse 'symm. o symm.'
cI
suboptimal implementation of sparse 'symm. o symm.'
cm
suboptimal implementation of sparse 'symm. o symm.'
d6
exploding <diag> o <diag> into dense matrix
exploding <diag> o <diag> into dense matrix
dLrg
dm
Note: Method with signature "sparseMatrix#ldiMatrix" chosen for function "&",
target signature "ngCMatrix#ldiMatrix".
"nsparseMatrix#sparseMatrix", "nMatrix#lMatrix" would also be valid
ds
f
fx
ge
l10
suboptimal implementation of sparse 'symm. o symm.'
l3
lx.x
suboptimal implementation of sparse 'symm. o symm.'
m.
m1
.TM.repl.i.2col(): drop 'matrix' case ...
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
m1g
m2
mC
mT
mc
mc0
mn
mo
ms
mt
mt0
ne.
ne1
nsel
s
.TM.repl.i.2col(): drop 'matrix' case ...
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.32
s2
s2.32
sel
sm
suboptimal implementation of sparse 'symm. o symm.'
spCol
ssel
st
sy3
t2
tI
suboptimal implementation of sparse 'symm. o symm.'
trH
uTr
.TM.repl.i.2col(): drop 'matrix' case ...
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
x.x
suboptimal implementation of sparse 'symm. o symm.'
x.x.
suboptimal implementation of sparse 'symm. o symm.'
> cat('Time elapsed: ', proc.time() - .pt,'\n') # "stats"
Time elapsed: 13.831 0.68 14.511 0 0
>
> if(!interactive()) warnings()
NULL
>
>
