REV.R 8.05 KB
Newer Older
Gilles Kratzer's avatar
Gilles Kratzer committed
1
REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
2

Gilles Kratzer's avatar
Gilles Kratzer committed
3
4
    rejection <- 1
    A <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
5

Gilles Kratzer's avatar
Gilles Kratzer committed
6
    # stor number of edges
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
7

Gilles Kratzer's avatar
Gilles Kratzer committed
8
    n.edges <- sum(dag.tmp)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
9

Gilles Kratzer's avatar
Gilles Kratzer committed
10
    # randomly select one
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
11
12

    if (sum(dag.tmp) != 0)
Gilles Kratzer's avatar
Gilles Kratzer committed
13
        {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
14

Gilles Kratzer's avatar
Gilles Kratzer committed
15
            selected.edge <- which(x = dag.tmp == 1, arr.ind = TRUE)[sample(x = 1:n.edges, size = 1), ]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
16

Gilles Kratzer's avatar
Gilles Kratzer committed
17
            # store current parent set
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
18

Gilles Kratzer's avatar
Gilles Kratzer committed
19
            parent.set <- dag.tmp[selected.edge[1], ]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
20

Gilles Kratzer's avatar
Gilles Kratzer committed
21
22
23
            # remove parent set
            dag.M.dot <- dag.tmp
            dag.M.dot[selected.edge[1], ] <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
24

Gilles Kratzer's avatar
Gilles Kratzer committed
25
            # store descendents j row i col
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
26

Gilles Kratzer's avatar
Gilles Kratzer committed
27
28
            descendents.M.dot.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.dot)
            descendents.M.dot.j <- descendents(nodes = unname(selected.edge[1]), dag = dag.M.dot)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
29

Gilles Kratzer's avatar
Gilles Kratzer committed
30
            # mark all parent sets of i score that do not include j
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
31

Gilles Kratzer's avatar
Gilles Kratzer committed
32
33
            sc.tmp <- sc[score.cache$children == selected.edge[2], ]
            sc.tmp <- sc.tmp[sc.tmp[, selected.edge[1]] == 1, ]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
34

Gilles Kratzer's avatar
Gilles Kratzer committed
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
            # remove descendents of i
            if (is.null(ncol(sc.tmp))) {
                if (!is.null(descendents.M.dot.i)) {
                  if (length(descendents.M.dot.i) == 1) {
                    sc.tmp <- sc.tmp[sc.tmp[descendents.M.dot.i] == 0]
                  } else {
                    sc.tmp <- sc.tmp[sum(sc.tmp[descendents.M.dot.i] == 0) == length(descendents.M.dot.i)]
                  }
                }
            } else {
                if (!is.null(descendents.M.dot.i)) {
                  if (length(descendents.M.dot.i) == 1) {
                    sc.tmp <- sc.tmp[sc.tmp[, descendents.M.dot.i] == 0, ]
                  } else {
                    sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.i] == 0) == length(descendents.M.dot.i), ]
                  }
                }
            }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
53
54
55



Gilles Kratzer's avatar
Gilles Kratzer committed
56
            # sample new parent set
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
57

Gilles Kratzer's avatar
Gilles Kratzer committed
58
            if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
59

Gilles Kratzer's avatar
Gilles Kratzer committed
60
61
62
63
64
65
66
                if (is.vector(sc.tmp)) {
                  new.parent.i <- sc.tmp
                  new.parent.i <- new.parent.i[-length(new.parent.i)]
                  # store partition function
                  z.star.x.i.M.dot <- sc.tmp[length(sc.tmp)]
                } else {
                  # !!!exp missing
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
67
                  new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
Gilles Kratzer's avatar
Gilles Kratzer committed
68
69
70
71
72
                    sum(sc.tmp[, ncol(sc.tmp)]))), ]
                  new.parent.i <- new.parent.i[-length(new.parent.i)]
                  # store partition function
                  z.star.x.i.M.dot <- (sum(sc.tmp[, ncol(sc.tmp)]))
                }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
