[Dream-commits] r22 - in pkg: R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 19 03:30:28 CET 2010
Author: josephguillaume
Date: 2010-03-19 03:30:25 +0100 (Fri, 19 Mar 2010)
New Revision: 22
Added:
pkg/demo/FME.nonlinear.model_parallelisation.R
Modified:
pkg/R/CompDensity.R
pkg/R/dream.R
pkg/R/plot.dream.R
pkg/R/summary.dream.R
pkg/demo/00Index
pkg/demo/FME.nonlinear.model.R
pkg/man/coef.dream.Rd
pkg/man/dream.Rd
Log:
Change plot to allow non-screen output by printing lattice objects.
Foreach now correctly uses %dopar%
dream checks that FUN returns an object of class numeric
syntax changes to work with a single parameter
Summary: Pre-compute coda output, so all results then appear at once
Clarified compdensity multicore class error
Changed default back to multicore - requires no other setup
foreach incurs substantial overhead. Implemented snow loop.
Made changes in CompDensity to avoid parameter clashes, allow variables on worker nodes.
Adapted FME as new parallelisation example, with e.g. output times
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/R/CompDensity.R 2010-03-19 02:30:25 UTC (rev 22)
@@ -21,7 +21,7 @@
## TODO: more appropriate naming of options?
## TODO: allow shortenings of option?
CompDensity <- function(pars,control,FUN,func.type,
- measurement=NULL,...){
+ measurement=NULL,FUN.pars=NULL){
## Should be guaranteed by dream
## stopifnot(!is.null(measurement) || func.type%in% c("posterior.density","logposterior.density"))
@@ -35,11 +35,11 @@
## SSR scalar
## temp. list of length nseq with elements of length 2: p and logp
- do.calc <- function (ii){
-
+ do.calc <- function (pp,control,MFUN,func.type,measurement,FUN.pars){
## Call model to generate simulated data
- ## TODO: correct use of optional pars?
- modpred <- FUN(pars[ii,],...)
+ FUN.pars[[names(formals(FUN))[1]]] <- pp
+ modpred <- do.call(MFUN,FUN.pars)
+ ##stopifnot(inherits(modpred,"numeric"))
switch(func.type,
## Model directly computes posterior density
@@ -61,7 +61,6 @@
## Model computes output simulation
## TODO: may need as.numeric
calc.rmse={
-
err <- as.numeric(measurement$data-modpred)
## Derive the sum of squared error
SSR <- sum(abs(err)^(2/(1+control$gamma)))
@@ -88,15 +87,32 @@
}) ##switch
c(p,logp)
}
+
+ pars <- lapply(1:nrow(pars),function(x) pars[x,])
+ switch(control$parallel,
+ "multicore"={
+ temp <- mclapply(pars,do.calc,control=control,MFUN=FUN,func.type=func.type,
+ measurement=measurement,FUN.pars=FUN.pars,
+ mc.preschedule=FALSE)
+ },
+ "foreach"={
+ ## foreach(pp=iter(pars,by="row")) %dopar%
+ temp <- foreach(pp=pars) %dopar% do.calc(pp,control=control,MFUN=FUN,func.type=func.type,
+ measurement=measurement,FUN.pars=FUN.pars)
+ },
+ "snow"={
+ temp <- clusterApplyLB(cl=cl,x=pars,fun=do.calc,
+ control=control,MFUN=FUN,func.type=func.type,
+ measurement=measurement,FUN.pars=FUN.pars)
+ },
+ temp <- lapply(pars,FUN=do.calc,control=control,MFUN=FUN,func.type=func.type,
+ measurement=measurement,FUN.pars=FUN.pars)
+ )##switch parallel
- if (control$use.multicore) temp <- mclapply(1:nrow(pars),do.calc,mc.preschedule=FALSE)
- else if (control$use.foreach) temp <- foreach(ii=1:nrow(pars)) %do% do.calc(ii)
- else temp <- lapply(1:nrow(pars),do.calc)
-
p <- sapply(temp,function(x) x[1])
logp <- sapply(temp,function(x) x[2])
- if(class(p)!="numeric") stop("Error with multicore, set control$use.multicore=FALSE to turn it off")
+ if(class(p)!="numeric") stop(sprintf("Expected class numeric, got class %s. Error with multicore? Set control$use.multicore=FALSE to turn it off",class(p)))
stopifnot(!any(is.na(p)))
##stopifnot(!any(is.na(logp))) ##Not used anyway
return(list(p=p,logp=logp))
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/R/dream.R 2010-03-19 02:30:25 UTC (rev 22)
@@ -18,8 +18,7 @@
thin.t=NA, ## parameter for reduced sample collection
### Efficiency improvements
REPORT = 1000, ## approximate number of function evaluations between reports. >0. 0=none
- use.multicore=FALSE,
- use.foreach=TRUE,
+ parallel=c("multicore","snow","foreach"), ##packages to use for parallel in order of preference: multicore,snow,foreach
### Parameters with auto-set values
ndim=NA, ## number of parameters (automatically set from length of pars)
DEpairs = NA, ## Number of DEpairs. defaults to max val floor((nseq-1)/2)
@@ -101,16 +100,19 @@
## Check INIT and FUN have required extra parameters in INIT.pars & FUN.pars
req.args.init <- names(formals(INIT))
req.args.FUN <- names(formals(FUN))
- req.args.FUN <- req.args.FUN[2:length(req.args.FUN)] ##optional pars only
-
- if(!all(req.args.init %in% c("pars","nseq",names(INIT.pars))))
- stop(paste(c("INIT Missing extra arguments:",
- req.args.init[!req.args.init %in% c("pars","nseq",names(INIT.pars))]),
- sep=" ",collapse=" "))
- if(!all(req.args.FUN %in% names(FUN.pars))) stop(paste(c("FUN Missing extra arguments:",
- req.args.FUN[!req.args.FUN %in% c(names(FUN.pars))]),
- collapse=" "))
+ if(!all(req.args.init %in% c("pars","nseq",names(INIT.pars))))
+ stop(paste(c("INIT Missing extra arguments:",
+ req.args.init[!req.args.init %in% c("pars","nseq",names(INIT.pars))]),
+ sep=" ",collapse=" "))
+
+ if (length(req.args.FUN)<length(FUN.pars)+1) stop("Some FUN.pars are not required by FUN")
+ if (length(req.args.FUN)>1){
+ req.args.FUN <- req.args.FUN[2:length(req.args.FUN)] ##optional pars only
+ if(!all(req.args.FUN %in% names(FUN.pars))) stop(paste(c("FUN Missing extra arguments:",
+ req.args.FUN[!req.args.FUN %in% c(names(FUN.pars))]),
+ collapse=" "))
+ }
## Update default settings with supplied settings
@@ -135,8 +137,19 @@
if (control$burnin.length<1) control$burnin.length <- control$burnin.length*control$ndraw
- if (control$use.multicore) control$use.multicore <- require(multicore)
- if (control$use.foreach) control$use.foreach <- require(foreach)
+ ##Choice of parallel backend
+ if (control$parallel!="none"){
+ parallel <- "none"
+ for (p in control$parallel) {
+ if (require(p,character.only=TRUE)) {
+ parallel <- p
+ break
+ }
+ }
+ if (parallel=="none") warning(sprintf("Requested parallel backends not available (%s)",paste(control$parallel,collapse=",")))
+ control$parallel <- parallel
+ }
+
## Check validity of settings
if (control$DEpairs==0) stop("control$DEpairs set to 0. Increase nseq?")
@@ -252,12 +265,18 @@
x <- do.call(INIT,modifyList(INIT.pars,list(pars=pars,nseq=NSEQ)))
+ ## Test that FUN returns numeric
+ test.pars <- FUN.pars
+ test.pars[[names(formals(FUN))[1]]] <- x[1,]
+ modpred <- do.call(FUN,test.pars)
+ if (!inherits(modpred,"numeric")) stop(sprintf("Result of FUN should be of class numeric, not %s",class(modpred)))
+
## make each element of pars a list and extract lower / upper
lower <- sapply(pars, function(x) min(x[[1]]))
upper <- sapply(pars, function(x) max(x[[1]]))
##Step 2: Calculate posterior density associated with each value in x
- tmp<-do.call(CompDensity,modifyList(FUN.pars,list(pars=x,control=control,FUN=FUN,func.type=func.type,measurement=measurement)))
+ tmp<-do.call(CompDensity,list(pars=x,control=control,FUN=FUN,func.type=func.type,measurement=measurement,FUN.pars=FUN.pars))
##Save the initial population, density and log density in one list X
X<-cbind(x=x,p=tmp$p,logp=tmp$logp)
@@ -288,7 +307,7 @@
counter.thin <- counter.thin + 1
## Define the current locations and associated posterior densities
- x.old <- X[,1:NDIM]
+ x.old <- X[,1:NDIM,drop=FALSE]
p.old <- X[,NDIM+1]
logp.old <- X[,NDIM+2]
@@ -303,7 +322,7 @@
CR[,gen.number] <- tmp$CR
## Now compute the likelihood of the new points
- tmp<-do.call(CompDensity,modifyList(FUN.pars,list(pars=x.new,control=control,FUN=FUN,func.type=func.type,measurement=measurement)))
+ tmp<-do.call(CompDensity,list(pars=x.new,control=control,FUN=FUN,func.type=func.type,measurement=measurement,FUN.pars=FUN.pars))
p.new <- tmp$p
logp.new <- tmp$logp
@@ -337,9 +356,9 @@
## Calculate the standard deviation of each dimension of X
## TODO: matlab syntax is unclear - seems to be columnwise
## element-wise: sd(c(X[,1:NDIM]))
- r <- apply(X[,1:NDIM],2,sd)
+ r <- apply(X[,1:NDIM,drop=FALSE],2,sd)
## Compute the Euclidean distance between new X and old X
- delta.normX <- rowSums(((x.old-X[,1:NDIM])/r)^2)
+ delta.normX <- rowSums(((x.old-X[,1:NDIM,drop=FALSE])/r)^2)
## Use this information to update sum_p2 to update N_CR
delta.tot <- CalcDelta(NCR,delta.tot,delta.normX,CR[,gen.number])
}
Modified: pkg/R/plot.dream.R
===================================================================
--- pkg/R/plot.dream.R 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/R/plot.dream.R 2010-03-19 02:30:25 UTC (rev 22)
@@ -13,8 +13,8 @@
plot(ss)
- xyplot(ss)
- densityplot(ss)
+ print(xyplot(ss))
+ print(densityplot(ss))
## Acceptance rate
plot(table(x$AR[,2]),main="Distribution of % acceptance rate")
@@ -29,7 +29,7 @@
## Multi-variate density for first chain
- splom(as.data.frame(x$Sequences[[1]]),
+ print(splom(as.data.frame(x$Sequences[[1]]),
upper.panel = panel.smoothScatter, nrpoints = 0,
lower.panel = function(x, y, ...) {
panel.grid(-1, -1)
@@ -37,6 +37,6 @@
panel.loess(x, y, span = 2/3, lwd = 2)
grid::grid.text(paste("cor =", round(cor(x, y),2)),
y = 0.1)
- })
+ }))
}##plot.dream
Modified: pkg/R/summary.dream.R
===================================================================
--- pkg/R/summary.dream.R 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/R/summary.dream.R 2010-03-19 02:30:25 UTC (rev 22)
@@ -1,6 +1,8 @@
summary.dream <- function(object,...){
+ coda.sum <- summary(window(object$Sequences, start = end(object$Sequences)/2 + 1))
+
cat(sprintf("
Exit message: %s
Num fun evals: %d
@@ -28,11 +30,12 @@
cat("
CODA summary for last 50% of MCMC chains:
")
- print(summary(window(object$Sequences, start = end(object$Sequences)/2 + 1)))
+ print(coda.sum)
+
cat("
Acceptance Rate
")
- summary(object$AR[,2])
-
+ print(summary(object$AR[,2]))
+
} ##summary.dream
Modified: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/demo/00Index 2010-03-19 02:30:25 UTC (rev 22)
@@ -1,3 +1,4 @@
FME.nonlinear.model Nonlinear model calibration as with FME vignette.
+FME.nonlinear.model_parallelisation Example of parallelisation options
run_example1 Function call for n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
example1 Setup for n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/demo/FME.nonlinear.model.R 2010-03-19 02:30:25 UTC (rev 22)
@@ -44,13 +44,13 @@
unloadNamespace("dream")
library(dream)
+Model.y <- function(p,x) p[1]*x/(x+p[2])
-Model.y <- function(p,x) as.ts(p[1]*x/(x+p[2]))
-
set.seed(456)
control <- list(
- nseq=4
+ nseq=4,
+ parallel="none"
## REPORT=0
## ndraw=1000
## Rthres=1+1e-3
@@ -58,16 +58,6 @@
pars <- list(p1=c(0,1),p2=c(0,100))
-
-## Optional parallelisation of function evaluations
-## Uses foreach, which allows Rmpi, SNOW or multicore to be registered,
-## or just normal sequential evaluation
-## Example here is for SNOW with a socket cluster, which doesn't require any other software or setup
-if (require(doSNOW)){
- cl <- makeCluster(c("localhost","localhost"), type = "SOCK")
- registerDoSNOW(cl)
-}
-
dd <- dream(
FUN=Model.y, func.type="calc.rmse",
pars = pars,
@@ -79,8 +69,6 @@
control = control
)
-if (require(doSNOW)) stopCluster(cl)
-
dd
summary(dd)
Added: pkg/demo/FME.nonlinear.model_parallelisation.R
===================================================================
--- pkg/demo/FME.nonlinear.model_parallelisation.R (rev 0)
+++ pkg/demo/FME.nonlinear.model_parallelisation.R 2010-03-19 02:30:25 UTC (rev 22)
@@ -0,0 +1,66 @@
+
+## Example from FME vignette, p21, section 7. Soetaert & Laine
+## Nonlinear model parameter estimation
+## TODO: document more clearly, better outputs
+##
+
+## http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/FME/inst/doc/FMEmcmc.Rnw?rev=96&root=fme&view=markup
+
+Obs <- data.frame(x=c( 28, 55, 83, 110, 138, 225, 375), # mg COD/l
+ y=c(0.053,0.06,0.112,0.105,0.099,0.122,0.125)) # 1/hour
+Obs2<- data.frame(x=c( 20, 55, 83, 110, 138, 240, 325), # mg COD/l
+ y=c(0.05,0.07,0.09,0.10,0.11,0.122,0.125)) # 1/hour
+obs.all <- rbind(Obs,Obs2)
+
+##########################
+## DREAM results
+
+unloadNamespace("dream")
+library(dream)
+
+Model.y <- function(p,x) p[1]*x/(x+p[2])
+
+control <- list(
+ nseq=4
+ )
+
+pars <- list(p1=c(0,1),p2=c(0,100))
+
+## Optional parallelisation of function evaluations
+## Uses foreach, which allows Rmpi, SNOW or multicore to be registered,
+## or just normal sequential evaluation
+## Example here is for SNOW with a socket cluster, which doesn't require any other software or setup
+## Parallelisation incurs an overhead and is not worthwhile for fast functions
+
+
+if (require(doSNOW)){
+ cl <- makeCluster(2, type = "SOCK")
+ registerDoSNOW(cl)
+}
+
+for (p in c("none","snow","foreach")){
+
+ set.seed(456)
+ control$parallel <- p
+
+ dd <- dream(
+ FUN=Model.y, func.type="calc.rmse",
+ pars = pars,
+ FUN.pars=list(
+ x=obs.all$x
+ ),
+ INIT = LHSInit,
+ measurement=list(data=obs.all$y),
+ control = control
+ )
+
+ print(dd$control$parallel)
+ print(coef(dd))
+ print(dd$time)
+}
+if (require(doSNOW)) stopCluster(cl)
+
+## Seconds taken:
+## none: 6.234, 6.219
+## snow: 12.234, 12.125
+## foreach: 107.999, 107.893
Property changes on: pkg/demo/FME.nonlinear.model_parallelisation.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/man/coef.dream.Rd
===================================================================
--- pkg/man/coef.dream.Rd 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/man/coef.dream.Rd 2010-03-19 02:30:25 UTC (rev 22)
@@ -16,3 +16,9 @@
vector=f(dream object)}
\item{...}{Unused. Matches generic}
}
+\details{
+ e.g. of using arbitrary function for method: 20\% quantile.
+ \code{
+ coef(object,method=function(sss) apply(as.matrix(sss),2,quantile,0.2)
+ }
+}
\ No newline at end of file
Modified: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd 2010-03-12 07:37:35 UTC (rev 21)
+++ pkg/man/dream.Rd 2010-03-19 02:30:25 UTC (rev 22)
@@ -64,12 +64,11 @@
Frequent calculation is likely to slow performance. The summary
function can use \code{\link{gelman.diag}} to calculate convergence statistic after completion
\cr
- \code{use.foreach} \tab \code{TRUE} \tab Whether to use foreach
- package if available. For parallelisation of function evaluations, a
- parallel backend must be registered from one of the doMC, doMPI,
+ \code{parallel} \tab multicore, snow, foreach \tab Character vector
+ of packages to
+ use for parallelisation of function evaluations, in order of
+ preference. Set to \code{"none"} to not parallelise. For foreach, a backend must be registered from one of the doMC, doMPI,
doSNOW packages. \cr
- \code{use.multicore} \tab \code{FALSE} \tab Whether to use multicore
- package if available. Parallelises function evaluations. \cr
\code{ndraw} \tab 1e5 \tab Maximum number of function evaluations. May
terminate before convergence. \cr
\code{maxtime} \tab Inf \tab Maximum duration of optimization in seconds. May
More information about the Dream-commits
mailing list