SCM

SCM Repository

[matrix] Annotation of /pkg/Matrix/tests/group-methods.R
ViewVC logotype

Annotation of /pkg/Matrix/tests/group-methods.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2813 - (view) (download)

1 : maechler 2011 ### Testing the group methods --- some also happens in ./Class+Meth.R
2 : maechler 573
3 :     library(Matrix)
4 : mmaechler 2490 source(system.file("test-tools.R", package = "Matrix"))# identical3() etc
5 :    
6 : maechler 573 set.seed(2001)
7 :    
8 :     mm <- Matrix(rnorm(50 * 7), nc = 7)
9 :     xpx <- crossprod(mm)# -> "factors" in mm !
10 :     round(xpx, 3) # works via "Math2"
11 :    
12 :     y <- rnorm(nrow(mm))
13 :     xpy <- crossprod(mm, y)
14 :     res <- solve(xpx, xpy)
15 :     signif(res, 4) # 7 x 1 Matrix
16 :    
17 : maechler 2011 stopifnot(all(signif(res) == signif(res, 6)),
18 :     all(round (xpx) == round (xpx, 0)))
19 :    
20 : maechler 573 ## exp(): component wise
21 :     signif(dd <- (expm(xpx) - exp(xpx)) / 1e34, 3)# 7 x 7
22 :    
23 :     stopifnot(validObject(xpx),
24 :     validObject(xpy),
25 :     validObject(dd))
26 :    
27 : maechler 925 ## "Math" also, for log() and [l]gamma() which need special treatment
28 :     stopifnot(identical(exp(res)@x, exp(res@x)),
29 :     identical(log(abs(res))@x, log(abs((res@x)))),
30 :     identical(lgamma(res)@x, lgamma(res@x)))
31 : maechler 573
32 :    
33 : maechler 925 ###--- sparse matrices ---------
34 :    
35 :     m <- Matrix(c(0,0,2:0), 3,5)
36 :     (mC <- as(m, "dgCMatrix"))
37 :     sm <- sin(mC)
38 : maechler 1226 stopifnot(class(sm) == class(mC), class(mC) == class(mC^2),
39 : maechler 925 dim(sm) == dim(mC),
40 :     class(0 + 100*mC) == class(mC),
41 :     all.equal(0.1 * ((0 + 100*mC)/10), mC),
42 : maechler 1226 all.equal(sqrt(mC ^ 2), mC),
43 :     all.equal(m^m, mC^mC),
44 :     identical(mC^2, mC * mC),
45 :     identical(mC*2, mC + mC)
46 :     )
47 : maechler 925
48 : maechler 956 x <- Matrix(rbind(0,cbind(0, 0:3,0,0,-1:2,0),0))
49 : maechler 1226 x # sparse
50 : maechler 1467 (x2 <- x + 10*t(x))
51 :     stopifnot(is(x2, "sparseMatrix"),
52 :     identical(x2, t(x*10 + t(x))),
53 :     identical(x, as((x + 10) - 10, class(x))))
54 :    
55 : maechler 956 (px <- Matrix(x^x - 1))#-> sparse again
56 :     stopifnot(px@i == c(3,4,1,4),
57 :     px@x == c(3,26,-2,3))
58 :    
59 : mmaechler 2767 ## From: "Florent D." .. Thu, 23 Feb 2012 -- bug report
60 :     ##---> MM: Make a regression test:
61 :     tst <- function(n, i = 1) {
62 :     stopifnot(i >= 1, n >= i)
63 :     D <- .sparseDiagonal(n)
64 :     ee <- numeric(n) ; ee[i] <- 1
65 :     stopifnot(all(D - ee == diag(n) - ee),
66 :     all(D * ee == diag(n) * ee),
67 :     all(ee - D == ee - diag(n)),
68 :     {C <- (ee / D == ee / diag(n)); all(is.na(C) | C)},
69 :     TRUE)
70 :     }
71 : mmaechler 2768 tmp <- sapply(1:16, tst)
72 : mmaechler 2767 i <- sapply(1:16, function(i) sample(i,1))
73 : mmaechler 2768 tmp <- mapply(tst, n= 1:16, i= i)
74 : mmaechler 2767
75 :    
76 : maechler 1571 ###----- Compare methods ---> logical Matrices ------------
77 :     l3 <- upper.tri(matrix(, 3, 3))
78 :     (ll3 <- Matrix(l3))
79 : mmaechler 2493 dt3 <- (99* Diagonal(3) + (10 * ll3 + Diagonal(3)))/10
80 : maechler 1571 (dsc <- crossprod(ll3))
81 :     stopifnot(validObject(ll3), validObject(dsc),
82 :     identical(ll3, t(t(ll3))),
83 : mmaechler 2493 identical(dsc, t(t(dsc))),
84 :     isValid(dsc + 3 * Diagonal(nrow(dsc)), "dsCMatrix"),
85 :     isValid(dt3, "triangularMatrix"), # remained triangular
86 :     isValid(dt3 > 0, "triangularMatrix")# ditto
87 : maechler 1571 )
88 : maechler 956
89 : maechler 1571 (lm1 <- dsc >= 1) # now ok
90 :     (lm2 <- dsc == 1) # now ok
91 :     nm1 <- as(lm1, "nMatrix")
92 :     (nm2 <- as(lm2, "nMatrix"))
93 :    
94 :     stopifnot(validObject(lm1), validObject(lm2),
95 :     validObject(nm1), validObject(nm2),
96 :     identical(dsc, as(dsc * as(lm1, "dMatrix"), "dsCMatrix")))
97 :    
98 : bates 1577 crossprod(lm1) # lm1: "lsC*"
99 : mmaechler 2490 cnm1 <- crossprod(nm1)
100 :     stopifnot(is(cnm1, "symmetricMatrix"), ## whereas the %*% is not:
101 :     Q.eq(cnm1, nm1 %*% nm1))
102 : mmaechler 2476 dn1 <- as(nm1, "denseMatrix")
103 :     stopifnot(all(dn1 == nm1))
104 : maechler 1571
105 : mmaechler 2503 dsc[2,3] <- NA ## now has an NA (and no longer is symmetric)
106 :     ## ----- and "everything" is different
107 : mmaechler 2556 ## also add "non-structural 0":
108 :     dsc@x[1] <- 0
109 : maechler 1571 dsc
110 :     dsc/ 5
111 :     dsc + dsc
112 :     dsc - dsc
113 :     dsc + 1 # -> no longer sparse
114 : mmaechler 2476 Tsc <- as(dsc, "TsparseMatrix")
115 : mmaechler 2556 dsc. <- drop0(dsc)
116 :     stopifnot(identical(dsc., Matrix((dsc + 1) -1)),
117 : mmaechler 2764 identical(as(-Tsc,"CsparseMatrix"), (-1) * Tsc),
118 :     identical(-dsc., (-1) * dsc.),
119 :     identical3(-Diagonal(3), Diagonal(3, -1), (-1) * Diagonal(3)),
120 :     identical(dsc., Matrix((Tsc + 1) -1)), # ok (exact arithmetic)
121 : mmaechler 2556 Q.eq(0 != dsc, dsc != Matrix(0, 3, 3)),
122 :     Q.eq(0 != dsc, dsc != c(0,0)) # with a warning ("not multiple ..")
123 : mmaechler 2764 )
124 : maechler 1571 str(lm1 <- dsc >= 1) # now ok (NA in proper place, however:
125 : maechler 1573 lm1 ## NA used to print as ' ' , now 'N'
126 : maechler 1582 (lm2 <- dsc == 1)# ditto
127 : maechler 1571
128 : mmaechler 2503 ddsc <- kronecker(Diagonal(7), dsc)
129 :     isValid(ddv <- rowSums(ddsc, sparse=TRUE), "sparseVector")
130 :     sv <- colSums(kC <- kronecker(mC,kronecker(mC,mC)), sparse=TRUE)
131 :     EQ <- ddv == rowSums(ddsc)
132 :     na.ddv <- is.na(ddv)
133 : mmaechler 2505 sM <- Matrix(pmax(0, round(rnorm(50*15, -1.5), 2)), 50,15)
134 : mmaechler 2503 stopifnot(sv == colSums(kC), is.na(as.vector(ddv)) == na.ddv,
135 : mmaechler 2505 isValid(sM/(-7:7), "CsparseMatrix"),
136 : mmaechler 2503 all(EQ | na.ddv))
137 :    
138 : maechler 1573 ## Just for print "show":
139 :     z <- round(rnorm(77), 2)
140 :     z[sample(77,10)] <- NA
141 :     (D <- Matrix(z, 7)) # dense
142 :     z[sample(77,15)] <- 0
143 :     (D <- Matrix(z, 7)) # sparse
144 :     abs(D) >= 0.5 # logical sparse
145 :    
146 : maechler 1582 stopifnot(identical(crossprod(lm1),# "lgC": here works!
147 : maechler 1571 crossprod(as(lm1, "dMatrix"))
148 :     ))
149 :    
150 : mmaechler 2813 D3 <- Diagonal(x=4:2) # for the checks below, also want a *diagonal*
151 :    
152 :    
153 : mmaechler 2476 ## Systematically look at all "Ops" group generics for "all" Matrix classes
154 :     ## -------------- Main issue: Detect infinite recursion problems
155 : mmaechler 2503 cl <- sapply(ls(), function(.) class(get(.)))
156 :     Mcl <- c(grep("Matrix$", cl, value=TRUE),
157 :     grep("sparseVector", cl, value=TRUE))
158 : mmaechler 2476 table(Mcl)
159 : mmaechler 2813 ## choose *one* of each class:
160 : mmaechler 2476 M.objs <- names(Mcl[!duplicated(Mcl)])
161 : mmaechler 2813 Mat.objs <- M.objs[vapply(M.objs, function(nm) is(get(nm), "Matrix"), NA)]
162 :     (MatDims <- t(vapply(Mat.objs, function(nm) dim(get(nm)), 0:1)))
163 :     mDims <- MatDims %*% (d.sig <- c(1, 1000)) # "dim-signature" to match against
164 : mmaechler 2476
165 : mmaechler 2813 m2num <- function(m) { if(is.integer(m)) storage.mode(m) <- "double" ; m }
166 : mmaechler 2476 cat("Checking all group generics for a set of arguments:\n",
167 :     "---------------------------------------------------\n", sep='')
168 :     for(gr in getGroupMembers("Ops")) {
169 :     cat(gr,"\n",paste(rep.int("=",nchar(gr)),collapse=""),"\n", sep='')
170 :     for(f in getGroupMembers(gr)) {
171 :     cat(sprintf("%9s : ", dQuote(f)))
172 :     for(nm in M.objs) {
173 :     M <- get(nm, inherits=FALSE)
174 : mmaechler 2503 nm <- NROW(M)
175 : mmaechler 2476 cat("o")
176 : mmaechler 2503 for(x in list(TRUE, -3.2, 0L, seq_len(nm))) {
177 : mmaechler 2476 cat(".")
178 :     validObject(r1 <- do.call(f, list(M,x)))
179 :     validObject(r2 <- do.call(f, list(x,M)))
180 :     stopifnot(dim(r1) == dim(M), dim(r2) == dim(M))
181 :     }
182 : mmaechler 2503 ## M o <sparseVector>
183 :     x <- numeric(nm)
184 :     x[c(1,length(x))] <- 1:2
185 :     sv <- as(x, "sparseVector")
186 :     cat("s.")
187 :     validObject(r3 <- do.call(f, list(M, sv)))
188 :     stopifnot(dim(r3) == dim(M))
189 : mmaechler 2813 if(is(M, "Matrix")) { ## M o <Matrix>
190 :     d <- dim(M)
191 :     ds <- sum(d * d.sig)
192 :     if(any(match. <- ds == mDims)) {
193 :     cat("\nM o M:")
194 :     for(oM in Mat.objs[match.]) {
195 :     M2 <- get(oM)
196 :     validObject(R4 <- do.call(f, list(M, M2)))
197 :     cat(".")
198 :     r4 <- m2num(do.call(f, list(as.mat(M), as.mat(M2))))
199 :     ## stopifnot(identical(r4, as.mat(R4)))
200 :     if(!identical(r4, as.mat(R4))) {
201 :     cat("\n ** not identical: r4 \\ R4 :\n")
202 :     print(r4); print(R4)
203 :     }
204 :     cat("i")
205 :     }
206 :     }
207 :     }
208 : mmaechler 2476 }
209 :     cat("\n")
210 :     }
211 :     }
212 :    
213 : mmaechler 2813 stopifnot(identical(lm2, lm1 & lm2),
214 :     identical(lm1, lm1 | lm2))
215 : maechler 1571
216 :    
217 : maechler 1216 cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''

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