+
+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)
+})