[Dream-commits] r34 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 10 07:44:31 CET 2012
Author: josephguillaume
Date: 2012-01-10 07:44:31 +0100 (Tue, 10 Jan 2012)
New Revision: 34
Added:
pkg/demo/parallelisation_chain_id.R
pkg/man/compareToMatlab.Rd
pkg/man/getMatlabControl.Rd
pkg/man/getMatlabSeq.Rd
pkg/man/likelihood_functions.Rd
pkg/man/plotMCMCQQ.Rd
pkg/man/snow.chains.Rd
pkg/man/writeMatlabDREAMSettings.Rd
Removed:
pkg/R/maxLikPars.R
pkg/man/calc.loglik.Rd
pkg/man/calc.rmse.Rd
pkg/man/calc.weighted.rmse.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/CompDensity.R
pkg/R/dream.R
pkg/R/dreamCalibrate.R
pkg/R/window.dream.R
pkg/demo/00Index
pkg/demo/FME.nonlinear.model_parallelisation.R
pkg/man/coef.dream.Rd
pkg/man/dream.Rd
pkg/man/plot.dream.Rd
pkg/man/predict.dream.Rd
pkg/man/simulate.dream.Rd
pkg/man/window.dream.Rd
Log:
Fixed Cb,Wb overwriting user input.
Added parallelisation="snow.chains" and demo
Passed R CMD check
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/DESCRIPTION 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,6 +1,6 @@
Package: dream
-Version: 0.3-3
-Date: 2010-12-14
+Version: 0.4
+Date: 2011-01-10
Title: DiffeRential Evolution Adaptive Metropolis
Author: Joseph Guillaume and Felix Andrews, based on MATLAB code by Jasper Vrugt
Maintainer: Joseph Guillaume <joseph.guillaume at anu.edu.au>
@@ -8,7 +8,8 @@
Based on the original MATLAB code written by Jasper Vrugt.
Methods for calibration and prediction using the DREAM algorithm
Depends: coda
-Suggests: multicore, foreach, SNOW, doSNOW, doMPI, doMC
+Suggests: snow, doSNOW, R.matlab, lattice, reshape, mgcv
+Enhances: multicore, foreach, doMPI, doMC
Imports: stats, utils
License: file LICENSE
URL: http://dream.r-forge.r-project.org/
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/NAMESPACE 2012-01-10 06:44:31 UTC (rev 34)
@@ -12,7 +12,6 @@
S3method(print,dream)
S3method(simulate,dream)
-export(maxLikPars)
export(maxLikCoda)
export(dreamCalibrate)
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/R/CompDensity.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -25,9 +25,21 @@
## Should be guaranteed by dream
## stopifnot(!is.null(measurement) || func.type%in% c("posterior.density","logposterior.density"))
-
+
stopifnot(!any(is.na(pars)))
-
+
+ if (control$parallel=="snow.chains"){
+ ## Custom code for John Joseph 6/8/2011
+ ## FUN should be logp=f(cluster instance identifier, list of pars)
+ ## func.type,measurement,FUN.pars should all be NULL
+ logp <- as.numeric(clusterApply(cl=cl,x=1:nrow(pars),fun=FUN,pars=pars))
+ if (any(is.na(logp))) {
+ stop("likelihood function produced invalid probabilities (NA/NaN)")
+ }
+ ##stopifnot(!any(is.na(logp))) ##Not used anyway
+ return(list(p=logp,logp=logp))
+ }
+
## dimensions:
## i. iter 1:nseq
## modpred. scalar or vector commensurate to measurement$data
@@ -48,10 +60,10 @@
if (any(modpred<0)) stop("Posterior density returned by FUN should be positive. Otherwise use logposterior.density?")
logp <- log(modpred)
},
- ## Model computes output simulation
+ ## Model computes output simulation
calc.loglik={
err <- as.numeric(measurement$data-modpred)
-
+
## Derive the log likelihood
logp <- measurement$N*log(control$Wb/measurement$sigma)-
control$Cb*(sum((abs(err/measurement$sigma))^(2/(1+control$gamma))))
@@ -67,7 +79,7 @@
## And retain in memory
p <- -SSR
logp <- -0.5*SSR
-
+
},
## Model directly computes log posterior density
logposterior.density={
@@ -98,8 +110,8 @@
},
"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)
+ 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,
@@ -109,7 +121,7 @@
temp <- lapply(pars,FUN=do.calc,control=control,MFUN=FUN,func.type=func.type,
measurement=measurement,FUN.pars=FUN.pars)
)##switch parallel
-
+
p <- sapply(temp,function(x) x[1])
logp <- sapply(temp,function(x) x[2])
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/R/dream.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,40 +1,40 @@
-## Copyright (c) 2008, Los Alamos National Security, LLC
-## All rights reserved.
-##
-## Copyright 2008. Los Alamos National Security, LLC. This software was produced under U.S.
-## Government contract DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is
-## operated by Los Alamos National Security, LLC for the U.S. Department of Energy. The U.S.
-## Government has rights to use, reproduce, and distribute this software.
-##
-## NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES A NY WARRANTY, EXPRESS OR
-## IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is modified to
-## produce derivative works, such modified software should be clearly marked, so as not to
-## confuse it with the version available from LANL.
-##
-## Additionally, redistribution and use in source and binary forms, with or without
-## modification, are permitted provided that the following conditions are met:
-## * Redistributions of source code must retain the above copyright notice, this list of
-## conditions and the following disclaimer.
-## * Redistributions in binary form must reproduce the above copyright notice, this list of
-## conditions and the following disclaimer in the documentation and/or other materials
-## provided with the distribution.
+## Copyright (c) 2008, Los Alamos National Security, LLC
+## All rights reserved.
+##
+## Copyright 2008. Los Alamos National Security, LLC. This software was produced under U.S.
+## Government contract DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is
+## operated by Los Alamos National Security, LLC for the U.S. Department of Energy. The U.S.
+## Government has rights to use, reproduce, and distribute this software.
+##
+## NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES A NY WARRANTY, EXPRESS OR
+## IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is modified to
+## produce derivative works, such modified software should be clearly marked, so as not to
+## confuse it with the version available from LANL.
+##
+## Additionally, redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+## * Redistributions of source code must retain the above copyright notice, this list of
+## conditions and the following disclaimer.
+## * Redistributions in binary form must reproduce the above copyright notice, this list of
+## conditions and the following disclaimer in the documentation and/or other materials
+## provided with the distribution.
## * Neither the name of Los Alamos National Security, LLC, Los Alamos National Laboratory, LANL
-## the U.S. Government, nor the names of its contributors may be used to endorse or promote
-## products derived from this software without specific prior written permission.
-##
-## THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS "AS IS" AND
-## ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+## the U.S. Government, nor the names of its contributors may be used to endorse or promote
+## products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS "AS IS" AND
+## ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
## OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS
## ALAMOS NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-## SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+## SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
## HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
-## EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-##
-## MATLAB code written by Jasper A. Vrugt, Center for NonLinear Studies (CNLS)
-##
-## Written by Jasper A. Vrugt: vrugt at lanl.gov
+## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
+## EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## MATLAB code written by Jasper A. Vrugt, Center for NonLinear Studies (CNLS)
+##
+## Written by Jasper A. Vrugt: vrugt at lanl.gov
################################################################################################
@@ -109,7 +109,7 @@
)
{
-
+
## dimensions
## hist.logp matrix. ndraw/nseq x nseq. length nearly ndraw.
## TODO: removed counter.fun.evals for simplicity. should have been kept?
@@ -143,14 +143,14 @@
stopifnot(func.type %in% c("calc.rmse","calc.loglik","calc.weighted.rmse","posterior.density","logposterior.density"))
stopifnot(!is.null(measurement) || func.type %in% c("posterior.density","logposterior.density"))
stopifnot(!func.type %in% c("calc.rmse","calc.loglik","calc.weighted.rmse") || "data" %in% names(measurement))
-
+
## Check INIT and FUN have required extra parameters in INIT.pars & FUN.pars
req.args.init <- names(formals(INIT))
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=" "))
-
+
req.args.FUN <- names(formals(FUN))
## if (length(req.args.FUN)<length(FUN.pars)+1) stop("Some FUN.pars are not required by FUN")
## if (length(req.args.FUN)>1){
@@ -159,7 +159,7 @@
## req.args.FUN[!req.args.FUN %in% c(names(FUN.pars))]),
## collapse=" "))
## }
-
+
## Update default settings with supplied settings
control <- modifyList(dreamDefaults(), control)
@@ -185,7 +185,7 @@
if (identical(tolower(control$outlierTest),'none')) control$burnin.length <- 0
##Choice of parallel backend
- if (control$parallel!="none"){
+ if (!control$parallel %in% c("none","snow.chains")){
parallel <- "none"
for (p in control$parallel) {
if (require(p,character.only=TRUE)) {
@@ -197,22 +197,27 @@
control$parallel <- parallel
}
-
## Check validity of settings
if (control$DEpairs==0) stop("control$DEpairs set to 0. Increase nseq?")
stopifnot(control$DEpairs<=(control$nseq-1)/2) ## Requirement of offde
- stopifnot(control$boundHandling %in% c("reflect", "bound", "fold", "none"))
+ stopifnot(control$boundHandling %in% c("reflect", "bound", "fold", "none"))
if (control$boundHandling == 'none') warning("No bound handling in use, parameters may cause errors elsewhere")
stopifnot(control$REPORT>=0)
+ if (control$parallel=="snow.chains"){
+ library(snow)
+ if (func.type != "logposterior.density") stop("control$parallel=snow.chains only supports func.type=logposterior.density")
+ if (!is.null(control$measurement)) stop("control$measurement is ignored for control$parallel=snow.chains")
+ if (!identical(FUN.pars,list())) stop("FUN.pars is ignored for control$parallel=snow.chains")
-
+ }
+
############################
## Initialize variables
NDIM <- control$ndim
NCR <- control$nCR
NSEQ <- control$nseq
-
+
## Counters
counter.fun.evals <- NSEQ
counter <- 2
@@ -224,24 +229,24 @@
max.counter.outloop <- ceiling((control$ndraw+control$steps*NSEQ)/NSEQ/control$steps)+1
if (!is.na(control$thin.t) && floor(max.counter/control$thin.t)<2) stop(sprintf("Thin parameter should be much smaller than number of iterations (currently %d,%d)",control$thin.t,max.counter))
-
+
## Calculate the parameters in the exponential power density function of Box and Tiao (1973)
cbwb <- CalcCbWb(control$gamma)
- control$Cb <- cbwb$Cb
- control$Wb <- cbwb$Wb
+ if (is.na(control$Cb)) control$Cb <- cbwb$Cb
+ if (is.na(control$Wb)) control$Wb <- cbwb$Wb
## Generate the Table with JumpRates (dependent on number of dimensions and number of pairs)
Table.JumpRate<-matrix(NA,NDIM,control$DEpairs)
for (zz in 1:control$DEpairs) Table.JumpRate[,zz] <- 2.38/sqrt(2*zz*1:NDIM)
-
+
## Initialize the array that contains the history of the log_density of each chain
hist.logp <- matrix(NA_real_,max.counter,NSEQ)
real.hist.logp <- matrix(NA_real_,max.counter,NSEQ)
-
+
if (control$pCR.Update){
## Calculate multinomial probabilities of each of the nCR CR values
pCR <- rep(1/NCR,NCR)
-
+
## Calculate the actual CR values based on p
CR <- GenCR(pCR,control)
lCR <- rep(0,NCR)
@@ -265,7 +270,7 @@
obj$control <- control
obj$in.burnin <- TRUE
-
+
obj$EXITFLAG <- NA
obj$EXITMSG <- NULL
@@ -273,7 +278,7 @@
obj$AR<-matrix(NA,max.counter,2)
obj$AR[1,2]<-NSEQ-1 ##Number if only one rejected
colnames(obj$AR) <- c("fun.evals", "AR")
-
+
##counter.fun.evals + R statistic for each variable at each step
## TODO: now using counter.report
obj$R.stat<-matrix(NA,max.counter.outloop,NDIM+1)
@@ -286,7 +291,7 @@
colnames(obj$CR) <- c("fun.evals",paste("CR",1:length(pCR),sep=""))
obj$outlier<-NULL
-
+
Sequences <- array(NA_real_, c(max.counter,NDIM+2,NSEQ))
colnames(Sequences) <- c(names(pars), "p", "logp")
## Sequences[1,] <- sapply(pars, mean) ## TODO: include?
@@ -298,28 +303,29 @@
} else Reduced.Seq <- NULL
############################
-
+
## Change MCMCPar.steps to make sure to get nice iteration numbers in first loop
control$steps<-control$steps-1
-
+
## initialize timer
tic <- as.numeric(Sys.time())
toc <- 0
counter.report <- 1
################################
-
+
## Step 1: Sample s points in the parameter space
x <- do.call(INIT,modifyList(INIT.pars,list(pars=pars,nseq=NSEQ)))
## Test that FUN returns numeric
+if (control$parallel!="snow.chains"){ ## TODO: something equivalent with snow.chains?
test.pars <- FUN.pars
test.pars[[names(formals(FUN))[1]]] <- x[1,]
modpred <- do.call(FUN,test.pars)
if (!is.numeric(modpred))
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]]))
@@ -331,7 +337,7 @@
##Save the initial population, density and log density in one list X
X <- cbind(x = x, p = tmp$p, logp = tmp$logp)
colnames(X) <- c(names(pars), "p", "logp")
-
+
##Initialise the sequences
for (qq in 1:NSEQ){
Sequences[1,,qq] <- X[qq,]
@@ -340,7 +346,7 @@
##Save pCR in memory and initialize delta.tot
obj$CR[1,] <- c(counter.fun.evals,pCR)
delta.tot <- rep(0,NCR)
-
+
##Save history log density of individual chains
hist.logp[1,] <- X[,"logp"]
real.hist.logp[1,] <- X[,"logp"]
@@ -354,7 +360,7 @@
for (gen.number in 1:control$steps) {
## Initialize DR properties
- counter.thin <- counter.thin + 1
+ counter.thin <- counter.thin + 1
## Define the current locations and associated posterior densities
x.old <- X[,1:NDIM,drop=FALSE]
@@ -402,7 +408,7 @@
## Reduced sample collection
Reduced.Seq[counter.redseq,,] <- t(X)
}
-
+
if (control$pCR.Update) {
## Calculate the standard deviation of each dimension of X
## TODO: matlab syntax is unclear - seems to be columnwise
@@ -417,7 +423,7 @@
## Update hist.logp
hist.logp[counter,] <- X[,NDIM+2]
real.hist.logp[counter,] <- X[,NDIM+2]
-
+
## Save Acceptance Rate
obj$AR[counter,] <- c(counter.fun.evals,100 * sum(accept) / NSEQ)
@@ -435,19 +441,19 @@
if (counter.outloop == 2) control$steps <- control$steps + 1
if (control$burnin.length!=0) outliers <- RemOutlierChains(X,hist.logp[1:(counter-1),],control)
-
+
if (counter.fun.evals <= end.burnin) {
- ## Check whether to update individual pCR values
+ ## Check whether to update individual pCR values
if (control$pCR.Update) {
## Update pCR values
tmp <- AdaptpCR(CR, delta.tot, lCR, control)
pCR <- tmp$pCR
lCR <- tmp$lCR
}
-
+
## Change any outlier chains to current best value of X
## TODO: matlab code didn't match paper. Outlier removal should be within burnin period
-
+
## Loop over each outlier chain (if length>0)
for (out.id in outliers$chain.id){
## Draw random other chain -- cannot be the same as current chain
@@ -471,7 +477,7 @@
}
} ##in burn in.
-
+
if (control$pCR.Update) {
## Generate CR values based on current pCR values
CR <- GenCR(pCR, control)
@@ -485,7 +491,7 @@
if (control$REPORT>0 && counter.fun.evals %% control$REPORT==0) {
counter.report <- counter.report+1
-
+
try({
obj$R.stat[counter.report,] <-
c(counter.fun.evals,
@@ -498,7 +504,7 @@
}
message(format(obj$R.stat[counter.report,], width = 10, digits = 4))
})
-
+
if (all(!is.na(obj$R.stat[counter.report,])) &&
all(obj$R.stat[counter.report,-1]<control$Rthres)) {
obj$EXITMSG <- 'Convergence criteria reached'
@@ -507,10 +513,10 @@
}
} ##counter.report
-
+
## Update the counter.outloop
counter.outloop = counter.outloop + 1
-
+
## break if maximum time exceeded
toc <- as.numeric(Sys.time()) - tic
if (toc > control$maxtime) {
@@ -547,7 +553,7 @@
obj$hist.logp <- real.hist.logp[1:(counter-1),,drop=FALSE]
obj$AR <- obj$AR[1:(counter-1),,drop=FALSE]
obj$CR <- obj$CR[1:(counter.outloop-1),,drop=FALSE]
-
+
## store number of iterations
obj$iterations <- counter.outloop
## store number of function evaluations
Modified: pkg/R/dreamCalibrate.R
===================================================================
--- pkg/R/dreamCalibrate.R 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/R/dreamCalibrate.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -6,9 +6,9 @@
stop(sprintf("Missing arguments to control: %s",
paste(req.control[!req.control %in% names(control)],collapse=", ")
))
-
+
err <- as.numeric(observed-predicted)
-
+
SSR <- sum(abs(err)^(2/(1+control$gamma)))
(-control$N*(1+control$gamma)/2)*0.5*SSR
}## calc.rmse
@@ -20,7 +20,7 @@
stop(sprintf("Missing arguments to control: %s",
paste(req.control[!req.control %in% names(control)],collapse=", ")
))
-
+
err <- as.numeric(observed-predicted)
SSR <- sum(abs(err)^(2/(1+control$gamma)))
-0.5*SSR/control$sigma^2
@@ -34,12 +34,12 @@
stop(sprintf("Missing arguments to control: %s",
paste(req.control[!req.control %in% names(control)],collapse=", ")
))
-
+
if (!all(c("Cb","Wb") %in% names(control))) control <- modifyList(control,CalcCbWb(control$gamma))
if (! "N" %in% names(control)) control$N <- length(observed)
err <- as.numeric(observed-predicted)
-
+
## Derive the log likelihood
with(control,N*log(Wb/sigma)-Cb*(sum((abs(err/sigma))^(2/(1+gamma)))))
}## calc.loglik
@@ -58,6 +58,16 @@
if (is.null(lik.control)) lik.control <- eval(formals(lik.fun)$control)
## TODO: remove ... from here and move do.call processing from dream to here.
wrap.lik.fun <- function(pars,...) lik.fun(FUN(pars,...),obs,lik.control)
+
+ ## TODO: need effective and efficient way of making lik.fun etc available inside function
+ ## wrap.lik.fun <- eval(parse(text=paste("function(id,pars,...){",
+ ## "LOClik.fun <- ",paste(deparse(lik.fun),collapse="\n"),
+ ## "LOCFUN <-",paste(deparse(FUN),collapse="\n"),
+ ## "LOCobs <-",paste(deparse(obs),collapse="\n"),
+ ## "LOClik.control <-",paste(deparse(lik.control),collapse="\n"),
+ ## "LOClik.fun(LOCFUN(id,pars,...),LOCobs,LOClik.control)}",
+ ## collapse="\n",sep="\n"
+ ## )))
dd <- dream(FUN=wrap.lik.fun,
pars=pars,
func.type="logposterior.density",
Deleted: pkg/R/maxLikPars.R
===================================================================
--- pkg/R/maxLikPars.R 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/R/maxLikPars.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,20 +0,0 @@
-##' Naive maximum density parameter selection
-##' @param ss MCMC chains. mcmc or mcmclist object
-##' @param ... extra parameters to pass to density
-##' @return named character vector of parameter values
-##'
-##' Uses which.max and density function
-maxLikPars <- function(ss,...){
- xx <- as.matrix(ss)
- pars.maxp <- rep(NA,ncol(xx))
- names(pars.maxp) <- colnames(xx)
- maxp.res <- rep(NA,ncol(xx))
- names(maxp.res) <- colnames(xx)
- for (n in colnames(xx)){
- den <- density(xx[,n],...)
- ii <- which.max(den$y)
- pars.maxp[n] <- den$x[ii]
- maxp.res[n] <- mean(diff(den$x))
- }
- return(pars.maxp)
-} ##maxLikPars
Modified: pkg/R/window.dream.R
===================================================================
--- pkg/R/window.dream.R 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/R/window.dream.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,10 +1,10 @@
window.dream <-
- function(object,
- start = 1+(end(object$Sequences)-1)*(1-fraction),
+ function(x,
+ start = 1+(end(x$Sequences)-1)*(1-fraction),
fraction = 0.5,
...)
{
- window(object$Sequences, start = start, ...)
+ window(x$Sequences, start = start, ...)
}
Modified: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/demo/00Index 2012-01-10 06:44:31 UTC (rev 34)
@@ -2,3 +2,4 @@
FME.nonlinear.model_parallelisation Example of parallelisation options
example1 n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
example2 n-dimensional Gaussian distribution from Vrugt's matlab code
+parallelisation_chain_id Example of parallelisation options passing id
Modified: pkg/demo/FME.nonlinear.model_parallelisation.R
===================================================================
--- pkg/demo/FME.nonlinear.model_parallelisation.R 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/demo/FME.nonlinear.model_parallelisation.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -39,8 +39,8 @@
print(p)
control$parallel <- p
-
- set.seed(456)
+
+ set.seed(456)
dd <- dreamCalibrate(
FUN=Model.y,
pars = pars,
@@ -65,16 +65,16 @@
## For function with 1e-3 sleep
## [1] "none"
-## p1 p2
-## 0.1511961 50.7765344
+## p1 p2
+## 0.1511961 50.7765344
## [1] 109.483
## [1] "snow"
-## p1 p2
-## 0.1511961 50.7765344
+## p1 p2
+## 0.1511961 50.7765344
## [1] 54.858
## [1] "foreach"
-## p1 p2
-## 0.1511961 50.7765344
+## p1 p2
+## 0.1511961 50.7765344
## [1] 109.593
@@ -82,7 +82,7 @@
## Test of number of function evaluations in fixed time for a time-expensive function
Model.y <- function(p,x) {
- Sys.sleep(1e-3)
+ Sys.sleep(0.1)
p[1]*x/(x+p[2])
}
@@ -98,7 +98,7 @@
set.seed(456)
control$parallel <- p
control$maxtime <- 20
-
+
dd <- dreamCalibrate(
FUN=Model.y,
pars = pars,
@@ -113,11 +113,11 @@
if (require(doSNOW)) stopCluster(cl)
-## Number of function evaluations
-## TODO: appears to be an error in foreach evaluation
+## Number of function evaluations 10/01/2012
+## Note foreach has significant overhead
## [1] "none"
-## [1] 1280
+## [1] "Number of function evaluations: 200.000000"
## [1] "snow"
-## [1] 2560
+## [1] "Number of function evaluations: 400.000000"
## [1] "foreach"
-## [1] 1280
+## [1] "Number of function evaluations: 280.000000"
Added: pkg/demo/parallelisation_chain_id.R
===================================================================
--- pkg/demo/parallelisation_chain_id.R (rev 0)
+++ pkg/demo/parallelisation_chain_id.R 2012-01-10 06:44:31 UTC (rev 34)
@@ -0,0 +1,105 @@
+## Example of parallelisation
+## Uses FME data (see other example)
+
+obs.all <- 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
+
+
+control <- list(
+ nseq=4
+ )
+pars <- list(p1=c(0,1),p2=c(0,100))
+
+
+##########################
+## Demonstration of parallelisation packages
+## Comparison of run times for known fixed termination problem
+
+lik <- function(Qobs,Qs,...){
+ res<-Qobs-Qs
+ sd.res<-sd(res)
+ var.res<-var(res)
+ logp<-sum(sapply(res,function(res_k)
+ log(1/sqrt(2*pi)*sd.res)*exp(-res_k^2/(2*var.res))))
+}
+
+Model.split <- function(id,pars){
+
+ ##Make sure that you have as many batch scripts as you have parameter chains
+ ## placeholder for this example
+
+ ##Select out this parameter set
+ this.par.set=pars[id,]
+
+ ##Do something to set the parameters in SWAT for this instance
+ ## Because they're running on other nodes, they don't have access to data or functions in R
+ x <- c( 28, 55, 83, 110, 138, 225, 375)
+ Model.y <- function(p,x) {
+ p[1]*x/(x+p[2])
+ }
+
+ ##Run the model
+ ans <- Model.y(this.par.set,x)
+
+ ##Do something to collect the answer, in appropriate form for lik.fun in dreamCalibrate
+ return(ans)
+}
+
+Model.fit <- function(id,pars){
+ ##Make sure that you have as many batch scripts as you have parameter chains
+ ## placeholder for this example
+
+ ##Select out this parameter set
+ this.par.set=pars[id,]
+
+ ##Do something to set the parameters in SWAT for this instance
+ ## Because they're running on other nodes, they don't have access to data or functions in R
+ x <- c( 28, 55, 83, 110, 138, 225, 375)
+ Model.y <- function(p,x) {
+ p[1]*x/(x+p[2])
+ }
+
+ ##Run the model
+ ans <- Model.y(this.par.set,x)
+
+ ##Do something to collect the answer as logp
+ Qs <- ans
+ Qobs <- c(0.053,0.06,0.112,0.105,0.099,0.122,0.125)
+
+ res<-Qobs-Qs
+ sd.res<-sd(res)
+ var.res<-var(res)
+ logp<-sum(sapply(res,function(res_k)
+ log(1/sqrt(2*pi)*sd.res)*exp(-res_k^2/(2*var.res))))
+ return(logp)
+}
+
+## 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
+library(doSNOW)
+cl <- makeCluster(2, type = "SOCK")
+registerDoSNOW(cl)
+
+## Using dream directly, with Model.fit, which returns a likelihood
+set.seed(456)
+control$parallel <- "snow.chains"
+dd <- dream(FUN = Model.fit,
+ pars = pars,
+ func.type = "logposterior.density",
+ control=control
+ )
+summary(dd)
+
+## TODO: dreamCalibrate cannot be used because of the way it calls FUN
+## ## Using dreamCalibrate, which allows the likelihood function to be changed separately from the model function
+## set.seed(456)
+## dd <- dreamCalibrate(
+## FUN=Model.split,
+## pars = pars,
+## obs=obs.all$y,
+## lik.fun=lik,
+## control = control
+## )
+## summary(dd)
+
+stopCluster(cl)
Property changes on: pkg/demo/parallelisation_chain_id.R
___________________________________________________________________
Added: svn:executable
+ *
Deleted: pkg/man/calc.loglik.Rd
===================================================================
--- pkg/man/calc.loglik.Rd 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/man/calc.loglik.Rd 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,3 +0,0 @@
-\name{calc.loglik}
-\title{Log likelihood function}
-\description{Calculate log likelihood of predicted values matching observed}
\ No newline at end of file
Deleted: pkg/man/calc.rmse.Rd
===================================================================
--- pkg/man/calc.rmse.Rd 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/man/calc.rmse.Rd 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,4 +0,0 @@
-\name{calc.rmse}
-\title{RMSE log likelihood function}
-\description{Calculate log RMSE using formulas used by Vrugt}
-\details{This is used as the default, but it is recommended to use an objective function better suited to the particular problem}
\ No newline at end of file
Deleted: pkg/man/calc.weighted.rmse.Rd
===================================================================
--- pkg/man/calc.weighted.rmse.Rd 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/man/calc.weighted.rmse.Rd 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,2 +0,0 @@
-\name{calc.weighted.rmse}
-\title{Sigma-weighted RMSE log likelihood function}
\ No newline at end of file
Modified: pkg/man/coef.dream.Rd
===================================================================
--- pkg/man/coef.dream.Rd 2010-12-16 05:53:07 UTC (rev 33)
+++ pkg/man/coef.dream.Rd 2012-01-10 06:44:31 UTC (rev 34)
@@ -1,24 +1,29 @@
\name{coef.dream}
\alias{coef.dream}
\alias{coef.dream_model}
+\alias{maxLikCoda}
\title{Extract parameter values from dream or dream_model object}
\usage{
-\method{coef}{dream}(object, method = c("sample.ml","uni.mode", "mean", "median"), \dots)
+\method{coef}{dream}(object, method = c("sample.ml","uni.mode", "mean",
+"median"), \dots)
+maxLikCoda(x)
}
\description{Extract parameter values using a choice of methods (or an
arbitrary function)}
\value{named vector of parameter values}
\arguments{
+ \item{x}{for maxLikCoda, usually mcmc or mcmc.list object, or any representation of
+ MCMC chains that can be converted with \code{as.matrix}}
\item{object}{dream object}
\item{method}{method for extracting a parameter set from the MCMC
chains. One of:
\describe{
- \item{\code{"uni.mode"}}{
+ \item{\code{"uni.mode"}}{ using maxLikCoda,
mode of the univariate density estimate
for each parameter, using settings as in
\code{\link{densityplot.mcmc}}.
May not find the optimal parameter combination of the
- distribution is multi-modal.
+ distribution is multi-modal.
}
\item{\code{"mean"}}{
mean of each univariate parameter distribution.
@@ -38,7 +43,8 @@
\item{...}{Passed to \code{\link{window.dream}}}
}
\details{
-
+ maxLikCoda re-uses code from \code{\link{densityplot.mcmc}}
+
e.g. of using arbitrary function for method: 20\% quantile.
\code{
coef(object,method=function(sss) apply(as.matrix(sss),2,quantile,0.2)
Added: pkg/man/compareToMatlab.Rd
===================================================================
--- pkg/man/compareToMatlab.Rd (rev 0)
+++ pkg/man/compareToMatlab.Rd 2012-01-10 06:44:31 UTC (rev 34)
@@ -0,0 +1,29 @@
+\name{compareToMatlab}
+\alias{compareToMatlab}
+\title{
+Compare dream and matlab results
+}
+\description{
+ Compare results from dream run to results from matlab
+}
+\usage{
+compareToMatlab(mat.file, dream.obj)
+}
+\arguments{
+ \item{mat.file}{
+ Path to matlab .mat file
+}
+ \item{dream.obj}{
+ object returned by dream
+}
+}
+\details{
+ Used in demos
+}
+\value{
+ Prints diagnostics. \cr
+ Did matlab and R use same inputs? \cr
+ Do outputs have same dimensions? \cr
+ ks.test that samples are from same distribution for each variable \cr
+ returns result of \code{\link{plotMCMCQQ}}
+}
Property changes on: pkg/man/compareToMatlab.Rd
___________________________________________________________________
Added: svn:executable
+ *
Modified: pkg/man/dream.Rd
===================================================================
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/dream -r 34
More information about the Dream-commits
mailing list