mcmcabn-utiles.R 9.7 KB
Newer Older
Gilles Kratzer's avatar
Gilles Kratzer committed
1
2
3
4
5
6
############################################################################### mcmcabn-utiles.R --- Author : Gilles Kratzer Document created : 13/02/2018 Last modification :

##-------------------------------------------------------------------------
## Internal function to test input of mcmcabn()
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
7
.tests.mcmcabn <- function(score.cache, data.dists, max.parents, mcmc.scheme, seed, verbose, start.dag, prior.dag, prior.lambda,
8
    prob.rev, prob.mbr, prior.choice,scaling) {
Gilles Kratzer's avatar
Gilles Kratzer committed
9
    # start tests
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
10
    if (is.null(score.cache))
Gilles Kratzer's avatar
website    
Gilles Kratzer committed
11
        stop("A cache of score should be provided. You can produce it using the R package abn.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
12

Gilles Kratzer's avatar
Gilles Kratzer committed
13
    if (max(rowSums(score.cache$node.defn)) > (max.parents+1))
Gilles Kratzer's avatar
Gilles Kratzer committed
14
        stop("Check max.parents. It should be the same as the one used in abn::buildscorecache() R function")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
15
16

    if (length(mcmc.scheme) != 3)
Gilles Kratzer's avatar
Gilles Kratzer committed
17
        stop("An MCMC scheme have to be provided. It should be such that c(returned,thinned,burned) made of non negative integers.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
18
19

    if (!is.numeric(mcmc.scheme[1]) | !is.numeric(mcmc.scheme[2]) | !is.numeric(mcmc.scheme[3]) | mcmc.scheme[1] < 0 |
Gilles Kratzer's avatar
Gilles Kratzer committed
20
21
22
        mcmc.scheme[2] < 0 | mcmc.scheme[3] < 0) {
        stop("An MCMC scheme have to be provided. It should be such that c(returned,thinned,burned) made of non negative integers.")
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
23
24

    if (max.parents < 1 || max.parents > length(data.dists))
Gilles Kratzer's avatar
Gilles Kratzer committed
25
        stop("max.parents makes no sense.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
26
27

    if (is.numeric(prob.rev) && prob.rev > 1 && prob.rev < 0)
Gilles Kratzer's avatar
Gilles Kratzer committed
28
        stop("prob.rev should be a probability.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
29
30

    if (is.numeric(prob.mbr) && prob.mbr > 1 && prob.mbr < 0)
Gilles Kratzer's avatar
Gilles Kratzer committed
31
        stop("prob.mbr should be a probability.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
32

Gilles Kratzer's avatar
Gilles Kratzer committed
33
    if (is.matrix(start.dag) && (dim(start.dag)[1] != length(data.dists)))
Gilles Kratzer's avatar
Gilles Kratzer committed
34
        stop("start.dag should be a squared matrix with dimension equal to the number of variables.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
35

Gilles Kratzer's avatar
Gilles Kratzer committed
36
37
38
    if (is.matrix(start.dag) && is.null(topoSortMAT((start.dag)))) {
        stop("start.dag should be a named DAG.")
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
39

Gilles Kratzer's avatar
Gilles Kratzer committed
40
41
42
43
    if (is.matrix(start.dag) && max(rowSums(start.dag))>max.parents) {
      stop("start.dag should not have more parent than max.parent argument.")
    }

Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
44
    if (length(data.dists) != max(score.cache$children))
Gilles Kratzer's avatar
Gilles Kratzer committed
45
        stop("data.dists should be a named list of all variables used to build the cache of precomputed score.")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
46

Gilles Kratzer's avatar
Gilles Kratzer committed
47
48
49
    if (length(data.dists) == 0 || is.null(names(data.dists))) {
        stop("data.dists should be a named list of all variables used to build the cache of precomputed score.")
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
50

Gilles Kratzer's avatar
Gilles Kratzer committed
51
52
53
    if (prior.choice %!in% 1:3) {
        stop("prior.choice should be either 1,2 or 3.")
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
54

Gilles Kratzer's avatar
Gilles Kratzer committed
55
    if (is.matrix(prior.dag) && dim(prior.dag)[1] != length(data.dists))
Gilles Kratzer's avatar
Gilles Kratzer committed
56
        stop("prior.dag should be a squared matrix with dimension equal to the number of variables.")
57

58
  if(scaling<0)
59
60
    stop("heating parameter shoud be between zero and one. Zero corresponds to a lower probability of accepting a lower score proposal. One corresponds to a larger probability of accepting a lowering of the score (values higher than one are possible but very unlikely to produce meaningful (or helpful) results.).")

Gilles Kratzer's avatar
Gilles Kratzer committed
61
62
}

63

Gilles Kratzer's avatar
Gilles Kratzer committed
64
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
65
## Internal function that call multiple times strsplit() and remove space
Gilles Kratzer's avatar
Gilles Kratzer committed
66
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
67

Gilles Kratzer's avatar
Gilles Kratzer committed
68
69
70
71
72
73
strsplits <- function(x, splits, ...) {
    for (split in splits) {
        x <- unlist(strsplit(x, split, ...))
    }
    x <- gsub(" ", "", x, fixed = TRUE)  #remove space
    return(x[!x == ""])  # Remove empty values
Gilles Kratzer's avatar
Gilles Kratzer committed
74
}
Gilles Kratzer's avatar
Gilles Kratzer committed
75

Gilles Kratzer's avatar
Gilles Kratzer committed
76
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
77
78
79
80
81
82
## Internal function that produce a square matrix length(name) with {0,1} depending on f. f have to start with ~
## terms are entries of name terms are separated by + term1 | term2 indicates col(term1) row(term2) puts a 1 term1 |
## term2:term3: ... : is used as a sep . = all terms in name
##-------------------------------------------------------------------------

formula.mcmcabn <- function(f, name) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
83

Gilles Kratzer's avatar
Gilles Kratzer committed
84
    name_orignial <- name
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
85

Gilles Kratzer's avatar
Gilles Kratzer committed
86
    f <- as.character(f)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
87

Gilles Kratzer's avatar
Gilles Kratzer committed
88
    ## tests for consistence ----------------------------------------------------------------------
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
89

Gilles Kratzer's avatar
Gilles Kratzer committed
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    ## transformation name + or | or : or . or name to name_name
    if (sum((c("+", "|", ":", ".") %in% unlist(strsplit(name, split = c(""))))) != 0) {
        for (i in 1:length(name)) {
            if (sum(unlist(strsplit(name[i], split = c(""))) %in% c("+")) != 0) {
                f[[2]] <- gsub(name[i], gsub("+", "_", name[i], fixed = TRUE), f[[2]], fixed = TRUE)
                name[i] <- gsub("+", "_", name[i], fixed = TRUE)
            }
            if (sum(unlist(strsplit(name[i], split = c(""))) %in% c("|")) != 0) {
                f[[2]] <- gsub(name[i], gsub("|", "_", name[i], fixed = TRUE), f[[2]], fixed = TRUE)
                name[i] <- gsub("|", "_", name[i], fixed = TRUE)
            }
            if (sum(unlist(strsplit(name[i], split = c(""))) %in% c(":")) != 0) {
                f[[2]] <- gsub(name[i], gsub(":", "_", name[i], fixed = TRUE), f[[2]], fixed = TRUE)
                name[i] <- gsub(":", "_", name[i], fixed = TRUE)
            }
            if (sum(unlist(strsplit(name[i], split = c(""))) %in% c(".")) != 0) {
                f[[2]] <- gsub(name[i], gsub(".", "_", name[i], fixed = TRUE), f[[2]], fixed = TRUE)
                name[i] <- gsub(".", "_", name[i], fixed = TRUE)
            }
        }
Gilles Kratzer's avatar
Gilles Kratzer committed
110
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
111

Gilles Kratzer's avatar
Gilles Kratzer committed
112
113
114
115
    ## collapse name
    name.c <- paste(name, collapse = ":")
    ## Split by terms
    f.p <- strsplit(x = f[[2]], split = "+", fixed = TRUE)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
116

Gilles Kratzer's avatar
Gilles Kratzer committed
117
118
119
120
121
122
    ## nothing more than name variable in the dag formula
    tmp.test <- strsplits(x = f[[2]], splits = c("+", "|", ":", "."), fixed = TRUE)
    if (sum(!(tmp.test %in% name)) != 0) {
        stop("Formula contains some variables not used in the mcmcabn search")
    }
    ## End of tests for consistence ----------------------------------------------------------------
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
123

Gilles Kratzer's avatar
Gilles Kratzer committed
124
125
    ## creat the void matrix
    out <- matrix(data = 0, nrow = length(name), ncol = length(name))
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
126

Gilles Kratzer's avatar
Gilles Kratzer committed
127
128
    ## delete all spaces
    f.p <- gsub(" ", "", f.p[[1]], fixed = TRUE)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
129

Gilles Kratzer's avatar
Gilles Kratzer committed
130
131
    ## replace '.' by all names
    f.p.completed <- gsub(".", name.c, f.p, fixed = TRUE)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
132

Gilles Kratzer's avatar
Gilles Kratzer committed
133
    ## atomization of left term
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
134
135


Gilles Kratzer's avatar
Gilles Kratzer committed
136
137
138
    ## contruction of the output matrix
    for (i in 1:length(f.p)) {
        tmp <- f.p.completed[i]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
139

Gilles Kratzer's avatar
Gilles Kratzer committed
140
141
        ## forget unique terms -> test for |
        if (grepl("|", tmp, fixed = TRUE)) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
142

Gilles Kratzer's avatar
Gilles Kratzer committed
143
144
            ## split wrt |
            tmp.p <- strsplit(x = tmp, split = "|", fixed = TRUE)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
145

Gilles Kratzer's avatar
Gilles Kratzer committed
146
147
148
149
150
151
152
            ## test for multiple terms and contruction of the list first term
            if (grepl(":", tmp.p[[1]][1])) {
                tmp.p.p.1 <- strsplit(x = tmp.p[[1]][1], split = ":", fixed = TRUE)
            }
            if (!grepl(":", tmp.p[[1]][1])) {
                tmp.p.p.1 <- tmp.p[[1]][1]
            }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
153

Gilles Kratzer's avatar
Gilles Kratzer committed
154
155
156
157
158
159
160
            ## test for multiple terms and contruction of the list second term
            if (grepl(":", tmp.p[[1]][2])) {
                tmp.p.p.2 <- strsplit(x = tmp.p[[1]][2], split = ":", fixed = TRUE)
            }
            if (!grepl(":", tmp.p[[1]][2])) {
                tmp.p.p.2 <- tmp.p[[1]][2]
            }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
161

Gilles Kratzer's avatar
Gilles Kratzer committed
162
163
164
165
166
            ## loop over the
            for (j in 1:length(tmp.p.p.1[[1]])) {
                for (k in 1:length(tmp.p.p.2[[1]])) {
                  ## update of matrix
                  out[grep(tmp.p.p.1[[1]][j], name), grep(tmp.p.p.2[[1]][k], name)] <- 1
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
167

Gilles Kratzer's avatar
Gilles Kratzer committed
168
169
                }
            }
Gilles Kratzer's avatar
Gilles Kratzer committed
170
        }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
171

Gilles Kratzer's avatar
Gilles Kratzer committed
172
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
173

Gilles Kratzer's avatar
Gilles Kratzer committed
174
175
    ## avoid auto dependance
    diag(out) <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
176

Gilles Kratzer's avatar
Gilles Kratzer committed
177
178
    ## only 0 and 1
    out[out > 1] <- 1
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
179

Gilles Kratzer's avatar
Gilles Kratzer committed
180
181
182
183
184
    ## naming
    colnames(out) <- name_orignial
    rownames(out) <- name_orignial
    ## output
    return(out)
Gilles Kratzer's avatar
Gilles Kratzer committed
185
}
Gilles Kratzer's avatar
Gilles Kratzer committed
186
187


Gilles Kratzer's avatar
updates    
Gilles Kratzer committed
188
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
189
## Ancestor function
Gilles Kratzer's avatar
updates    
Gilles Kratzer committed
190
191
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
Gilles Kratzer committed
192
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
193
## descendent function
Gilles Kratzer's avatar
Gilles Kratzer committed
194
195
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
Gilles Kratzer committed
196
197
descendents <- function(nodes, dag) {
    if (!identical(topoSortMAT(amat = dag, index = FALSE), character(0))) {
Gilles Kratzer's avatar
Gilles Kratzer committed
198
        if (!is.integer0(which(x = dag[, nodes] == 1, arr.ind = TRUE)) && sum(diag(dag)) == 0) {
Gilles Kratzer's avatar
Gilles Kratzer committed
199
200
201
202
203
204
205
206
207
            if (length(nodes) == 1) {
                tmp <- unname(which(x = dag[, nodes] == 1, arr.ind = TRUE))
            } else {
                tmp <- t(unname(which(x = dag[, nodes] == 1, arr.ind = TRUE)))[1, ]
            }
            return(unique(unname(c(tmp, descendents(nodes = (tmp), dag = dag)))))
        } else {
            return(NULL)
        }
Gilles Kratzer's avatar
Gilles Kratzer committed
208
    }
Gilles Kratzer's avatar
Gilles Kratzer committed
209
210
    return("Not a DAG")
}
Gilles Kratzer's avatar
Gilles Kratzer committed
211

Gilles Kratzer's avatar
Gilles Kratzer committed
212
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
213
## range function
Gilles Kratzer's avatar
Gilles Kratzer committed
214
215
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
Gilles Kratzer committed
216
range01 <- function(x) {
Gilles Kratzer's avatar
update    
Gilles Kratzer committed
217
218
219
220
  if(length(x)==1){return(1)}
  if(length(x)==2){return(x-min(x))}else{
    return((x - min(x))/(max(x) - min(x)))
}}
Gilles Kratzer's avatar
updates    
Gilles Kratzer committed
221

Gilles Kratzer's avatar
Gilles Kratzer committed
222
##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
223
## not in function
Gilles Kratzer's avatar
Gilles Kratzer committed
224
225
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
Gilles Kratzer committed
226
"%!in%" <- function(x, y) !(x %in% y)
Gilles Kratzer's avatar
Gilles Kratzer committed
227
228

##-------------------------------------------------------------------------
Gilles Kratzer's avatar
Gilles Kratzer committed
229
## is.integer0
Gilles Kratzer's avatar
Gilles Kratzer committed
230
231
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
Gilles Kratzer committed
232
233
is.integer0 <- function(x) {
    is.integer(x) && length(x) == 0L
Gilles Kratzer's avatar
Gilles Kratzer committed
234
}
Gilles Kratzer's avatar
Gilles Kratzer committed
235
236
237
238
239
240
241
242
243

##-------------------------------------------------------------------------
## scoring DAGs
##-------------------------------------------------------------------------

score.dag <- function(dag,bsc.score,sc){
  n.var <- dim(dag)[1]
  score <- 0
  for (a in 1:n.var) {
Gilles Kratzer's avatar
update    
Gilles Kratzer committed
244
    sc.tmp <- sc[bsc.score$children == a, ,drop=FALSE]
Gilles Kratzer's avatar
Gilles Kratzer committed
245
246
247
248
    score <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag[a,])))), n.var + 1]), score)
  }
  return(score)
}