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