-
-test_that("W computed in C and in R are the same",
-{
- # TODO: provide data X,Y + parameters theta
- dd <- d + d^2 + d^3
- p <- θ$p
- β <- θ$β
- λ <- sqrt(colSums(β^2))
- b <- θ$b
- β2 <- apply(β, 2, function(col) col %o% col)
- β3 <- apply(β, 2, function(col) col %o% col %o% col)
- M <- c(
- β %*% (p * .G(li,1,λ,b)),
- β2 %*% (p * .G(li,2,λ,b)),
- β3 %*% (p * .G(li,3,λ,b)))
- Id <- as.double(diag(d))
- E <- diag(d)
- v1 <- Y * X
- v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )
- v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
- - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
- - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
- - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) )
- Omega1 <- matrix(0, nrow=dd, ncol=dd)
- for (i in 1:n)
- {
- gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M))
- Omega1 <- Omega1 + t(gi) %*% gi / n
- }
- Omega2 <- matrix( .C("Compute_Omega",
- X=as.double(X), Y=as.double(Y), M=as.double(M),
- pn=as.integer(n), pd=as.integer(d),
- W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
- rg <- range(Omega1 - Omega2)
- expect_that(rg[2] - rg[1] <= 1e8)
-})