-
Notifications
You must be signed in to change notification settings - Fork 1
/
gennetwork_test.R
59 lines (43 loc) · 1.07 KB
/
gennetwork_test.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# Copyright (C) 2011-2012 Gad Abraham and National ICT Australia (NICTA).
# All rights reserved.
#
set.seed(2879)
Y <- round(matrix(rnorm(20 * 5), 20, 5), 2)
cortype <- 1
corthresh <- 0
R <- cor(Y)
R[abs(R) < corthresh] <- 0
K <- ncol(Y)
nV <- K
UR <- R
UR[lower.tri(UR)] <- 0
diag(UR) <- 0
W <- if(cortype == 0) {
(abs(UR) > corthresh) + 0
} else if(cortype == 1) {
abs(UR)
} else if(cortype == 2) {
UR^2
} else {
stop("unknown cortype:", cortype)
}
nzUR <- which(UR != 0)
E <- which(UR != 0, arr.ind=TRUE)
Ecoef <- W[nzUR]
Esign <- sign(R[nzUR])
#nE <- nrow(E)
nE <- K * (K - 1) / 2
C_I <- c(1:nE, 1:nE)
C_J <- as.numeric(E)
C_S <- cbind(Ecoef, -Ecoef * Esign)
C <- matrix(0, nE, nV)
C[cbind(C_I, C_J)] <- C_S
system("./gennetwork_test")
C2 <- matrix(scan("C.txt", sep=","), nE, nV, byrow=TRUE)
mean((C2 - C)^2)