pkg/tests/testthat/test-optimParams.R

   1 context("OptimParams")
   2
   3 naive_f <- function(link, M1,M2,M3, p,β,b)
   4 {
   5   d <- length(M1)
   6   K <- length(p)
   7   λ <- sqrt(colSums(β^2))
   8
   9   # Compute β x2,3 (self) tensorial products
  10   β2 <- array(0, dim=c(d,d,K))
  11   β3 <- array(0, dim=c(d,d,d,K))
  12   for (k in 1:K)
  13   {
  14     for (i in 1:d)
  15     {
  16       for (j in 1:d)
  17       {
  18         β2[i,j,k] = β[i,k]*β[j,k]
  19         for (l in 1:d)
  20           β3[i,j,l,k] = β[i,k]*β[j,k]*β[l,k]
  21       }
  22     }
  23   }
  24
  25   res <- 0
  26   for (i in 1:d)
  27   {
  28     term <- 0
  29     for (k in 1:K)
  30       term <- term + p[k]*.G(link,1,λ[k],b[k])*β[i,k]
  31     res <- res + (term - M1[i])^2
  32     for (j in 1:d)
  33     {
  34       term <- 0
  35       for (k in 1:K)
  36         term <- term + p[k]*.G(link,2,λ[k],b[k])*β2[i,j,k]
  37       res <- res + (term - M2[i,j])^2
  38       for (l in 1:d)
  39       {
  40         term <- 0
  41         for (k in 1:K)
  42           term <- term + p[k]*.G(link,3,λ[k],b[k])*β3[i,j,l,k]
  43         res <- res + (term - M3[i,j,l])^2
  44       }
  45     }
  46   }
  47   res
  48 }
  49
  50 # TODO: understand why it fails and reactivate this test
  51 #test_that("naive computation provides the same result as vectorized computations",
  52 #{
  53 #  h <- 1e-7 #for finite-difference tests
  54 #  tol <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK)
  55 #  n <- 10
  56 #  for (dK in list( c(2,2), c(5,3)))
  57 #  {
  58 #    d <- dK[1]
  59 #    K <- dK[2]
  60 #
  61 #    M1 <- runif(d, -1, 1)
  62 #    M2 <- matrix(runif(d^2, -1, 1), ncol=d)
  63 #    M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d))
  64 #
  65 #    for (link in c("logit","probit"))
  66 #    {
  67 #      # X and Y are unused here (W not re-computed)
  68 #      op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5),
  69 #        K, link, M=list(M1,M2,M3))
  70 #      op$W <- diag(d + d^2 + d^3)
  71 #
  72 #      for (var in seq_len((2+d)*K-1))
  73 #      {
  74 #        p <- runif(K, 0, 1)
  75 #        p <- p / sum(p)
  76 #        β <- matrix(runif(d*K,-5,5),ncol=K)
  77 #        b <- runif(K, -5, 5)
  78 #        x <- c(p[1:(K-1)],as.double(β),b)
  79 #
  80 #        # Test functions values
  81 #        expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) )
  82 #
  83 #        # Test finite differences ~= gradient values
  84 #        dir_h <- rep(0, (2+d)*K-1)
  85 #        dir_h[var] = h
  86 #        expect_equal(op$grad_f(x)[var], (op$f(x+dir_h) - op$f(x)) / h, tol)
  87 #      }
  88 #    }
  89 #  }
  90 #})
  91
  92 test_that("W computed in C and in R are the same",
  93 {
  94   tol <- 1e-8
  95   n <- 500
  96   for (dK in list( c(2,2), c(5,3)))
  97   {
  98     d <- dK[1]
  99     K <- dK[2]
 100     link <- ifelse(d==2, "logit", "probit")
 101     θ <- list(
 102       p=rep(1/K,K),
 103       β=matrix(runif(d*K),ncol=K),
 104       b=rep(0,K))
 105     io <- generateSampleIO(n, θ$p, θ$β, θ$b, link)
 106     X <- io$X
 107     Y <- io$Y
 108     dd <- d + d^2 + d^3
 109     p <- θ$p
 110     β <- θ$β
 111     λ <- sqrt(colSums(β^2))
 112     b <- θ$b
 113     β2 <- apply(β, 2, function(col) col %o% col)
 114     β3 <- apply(β, 2, function(col) col %o% col %o% col)
 115     M <- c(
 116       β  %*% (p * .G(link,1,λ,b)),
 117       β2 %*% (p * .G(link,2,λ,b)),
 118       β3 %*% (p * .G(link,3,λ,b)))
 119     Id <- as.double(diag(d))
 120     E <- diag(d)
 121     v1 <- Y * X
 122     v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )
 123     v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
 124       - Reduce('+', lapply(1:d, function(j)
 125         as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
 126       - Reduce('+', lapply(1:d, function(j)
 127         as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
 128       - Reduce('+', lapply(1:d, function(j)
 129         as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) )
 130     Omega1 <- matrix(0, nrow=dd, ncol=dd)
 131     for (i in 1:n)
 132     {
 133       gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M))
 134       Omega1 <- Omega1 + t(gi) %*% gi / n
 135     }
 136     W <- matrix(0, nrow=dd, ncol=dd)
 137     Omega2 <- matrix( .C("Compute_Omega",
 138       X=as.double(X), Y=as.integer(Y), M=as.double(M),
 139       pn=as.integer(n), pd=as.integer(d),
 140       W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
 141     rg <- range(Omega1 - Omega2)
 142     expect_equal(rg[1], rg[2], tol)
 143   }
 144 })