Commit | Line | Data |
---|---|---|
cbd88fe5 BA |
1 | context("OptimParams") |
2 | ||
3 | naive_f = function(link, M1,M2,M3, p,β,b) | |
4 | { | |
6dd5c2ac BA |
5 | d = length(M1) |
6 | K = length(p) | |
7 | λ <- sqrt(colSums(β^2)) | |
cbd88fe5 | 8 | |
6dd5c2ac BA |
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 | } | |
cbd88fe5 | 24 | |
6dd5c2ac BA |
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 | |
cbd88fe5 BA |
48 | } |
49 | ||
50 | test_that("naive computation provides the same result as vectorized computations", | |
51 | { | |
6dd5c2ac BA |
52 | h <- 1e-7 #for finite-difference tests |
53 | tol <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK) | |
54 | for (dK in list( c(2,2), c(5,3))) | |
55 | { | |
56 | d = dK[1] | |
57 | K = dK[2] | |
cbd88fe5 | 58 | |
6dd5c2ac BA |
59 | M1 = runif(d, -1, 1) |
60 | M2 = matrix(runif(d*d,-1,1), ncol=d) | |
61 | M3 = array(runif(d*d*d,-1,1), dim=c(d,d,d)) | |
cbd88fe5 | 62 | |
6dd5c2ac BA |
63 | for (link in c("logit","probit")) |
64 | { | |
65 | op = new("OptimParams", "li"=link, "M1"=as.double(M1), | |
66 | "M2"=as.double(M2), "M3"=as.double(M3), "K"=as.integer(K)) | |
cbd88fe5 | 67 | |
6dd5c2ac BA |
68 | for (var in seq_len((2+d)*K-1)) |
69 | { | |
70 | p = runif(K, 0, 1) | |
71 | p = p / sum(p) | |
72 | β <- matrix(runif(d*K,-5,5),ncol=K) | |
73 | b = runif(K, -5, 5) | |
74 | x <- c(p[1:(K-1)],as.double(β),b) | |
cbd88fe5 | 75 | |
6dd5c2ac BA |
76 | # Test functions values |
77 | expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) ) | |
cbd88fe5 | 78 | |
6dd5c2ac BA |
79 | # Test finite differences ~= gradient values |
80 | dir_h <- rep(0, (2+d)*K-1) | |
81 | dir_h[var] = h | |
cbd88fe5 | 82 | |
6dd5c2ac BA |
83 | expect_equal( op$grad_f(x)[var], ( op$f(x+dir_h) - op$f(x) ) / h, tol ) |
84 | } | |
85 | } | |
86 | } | |
cbd88fe5 | 87 | }) |
0f5fbd13 | 88 | |
44559add BA |
89 | test_that("W computed in C and in R are the same", |
90 | { | |
91 | # TODO: provide data X,Y + parameters theta | |
92 | dd <- d + d^2 + d^3 | |
93 | p <- θ$p | |
94 | β <- θ$β | |
95 | λ <- sqrt(colSums(β^2)) | |
96 | b <- θ$b | |
97 | β2 <- apply(β, 2, function(col) col %o% col) | |
98 | β3 <- apply(β, 2, function(col) col %o% col %o% col) | |
99 | M <- c( | |
100 | β %*% (p * .G(li,1,λ,b)), | |
101 | β2 %*% (p * .G(li,2,λ,b)), | |
102 | β3 %*% (p * .G(li,3,λ,b))) | |
103 | Id <- as.double(diag(d)) | |
104 | E <- diag(d) | |
105 | v1 <- Y * X | |
106 | v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) ) | |
107 | v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi | |
108 | - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d)) | |
109 | - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d)) | |
110 | - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) ) | |
111 | Omega1 <- matrix(0, nrow=dd, ncol=dd) | |
112 | for (i in 1:n) | |
113 | { | |
114 | gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M)) | |
115 | Omega1 <- Omega1 + t(gi) %*% gi / n | |
116 | } | |
117 | Omega2 <- matrix( .C("Compute_Omega", | |
118 | X=as.double(X), Y=as.double(Y), M=as.double(M), | |
119 | pn=as.integer(n), pd=as.integer(d), | |
120 | W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) | |
121 | rg <- range(Omega1 - Omega2) | |
122 | expect_that(rg[2] - rg[1] <= 1e8) | |
123 | }) |