SCM

SCM Repository

[matrix] Annotation of /pkg/tests/simple.R
ViewVC logotype

Annotation of /pkg/tests/simple.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1548 - (view) (download)

1 : maechler 948 #### Currently a collection of simple tests
2 :     ## (since 'Matrix' takes long to load, rather have fewer source files!)
3 :    
4 : bates 346 library(Matrix)
5 : maechler 509
6 : maechler 948 source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
7 :    
8 : maechler 1174 ### Matrix() ''smartness''
9 :     (d4 <- Matrix(diag(4)))
10 :     (z4 <- Matrix(0*diag(4)))
11 :     (o4 <- Matrix(1+diag(4)))
12 :     (m4 <- Matrix(cbind(0,rbind(6*diag(3),0))))
13 :     dm4 <- Matrix(m4, sparse = FALSE)
14 : maechler 1467 class(mN <- Matrix(NA, 3,4)) # NA *is* logical
15 : maechler 1174 stopifnot(validObject(d4), validObject(z4), validObject(o4),
16 : maechler 1467 validObject(m4), validObject(dm4), validObject(mN))
17 : maechler 1174 assert.EQ.mat(dm4, as(m4, "matrix"))
18 : maechler 1467 if(FALSE)# assert.EQ.. fails for all NA
19 :     assert.EQ.mat(mN, matrix(NA, 3,4))
20 : maechler 1548 sL <- Matrix(, 3,4, sparse=TRUE)# -> "lgC
21 :     stopifnot(##length(sN@i) == 0, # all "FALSE"
22 : maechler 1467 validObject(Matrix(c(NA,0), 4, 3, byrow = TRUE)),
23 :     validObject(Matrix(c(NA,0), 4, 4)),
24 :     is(Matrix(c(NA,0,0,0), 4, 4), "sparseMatrix"))
25 :    
26 : maechler 1315 ## large sparse ones: these now directly "go sparse":
27 :     str(m0 <- Matrix(0, nrow=100, ncol = 1000))
28 :     str(l0 <- Matrix(FALSE, nrow=100, ncol = 200))
29 : maechler 1174
30 : maechler 948 ###-- Sparse Triangular :
31 :    
32 :     (t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2,
33 :     Dim= as.integer(c(4,4))))
34 :     stopifnot(validObject(t1),
35 :     validObject(t1c <- as(t1, "dtCMatrix")))
36 :     assert.EQ.mat(t1, as(t1c, "matrix"))
37 :    
38 :     ## from 0-diagonal to unit-diagonal {low-level step}:
39 :     tu <- t1 ; tu@diag <- "U"
40 :     tu
41 :     stopifnot(validObject(cu <- as(tu, "dtCMatrix")),
42 : maechler 1253 validObject(tu. <- as(cu, "dtTMatrix")),
43 :     ## NOT: identical(tu, tu.), # since T* is not unique!
44 :     identical(cu, as(tu., "dtCMatrix")),
45 :     validObject(t(cu)),
46 :     validObject(t(tu)))
47 : maechler 948 assert.EQ.mat(cu, as(tu,"matrix"), tol=0)
48 :    
49 :    
50 :     ###-- Numeric Dense: Crossprod & Solve
51 :    
52 : maechler 509 set.seed(123)
53 : maechler 1207 mm. <- mm <- Matrix(rnorm(500 * 150), nc = 150)
54 : maechler 509 stopifnot(validObject(mm))
55 : maechler 1207 xpx <- crossprod(mm)
56 :     stopifnot(identical(mm, mm.),# once upon a time, mm was altered by crossprod()
57 : maechler 873 validObject(xpx))
58 :     str(mm) # 'dge*"
59 :     str(xpx)# 'dpo*"
60 : bates 346 xpy <- crossprod(mm, rnorm(500))
61 :     res <- solve(xpx, xpy)
62 : maechler 873 str(xpx)# now with Cholesky factor
63 :     stopifnot(validObject(xpx),
64 :     validObject(xpy),
65 :     validObject(res))
66 :     stopifnot(all.equal(xpx %*% res, xpy, tol= 1e-12))
67 : maechler 509
68 : maechler 948 ###-- more solve() methods {was ./solve.R }
69 :    
70 :     ## first for "dgeMatrix" and all kinds of RHS :
71 :     (m6 <- 1 + as(diag(0:5), "dgeMatrix"))
72 :     rcond(m6)
73 :     I6 <- as(diag(6), "dgeMatrix")
74 :     stopifnot(all.equal(I6, m6 %*% solve(m6)),
75 :     all.equal(I6, solve(m6) %*% m6) )
76 :    
77 :     (i6 <- solve(m6, Matrix(1:6)))
78 :     stopifnot(identical(i6, as(cbind(c(-4, rep(1,5))), "dgeMatrix")),
79 : maechler 949 identical(i6, solve(m6, 1:6)),
80 : maechler 948 identical(i6, solve(m6, matrix(1:6))),
81 : maechler 949 identical(i6, solve(m6, matrix(c(1,2,3,4,5,6))))
82 : maechler 948 )
83 :    
84 :     ###-- row- and column operations {was ./rowcolOps.R }
85 :    
86 :     set.seed(321)
87 : maechler 1467 (m1 <- round(Matrix(rnorm(25), 5), 2))
88 : maechler 956 m1k <- Matrix(round(rnorm(1000), 2), 50, 20)
89 :     m.m <- as(m1k, "matrix")
90 :     stopifnot(all.equal(colMeans(m1k), colMeans(m.m)),
91 :     all.equal(colSums (m1k), colSums (m.m)),
92 :     all.equal(rowMeans(m1k), rowMeans(m.m)),
93 :     all.equal(rowSums (m1k), rowSums (m.m)))
94 : maechler 948
95 : maechler 956 ###-- kronecker for nonsparse uses Matrix(.):
96 :     stopifnot(is(kr <- kronecker(m1, m6), "Matrix"))
97 :     assert.EQ.mat(kr,
98 :     kronecker(as(m1, "matrix"),
99 :     as(m6, "matrix")),
100 :     tol = 0)
101 :     ## sparse:
102 :     (kt1 <- kronecker(t1, tu))
103 :     kt2 <- kronecker(t1c, cu)
104 :     ktf <- kronecker(as.matrix(t1), as.matrix(tu))
105 :    
106 :     assert.EQ.mat(kt1, ktf, tol= 0)
107 :     assert.EQ.mat(kt2, ktf, tol= 0)
108 :     ## but kt1 and kt2, both "dgT" are different since entries are not ordered!
109 :    
110 : maechler 1331 ## coercion from "dpo" or "dsy"
111 :     xx <- as(xpx, "dsyMatrix")
112 :     stopifnot(isSymmetric(xxS <- as(xx, "sparseMatrix")),
113 :     isSymmetric(xpxS <- as(xpx, "sparseMatrix")))
114 :    
115 : maechler 1238 tm <- matrix(0, 8,8)
116 :     tm[cbind(c(1,1,2,7,8),
117 :     c(3,6,4,8,8))] <- c(2,-30,15,20,80)
118 :     (tM <- Matrix(tm)) ## dtC
119 :     (mM <- Matrix(m <- (tm + t(tm)))) ## dsC
120 :     mT <- as(mM, "dsTMatrix")
121 : maechler 1319 gC <- as(as(mT, "dgTMatrix"), "dgCMatrix")
122 :     ## Check that 'mT' and gC print properly :
123 :     pr.mT <- capture.output(mT)
124 :     nn <- unlist(strsplit(gsub(" +\\.", "", sub("^....", "", pr.mT[-(1:2)])), " "))
125 :     stopifnot(as.numeric(nn[nn != ""]) == m[m != 0],
126 :     capture.output(gC)[-1] == pr.mT[-1])
127 : maechler 1238 assert.EQ.mat(tM, tm, tol=0)
128 : maechler 1319 assert.EQ.mat(gC, m, tol=0)
129 : maechler 1238 assert.EQ.mat(mT, m, tol=0)
130 :     stopifnot(is(mM, "dsCMatrix"), is(tM, "dtCMatrix"),
131 :     ## coercions general <-> symmetric
132 :     identical(as(as(mM, "dgCMatrix"), "dsCMatrix"), mM),
133 :     identical(as(as(mM, "dgTMatrix"), "dsTMatrix"), mT),
134 :     identical(as(as(tM, "dgCMatrix"), "dtCMatrix"), tM)
135 :     )
136 :     eM <- eigen(mM) # works thanks to base::as.matrix hack in ../R/zzz.R
137 :     stopifnot(all.equal(eM$values,
138 :     { v <- c(162.462112512353, 30.0665927567458)
139 :     c(v, 15, 0, 0, 160-v[1], -15, -v[2])}, tol=1e-14))
140 :    
141 : maechler 974 ##--- symmetric -> pos.def. needs valid test:
142 : maechler 1238 m5 <- Matrix(diag(5) - 1)
143 :     if(FALSE) # FIXME: this happily "works" but MM thinks it shouldn't:
144 :     assertError(as(m5, "dpoMatrix"))
145 : maechler 956
146 : maechler 974
147 : maechler 1548 ###-- sparse nonzero pattern : ----------
148 : maechler 956
149 : maechler 1548 (nkt <- as(as(kt1, "dgCMatrix"), "ngCMatrix"))# ok
150 :     (clt <- crossprod(nkt))
151 : maechler 1238 if(FALSE) ## FIXME!
152 : maechler 956 crossprod(clt)
153 :     ## CHOLMOD error: matrix cannot be symmetric
154 :    
155 : maechler 961
156 : maechler 956 ### "d" <-> "l" for (symmetric) sparse :
157 : maechler 1228 data(KNex)
158 :     mm <- KNex$mm
159 : maechler 956 xpx <- crossprod(mm)
160 : maechler 1548 ## extract nonzero pattern
161 :     if(FALSE) ## FIXME -- {used to work..
162 :     nxpx <- as(xpx, "nsCMatrix")
163 : maechler 956 if(FALSE)
164 : maechler 1548 show(nxpx) ## gives error about "lsC" -> "lgT" coercion ..
165 : maechler 956 ## The bug is actually from *subsetting* the large matrix:
166 :     if(FALSE) ## FIXME
167 : maechler 1548 r <- nxpx[1:2,]
168 : maechler 956
169 :     lmm <- as(mm, "lgCMatrix")
170 : maechler 1548 nmm <- as(lmm, "nMatrix")
171 : maechler 956 xlx <- crossprod(lmm)
172 : maechler 1548 x.x <- crossprod(nmm)
173 : maechler 961 ## now A = lxpx and B = xlx should be close, but not quite the same
174 :     ## since <x,y> = 0 is well possible when x!=0 and y!=0 .
175 :     ## However, A[i,j] != 0 ==> B[i,j] != 0:
176 : maechler 1548 if(FALSE) { ## FIXME : nxpx above
177 :     A <- as(as(nxpx, "lgCMatrix"), "lgTMatrix")
178 : maechler 961 B <- as(as(xlx, "lgCMatrix"), "lgTMatrix")
179 :     ij <- function(a) a@i + ncol(a) * a@j
180 :     stopifnot(all(ij(A) %in% ij(B)))
181 : maechler 1548 }
182 : maechler 956
183 : maechler 974
184 :    
185 :     cat('Time elapsed: ', proc.time(),'\n') # "stats"

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge