Fix OptimParams + clean test file

[morpheus.git] / pkg / tests / testthat / test-optimParams.R

diff --git a/pkg/tests/testthat/test-optimParams.R b/pkg/tests/testthat/test-optimParams.R
index 8f65e46..305c36f 100644 (file)
--- a/pkg/tests/testthat/test-optimParams.R
+++ b/pkg/tests/testthat/test-optimParams.R
@@ -85,3 +85,39 @@ test_that("naive computation provides the same result as vectorized computations
     }
   }
 })
+
+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)
+})