Small plotting change for black and white printing
[morpheus.git] / pkg / R / multiRun.R
CommitLineData
cbd88fe5
BA
1#' multiRun
2#'
3#' Estimate N times some parameters, outputs of some list of functions. This method is
4#' thus very generic, allowing typically bootstrap or Monte-Carlo estimations of matrices
5#' μ or β. Passing a list of functions opens the possibility to compare them on a fair
6#' basis (exact same inputs). It's even possible to compare methods on some deterministic
7#' design of experiments.
8#'
9#' @param fargs List of arguments for the estimation functions
10#' @param estimParams List of nf function(s) to apply on fargs - shared signature
11#' @param prepareArgs Prepare arguments for the functions inside estimParams
12#' @param N Number of runs
13#' @param ncores Number of cores for parallel runs (<=1: sequential)
cf673dee 14#' @param agg Aggregation method (default: lapply)
cbd88fe5
BA
15#' @param verbose TRUE to indicate runs + methods numbers
16#'
17#' @return A list of nf aggregates of N results (matrices).
18#'
19#' @examples
1b53e3a5 20#' \donttest{
cbd88fe5
BA
21#' β <- matrix(c(1,-2,3,1),ncol=2)
22#'
23#' # Bootstrap + computeMu, morpheus VS flexmix ; assumes fargs first 3 elts X,Y,K
24#' io <- generateSampleIO(n=1000, p=1/2, β=β, b=c(0,0), "logit")
25#' μ <- normalize(β)
26#' res <- multiRun(list(X=io$X,Y=io$Y,optargs=list(K=2,jd_nvects=0)), list(
27#' # morpheus
28#' function(fargs) {
29#' library(morpheus)
30#' ind <- fargs$ind
31#' computeMu(fargs$X[ind,],fargs$Y[ind],fargs$optargs)
32#' },
33#' # flexmix
34#' function(fargs) {
35#' library(flexmix)
36#' ind <- fargs$ind
37#' K <- fargs$optargs$K
38#' dat = as.data.frame( cbind(fargs$Y[ind],fargs$X[ind,]) )
39#' out = refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K,
40#' model=FLXMRglm(family="binomial") ) )
41#' normalize( matrix(out@@coef[1:(ncol(fargs$X)*K)], ncol=K) )
42#' } ),
5fc1b9d9
BA
43#' prepareArgs = function(fargs,index) {
44#' if (index == 1)
45#' fargs$ind <- 1:nrow(fargs$X)
46#' else
47#' fargs$ind <- sample(1:nrow(fargs$X),replace=TRUE)
cbd88fe5
BA
48#' fargs
49#' }, N=10, ncores=3)
50#' for (i in 1:2)
51#' res[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact")
52#'
53#' # Monte-Carlo + optimParams from X,Y, morpheus VS flexmix ; first args n,p,β,b
54#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),optargs=list(link="logit")),list(
55#' # morpheus
56#' function(fargs) {
57#' library(morpheus)
58#' K <- fargs$optargs$K
59#' μ <- computeMu(fargs$X, fargs$Y, fargs$optargs)
60#' V <- list( p=rep(1/K,K-1), β=μ, b=c(0,0) )
61#' optimParams(V,fargs$optargs)$β
62#' },
63#' # flexmix
64#' function(fargs) {
65#' library(flexmix)
66#' K <- fargs$optargs$K
67#' dat <- as.data.frame( cbind(fargs$Y,fargs$X) )
68#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K,
69#' model=FLXMRglm(family="binomial") ) )
70#' sapply( seq_len(K), function(i) as.double( out@@components[[1]][[i]][,1] ) )
71#' } ),
5fc1b9d9 72#' prepareArgs = function(fargs,index) {
cbd88fe5
BA
73#' library(morpheus)
74#' io = generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$optargs$link)
75#' fargs$X = io$X
76#' fargs$Y = io$Y
77#' fargs$optargs$K = ncol(fargs$β)
78#' fargs$optargs$M = computeMoments(io$X,io$Y)
79#' fargs
80#' }, N=10, ncores=3)
81#' for (i in 1:2)
82#' res[[i]] <- alignMatrices(res[[i]], ref=β, ls_mode="exact")}
83#' @export
84multiRun <- function(fargs, estimParams,
5fc1b9d9 85 prepareArgs = function(x,i) x, N=10, ncores=3, agg=lapply, verbose=FALSE)
cbd88fe5
BA
86{
87 if (!is.list(fargs))
88 stop("fargs: list")
89 # No checks on fargs: supposedly done in estimParams[[i]]()
90 if (!is.list(estimParams))
91 estimParams = list(estimParams)
92 # Verify that the provided parameters estimations are indeed functions
93 lapply(seq_along(estimParams), function(i) {
94 if (!is.function(estimParams[[i]]))
95 stop("estimParams: list of function(fargs)")
96 })
97 if (!is.numeric(N) || N < 1)
98 stop("N: positive integer")
99
100 estimParamAtIndex <- function(index)
101 {
5fc1b9d9 102 fargs <- prepareArgs(fargs, index)
cbd88fe5
BA
103 if (verbose)
104 cat("Run ",index,"\n")
105 lapply(seq_along(estimParams), function(i) {
106 if (verbose)
107 cat(" Method ",i,"\n")
108 out <- estimParams[[i]](fargs)
109 if (is.list(out))
110 do.call(rbind, out)
111 else
112 out
113 })
114 }
115
116 if (ncores > 1)
117 {
118 cl = parallel::makeCluster(ncores, outfile="")
119 parallel::clusterExport(cl, c("fargs","verbose"), environment())
120 list_res = parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex)
121 parallel::stopCluster(cl)
122 }
123 else
124 list_res = lapply(1:N, estimParamAtIndex)
125
126 # De-interlace results: output one list per function
127 nf <- length(estimParams)
128 lapply( seq_len(nf), function(i) lapply(seq_len(N), function(j) list_res[[j]][[i]]) )
129}