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 2185 - (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 : mmaechler 2175 if(interactive()) {
9 :     options(error = recover)
10 :     } else options(verbose = TRUE)# to show message()s
11 : maechler 2113
12 : maechler 1174 ### Matrix() ''smartness''
13 :     (d4 <- Matrix(diag(4)))
14 :     (z4 <- Matrix(0*diag(4)))
15 :     (o4 <- Matrix(1+diag(4)))
16 : maechler 1825 (tr <- Matrix(cbind(1,0:1)))
17 : maechler 1174 (m4 <- Matrix(cbind(0,rbind(6*diag(3),0))))
18 :     dm4 <- Matrix(m4, sparse = FALSE)
19 : maechler 1467 class(mN <- Matrix(NA, 3,4)) # NA *is* logical
20 : maechler 1174 stopifnot(validObject(d4), validObject(z4), validObject(o4),
21 : maechler 1467 validObject(m4), validObject(dm4), validObject(mN))
22 : maechler 1174 assert.EQ.mat(dm4, as(m4, "matrix"))
23 : maechler 1467 assert.EQ.mat(mN, matrix(NA, 3,4))
24 : maechler 1825 sL <- Matrix(, 3,4, sparse=TRUE)# -> "lgC"
25 :     trS <- Matrix(tr, sparse=TRUE)# failed in 0.9975-11
26 : maechler 2157 stopifnot(isValid(d4, "diagonalMatrix"), isValid(z4, "diagonalMatrix"),
27 :     isValid(tr, "triangularMatrix"), isValid(trS, "triangularMatrix"),
28 : maechler 1825 all(is.na(sL@x)), ## not yet: all(is.na(sL)),
29 : maechler 1883 !any(sL, na.rm=TRUE), all(!sL, na.rm=TRUE),
30 : maechler 1467 validObject(Matrix(c(NA,0), 4, 3, byrow = TRUE)),
31 :     validObject(Matrix(c(NA,0), 4, 4)),
32 : maechler 2157 isValid(Matrix(c(NA,0,0,0), 4, 4), "sparseMatrix"))
33 : maechler 1467
34 : maechler 1315 ## large sparse ones: these now directly "go sparse":
35 :     str(m0 <- Matrix(0, nrow=100, ncol = 1000))
36 :     str(l0 <- Matrix(FALSE, nrow=100, ncol = 200))
37 : maechler 1705 stopifnot(all(!l0),
38 :     identical(FALSE, any(l0)))
39 : maechler 1174
40 : maechler 1660 ## with dimnames:
41 :     m. <- matrix(c(0, 0, 2:0), 3, 5)
42 :     dimnames(m.) <- list(LETTERS[1:3], letters[1:5])
43 : maechler 1904 (m0 <- m <- Matrix(m.))
44 : maechler 1660 m@Dimnames[[2]] <- m@Dimnames[[1]]
45 :     ## not valid anymore:
46 : maechler 2158 (val <- validObject(m, test=TRUE)); stopifnot(is.character(val))
47 : maechler 2102 dm <- as(m0, "denseMatrix")
48 : maechler 2158 stopifnot(identical(rcond(dm), rcond(m.)),
49 : maechler 2102 all.equal(rcond(dm), 0.4899474520656))
50 : maechler 1725 rm(m)
51 : maechler 1660
52 : maechler 948 ###-- Sparse Triangular :
53 :    
54 : maechler 2130 g5 <- new("dgCMatrix", Dim = c(5L, 5L),
55 :     x = c(10, 1, 3, 10, 1, 10, 1, 10, 10),
56 :     i = c(0L,2L,4L, 1L, 3L,2L,4L, 3L, 4L),
57 :     p = c(0L, 3L, 5L, 7:9))
58 :     t5 <- as(g5, "triangularMatrix") # works fine (but slowly) FIXME
59 :     stopifnot(class(t5) == "dtCMatrix",
60 :     identical(t5, tril(g5)))
61 :    
62 : maechler 948 (t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2,
63 :     Dim= as.integer(c(4,4))))
64 : maechler 2128 ## Diagonal o Sparse
65 :     I4 <- Diagonal(4)
66 : maechler 2144 D4 <- Diagonal(4, x=1:4)
67 : maechler 2128 validObject(t2 <- t1 + I4)
68 :     validObject(tt2 <- t(t1) + I4)
69 :     validObject(t1c <- as(t1, "CsparseMatrix"))
70 :     validObject(t2c <- as(t2, "CsparseMatrix"))
71 : maechler 948 stopifnot(validObject(t1),
72 : maechler 2157 isValid(2 * I4, "diagonalMatrix"),
73 :     isValid(D4 * 3, "diagonalMatrix"),
74 :     isValid(I4 / 5, "diagonalMatrix"),
75 :     isValid(D4 / 2, "diagonalMatrix"),
76 : maechler 2128 identical(t1, t(t(t1))),
77 :     identical(t1c, t(t(t1c))),
78 : maechler 2157 isValid(t1c + I4,"triangularMatrix"), isValid(t2c + I4,"triangularMatrix"),
79 : maechler 2128 c(class(t2), class(t1c), class(t2c), class(tt2)) == "dtCMatrix",
80 :     identical(t(tt2), t2))
81 : maechler 948 assert.EQ.mat(t1, as(t1c, "matrix"))
82 :    
83 : maechler 2128 ## as(<diag>, <anything>) :
84 :     str(cls <- names(getClass("Matrix")@subclasses))# all Matrix classes
85 :     for(cl in cls)
86 :     if(canCoerce(I4, cl)) {
87 :     cat(cl,":")
88 :     M <- as(I4, cl)
89 :     M. <- as(D4, cl)
90 :     stopifnot(diag(4) == as(M,"matrix"),
91 :     if(is(cl,"dMatrix")) diag(x=1:4) == as(M.,"matrix") else TRUE)
92 :     cat(" [Ok]\n")
93 :     }
94 : mmaechler 2175 s4 <- as(D4,"sparseMatrix")
95 :     v <- c(11,2,2,12); s4[2:3,2:3] <- v; validObject(s4)
96 :     s4. <- D4; s4.[2:3,2:3] <- v; validObject(s4.)
97 :     stopifnot(all(s4 == s4.))
98 :     ## now assign symmetrically to symmetricMatrix
99 :     s4 <- as(as(D4,"sparseMatrix"),"symmetricMatrix")
100 :     s4[2:3,2:3] <- v
101 :     validObject(s4)
102 :     stopifnot(is(s4,"symmetricMatrix"))
103 :     assert.EQ.mat(s4, as(s4.,"matrix"),tol=0)
104 : maechler 2128
105 : mmaechler 2175 ## lower-triangular unit-diagonal
106 :     L <- new("dtCMatrix", i = 1L, p = c(0:1, 1L), Dim = c(2L, 2L),
107 :     x = 0.5, uplo = "L", diag = "U")
108 :     stopifnot(range(L) == 0:1, all.equal(mean(L), 5/8))
109 :    
110 : maechler 2144 ## from 0-diagonal to unit-diagonal triangular {low-level step}:
111 : maechler 948 tu <- t1 ; tu@diag <- "U"
112 :     tu
113 : mmaechler 2175 validObject(cu <- as(tu, "dtCMatrix"))
114 : maechler 2125 validObject(cnu <- Matrix:::diagU2N(cu))# <- testing diagU2N
115 : mmaechler 2175 validObject(tu. <- as(cu, "dtTMatrix"))
116 :     validObject(tt <- as(cu, "TsparseMatrix"))
117 :     stopifnot(## NOT: identical(tu, tu.), # since T* is not unique!
118 : maechler 1253 identical(cu, as(tu., "dtCMatrix")),
119 : maechler 2125 length(cnu@i) == length(cu@i) + nrow(cu),
120 :     identical(cu, Matrix:::diagN2U(cnu)),# <- testing diagN2U
121 : maechler 1883 all(cu >= 0, na.rm = TRUE), all(cu >= 0),
122 : mmaechler 2175 any(cu >= 7))
123 :     validObject(tcu <- t(cu))
124 :     validObject(ttu <- t(tu))
125 :     validObject(ltu <- as(ttu, "lMatrix"))
126 :     validObject(ldtu <- as(ltu, "denseMatrix"))
127 :     validObject(Cltu <- as(ltu, "CsparseMatrix"))
128 :     stopifnot(identical(asCsp(ttu > 0), asCsp(ltu)),
129 :     all(ltu == as(ttu > 0,"denseMatrix")))
130 :     ltu - (ttu > 0) # failed
131 : maechler 2125
132 : mmaechler 2175
133 :     lcu <- new("ltCMatrix", Dim = c(4L, 4L), i = c(0:1, 0L), p = c(0L, 0:3),
134 :     x = c(TRUE, FALSE, FALSE), uplo = "U", diag = "U")
135 : mmaechler 2185 stopifnot(identical(rowSums(lcu), rowSums(drop0(lcu))))
136 :     (ncu <- as(lcu, "nMatrix"))# -- gives the "pattern" of lcu, i.e. FALSE are *there*
137 :     stopifnot(identical(ncu, as(lcu,"nsparseMatrix")),
138 :     identical(rowSums(ncu), c(3:1, 1L)))
139 : mmaechler 2175
140 : mmaechler 2185
141 : maechler 948 assert.EQ.mat(cu, as(tu,"matrix"), tol=0)
142 : maechler 2125 assert.EQ.mat(cnu, as(tu,"matrix"), tol=0)
143 : maechler 1747
144 :     ## <sparse> o <numeric> (of length > 1):
145 : maechler 2157 stopifnot(isValid(tm <- tu * 1:8, "sparseMatrix"),
146 : maechler 1747 identical4(tm, cu * 1:8, 1:8 * cu, 1:8 * tu))
147 :    
148 : maechler 1593 cu[1,2] <- tu[1,2] <- NA
149 : maechler 1751 mu <- as(tu,"matrix")
150 : maechler 2157 stopifnot(isValid(cu, "CsparseMatrix"), isValid(cu, "triangularMatrix"),
151 :     isValid(tu, "TsparseMatrix"), isValid(tu, "triangularMatrix"),
152 : maechler 1751 identical(cu * 1:8, tu * 1:8), # but are no longer triangular
153 : mmaechler 2175 all(cu >= 0, na.rm=TRUE), !all(cu >= 1), is.na(all(tu >= 0)),
154 :     ## Csparse_drop: preserves triangularity incl diag="U"
155 :     identical(cu, .Call(Matrix:::Csparse_drop, cu, 0.))
156 :     )
157 : maechler 1751 assert.EQ.mat(cu * 1:8, mu * 1:8)
158 : maechler 1747
159 : mmaechler 2175 ina <- is.na(as(cu,"matrix"))
160 :     ## These 3 were each different (2008-03) !!
161 :     stopifnot(all(ina == is.na(cu)),
162 :     all(ina == is.na(as(cu,"generalMatrix"))),
163 :     all(ina == as(is.na(as(cu,"matrix")),"nMatrix")))
164 :    
165 :    
166 : maechler 1738 ## tu. is diag "U", but tu2 not:
167 : maechler 2115 tu2 <- as(as(tu., "generalMatrix"), "triangularMatrix")
168 : maechler 1593 assert.EQ.mat(cu, mu, tol=0)
169 :     stopifnot(identical3(cu[cu > 1], tu [tu > 1], mu [mu > 1]),
170 : maechler 1747 identical3(cu <= 1, tu <= 1, as(mu <= 1, "lMatrix")),# all lgeMatrix
171 : maechler 1738 identical3(cu[cu <= 1], tu[tu <= 1], mu[mu <= 1]),
172 :     identical3(cu , triu(cu ), t(t(cu))),
173 :     identical3(tu , triu(tu ), t(t(tu))),
174 :     identical3(tu., triu(tu.), t(t(tu.))),
175 :     identical(tu2, triu(tu2)),
176 :     identical(tcu , tril(tcu)),
177 :     identical(ttu , tril(ttu)),
178 :     identical(t(tu), tril(t(tu)))
179 :     )
180 : maechler 2113 assert.EQ.mat(triu(cu), as.matrix(triu(as.matrix(cu))))
181 :     for(k in -1:1)
182 :     assert.EQ.mat(tril(cu,k), as.matrix(tril(as.matrix(cu),k)))
183 : maechler 948
184 : maechler 2128 (dtr <- Matrix(local({m <- diag(2); m[1,2] <- 3;m})))
185 :     identical(dtr, triu(dtr))
186 :     assert.EQ.mat(tril(dtr), diag(2))
187 : maechler 2113
188 :    
189 : maechler 1767 (t4 <- new("dgTMatrix", i = 3:0, j = 0:3, x = rep(1,4), Dim = as.integer(c(4,4))))
190 :     c4 <- as(t4, "CsparseMatrix")
191 :     ## the same but "dsT" (symmetric)
192 :     M <- Matrix(c(0, rep(c(0,0:1),4)), 4,4)#ok warning
193 :     tt <- as(M, "TsparseMatrix")
194 :     stopifnot(all.equal(triu(t4) + tril(t4), c4),
195 :     all.equal(triu(tt) + tril(tt), c4))
196 : maechler 1738
197 : maechler 1767
198 : maechler 948 ###-- Numeric Dense: Crossprod & Solve
199 :    
200 : maechler 509 set.seed(123)
201 : maechler 1207 mm. <- mm <- Matrix(rnorm(500 * 150), nc = 150)
202 : maechler 509 stopifnot(validObject(mm))
203 : maechler 1207 xpx <- crossprod(mm)
204 :     stopifnot(identical(mm, mm.),# once upon a time, mm was altered by crossprod()
205 : maechler 873 validObject(xpx))
206 :     str(mm) # 'dge*"
207 :     str(xpx)# 'dpo*"
208 : bates 346 xpy <- crossprod(mm, rnorm(500))
209 :     res <- solve(xpx, xpy)
210 : maechler 873 str(xpx)# now with Cholesky factor
211 :     stopifnot(validObject(xpx),
212 :     validObject(xpy),
213 :     validObject(res))
214 :     stopifnot(all.equal(xpx %*% res, xpy, tol= 1e-12))
215 : maechler 1705 lp <- xpx >= 1
216 :     slp <- as(lp, "sparseMatrix")
217 : maechler 1749
218 : maechler 1921 ltlp <- lp[ lower.tri(lp) ]
219 :     sltlp <- slp[ lower.tri(slp) ]
220 : maechler 2113 dim(ij <- which(lower.tri(lp), arr.ind = TRUE))
221 :     ss <- slp[ij] # now fast (!)
222 :     stopifnot(identical4(lp[ij], ltlp, sltlp, as(lp, "matrix")[ij]),
223 :     identical(ss, sltlp),
224 : maechler 2157 isValid(lp, "lsyMatrix"), lp@uplo == "U")
225 : maechler 509
226 : maechler 948 ###-- more solve() methods {was ./solve.R }
227 :    
228 :     ## first for "dgeMatrix" and all kinds of RHS :
229 :     (m6 <- 1 + as(diag(0:5), "dgeMatrix"))
230 :     rcond(m6)
231 :     I6 <- as(diag(6), "dgeMatrix")
232 :     stopifnot(all.equal(I6, m6 %*% solve(m6)),
233 :     all.equal(I6, solve(m6) %*% m6) )
234 :    
235 :     (i6 <- solve(m6, Matrix(1:6)))
236 :     stopifnot(identical(i6, as(cbind(c(-4, rep(1,5))), "dgeMatrix")),
237 : maechler 949 identical(i6, solve(m6, 1:6)),
238 : maechler 948 identical(i6, solve(m6, matrix(1:6))),
239 : maechler 949 identical(i6, solve(m6, matrix(c(1,2,3,4,5,6))))
240 : maechler 948 )
241 :    
242 : maechler 1655 ## solve(<sparse>)
243 :     (m <- t1+ t(t1) + Diagonal(4))
244 :     i.m <- solve(as.mat(m))
245 : maechler 1656 I1 <- m %*% i.m
246 :     o4 <- diag(I1)
247 : maechler 1655 im <- solve(m)
248 :     (I2 <- m %*% im)
249 : maechler 2115 (ms <- as(m, "symmetricMatrix"))
250 : maechler 1665 ## solve(<sparse>, <sparse>):
251 :     s.mm <- solve(m,m)
252 :     s.mms <- solve(m, ms)
253 :     ## these now work "fully-sparse"
254 :     s.ms2 <- solve(ms, ms)
255 :     s.msm <- solve(ms, m)
256 : maechler 2115 I4c <- as(Matrix(diag(4),sparse=TRUE), "generalMatrix")
257 : maechler 2157 stopifnot(isValid(im, "Matrix"), isValid(I2, "Matrix"), class(I4c) == "dgCMatrix",
258 : maechler 1655 all.equal(I1, I2, tol = 1e-14),
259 :     all.equal(diag(4), as.mat(I2), tol = 1e-12),
260 : maechler 1665 all.equal(s.mm, I2, tol = 1e-14),
261 :     all.equal(s.mms, I2, tol = 1e-14),
262 :     all.equal(s.ms2, s.msm, tol = 4e-15),
263 :     all.equal(s.ms2, I4c , tol = 4e-15),
264 : maechler 1655 abs(o4 - 1) < 1e-14)
265 :    
266 : maechler 948 ###-- row- and column operations {was ./rowcolOps.R }
267 :    
268 :     set.seed(321)
269 : maechler 1467 (m1 <- round(Matrix(rnorm(25), 5), 2))
270 : maechler 956 m1k <- Matrix(round(rnorm(1000), 2), 50, 20)
271 :     m.m <- as(m1k, "matrix")
272 :     stopifnot(all.equal(colMeans(m1k), colMeans(m.m)),
273 :     all.equal(colSums (m1k), colSums (m.m)),
274 :     all.equal(rowMeans(m1k), rowMeans(m.m)),
275 :     all.equal(rowSums (m1k), rowSums (m.m)))
276 : maechler 948
277 : maechler 956 ###-- kronecker for nonsparse uses Matrix(.):
278 : maechler 2157 stopifnot(isValid(kr <- kronecker(m1, m6), "Matrix"))
279 : maechler 956 assert.EQ.mat(kr,
280 :     kronecker(as(m1, "matrix"),
281 : maechler 1921 as(m6, "matrix")), tol = 0)
282 :    
283 : maechler 956 ## sparse:
284 :     (kt1 <- kronecker(t1, tu))
285 :     kt2 <- kronecker(t1c, cu)
286 : maechler 1883 stopifnot(identical(Matrix:::uniq(kt1), Matrix:::uniq(kt2)))
287 :     ## but kt1 and kt2, both "dgT" are different since entries are not ordered!
288 : maechler 956 ktf <- kronecker(as.matrix(t1), as.matrix(tu))
289 : maechler 1593 if(FALSE) # FIXME? our kronecker treats "0 * NA" as "0" for structural-0
290 : maechler 956 assert.EQ.mat(kt2, ktf, tol= 0)
291 : maechler 1887 (cs1 <- colSums(kt1))
292 :     NA.or.True <- function(x) is.na(x) | x
293 :     eq <- (cs1 == colSums(as(kt1, "matrix")))
294 :     stopifnot(NA.or.True(eq), identical(is.na(eq), is.na(cs1)))
295 :     nt1 <- as(kt1, "nMatrix") # no NA's anymore
296 : maechler 1911 (ng1 <- as(as(nt1, "generalMatrix"),"CsparseMatrix")) # ngC
297 :     dg1 <- as(ng1, "dMatrix")# dgC
298 : maechler 1921 lt1 <- kt1 > 5
299 :     nt1 <- as(lt1, "nMatrix")
300 :     (colSums(nt1, sparseResult = TRUE))
301 :     (colSums(kt1, sparseResult = TRUE)) # dsparse, with NA
302 :     (colSums(lt1, sparseResult = TRUE)) # isparse, with NA
303 :     (colSums(lt1, sparseResult = TRUE, na.rm = TRUE))
304 :     (colSums(nt1, sparseResult = TRUE)) # isparse, no NA
305 : maechler 1887 ## check correct sparseness of both:
306 : maechler 1921 for(M in list(kt1, nt1, ng1, dg1, lt1, nt1)) {
307 :     m <- as(M, "matrix")
308 :     for(na.rm in c(FALSE,TRUE)) {
309 :     cs <- colSums(M, na.rm = na.rm)
310 :     cs. <- colSums(M, na.rm = na.rm, sparseResult = TRUE)
311 :     rs <- rowSums(M, na.rm = na.rm)
312 :     rs. <- rowSums(M, na.rm = na.rm, sparseResult = TRUE)
313 : maechler 2157 stopifnot(isValid(cs., "sparseVector"), identical(cs, as(cs., "vector")),
314 :     isValid(rs., "sparseVector"), identical(rs, as(rs., "vector")),
315 : maechler 1921 {eq <- cs == colSums(m, na.rm = na.rm) ; ineq <- is.na(eq)
316 :     all(ineq | eq) && identical(ineq, is.na(cs)) },
317 :     {eq <- rs == rowSums(m, na.rm = na.rm) ; ineq <- is.na(eq)
318 : maechler 2157 all(ineq | eq) && identical(ineq, is.na(rs)) } )
319 : maechler 1921 }
320 :     }
321 : maechler 956
322 : mmaechler 2175 (N <- as(crossprod(kronecker(diag(2), Matrix(c(2:0,1),2))) > 0,
323 :     "nMatrix"))
324 :     (L. <- as(N,"lMatrix"))
325 :     stopifnot(identical(N, as(L.,"nMatrix")))
326 :    
327 : maechler 1331 ## coercion from "dpo" or "dsy"
328 :     xx <- as(xpx, "dsyMatrix")
329 :     stopifnot(isSymmetric(xxS <- as(xx, "sparseMatrix")),
330 :     isSymmetric(xpxS <- as(xpx, "sparseMatrix")))
331 :    
332 : maechler 1238 tm <- matrix(0, 8,8)
333 :     tm[cbind(c(1,1,2,7,8),
334 :     c(3,6,4,8,8))] <- c(2,-30,15,20,80)
335 :     (tM <- Matrix(tm)) ## dtC
336 :     (mM <- Matrix(m <- (tm + t(tm)))) ## dsC
337 :     mT <- as(mM, "dsTMatrix")
338 : maechler 1319 gC <- as(as(mT, "dgTMatrix"), "dgCMatrix")
339 :     ## Check that 'mT' and gC print properly :
340 :     pr.mT <- capture.output(mT)
341 :     nn <- unlist(strsplit(gsub(" +\\.", "", sub("^....", "", pr.mT[-(1:2)])), " "))
342 :     stopifnot(as.numeric(nn[nn != ""]) == m[m != 0],
343 :     capture.output(gC)[-1] == pr.mT[-1])
344 : maechler 1238 assert.EQ.mat(tM, tm, tol=0)
345 : maechler 1319 assert.EQ.mat(gC, m, tol=0)
346 : maechler 1238 assert.EQ.mat(mT, m, tol=0)
347 : maechler 2157 stopifnot(isValid(mM, "dsCMatrix"), isValid(tM, "dtCMatrix")
348 : maechler 1825 , identical(mT, as(mM, "TsparseMatrix"))
349 :     , identical(gC, as(mM, "generalMatrix"))
350 :     ## coercions general <-> symmetric
351 : maechler 2115 , identical(as(as(mM, "generalMatrix"), "symmetricMatrix"), mM)
352 :     , identical(as(as(mM, "dgTMatrix"), "symmetricMatrix"), mT)
353 :     , identical(as(as(tM, "generalMatrix"),"triangularMatrix"), tM)
354 : maechler 2113 , identical(tM + Diagonal(8), tMD <- Diagonal(8) + tM)
355 : maechler 2157 , isValid(tMD, "dtCMatrix")
356 : maechler 1825 )
357 : maechler 1238 eM <- eigen(mM) # works thanks to base::as.matrix hack in ../R/zzz.R
358 :     stopifnot(all.equal(eM$values,
359 :     { v <- c(162.462112512353, 30.0665927567458)
360 :     c(v, 15, 0, 0, 160-v[1], -15, -v[2])}, tol=1e-14))
361 :    
362 : maechler 974 ##--- symmetric -> pos.def. needs valid test:
363 : maechler 1238 m5 <- Matrix(diag(5) - 1)
364 : mmaechler 2183 assertError(as(m5, "dpoMatrix"))# not pos.definite!
365 :     pm5 <- as(m5, "dspMatrix") # packed
366 :     assertError(as(pm5, "dppMatrix"))# not pos.definite!
367 : maechler 956
368 : maechler 1738 ###-- dense nonzero pattern:
369 :     class(m <- Matrix(TRUE,2,2)) # lsy
370 :     (n <- as(m, "nMatrix")) # nsy
371 :     validObject(n)
372 :    
373 :     ## 1)
374 :     as(n,"CsparseMatrix") # used to give CHOLMOD error: invalid xtype...
375 :     ls2 <- as(m, "CsparseMatrix") # works fine
376 :     ## and really 'm' and 'n' are interally slot identical (!!!)
377 :    
378 :     as(n,"sparseMatrix")
379 :     as(m,"sparseMatrix")
380 :    
381 :     ### -- now when starting with nsparse :
382 :     nT <- new("ngTMatrix",
383 :     i = as.integer(c(0, 1, 0)),
384 :     j = as.integer(c(0, 0, 1)), Dim = as.integer(c(2,2)),
385 :     Dimnames = list(NULL, NULL))
386 :     (nC <- as(nT, "ngCMatrix"))
387 :     str(nC)# of course, no 'x' slot
388 :    
389 : maechler 1749 stopifnot(identical(tt <- as(nT,"denseMatrix"), # lge
390 :     as(as(nT, "lMatrix"),"denseMatrix")))
391 :     tt
392 : maechler 1738 as(nC,"denseMatrix")
393 :    
394 :    
395 : maechler 1548 ###-- sparse nonzero pattern : ----------
396 : maechler 956
397 : maechler 1825 (nkt <- as(as(as(kt1, "generalMatrix"), "CsparseMatrix"), "ngCMatrix"))# ok
398 : maechler 1738 dkt <- as(nkt, "denseMatrix")
399 : maechler 1548 (clt <- crossprod(nkt))
400 : maechler 2157 stopifnot(isValid(nkt, "ngCMatrix"),
401 :     isValid(clt, "nsCMatrix"))
402 : maechler 1593 crossprod(clt) ## a warning: crossprod() of symmetric
403 : maechler 956
404 : maechler 1893 ## a Csparse with *repeated* entry is not valid!
405 : maechler 1894 assertError(new("ngCMatrix", p = c(0L,2L), i = c(0L,0L), Dim = 2:1))
406 : maechler 961
407 : maechler 1893
408 : maechler 956 ### "d" <-> "l" for (symmetric) sparse :
409 : maechler 1228 data(KNex)
410 :     mm <- KNex$mm
411 : maechler 956 xpx <- crossprod(mm)
412 : maechler 1548 ## extract nonzero pattern
413 :     nxpx <- as(xpx, "nsCMatrix")
414 : maechler 2113 show(nxpx) ## now ok, since subsetting works
415 :     r <- nxpx[1:2,]
416 : maechler 1883 lmm <- as(mm, "lgCMatrix")
417 : maechler 1548 nmm <- as(lmm, "nMatrix")
418 : maechler 1883 xlx <- crossprod(lmm)
419 :     x.x <- crossprod(nmm)
420 : maechler 961 ## now A = lxpx and B = xlx should be close, but not quite the same
421 :     ## since <x,y> = 0 is well possible when x!=0 and y!=0 .
422 :     ## However, A[i,j] != 0 ==> B[i,j] != 0:
423 : maechler 1593 A <- as(as(nxpx, "lMatrix"), "TsparseMatrix")
424 :     B <- as(as(xlx, "lMatrix"), "TsparseMatrix")
425 : maechler 961 ij <- function(a) a@i + ncol(a) * a@j
426 :     stopifnot(all(ij(A) %in% ij(B)))
427 : maechler 956
428 : maechler 1618 l3 <- upper.tri(matrix(,3,3))
429 : maechler 2157 stopifnot(isValid(c3 <- as(l3, "CsparseMatrix"), "CsparseMatrix"),# lgC
430 :     is(c3, "lMatrix"))
431 :     (M <- Matrix(l3))
432 :     stopifnot(isValid(M, "ltCMatrix"),
433 :     isValid(M2 <- M %x% M, "triangularMatrix"), # is "dtT" (why not "dtC" ?)
434 :     dim(M2) == c(9,9), identical(M2, kronecker(M,M)))
435 : maechler 1904 M3 <- M %x% M2 #ok
436 :     (cM3 <- colSums(M3, sparse=TRUE))
437 :     identical(as.vector(cM3),
438 :     as(rev(rowSums(M3, sparse=TRUE)), "vector"))
439 :     M. <- M2 %x% M # gave infinite recursion
440 : maechler 974
441 : maechler 1654 ## diagonal, sparse & interactions
442 : maechler 2157 stopifnot(isValid(as(Diagonal(3), "TsparseMatrix"), "TsparseMatrix"),
443 :     isValid(X <- Diagonal(7) + 1.5 * tM[1:7,1:7], "sparseMatrix"),
444 :     isValid(X, "triangularMatrix"),
445 :     isValid(XX <- X - chol(crossprod(X)), "triangularMatrix"))
446 : maechler 1654 X
447 : maechler 2115 XX
448 : maechler 2113 XX <- as(drop0(XX), "dsCMatrix")
449 : maechler 1654 stopifnot(identical(XX, Matrix(0, nrow(X), ncol(X))))
450 :    
451 : maechler 1725 M <- Matrix(m., sparse = FALSE)
452 :     (sM <- Matrix(m.))
453 :     class(dlM <- M >= 1)
454 :     stopifnot(identical(dlM, !(M < 1)),
455 : maechler 2157 isValid(sM, "sparseMatrix"),
456 :     isValid(dlM, "denseMatrix"))
457 : maechler 1725 (lM <- as(dlM, "sparseMatrix"))
458 :     lM2 <- as(dlM, "CsparseMatrix") #-> now ok
459 :     lM0 <- Matrix:::as_Csparse(dlM)
460 :     stopifnot(identical3(lM, lM2, lM0))
461 : maechler 1654
462 : maechler 1725 selectMethod("coerce", c("lgeMatrix", "CsparseMatrix"),
463 :     useInherited = c(from = TRUE, to = FALSE))
464 :    
465 : maechler 1747 ms0 <- Matrix(c(0,1,1,0), 2,2)
466 : maechler 2157 ms <- as(ms0, "TsparseMatrix")
467 : maechler 1767 cs <- as(ms, "CsparseMatrix")
468 : maechler 1747 ll <- as(ms, "lMatrix")
469 :     lt <- as(ll, "lgTMatrix")
470 : maechler 1767 nn <- as(cs, "nsparseMatrix")
471 :     l2 <- as(cs, "lsparseMatrix")
472 :     nt <- triu(nn)
473 :     n3 <- as(nt, "lsparseMatrix")
474 :     da <- nt + t(nt)
475 :     dm <- nt * t(nt) + da
476 : maechler 2157 stopifnot(isValid(ms, "dsTMatrix"),
477 :     as(ms0,"matrix") == as(ll, "matrix"), # coercing num |-> log
478 : maechler 1767 as(lt, "matrix") == as(ll, "matrix"),
479 : maechler 2157 identical(ms, as(ll, "dMatrix")),
480 : maechler 1767 identical4(as(ll, "CsparseMatrix"), as(cs, "lMatrix"),# lsC*
481 :     as(nn, "lsparseMatrix"), l2),
482 :     identical3(da, dm, as(cs, "generalMatrix")), # dgC*
483 :     identical(as(da, "lMatrix"), as(lt, "CsparseMatrix")) # lgC*
484 :     )
485 : mmaechler 2175 ## Dense *packed* ones:
486 :     s4 <- as(D4, "symmetricMatrix")
487 :     sp <- as(as(as(D4, "symmetricMatrix"),"denseMatrix"),"dspMatrix")
488 :     tp <- as(triu(sp),"dtpMatrix")
489 :     tpL <- as(tril(sp),"dtpMatrix")
490 :     (spL <- t(sp))
491 :     stopifnot(sp @uplo=="U", tp @uplo=="U",
492 :     spL@uplo=="L", tpL@uplo=="L")
493 : maechler 1747
494 : mmaechler 2175 D. <- Diagonal(x= c(-2,3:4)); D.[lower.tri(D.)] <- 1:3 ; D.
495 :     D0 <- Diagonal(x= 0:3); D0[upper.tri(D0)] <- 1:6 ; D0
496 :     stopifnot(all.equal(list(modulus = structure(24, logarithm = FALSE), sign = -1L),
497 :     unclass(determinant(D.,FALSE)), tol=1e-15),
498 :     all.equal(list(modulus = structure(0, logarithm = FALSE), sign = 1L),
499 :     unclass(determinant(D0,FALSE)), tol=0)
500 :     )
501 : maechler 1767
502 : maechler 1747
503 : mmaechler 2175 cat('Time elapsed: ', (.pt <- proc.time()),'\n') # "stats"
504 :     ##
505 :     cat("checkMatrix() of all: \n---------\n")
506 :     Sys.setlocale("LC_COLLATE", "C")# to keep ls() reproducible
507 :     for(nm in ls()) if(is(.m <- get(nm), "Matrix")) {
508 :     cat("\n", rep("-",nchar(nm)),"\n",nm, ":\n", sep='')
509 :     checkMatrix(.m)
510 :     }
511 :     cat('Time elapsed: ', proc.time() - .pt,'\n') # "stats"
512 :    
513 : maechler 1747 if(!interactive()) warnings()

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