reconhub / outbreaker2

@@ -33,7 +33,7 @@
Loading
33 33
34 34
  ## remove the contact likelihood from outlier detection
35 35
  tmp_likelihoods <- likelihoods
36 -
  tmp_likelihoods$contact <- function(data, param) return(0)
36 +
  tmp_likelihoods$contact <- list(function(data, param) return(0), 0)
37 37
38 38
  ## RUN MCMC ##
39 39
  for (i in seq.int(2, config$n_iter_import, 1)) {

@@ -56,7 +56,11 @@
Loading
56 56
  size_t N = static_cast<size_t>(data["N"]);
57 57
  if (N < 2) return 0.0;
58 58
59 -
  if (custom_function == R_NilValue) {
59 +
  Rcpp::List l;
60 +
  if (custom_function != R_NilValue) {
61 +
    l = Rcpp::as<Rcpp::List>(custom_function);
62 +
  }
63 +
  if (custom_function == R_NilValue || (l.size() > 0 && l[0] == R_NilValue)) {
60 64
61 65
    // Variables from the data & param
62 66
    Rcpp::NumericMatrix w_dens = data["log_w_dens"];
@@ -169,8 +173,9 @@
Loading
169 173
    return(out);
170 174
    
171 175
  } else { // use of a customized likelihood function
172 -
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(custom_function);
173 -
176 +
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(l[0]);
177 +
    int arity = l[1];
178 +
    if (arity == 3) return Rcpp::as<double>(f(data, param, i));
174 179
    return Rcpp::as<double>(f(data, param));
175 180
  }
176 181
}
@@ -198,7 +203,11 @@
Loading
198 203
  size_t N = static_cast<size_t>(data["N"]);
199 204
  if(N < 2) return 0.0;
200 205
201 -
  if (custom_function == R_NilValue) {
206 +
  Rcpp::List l;
207 +
  if (custom_function != R_NilValue) {
208 +
    l = Rcpp::as<Rcpp::List>(custom_function);
209 +
  }
210 +
  if (custom_function == R_NilValue || (l.size() > 0 && l[0] == R_NilValue)) {
202 211
203 212
    Rcpp::IntegerVector alpha = param["alpha"];
204 213
    Rcpp::IntegerVector t_inf = param["t_inf"];
@@ -246,8 +255,9 @@
Loading
246 255
247 256
    return out;
248 257
  } else { // use of a customized likelihood function
249 -
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(custom_function);
250 -
258 +
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(l[0]);
259 +
    int arity = l[1];
260 +
    if (arity == 3) return Rcpp::as<double>(f(data, param, i));
251 261
    return Rcpp::as<double>(f(data, param));
252 262
  }
253 263
}
@@ -275,7 +285,11 @@
Loading
275 285
  size_t N = static_cast<size_t>(data["N"]);
276 286
  if(N < 2) return 0.0;
277 287
278 -
  if (custom_function == R_NilValue) {
288 +
  Rcpp::List l;
289 +
  if (custom_function != R_NilValue) {
290 +
    l = Rcpp::as<Rcpp::List>(custom_function);
291 +
  }
292 +
  if (custom_function == R_NilValue || (l.size() > 0 && l[0] == R_NilValue)) {
279 293
280 294
    Rcpp::IntegerVector dates = data["dates"];
281 295
    Rcpp::IntegerVector t_inf = param["t_inf"];
@@ -308,8 +322,9 @@
Loading
308 322
309 323
    return out;
310 324
  }  else { // use of a customized likelihood function
311 -
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(custom_function);
312 -
325 +
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(l[0]);
326 +
    int arity = l[1];
327 +
    if (arity == 3) return Rcpp::as<double>(f(data, param, i));
313 328
    return Rcpp::as<double>(f(data, param));
314 329
  }
315 330
}
@@ -354,7 +369,11 @@
Loading
354 369
    return R_NegInf;
355 370
  }
356 371
357 -
  if (custom_function == R_NilValue) {
372 +
  Rcpp::List l;
373 +
  if (custom_function != R_NilValue) {
374 +
    l = Rcpp::as<Rcpp::List>(custom_function);
375 +
  }
376 +
  if (custom_function == R_NilValue || (l.size() > 0 && l[0] == R_NilValue)) {
358 377
359 378
    double out = 0.0;
360 379
@@ -385,8 +404,9 @@
Loading
385 404
386 405
    return out;
387 406
  } else { // use of a customized likelihood function
388 -
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(custom_function);
389 -
407 +
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(l[0]);
408 +
    int arity = l[1];
409 +
    if (arity == 3) return Rcpp::as<double>(f(data, param, i));
390 410
    return Rcpp::as<double>(f(data, param));
391 411
  }
392 412
}
@@ -435,7 +455,11 @@
Loading
435 455
  size_t N = static_cast<size_t>(data["N"]);
436 456
  if (N < 2) return 0.0;
437 457
438 -
  if (custom_function == R_NilValue) {
458 +
  Rcpp::List l;
459 +
  if (custom_function != R_NilValue) {
460 +
    l = Rcpp::as<Rcpp::List>(custom_function);
461 +
  }
462 +
  if (custom_function == R_NilValue || (l.size() > 0 && l[0] == R_NilValue)) {
439 463
440 464
    double out;
441 465
    double eps = Rcpp::as<double>(param["eps"]);
@@ -497,8 +521,9 @@
Loading
497 521
    return out;
498 522
    
499 523
  } else { //use of a customized likelihood function
500 -
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(custom_function);
501 -
524 +
    Rcpp::Function f = Rcpp::as<Rcpp::Function>(l[0]);
525 +
    int arity = l[1];
526 +
    if (arity == 3) return Rcpp::as<double>(f(data, param, i));
502 527
    return Rcpp::as<double>(f(data, param));
503 528
  }
504 529
}

