X-Git-Url: https://git.auder.net/?a=blobdiff_plain;ds=sidebyside;f=pkg%2Ftests%2Ftestthat%2Ftest-optimParams.R;h=305c36fd04bc1ab21cbfe6c82cf6c357ea9435f6;hb=44559add0e38058d9ce539c4b91246e4a088f67a;hp=8f65e46c375b2d6c7eee619510aee7af69efd77c;hpb=6dd5c2acccd10635449230faa824b7e8906911bf;p=morpheus.git

diff --git a/pkg/tests/testthat/test-optimParams.R b/pkg/tests/testthat/test-optimParams.R
index 8f65e46..305c36f 100644
--- 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)
+})