73

Gilles Kratzer's avatar
Gilles Kratzer committed
74
75
76
                # M direct sum
                dag.M.cross <- dag.M.dot
                dag.M.cross[selected.edge[2], ] <- new.parent.i
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
77

Gilles Kratzer's avatar
Gilles Kratzer committed
78
79
80
                # descendant of i
                descendents.M.cross.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.cross)
                descendents.M.cross.j <- descendents(nodes = unname(selected.edge[1]), dag = dag.M.cross)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
81

Gilles Kratzer's avatar
Gilles Kratzer committed
82
83
                # score node j not descendant of M
                sc.tmp <- sc[score.cache$children == selected.edge[1], ]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
84

Gilles Kratzer's avatar
Gilles Kratzer committed
85
86
87
88
89
                # remove descendents of i
                if (!is.null(descendents.M.cross.j)) {
                  if (length(descendents.M.cross.j) == 1) {
                    sc.tmp <- sc.tmp[sc.tmp[, descendents.M.cross.j] == 0, ]
                  } else {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
90
                    sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.cross.j] == 0) == length(descendents.M.cross.j),
Gilles Kratzer's avatar
Gilles Kratzer committed
91
92
93
                      ]
                  }
                }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
94

Gilles Kratzer's avatar
Gilles Kratzer committed
95
96
                # sample new parent set
                if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
97

Gilles Kratzer's avatar
Gilles Kratzer committed
98
99
100
101
102
103
104
105
                  if (is.vector(sc.tmp)) {
                    new.parent.j <- sc.tmp
                    # print(new.parent.i)
                    new.parent.j <- new.parent.j[-length(new.parent.j)]
                    # print(new.parent.i) store partition function
                    z.x.j.M.cross <- sc.tmp[length(sc.tmp)]
                  } else {
                    ## !!! exp missing
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
106
                    new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
Gilles Kratzer's avatar
Gilles Kratzer committed
107
108
109
110
111
                      sum(sc.tmp[, ncol(sc.tmp)]))), ]
                    new.parent.j <- new.parent.j[-length(new.parent.j)]
                    # store partition function
                    z.x.j.M.cross <- (sum(sc.tmp[, ncol(sc.tmp)]))
                  }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
112

Gilles Kratzer's avatar
Gilles Kratzer committed
113
                  # M tilde
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
114

Gilles Kratzer's avatar
Gilles Kratzer committed
115
116
                  dag.M.tilde <- dag.M.cross
                  dag.M.tilde[selected.edge[1], ] <- new.parent.j
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
117

Gilles Kratzer's avatar
Gilles Kratzer committed
118
                  n.edges.tilde <- sum(dag.M.tilde)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
119

Gilles Kratzer's avatar
Gilles Kratzer committed
120
                  ############################## computing acceptance proba
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
121

Gilles Kratzer's avatar
Gilles Kratzer committed
122
123
124
125
126
127
128
                  # score node j that do not include i
                  sc.tmp <- sc[score.cache$children == selected.edge[1], ]
                  sc.tmp <- sc.tmp[sc.tmp[, selected.edge[2]] == 1, ]
                  # remove descendents of i
                  if (is.vector(sc.tmp)) {
                    if (!is.null(descendents.M.dot.j)) {
                      if (length(descendents.M.dot.j) == 1) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
129

Gilles Kratzer's avatar
Gilles Kratzer committed
130
131
132
133
134
135
136
137
138
                        sc.tmp <- sc.tmp[sc.tmp[descendents.M.dot.j] == 0]  #
                      } else {
                        # print(class(sc.tmp)) print(descendents.M.dot.j)
                        sc.tmp <- sc.tmp[sum(sc.tmp[descendents.M.dot.j] == 0) == length(descendents.M.dot.j)]
                      }
                    }
                  } else {
                    if (!is.null(descendents.M.dot.j)) {
                      if (length(descendents.M.dot.j) == 1) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
139

Gilles Kratzer's avatar
Gilles Kratzer committed
140
141
142
                        sc.tmp <- sc.tmp[sc.tmp[, descendents.M.dot.j] == 0, ]
                      } else {
                        # print(class(sc.tmp)) print(descendents.M.dot.j)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
143
                        sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.j] == 0) == length(descendents.M.dot.j),
Gilles Kratzer's avatar
Gilles Kratzer committed
144
145
146
147
                          ]
                      }
                    }
                  }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
148

Gilles Kratzer's avatar
Gilles Kratzer committed
149
150
151
152
153
                  if (is.vector(sc.tmp)) {
                    z.star.x.j.M.dot <- sc.tmp[length(sc.tmp)]
                  } else {
                    z.star.x.j.M.dot <- (sum(sc.tmp[, ncol(sc.tmp)]))
                  }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
154
155


Gilles Kratzer's avatar
Gilles Kratzer committed
156
157
                  dag.M.tilde.cross <- dag.M.dot
                  dag.M.tilde.cross[selected.edge[1], ] <- parent.set
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
158

Gilles Kratzer's avatar
Gilles Kratzer committed
159
160
                  # descendant
                  descendents.M.tilde.cross.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.tilde.cross)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
161
162


Gilles Kratzer's avatar
Gilles Kratzer committed
163
164
165
166
167
                  sc.tmp <- sc[score.cache$children == selected.edge[2], ]
                  if (!is.null(descendents.M.tilde.cross.i)) {
                    if (length(descendents.M.tilde.cross.i) == 1) {
                      sc.tmp <- sc.tmp[sc.tmp[, descendents.M.tilde.cross.i] == 0, ]
                    } else {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
168
                      sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.tilde.cross.i] == 0) == length(descendents.M.tilde.cross.i),
Gilles Kratzer's avatar
Gilles Kratzer committed
169
170
171
                        ]
                    }
                  }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
172

Gilles Kratzer's avatar
Gilles Kratzer committed
173
174
175
176
177
                  if (is.vector(sc.tmp)) {
                    z.x.i.M.tilde.cross <- sc.tmp[length(sc.tmp)]
                  } else {
                    z.x.i.M.tilde.cross <- (sum(sc.tmp[, ncol(sc.tmp)]))
                  }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
178

Gilles Kratzer's avatar
Gilles Kratzer committed
179
                  ############################## Acceptance probability
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
180

Gilles Kratzer's avatar
Gilles Kratzer committed
181
182
                  score.A <- n.edges/n.edges.tilde * exp(z.star.x.i.M.dot - z.star.x.j.M.dot + z.x.j.M.cross - z.x.j.M.cross)
                  A <- min(1, score.A)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
183

Gilles Kratzer's avatar
Gilles Kratzer committed
184
185
186
187
188
                  if (rbinom(n = 1, size = 1, prob = A) == 1) {
                    rejection <- 0
                    dag.tmp <- dag.M.tilde
                    score.M.tilde <- 0
                    for (a in 1:n.var) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
189

Gilles Kratzer's avatar
Gilles Kratzer committed
190
                      sc.tmp <- sc[score.cache$children == a, ]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
191
                      score.M.tilde <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
Gilles Kratzer's avatar
Gilles Kratzer committed
192
193
194
195
196
197
198
                        ])))), n.var + 1]), score.M.tilde)
                    }
                    score <- score.M.tilde
                  }
                }
            }
        }  #EOFif
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
199
200


Gilles Kratzer's avatar
Gilles Kratzer committed
201
    ############################## Return
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
202

Gilles Kratzer's avatar
Gilles Kratzer committed
203
204
    return(list(dag.tmp = dag.tmp, score = score, alpha = A, rejection = rejection))
}  #EOF