@@ -39,10 +39,11 @@
Loading
39 39
#'
40 40
#'
41 41
#' @return
42 -
#' A named list of functions with the class \code{custom_likelihood}, each
43 -
#'     implementing a customised log-likelihood components of
44 -
#'     outbreaker. Functions which are not customised will result in a NULL
45 -
#'     component.
42 +
#' A named list of list(function, arity) pairs with the class
43 +
#'     \code{custom_likelihood}, each function implementing a customised
44 +
#'     log-likelihood component of outbreaker. Functions which are not
45 +
#'     customised will result in a list(NULL, 0) component. Any function with
46 +
#'     arity 3 must have the third parameter default to NULL
46 47
#'
47 48
#' @author Thibaut Jombart (\email{thibautjombart@@gmail.com})
48 49
#'
@@ -88,13 +89,17 @@
Loading
88 89
## Likelihood functions in outbreaker2 are implemented using Rcpp. However,
89 90
## these functions can also be replaced by customized functions. These can be
90 91
## specified by the user, through the '...' argument of
91 -
## 'custom_likelihoods'. These functions must have 2 arguments:
92 +
## 'custom_likelihoods'. These functions must have at least 2 arguments:
92 93
93 94
## - data: a valid 'outbreaker_data' list
94 95
95 96
## - param: a list containing current parameter states, as returned by
96 97
## - create_param
97 98
99 +
## - [i=NULL]: (optional) a list of the cases for which the loglikelihoods
100 +
## - should be calculated. Needs to default to `NULL` in which case the
101 +
## - loglikelihood of the entire tree is calculated.
102 +
98 103
custom_likelihoods <- function(...) {
99 104
100 105
    ll_functions <- list(...)
@@ -121,40 +126,49 @@
Loading
121 126
    function_or_null <- function(x) {
122 127
        is.null(x) || is.function(x)
123 128
    }
129 +
    list_function_or_null <- function(x) {
130 +
        is.list(x) && (is.null(x[[1]]) || is.function(x[[1]]))
131 +
    }
124 132
125 -
    is_ok <- vapply(likelihoods, function_or_null, logical(1))
133 +
    # Ensure that custom_likelihoods(l) == custom_likelihoods(custom_likelihoods(l))
134 +
    is_list_function_or_null <- vapply(likelihoods, list_function_or_null, logical(1))
135 +
    is_function_or_null <- vapply(likelihoods, function_or_null, logical(1))
126 136
127 -
    if (!all(is_ok)) {
128 -
        culprits <- likelihoods_names[!is_ok]
137 +
    if (!all(is_function_or_null) & !all(is_list_function_or_null)) {
138 +
        culprits <- likelihoods_names[!is_function_or_null]
129 139
        msg <- paste0("The following likelihoods are not functions: ",
130 140
                      paste(culprits, collapse = ", "))
131 141
        stop(msg)
132 142
    }
133 143
144 +
    # If the arity of the likelihood functions is three, the last argumen should
145 +
    # be the (1-based) indices of the cases we're currently perturbing. This
146 +
    # allows us to calculate the local likelihood delta, rather than having to
147 +
    # calculate the likelihood of the entire tree twice for every single
148 +
    # perturbation we make.
149 +
    if (!all(is_list_function_or_null)) {
150 +
        likelihoods <- lapply(likelihoods, function(x) { if (is.null(x)) return(list(x, 0)); list(x, length(methods::formalArgs(x))) })
151 +
    }
134 152
135 -
    ## check they all have a single argument
136 -
137 -
    with_two_args <- function(x) {
138 -
        if(is.function(x)) {
139 -
            return (length(methods::formalArgs(x)) == 2L)
153 +
    arity_two_or_three <- function(x) {
154 +
        if (is.function(x[[1]])) {
155 +
            return (x[[2]] == 2L | x[[2]] == 3L)
140 156
        }
141 -
142 -
        return(TRUE)
157 +
        return(T)
143 158
    }
144 159
145 -
    two_args <- vapply(likelihoods, with_two_args, logical(1))
160 +
    legal_arity <- vapply(likelihoods, arity_two_or_three, logical(1))
146 161
147 -
    if (!all(two_args)) {
148 -
        culprits <- likelihoods_names[!two_args]
149 -
        msg <- paste0("The following likelihoods dont' have two arguments: ",
150 -
                      paste(culprits, collapse = ", "))
162 +
    if (!all(legal_arity)) {
163 +
        culprits <- likelihood_names[!legal_arity]
164 +
        msg <- paste0("The following likelihoods do not have arity two or three: ",
165 +
                      paste(culprits, collapse=", "))
151 166
        stop(msg)
152 167
    }
153 168
154 -
169 +
    names(likelihoods) <- likelihoods_names
155 170
    class(likelihoods) <- c("custom_likelihoods", "list")
156 171
    return(likelihoods)
157 -
158 172
}
159 173
160 174
Files Coverage
R 60.12%
src 92.23%
Project Totals (25 files) 69.92%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading