Nothing
R/twDEMCBlocks.R
In twDEMC: parallel DEMC
twDEMCBlockInt <- function(
### Parallelized differential Evolution Markov Chain with blocked parameter update
pops = list( list( ##<< list of population infos for each population, each info is a list with components
##describe<<
parms=list() ##<< list of matrices (nState x nParm x nChain) initial states for each population see details and \code{\link{initZtwDEMCNormal}}.
, nGen = 12 ##<< number of generations, i.e. steps to proceed
, TProp=NULL ##<< numeric vector (nResComp) temperature proportions of result components.
## It can also be given as character vector with names of result components, however, be aware that htis fails if several logDen may return components of the same name
, X=NULL ##<< numeric matrix (nParm x nChainPop) initial state
, logDenCompX=NULL ##<< numeric matrix (nComp x nChain): logDen components of initial state X, see details
, spaceInd = 1 ##<< the space replicate that this population belongs to
, upperParBounds = numeric(0) ##<< named numeric vectors: giving upper parameter bounds: lowerBound < par <= upperBound.
## For exploring subspaces of the limiting distribution, see details
, lowerParBounds = numeric(0) ##<< similar to upperParBounds: sample > bound
, splits=numeric(0) ##<< named numeric vector of splitting points, used to remerge divided subspaces
))
##end<<
,dInfos = list( list( ##<< named list of used density functions. Each entry is a list with components
##describe<<
fLogDen=NULL ##<< \code{function(theta, ...)} calculates a vector of logDensities
## corresponding to different data streams of parameter vector theta
## , or \code{function(theta, logDenCompX, metropolisStepTemp, ...)}
, argsFLogDen=list() ##<< further arguments passed to fLogDen
, compPosDen #=dInfos[[1]]$compPosDen ##<< names or index of the parameters used in with this density, defaults to all parameters
, intResComp=vector("integer",0) ##<< integer or character vector: indices or names of results components of fLogDen
## that are used for internal Metropolis decisions
#, TFix = vector("numeric",0) ##<< named numeric vector with Temperature for result components that have fixed temperature
, maxLogDen=c() ##<< integer vector (nRespComp): maximum logDensity (=usually -1/2* number of observations) for each result component, used in control of overfitting
, fDiscrProp=NULL ##<< function applied to proposal, e.g. to round proposals to to discrete possible values function(theta,...)
, argsFDiscrProp=list() ##<< further arguments to fDiscrProp
))
##end<<
, blocks = list( list( ##<< list of parameter blocks, each block is a list with entries
##describe<<
compPos = dInfos[[1]]$compPosDen ##<< names or index of the parameter components to be updated
, fUpdateBlock = NULL ##<< function to update the parameters.
## It must return a list with first named components \code{accepted}, \code{xC}, and optionally \code{intermediate}
## as described in \code{\link{updateBlockTwDEMC}}. If not specified (default), then \code{\link{updateBlockTwDEMC}} is used and \code{dInfoPos} must be specified
, argsFUpdate=list() ##<< further arguments passed to fUpdate
, intermediateId=character(0) ##<< string: identifier in list of intermediate results that are shared between block udpate functions. See details.
, dInfoPos=1 ##<< name or position to \code{fLogDenInfo} if using Metropolis update.
## Several blocks may share the same density but update different parameters
))
##end<<
, stateCurrent = list(
,xC=c() ##<< numeric vector: current accepted state
,logDenCompC=c() ##<< numeric vector: result components of all fLogDens for current position
,parUpdateDenC=c() ##<< integer vector: logDensity that recently updated parameter at given index
## to handle first steps in Multi-step Metropolis decision internally.
## See details.
)
, TSpec = matrix(1, ncol=2, nrow=1, dimnames=list(NULL, c("T0","TEnd")) )
### numeric matrix (nResComp x 2): specifying Initial and End Temperature of each fLogDen result component.
### If only one row is specified, the Temperature is taken for all result components
, m0 = numeric(0) ##<< scalar integer: number of samples in initial population, if length==0 will be calculated from number of chains and number of parameters
, nGen = integer(0) ##<< scalar integer: number of generations, if given overwrites \code{pops[[i]]$nGen}
, spacePop = 1:nPop ##<< the space replicate that each population belongs to. By default assume only one population per space, overridden by entry in pops
, controlTwDEMC = list() ##<< control parameters influencing the update and the convergens of the chains (see details)
, debugSequential=FALSE ##<< if TRUE apply is used instead of sfApply, for easier debugging
, remoteDumpfileBasename=NULL ##<< the basename of a dumpfile that is created on error on remote process (see example section)
, doRecordProposals=FALSE ##<< if TRUE then an array of each proposal together with the results of fLogDen are recorded and returned in component Y
, progressOutput="." ##<< output after each thinning interval
){
##details<<
## This is the central method for applying a Differential Evolution Markov Chain given a function of
## an unnormalized probability density.
## It is invoked usually by (\code{\link{twDEMCBlock.array}} or \code{\link{twDEMCBlock.twDEMCPops}})
##seealso<<
## \code{\link{calcDEMCTemp}}
## \code{\link{logDenGaussian}}
# # \code{\link{.generateXPropThin}}
# # \code{\link{.doDEMCSteps}}
# # \code{\link{.doDEMCStep}}
## #details<<
## # further elements in pop that are calculated are \itemize{
## # \item resCompPos: position of the result components in the concatenation across all density functions
## # }
##details<< \describe{ \item{recognized control parameters in \code{controlTwDEMC}: }{
##describe<<
ctrl = list(
thin = 4 ##<< thinning interval, Default: 4
,F = 2.38 ##<< related to multiplicative error (F2=F/sqrt(2*Npar), see eps.mult, Default: 2.38
,pSnooker= 0.1 ##<< probability of a snooker update (others parallel updates), Default 0.1
,pGamma1 = 0.1 ##<< probability of jumping to state of another chain (different modes), Default 0.1
,epsMult =0.2 ##<< >0 gives d-dimensional variation around gamma.
## It adds scaled uncorrelated noise to the proposal. Default: 0.2 ,
## Its advantage over eps.add is
## that its effect scales with the differences of vectors in the population whereas eps.add does not.
## if the variance of a dimensions is close to 0, eps.mult gives smaller changes. ,
## A uniformly distributed error, i.e. F2*runif(1+-epsMult*prop) multiplied to difference vector
## from parallel update
,epsAdd = 0 ##<< >0 is needed to ensure that all positions in the space can be reached.
## For targets without gaps, it can set small or even to 0. Default 0 ,
## sd of normally distributed error added to proposal by parallel or snooker update.
,pAcceptWindowWidth = 8 ##<< number of thinning intervals back over which the acceptance rate is calculated Default 8
,probUpDir=0.5 ##<< probability of direction between two states of increasing Density, Default 0.5
## , Increasing this during burin may accelerate convergence
,initialAcceptanceRate=0.25 ##<< numeric matrix (nBlock x nChains) initially assumed acceptance rate. Default 0.25
## , Used to calculate the number of generations backwards to sample from
,DRgamma=0 ##<< factor for reducing step length [0..1) in delayed rejection step, 0 means no DR step, Default 0
,minPCompAcceptTempDecr=0.15 ##<< if acceptance rate drops below minPCompAcceptTempDecr+0.02 times this level, employ delayed rejection (DR), Default 0.15
,pIndStep = 1.5 ##<< assume state to be independent, after on average about those number of accepted steps, Default 1.5
,nPastGen = 10 ##<< factor for determining the number of recent past states to sample during burnin. Default 10
## It is multiplied by the number of parameters. Past generations are calculated
## by deviding by the number of chains per population
#,useLoadBalancing=FALSE ##<< if set to TRUE, submits each step separaetely to a free node, else each node gets an entire chain process
#,freeMasterNode=FALSE ##<< if set to TRUE, no job is submitted to first node, so that this node can dispatch jobs without waiting,, see \code{\link[twSnowfall]{sfFArgsApplyDep}}
,returnIntermediates=TRUE ##<< set to FALSE if intermediate result is large and transferring it between slaves slows down calculation
#,useMultiT = FALSE ##<< whether to downscale Temperature of result components during burnin
#moved to block,TFix = vector("numeric",0) ##<< named numeric vector: result components for which temperature shoudl not change
,useConditionalProposal = FALSE ##<< if set to TRUE, poposal steps are generated conditional on the state of the other blocks (experimental)
#,controlOverfittingMinNObs = Inf ##<< scalar positive integer: minimum number of observations where bias is controlled. Need at least about 20 observations to work. Default Inf says no bias correction.
)
##details<< }}
ctrl[names(controlTwDEMC)] <- controlTwDEMC
#
#-- dimensions of pops
nPop <- length(pops)
iPops <- 1:nPop
pops <- if( 0 == length(nGen) ){
lapply( iPops, function(iPop){ .completePopDescription( pops[[iPop]], spaceInd=spacePop[iPop] )}) # fill in default values for missing entries
}else{
if( length(nGen)==1) nGen=rep(nGen[1],nPop)
if( length(nGen) != nPop ) stop("twDEMCBlockInt: lenght of nGen must equal the number of populations.")
lapply( iPops, function(iPop){ .completePopDescription( pops[[iPop]], nGen=nGen[iPop], spaceInd=spacePop[iPop])}) # fill in default values for missing entries
}
ZinitPops <- lapply(pops,"[[","parms")
parNames <- colnames(ZinitPops[[1]])
nParm <- ncol(ZinitPops[[1]])
nChainPop <- dim(ZinitPops[[1]])[3]
if( nChainPop < 2) stop(paste("Initialized with only ",nChainPop," per population, but must have at least two chains per population (recommened is at least 4)"))
nChain <- nChainPop*nPop
popChain <- rep(1:nPop,each=nChainPop) # population for chain at given index
chainsPop <- lapply( iPops, function(iPop){ (iPop-1)*nChainPop+(1:nChainPop)}) # chains for given population
#
# check that have enough values for initial populations of each population
if( 0 == length(m0) ) m0 = calcM0twDEMC(nParm,nChainPop) # number of samples in initial population
M0Pops <- sapply( ZinitPops, nrow )
if( any(M0Pops < 0.8*m0) ) warning(paste("twDEMCBlockInt: too few initial cases for populations",paste(which(M0Pops < m0),collapse=",")) )
#
# determine number of thinning intervals and create result arrays
nGenPops <- sapply( pops, "[[", "nGen")
nThinnedGenPops = sapply( nGenPops %/% ctrl$thin, max, 0 ) #number of thinning intervals
nGenPops = nThinnedGenPops * ctrl$thin #number of generations in this batch (do not need to move behind last thinning interval)
Nz <- M0Pops +(nThinnedGenPops) #number of rows in Z after batch
XPops <- lapply( pops, "[[", "X" )
XChains <- do.call(cbind,XPops)
logDenCompXPops <- lapply( pops, "[[", "logDenCompX" ) # see initial states for initializing missing entries
upperParBoundsPop <- lapply( pops, "[[", "upperParBounds" )
lowerParBoundsPop <- lapply( pops, "[[", "lowerParBounds" )
#
#-- dimensions of dInfo
nDen <- length(dInfos)
iDens <- 1:nDen
if( is.null(names(dInfos))) names(dInfos) <- paste("dInfo",iDens,sep="") # make up names if not given
posInfosEmptyString <- which(names(dInfos) == "")
names(dInfos)[posInfosEmptyString] <- paste("dInfo",posInfosEmptyString,sep="")
#compPosDens <- lapply(dInfos, "[[", "compPosDen")
#calculate each logDen once to obtain result component names required for proper dInfos initialization
#dInfo <- dInfos[[1]]
resCompNamesList <- lapply(iDens, function(iInfo){
dInfo <- dInfos[[iInfo]]
.resFLogDen <- .calcFLogDen(dInfo$fLogDen, xP=XChains[,1], compInt=c(), TiInt=c()
, argsFLogDen=dInfo$argsFLogDen, maxLogDen=c() )
resCompNames <- .getResFLogDenNames( .resFLogDen, logDenName=names(dInfos)[iInfo] )
})
for( i in seq_along(dInfos)){
dInfos[[i]] <- .completeDInfo( dInfos[[i]], parNames=parNames, resCompNames = resCompNamesList[[i]] )
}
#
#-- calculate positions in resCompNames of results of different blocks
resCompNamesU <- lapply( dInfos, "[[", "resCompNames" )
resCompNamesUFlat <- as.vector(do.call(c, resCompNamesU)) # concatenate result names
nResComp <- length(resCompNamesUFlat)
resCompPosPops <- { # list: for each dInfo: position of resCompNames in resCompNamesUFlat
tmp.length <- sapply(resCompNamesU, length)
tmp.end <- cumsum(tmp.length)
tmp.i <- vector("list",nDen)
for( iDen in iDens ){
tmp.i[[iDen]] <- dInfos[[iDen]]$resCompPos <- tmp.end[iDen]+1-(tmp.length[iDen]:1)
}
tmp.i
}
#
#-- dimensions of blocks
#mtrace(.checkBlock)
#recover()
blocks <- lapply( blocks, .completeBlockInfo, dInfos=dInfos, parNames=parNames)
nBlock <- length(blocks)
iBlocks <- 1:length(blocks)
#block <- blocks[[1]]
dInfoPosBlock <- sapply(blocks, "[[", "dInfoPos")
# following refers to position in parNames
compPosBlock <- lapply( blocks, "[[", "compPos" ) # already transformed to positions in .checkBlock
# following refers to position in dInfo$compPosDen
compPosInDenBlock <- lapply( blocks, "[[", "compPosInDen" )
# make sure that all parameters are updated at least once
if( 0==length(compPosBlock) || (length(sBlock <- sort(unique(do.call(c, compPosBlock)))) != nParm) || !all.equal( sBlock, 1:nParm, check.attributes = FALSE) )
stop(paste("each parameter (columns of parms) must be updated in at least one block. nParm=",nParm,"updatedVars=",paste(sBlock,collapse=",")) )
#
#-- initialize further parameters to parallel and snooker steps
# that depend on pops (nParm, nChainPop)
ctrl$Npar12 =(nParm - 1)/2 # factor for Metropolis ratio DE Snooker update
ctrl$F2 = ctrl$F/sqrt(2*nParm)
ctrl$F1 = 1
# terBraak report that may not sample the distribution, if not using the full past
# but together with decreasing temperature acceptance rate drops very low if using the full past
# hence constrain to recent past during burnin
ctrl$gInd <- ctrl$nPastGen*nParm/nChainPop #number of independent generations to sample from, similarly to M0 (usually ctrl$nPastGen=10)
#
#-- evaluate logDen of initial states and get names of the result components (dInfo$resCompNames)
##details<< \describe{ \item{Initial state: \code{X}}{
## If initial state X is not specified, the last column (generation) of Z is utilized.
## If in addition to X, logDenX is specified as a numeric vector, the fLogDen will not be avaluated for the initial state of the chains.
## All the results of fLogDen for the initial state must be finite.
## }}
# find those chains where no logDenCompX is given, or not all components are finite
popsMissingResX <- ( 0 == sapply( logDenCompXPops, length) )
iMissingPops <- which(popsMissingResX)
iNonMissingPops <- which(!popsMissingResX)
missingPopInfo <- do.call( rbind, lapply( iMissingPops, function(iPop){
cbind(iPop=iPop,iChainInPop=1:nChainPop)
} ))
missingPopInfo <- rbind( missingPopInfo, do.call( rbind, lapply( iNonMissingPops, function(iPop){
#all( is.finite(logDenCompXPops[[iPop]][,1]))
#chainIsAllFinite <- apply( logDenCompXPops[[iPop]], 2, all, is.finite )
chainIsAllFinite <- apply( logDenCompXPops[[iPop]], 2, function(cChain){ all(is.finite(cChain))} )
iChainInPop=which(!chainIsAllFinite)
if( 0 != length(iChainInPop) )
cbind(iPop=iPop,iChainInPop=iChainInPop)
else
cbind(iPop=iPop, iChainInPop=NA_integer_ )[FALSE,]
})))
if( 0 < nrow(missingPopInfo) ){
# evaluate logDen for those
missingPopInfo <- cbind(missingPopInfo, iChain=(missingPopInfo[,"iPop"]-1)*nChainPop + missingPopInfo[,"iChainInPop"])
XChainsMissing <- XChains[,missingPopInfo[,"iChain"] ,drop=FALSE] # initial states of the missings
logDenCompDen <- vector("list", nDen )
#XXTODO calculate using intermediate
#XXTODO on restart provide existing resLogDen
logDenRes <- twCalcLogDensPar(
### Invokes all fLogDens with proposal in a parallel way, taking care of intermediates between densities
dInfos ##<< list describing the logDensities, see \code{\link{twDEMCBlockInt}}
,xProp=t(XChainsMissing[,,drop=FALSE])
, debugSequential=debugSequential
, remoteDumpfileBasename=remoteDumpfileBasename
)
if( any(!is.finite(logDenRes$logDen)) )
stop(paste("twDEMCBlockInt: non-finite logDensity of starting value ",sep=""))
#iDen <- 1
for( iDen in seq_along(dInfos) ){
#calculated further down, for all dInfos, not just the missing ones: dInfos[[iDen]]$resCompPos <-
.resCompPosI <- which(logDenRes$logDenCompPos==iDen)
logDenCompDen[[iDen]] <- .logDenCompI <- logDenRes$logDenComp[,.resCompPosI]
}
logDenCompXMiss <- logDenRes$logDenComp
# replace statesin initial history for proposal generation with infinite logDen by
# ones with finite states from other chains/populations
#jPop <- 1
for( jPop in seq_along(iMissingPops) ){
iiPop <- iMissingPops[jPop]
iInfo <- which( missingPopInfo[,"iPop"] == iiPop )
logDenCompXPops[[ iMissingPops[jPop] ]] <-
t(logDenCompXMiss[iInfo,])
}
for( jPop in seq_along(iNonMissingPops) ){
iiPop <- iNonMissingPops[jPop]
iInfo <- which( missingPopInfo[,"iPop"] == iiPop )
logDenCompXPops[[iiPop]][,missingPopInfo$iChainInPop[iInfo] ] <-
t(logDenCompXMiss[iInfo,])
}
}
#-- do overfitting control on the intial logDenComp
maxLogDen <- structure( numeric(length(resCompNamesUFlat)), names=resCompNamesUFlat)
#dInfo <- dInfos[[1]]
for( dInfo in dInfos){
maxLogDen[dInfo$resCompNames] <- dInfo$maxLogDen
}
#logDenCompXPop <- logDenCompXPops[[2]]
tmp <- logDenCompXPops
logDenCompXPops <- lapply(logDenCompXPops, function(logDenCompXPop){
logDenCompXPop[] <- apply( logDenCompXPop, 2, function(logDenComp){ pmin( maxLogDen, logDenComp) })
logDenCompXPop # keep structure
})
#
logDenCompXChains <- do.call( cbind, logDenCompXPops )
# update the names in intial conditons to hold unique names
rownames(logDenCompXChains) <- resCompNamesUFlat
#
#-- get indices of internal components and TFix (dInfo$ intResCompPosWithin, intResCompPos, posTFix)
# depends on proper resCompNames (inferred from invoking fLogDen with initial states)
intResCompDen <- vector("list",nDen)
#.TFixL <- vector("list", nDen )
#.posTFixL <- vector("list", nDen )
for( iDen in iDens ){
dInfo = dInfos[[iDen]]
#.TFixL[[iDen]] <- dInfo$TFix
#.posTFixL[[iDen]] <- dInfos[[iDen]]$posTFix <- dInfo$resCompPos[ match( names(dInfo$TFix), dInfo$resCompNames ) ] # positions in concatenated result components
intResCompI <- dInfo$intResComp
iPosW <- if( is.character(intResCompI) ){
iPosW <- match(intResCompI, dInfo$resCompNames )
if( any(is.na(iPosW)))
stop(paste("not all names of intResCompNames (",paste(intResCompI,collapse=","),") in return of fLogDen: ",paste(dInfo$resCompNames,collapse=",")))
iPosW
} else intResCompI
dInfos[[iDen]]$intResCompPosWithin <- iPosW
intResCompDen[[iDen]] <- dInfos[[iDen]]$intResCompPos <- dInfo$resCompPos[iPosW] # position in concatenation of all result components
} #iDen
#TFixAll <- do.call(c, .TFixL )
#posTFixAll <- do.call(c, .posTFixL ) # position in concatenation of all resComp
#-- Temperature of result components, i.e. data streams
if( nrow(TSpec) == 1 ){
TSpec <- do.call( rbind, lapply(1:nResComp,function(i){TSpec[1,]}) )
rownames(TSpec) <- resCompNamesUFlat
}
if( nrow(TSpec) != nResComp ){
.missingResCompNames <- unique(resCompNamesUFlat[ !(resCompNamesUFlat %in% rownames(TSpec)) ])
if( length(.missingResCompNames) ){
stop(paste("twDEMCBlockInt: TSpec missing for result components",paste(.missingResCompNames,collapse=",")))
}
TSpec <- TSpec[ resCompNamesUFlat, ]
}
tempResCompPops <- vector("list", nPop )
for( iPop in iPops){
.tempResL <- lapply( 1:nResComp, function(iResComp){
.tempG <- pmax(1,c( TSpec[iResComp,"T0"],
calcDEMCTemp( TSpec[iResComp,"T0"], TSpec[iResComp,"TEnd"], nGenPops[iPop] )))
})
tempResCompPops[[iPop]] <- do.call( cbind, .tempResL)
colnames(tempResCompPops[[iPop]]) <- resCompNamesUFlat
}
#-- setup acceptance rate recording
# number of requried independent generations to choose from (10*d independent states: TerBraak06 TerBraak08) devided by number of chains
if( length(ctrl$initialAcceptanceRate) == 1)
ctrl$initialAcceptanceRate <- matrix(ctrl$initialAcceptanceRate, nrow=nBlock, ncol=nChain, dimnames=list(blocks=NULL, chains=NULL) )
if( !is.matrix(ctrl$initialAcceptanceRate) || dim(ctrl$initialAcceptanceRate) != c(nBlock,nChain))
stop("ctrl$initialAcceptanceRate must be matrix of dim nBlock x nChain")
pAcceptPops <- popMeansTwDEMC(ctrl$initialAcceptanceRate, nPop) # current acceptance rate of population (block x pop)
pAcceptPops[] <- pmax(1/ctrl$thin,pAcceptPops ) # lower bound
#nGenBack <- pmin(mZ,ceiling(ctrl$rIndToAcc * ctrl$gInd * pmax(1,ctrl$pIndStep/(ar * ctrl$thin)))) #number of genrations to select states from for each population, ctrl$gInd is multiplied by number of rows for one step depending on acceptance rate and thinning but at least one
pAcceptMinPops <- apply(pAcceptPops,2,min) # minimum of acceptance rate across blocks for each population
nGenBackPops <- pmin(M0Pops,ceiling(ctrl$gInd * pmax(1,ctrl$pIndStep/(pAcceptMinPops * ctrl$thin)))) #number of genrations to select states from for each population, ctrl$gInd is multiplied by number of rows for one step depending on acceptance rate and thinning but at least one
##details<< \describe{ \item{Acceptance rate}{
## The acceptance rate is tracked for each chain across ctrl$pAcceptWindowWidth thinning intervals. ,
## If acceptance rates drops to low values, this might be because of bad localization
## ,i.e the requirement to use states from more distant past.
## In order to improve localization, less parameters or more chains per population are required.
## }}
# acceptWindow holds the number of accepted steps for each thinning interval for each chain
# if its boundaries are filled up, the last nGenBack states are copied to the first interval
aThinWinW <- ctrl$pAcceptWindowWidth #%/% ctrl$thin # number of thinning intervals over which to calculate average acceptance rate
#ctrl$pAcceptWindowWidth <- aThinWinW * ctrl$thin
# acceptWindow records number of accepted steps in thinnging interval
# (stepslots x blocks x chains)
#acceptWindow <- matrix( rep(ctrl$thin*ar, each=2*aThinWinW*nChainPop), nrow=2*aThinWinW, ncol=d$chains ) #record number of accepted steps in thinning interval
# (width x block x chain )
acceptWindow <- array( rep(ctrl$thin*ctrl$initialAcceptanceRate, each=2*aThinWinW)
, dim=c(2*aThinWinW, dim(ctrl$initialAcceptanceRate) )
, dimnames=c(list(steps=NULL),dimnames(ctrl$initialAcceptanceRate)) ) #record number of accepted steps in thinning interval
#-- preallocate output of parameters, calculated LogDen, indices of accepted steps, and next proposal
# need lists for populations because Nz (number of samples including intial population ) and nGen (samlples without initial) and nThinnedGen is differ across populations in general
ZPops <- lapply( iPops, function(iPop){
dNames <- if( !is.null(dimnames(ZinitPops[[iPop]])) ) dimnames(ZinitPops[[iPop]]) else list( samples=NULL, parms=paste("p",1:nParm,sep="_"), chains=NULL )
Z <- array(NA_real_, dim=c(Nz[iPop],nParm,nChainPop), dimnames=dNames)
Z[1:M0Pops[iPop],,] <- ZinitPops[[iPop]]
Z
})
# several resLogDen per chain (dInfos are concatenated)
resLogDen <- lapply(iPops,function(iPop){ res <- array( NA_real_
, dim=c(1+nThinnedGenPops[iPop], nResComp, nChainPop)
, dimnames=list(steps=NULL, resComp=resCompNamesUFlat, chains=NULL))
res[1,,] <- logDenCompXPops[[iPop]]
res
})
# one rLogDen per chain and per dInfo
logDen <- lapply(iPops, function(iPop){ res <- array(NA_real_
, dim=c(1+nThinnedGenPops[iPop], nDen, nChainPop)
, dimnames=list(steps=NULL, den=names(dInfos), chains=NULL) )
for( iChainPop in 1:nChainPop)
for( iDen in iDens )
res[ 1,iDen,iChainPop] <- sum(logDenCompXPops[[iPop]][resCompPosPops[[iDen]],iChainPop ])
res
})
# one pAcceptance indication per chain and per block
pAccept <- lapply(iPops, function(iPop){ res <- array(NA_real_
, dim=c(1+nThinnedGenPops[iPop], nBlock, nChainPop)
, dimnames=list(steps=NULL, block=names(blocks), chains=NULL) )
})
for( iPop in iPops ){
pAccept[[iPop]][1,,] <- ctrl$initialAcceptanceRate[,chainsPop[[iPop]] ]
}
# record of proposals and fLogDen results, rows c(accepted, parNames, resCompNames, rLogDen)
# if doRecordProposals is FALSE record only the thinning intervals for the last 128 generations
# +1 for the initial state
nThinLastPops <- if(doRecordProposals) nThinnedGenPops else pmin(nThinnedGenPops, ceiling(128/ctrl$thin)) #number of thinning steps to record
nThinOmitRecordPops = nThinnedGenPops-nThinLastPops #the Thinning intervals with no recording of outputs
nGenOmitRecordPops = nThinOmitRecordPops*ctrl$thin
nGenYPops <- nThinLastPops*ctrl$thin #+1?
YPops <- lapply(iPops, function(iPop){ res <- array( NA_real_
, dim=c( nGenYPops[iPop], nParm+nBlock+nResComp, nChainPop)
, dimnames=list(steps=NULL, vars=c(rownames(XPops[[iPop]]),paste("accepted",iBlocks,sep=""),colnames(resLogDen[[iPop]])), chains=NULL ) )
# already initialized first row in nGen cycle
#if( nGenYPops[iPop] == nGenPops[iPop] ) # record first state as accepted
# res[1,,] <- rbind( XPops[[iPop]], matrix(TRUE, nrow=nBlock, ncol=nChainPop), logDenCompXPops[[iPop]] )
res
}) # records proposals, indication of acceptance per Block for each Chain
# here accepted is 0: not accepted, 1: full step accepted, between DR step accepted at DRGamma
#
##details<< \describe{ \item{Overfitting control}{
## Cost, i.e. -2*logDensity, below the number of observations indicates an overfitting to the data stream.
##
## By supplying a maximum logDensity for one result component with \code{dInfo$maxLogDen = -1/2*nObs} one can
## activate the control for this overfitting. This works only, however, if there are enough observations.
## When a result component is based on less than about 30 observartions, one should give a maximum Densitiy of 0.
##
## With overfitting control, there is an upper bound bound for the logDensity of one result compoent-
## Calculated higher logDensities indicate overfitting and are decreases to the maxLogDen value.
##}}
#
#
#-- arguments to fUpdate argument argsFUpdate that do not change within thinning interval
argsFUpdateDefault = list(
ctrl=ctrl
#,fCalcComponentTemp=calcComponentTemp
#,fDiscrProp=fDiscrProp
#,argsFDiscrProp=argsFDiscrProp
)
# construct an list for each block that includes compPos, and intResCompPos
#block <- blocks[[1]]
argsFUpdateBlocks <- lapply( blocks, function(block){
dInfo <- if( block$useMetropolisUpdate) dInfos[[block$dInfoPos]] else NULL
c( argsFUpdateDefault, dInfo=list(dInfo), block
#, list(maxLogDen=maxLogDenL[[block$dInfoPos]]) # defined in dInfo on .checkDInfoResCompNames
#, list(TPropPops=TPropPops[dInfo$compPosDen,])
)
})
remoteArgsFUpdateBlocksTwDEMC <- list(
remoteFun = .updateIntervalTwDEMCChain # this function is called remotely and invokes the blockupdate for all steps of a chain within a thinning interval.
,remoteDumpfileBasename = remoteDumpfileBasename #if null delete
,argsUpdateBlocksTwDEMC=list(
argsFUpdateBlocks = argsFUpdateBlocks
,nBlock=nBlock #, iBlocks=iBlocks
,upperParBoundsPop=upperParBoundsPop
,lowerParBoundsPop=lowerParBoundsPop
,returnIntermediate=TRUE # here stays on slave, so no overhead of returning it between block updates
#,temp=temp
)
)
#
if( !debugSequential && sfParallel() ){
sfExport("remoteArgsFUpdateBlocksTwDEMC") #evaluated in remote process, only passed one time
sfExport(".updateIntervalTwDEMCChain",".updateBlocksTwDEMC",".calcFLogDen", namespace="twDEMC") # used in remote process
suppressWarnings(sfLibrary(abind))
}
#else( assign("remoteArgsFUpdateBlocksTwDEMC", remoteArgsFUpdateBlocksTwDEMC, pos=1)) #export to current
# pass remoteArgsFUpdateBlocksTwDEMC as a name instead of each time complete object in a call
# by using sfRemoteWrapper it is obtained from exported object on the remote node
nameRemoteArgs <- if( !debugSequential && sfParallel() ) as.name("remoteArgsFUpdateBlocksTwDEMC") else remoteArgsFUpdateBlocksTwDEMC # eval.parent fails for remoteArgsFUpdateBlocksTwDEMC within sequential function
#
# arguments to demc that change between thinning intervals
# substract logDen components of logDenCompX from logDenXExt
chainState <- list(
X = XChains
,logDenCompX = logDenCompXChains
,intermediatesChains=lapply( 1:nChain, function(iChain){list()} ) # initially empty
,parUpdateDen = array(TRUE, dim=c(nDen,nParm,nChain), dimnames=list(dens=NULL,parms=parNames,chains=NULL)) ##<< for each parameter/density combination: is the density up to date
#,logDenX = logDenX
)
parUpdateDenLast = chainState$parUpdateDen # to store the last parUpdateDen
#
maxNThinned <- max(nThinnedGenPops)
isSamePopLength = all( nThinnedGenPops == maxNThinned)
isPops <- 1:nPop # populations that take part in this step (some may need less steps)
isChains <- 1:nChain # chains that take part in this step
#iThin0 <- maxNThinned-1
#
#--- main cycle of updates (each ctrl$thin generations)
for( iThin0 in (0:(maxNThinned-1)) ){
if( length(progressOutput)!=0 && nchar(progressOutput)!=0 ){ cat(progressOutput) }
iGen = iThin0*ctrl$thin+(1:ctrl$thin) # the generations within this thinning step
mZPops <- M0Pops + iThin0
iPopsOut <- which(nThinnedGenPops == iThin0) # those pops drop out
#if( length(iPopsOut)) recover()
if( 0 != length(iPopsOut) ){
iChainsOut <- do.call( c, lapply( iPopsOut, function(iPop){ chainsPop[[iPop]] }) )
iiPopsOut <- which(isPops %in% iPopsOut ) # index in set of current populations
iiChainsOut <- rep((iiPopsOut-1)*nChainPop, each=nChainPop) + (1:nChainPop)
isPops <- isPops[-iiPopsOut]
isChains <- isChains[-iiChainsOut]
# recored the last parUpdateDen for dropouts
parUpdateDenLast[,,iChainsOut] <- chainState$parUpdateDen[,,iiChainsOut]
# adapt chainState that refer to current populations only
chainState <- within(chainState,{
X <- X[,-iiChainsOut ,drop=FALSE]
intermediatesChains <- intermediatesChains[-iiChainsOut]
logDenCompX <- logDenCompX[,-iiChainsOut ,drop=FALSE]
parUpdateDen <- parUpdateDen[,,-iiChainsOut ,drop=FALSE]
})
#acceptance rates and nGenBack refer to all chains
}
#
# sample random vectors (note here parameters in rows)
# calculate proposed steps (differences not destinations) within next thinning interval
zxl <- if( (length(blocks)!=1) && ctrl$useConditionalProposal ){
lapply( iPops[isPops], function(iPop){
.sampleStatesConditional(Z=ZPops[[iPop]],mZ=mZPops[iPop],nGenBack=nGenBackPops[[iPop]], nSample=ctrl$thin, blocks=blocks )
})
} else {
lapply( iPops[isPops], function(iPop){
.sampleStates(Z=ZPops[[iPop]],mZ=mZPops[iPop],nGenBack=nGenBackPops[[iPop]], nSample=ctrl$thin )
})
}
zx <- abind( zxl, along= 2 ) # combine chains (and steps within each chain) of all populations
# second dim: nSample x nChainPop x nPop
genPropRes <- .generateXPropThin(zx,ctrl=ctrl, nChain=length(isPops)*nChainPop)
#
# check if proposals should be recorded
# for the first time
iPopsRecord <- which( (1:nPop %in% isPops) & (iThin0 == nThinOmitRecordPops))
for( iiPop in iPopsRecord ){
# record the initial state of the proposals record
iisPop <- match(iiPop,isPops) # index in iPop
YPops[[iiPop]][1,seq_along(parNames),] <- chainState$X[ ,chainsPop[[ iisPop ]] ]
YPops[[iiPop]][1,nParm+iBlocks,] <- 1 #TRUE
YPops[[iiPop]][1,nParm+nBlock+seq_along(resCompNamesUFlat),] <- chainState$logDenCompX[ ,chainsPop[[ iisPop ]] ]
}
# after the proposal has been made
isRecordProposalsPop = ((iThin0 >= nThinOmitRecordPops))[isPops] # only record the last proposals after nThinOmitRecord
#
# numeric matrix (nGenThin x nPop)
#tempGlobalThinStepsL <- lapply( iPops[isPops], function(iPop){ tempGlobalPops[[iPop]][ 1+iGen ] })
#tempGlobalThinSteps <- structure( do.call( cbind,tempGlobalThinStepsL ), dimnames=list(steps=NULL,pops=NULL) )
# numeric array (nGenThin x nResComp x nPop)
tempDenCompThinStepsL <- lapply( iPops[isPops], function(iPop){ tempResCompPops[[iPop]][1+iGen, ,drop=FALSE] })
tempDenCompThinSteps <- structure( abind( tempDenCompThinStepsL, along=3 ), dimnames=list(steps=NULL,resComp=resCompNamesUFlat, pops=NULL) )
#
# here may use code in .tmp.f.testStep
#
#-- do the steps of next thinning interval
resUpdate <- .updateIntervalTwDEMCParChains( X=chainState$X, logDenCompX=chainState$logDenCompX
, intermediatesChains=chainState$intermediatesChains
, parUpdateDen=chainState$parUpdateDen
, xStep=genPropRes$xStep
, rExtra=genPropRes$rExtra
#, tempGlobalSteps=tempGlobalThinSteps
, tempDenCompSteps=tempDenCompThinSteps
, nsPop=length(isPops)
, pAcceptPar=pAcceptPops[,isPops ,drop=FALSE]
,remoteFunArgs=nameRemoteArgs
,debugSequential=debugSequential
,isRecordProposalsPop=isRecordProposalsPop
,isPops=isPops
#,freeMasterNode=ctrl$freeMasterNode
)
# if( any(!is.finite(resUpdate$logDenCompX)) ){
# if( length(remoteDumpfileBasename) ) dump.frames(remoteDumpfileBasename, TRUE)
# stop("twDEMCBlocks (10): encountered non-finite resLogDen of accepted states")
# }
chainState <- resUpdate[ names(chainState) ]
#
#-- calculate accepteance rate
#row in acceptance Window to record acceptance, if exceeds window, copy second part to first (rewind)
acceptPos0 <- aThinWinW + (iThin0 %% aThinWinW)
if( acceptPos0 == aThinWinW ){ #interval exeeded or new, copy second half to first half
acceptWindow[ 1:aThinWinW,, ] <- acceptWindow[ aThinWinW+(1:aThinWinW),, ]
acceptWindow[ aThinWinW+(1:aThinWinW),, ] <- NA
}
acceptWindow[ acceptPos0+1,, isChains] <- resUpdate$accepted
curAcceptRows <- (acceptPos0+1)-(min(aThinWinW-1,max(iThin0,8)):0) # at the start average over at least 4 slots,
#acceptWindow[curAcceptRows,]
pAcceptChains <- apply(acceptWindow[curAcceptRows,, ,drop=FALSE],c(2,3),sum) / (length(curAcceptRows)*ctrl$thin)
#if( !all(is.finite(pAcceptChains[,isChains ,drop=FALSE])))
# stop("twDEMCBlocks: encountered non-finite pAcceptChains")
pAcceptPops <- popMeansTwDEMC(pAcceptChains, nPop) # current acceptance rate of population (block x pop)
pAcceptPops[] <- pmax(1/ctrl$thin,pAcceptPops ) # lower bound
pAcceptMinPops <- apply(pAcceptPops,2,min) # minimum of acceptance rate across blocks for each population
nGenBackPops <- pmin(M0Pops,ceiling(ctrl$gInd * pmax(1,ctrl$pIndStep/(pAcceptMinPops * ctrl$thin)))) #number of genrations to select states from for each population, ctrl$gInd is multiplied by number of rows for one step depending on acceptance rate and thinning but at least one
#
#-- record current thinning step for each population
for( iiPop in seq_along(isPops) ){
iPop <- isPops[iiPop] # translate index in current populations to all populations
mZ <- mZPops[ iPop ]
ipChains <- chainsPop[[iiPop]] # chains within current chains
igChains <- chainsPop[[iPop]] # chains within all chains
ZPops[[iPop]][mZ+1,,] <- chainState$X[,ipChains]
resLogDen[[iPop]][iThin0+2,,] <- chainState$logDenCompX[,ipChains] # first state is initial state
pAccept[[iPop]][iThin0+2,,] <- pAcceptChains[,igChains]
#if( iPop==13 && any(chainState$X["a",ipChains] < -0.6) ) recover()
}
#
#-- record poprosals
iYSteps0 <- (iThin0 - nThinOmitRecordPops)*ctrl$thin
for( iiPop in which(isRecordProposalsPop) ){
# recordProposals refers to current populations (some may have dropped out already)
iPop <- isPops[iiPop]
iiChains = chainsPop[[iiPop]]
YPops[[iPop]][iYSteps0[iPop]+(1:ctrl$thin),,] <- resUpdate$Y[,,iiChains]
}
}# for iThin0
if( length(progressOutput)!=0 && nchar(progressOutput)!=0 ){ cat("\n") }
# recored the last parUpdateDen for dropouts
iChainsOut <- do.call( c, lapply( isPops, function(iPop){ chainsPop[[iPop]] }) )
parUpdateDenLast[,,iChainsOut] <- chainState$parUpdateDen
#-- calculate the log-Density of last states that are not up to date, because they might be used for reinitialization
for( iDen in iDens ){
dInfo <- dInfos[[iDen]]
cd <- dInfo$compPosDen # components used by the current density function
udi <- adrop(parUpdateDenLast[iDen,, ,drop=FALSE],1) # (iDen= x par x chain
iChains <- which( apply(udi,2,function(ud){ !all(ud) }) )
# TODO parallelize those calculations
for( iChain in iChains ){
iPop <- popChain[iChain]
iChainInPop <- ((iChain-1) %% nChainPop) +1
Xc <- ZPops[[iPop]][(M0Pops+nThinnedGenPops)[iPop],cd,iChainInPop ]
#resCompDenC <- do.call( dInfo$fLogDen, c(list(Xc), dInfo$argsFLogDen) ) # evaluate logDen
resCompDenC <- .calcFLogDen( dInfo$fLogDen, xP=Xc, argsFLogDen=dInfo$argsFLogDen, compInt=c(), TiInt=c(), maxLogDen = dInfo$maxLogDen ) # evaluate logDen
resLogDen[[iPop]][ 1+nThinnedGenPops[iPop],dInfo$resCompPos,iChainInPop] <- resCompDenC
} # end need update
} # for iDen
#-- sum over resLogDen components to yield logDens for each given density
# (might be slightly off because not up to date, but work for inspecting the trend)
for( iDen in iDens ){
dInfo <- dInfos[[iDen]]
rcp <- dInfo$resCompPos # components used by the current density function
if( length(rcp) == 1){
for( iPop in iPops ){
logDen[[iPop]][,iDen,] <- resLogDen[[iPop]][,rcp,]
}
}else{
for( iPop in iPops ){
resLogDenI <- resLogDen[[iPop]][,rcp, ,drop=FALSE]
logDen[[iPop]][,iDen,] <- apply(resLogDenI,c(1,3),sum)
}
} # length(resComp) == 1
} # iDen
##details<< \describe{ \item{subspace localized sampling}{
## The covariance matix for new proposals might look very different in different locations of the parameter space.
## For example, overall CO2 efflux might be composed of two parts microbial respiration and inorganic CO2 production.
## that are relevant in different parts of the parameter space (governed by another parameter). The mirobial growth
## rate will be very sensitive in the respiration part and is allowed to to only small jumps, while it might vary large
## in the inorganic part. Hence, it might be beneficial to use different jumping proposal and different population
## samples in different subspaces.
##
## In order to support such supspace-localized sampling XX
##
## When using parallel execution on multiple cores, the transfer of very big intermediates can cause
## performance issues. Use argument \code{controlTwDEMC$returnIntermediate=FALSE} to ensure that
## intemediates are stored only within process boundaries. When executing the next interval or the next round of
## updates of all blocks, the intermediate will then be the empty list initially.
##}}
##details<< \describe{ \item{intermediate results}{
## Different density functions can share intermediate results, so that these do not need to be recomputed.
## For example two densities can be based on the same model predictions for given parameters, but compare different
## subsets of the predictions against different observations. If for a given parameter vector, the model output has been
## computed in one density, it does not need to be recomputed in the second density function.
##
## This is an advanced option. Care must be taken that the intermediate is really the same between
## densities. And care must be taken that intermediate is recalculated if the parameters change on
## which the intermediate is based on.
##
## Description of both block update in argument \code{blocks} (or argument \code{dInfos} if default Metropolis is used)
## should specify the same identifier string in argument \code{dInfo$intermediateId}.
## The densities or block update functions should return the model output in attribute "intermediate"
## with the result vector of logDensity
## components. \code{twDEMC} ensures that this result is provided with argument \code{intermediate}
## in another call to the density or block update function.
## With a new parameter proposal in Metropolis update, an empty list will be given with this argument.
##
## For an example see test case \code{ofMultiIntermediate} in file unitTests/runittwDEMC.R.
## Using densities \code{\link{denSparsePrior}} and \code{\link{denRichPrior}}.
##
## When using parallel execution on multiple cores, the transfer of very big intermediates can cause
## performance issues. Use argument \code{controlTwDEMC$returnIntermediate=FALSE} to ensure that
## intemediates are stored only within process boundaries. When executing the next interval or the next round of
## updates of all blocks, the intermediate will then be the empty list initially.
##}}
#-- return
for( iPop in iPops ){
if( any(!is.finite(resLogDen[[iPop]])) ){
if( length(remoteDumpfileBasename) ) dump.frames(remoteDumpfileBasename, TRUE)
stop("twDEMCBlocks (20): encountered non-finite resLogDen of accepted states")
}
}
spacesPop <- sapply(pops,"[[","spaceInd")
##value<< a list with entries of populations, each entry is a list
res <- list(
thin=ctrl$thin ##<< thinning interval that has been used
,dInfos=dInfos ##<< list of information on densities (argument \code{dInfos})
,blocks=blocks ##<< list of information on blocks (argument \code{blocks})
,YPos = list( ##<< list of column positions in pops[i]$Y, a list with entries \describe{
##describe<<
accepted = nParm+iBlocks ##<< integer vector of positions of acceptance indication of block at given index
,resLogDen0 = nParm+nBlock ##<< integer scalar: postion before first column of results of fLogDen
)
##end<<
,pops = lapply( iPops, function(iPop){ list( ##<< info on each population. A list with entries: \describe{
##describe<<
upperParBounds = upperParBoundsPop[[iPop]] ##<< upper parameter bounds for sampling
,lowerParBounds = lowerParBoundsPop[[iPop]] ##<< lower parameter bounds for sampling
,splits=pops[[iPop]]$splits ##<< named numeric vector: splitting history
,spaceInd = spacesPop[iPop] ##<< the space replicate that the population belongs to
,parms = ZPops[[iPop]][ M0Pops[iPop]:nrow(ZPops[[iPop]]),, ,drop=FALSE] ##<< numeric array (steps x parms x chains): collected states, including the initial states
,temp = tempResCompPops[[iPop]][seq(1,nGenPops[iPop]+1,by=ctrl$thin), ,drop=FALSE] ##<< numeric array (nSample+1, nResComp): temperature, i.e. cost reduction factor in each step for each datastream
,pAccept= pAccept[[iPop]] ##<< acceptance rate of chains (nStep x nChainPop)
,resLogDen = resLogDen[[iPop]] ##<< numeric array (steps x resComps x chains): results components of fLogDen of blocks
,logDen = logDen[[iPop]] ##<< numberic array (steps x iDen x chains): results summed over blocks
,Y = YPops[[iPop]] ##<< numeric matrix (steps x nParam+nBlock+nResComp x chains): parms, isAccepted for each block, and all fLogDen result components for each proposal
)})
##end<<
)
##end<<
#names(res$pops) <- names(pops)
class(res) <- c( class(res), "twDEMCPops" ) #monte carlo proposal list
res
}
attr(twDEMCBlockInt,"ex") <- function(){
# see unit test
# twUtestF("twDEMCPops","distinctLogDen")
data(twLinreg1); attach( twLinreg1 )
# tracing error in remote session:
# Provide argument \code{remoteDumpfileBasename="dumpFile"}
# Then a dumpfile is created on remote error by \code{\link{sfRemoteWrapper}}.
# In order to trace the density function, the following can be done
if( FALSE ){
.remoteDumpfileBasename="dumpfile"
.remoteDumpfile <- paste(.remoteDumpfileBasename,".rda",sep="")
load(.remoteDumpfile)
debugger(get(.remoteDumpfileBasename))
# choose last step (18)
#require(debug)
fDen <- remoteFunArgs$argsUpdateBlocksTwDEMC$argsFUpdateBlocks[[1]]$fLogDen
mtrace(fDen); remoteFunArgs$argsUpdateBlocksTwDEMC$argsFUpdateBlocks[[1]]$fLogDen <- fDen
do.call(remoteFun, c(remoteFunArgs, list(...)))
# go()
# qqq()
}
}
.tmp.f <- function(){ # testStep
# test one step for one chain
iStep <- 1
iChain <- 1
# those are arguments, generated updateBlocksTwDEMCPar
iPop <- popChain[iChain]
step <- genPropRes$xStep[,iChain,iStep]
rExtra <- genPropRes$rExtra[iChain,iStep]
Xc <- chainState$X[,iChain]
logDenCompC <- chainState$logDenCompX[,iChain]
parUpdateDenC <- chainState$parUpdateDen[,,iChain]
tempC <- tempThinSteps[iChain,iStep]
.updateBlocksTwDEMC(iPop = iPop,step=step,rExtra=rExtra,Xc=Xc,logDenCompC=logDenCompC
,parUpdateDenC=parUpdateDenC,pAccept=pAccept, temp=tempC
, argsUpdateBlocksTwDEMC=remoteArgsFUpdateBlocksTwDEMC$argsUpdateBlocksTwDEMC )
}
.tmp.f <- function(){
# before (code here) calculated temperature for each chain, but actually need it only by population
tempThinStepsL <- lapply( iPops, function(iPop){
matrix( temp[[iPop]][ 1+iGen ], nrow=length(iGen), ncol=nChainPop )
})
tempThinSteps <- t(do.call( cbind, tempThinStepsL)) # steps need to be in last dimension (so chains in first)
}
.getRescompNameBlock <- function(
### appending iBlock to the resCompNames to make the names unique among blocks
resCompNames
,iBlock
){
paste(resCompNames,iBlock,sep="_")
}
.completePopDescription <- function(
### filling default values in pop description
pop
,nGen=numeric(0) ##<< number of generations to take: overwrites settings given in pop
,spaceInd=numeric(0) ##<< the space replicate that this population belongs to, overwritten by setting in pop
){
##details<<
## If several entries are null, they are set to default values
if( 0==length(pop$parms) || !is.numeric(pop$parms) || length(dim(pop$parms)) != 3 )
stop(".checkPop: entry parms must be a numeric 3D array")
if( 0 != length(nGen) ) pop$nGen <- nGen
if( (1 != length(pop$nGen)) )
stop(".checkPop: scalar integer nGen must be given as entry nGen in population or in invocation of twDEMCBlockInt.")
if( 0 == length(pop$spaceInd) ) pop$spaceInd = spaceInd
if( (1 != length(pop$spaceInd)) )
stop(".checkPop: scalar integer spaceInd must be given as entry nGen in population or in invocation of twDEMCBlockInt.")
# if( 0==length(pop$T0) ) pop$T0 <- 1
# if( 0==length(pop$TEnd) ) pop$TEnd <- 1
if( 0==length(pop$X) ){
pop$X <- adrop( pop$parms[nrow(pop$parms),, ,drop=FALSE],1)
}else{
if( !all.equal( rownames(pop$X), colnames(pop$parms), check.attributes = FALSE) )
stop(".checkPop: rownames of X and colnames of parms must match")
}
if( 0 != length(pop$upperParBounds) && any(is.na(pop$upperParBounds)))
stop(".checkPop: encountered NA in upperParBounds")
if( 0 != length(pop$lowerParBounds) && any(is.na(pop$lowerParBounds)))
stop(".checkPop: encountered NA in lowerParBounds")
### pop with entries nGen, and X defined
pop
}
.completeDInfo <- function(
### filling default values in information on logDensities (dInfo)
dInfo ##<< list of information on logDensities to be checked for entries
,parNames ##<< string vector: names of all parameters. A subset (compPosDen) of these is required by the logDensity
,resCompNames ##<< string vector: list of names of result vector of logDensity function.
){
##details<<
## \code{dInfo} stores ample information on a log-Density function in order to
## help fast execution of calling the log-Density function and the metropolis update. However, the user
## is required to specify only minimal detail (see \code{\link{twDEMCBlockInt}}) and most of the information is
## filled in this function.
## Some of the information can only be completed after the names of logDensity results are
## know after calling the function. Therefore there is a second function \code{\link{.completeDInfoResCompNames}}.
#
if( 0 == length(dInfo$fLogDen) || !is.function(dInfo$fLogDen))
stop(".checkLogDenInfo: entry fLogDen (function of log of unnomalized density) must be specified")
#
##details<<
## \code{compPosDen} specifies the parameters provided to the given logDensity.
## If entry compPosDen is not given, then assume all parameter dimensions are required.
## If compPosDen is given as a string vector, then it is translated to positions within
## vector of all parameter names.
if( 0 == length(dInfo$compPosDen) ) dInfo$compPosDen <- seq_along(parNames)
compPosDenOrig <- dInfo$compPosDen
if( is.character(compPosDenOrig) ) dInfo$compPosDen <- match(compPosDenOrig,parNames)
if( any(is.na(dInfo$compPosDen)) ){
stop(".completeDInfo: some of compPosDen specified in dInfo are among parameter names (of initial values).")
}
#
if( 0==length(dInfo$intResComp) ) dInfo$intResComp <- vector("character",0)
#if( 0==length(dInfo$TFix) ) dInfo$TFix <- vector("numeric",0)
#
if( 0!=length(dInfo$TFix) ) warning("Usage of TFix is deprecated. Use arguemt TSpec of twDEMCBlockInt with T0 and TEnd of one for given components.")
if( 0!=length(dInfo$intermediate) || 0!=length(dInfo$intermediateId) ) warning("intermediateId now must be specified with blocks.")
if( length(i <- which(!(names( dInfo) %in% c("fLogDen","xC","logDenCompC","parUpdateDenC","compPosDen","intResComp","argsFLogDen"
,"fDiscrProp","argsFDiscrProp"
,"resCompNames","resCompPos","intResCompPos","intResCompPosWithin","nObs","maxLogDen")))) )
stop(paste(".completeDInfo: unknown entries in dInfo:",paste(names(dInfo)[i],collapse=",") ))
if( 0 != length(dInfo$nObs))
warning(".checkDInfoResCompNames: usage of nObs in dInfos is deprecated and not considered. Instead use: maxLogDen=-1/2*nObs")
#
dInfo$resCompNames <- resCompNames
#
if( 0 == length(dInfo$maxLogDen) ) dInfo$maxLogDen <- rep( 0, length(resCompNames) )
if( length(names(dInfo$maxLogDen)) ){
# check that all names matches the resCompNames, and sort along resCompNames
iM <- match( resCompNames, names(dInfo$maxLogDen) )
if( any(is.na(iM)) )
stop(paste(".checkLogDenInfo: maximum logDensitiy, maxLogDen, of observations missing for result components ",paste(resCompNames[which(is.na(iM))], collapse=","),sep="") )
dInfo$maxLogDen <- dInfo$maxLogDen[iM]
}
if( length(dInfo$maxLogDen) != length(resCompNames) )
stop(paste("dInfo maxLogDen (length=",length(dInfo$maxLogDen),") must have same lenght as number of result components=",length(resCompNames),sep=""))
### argument \code{dfo} with entries fLogDen, compPosDen, intResComp, resCompNames, and maxLogDen properly specified
dInfo
}
.completeBlockInfo <- function(
### filling default values in block description
block ##<< list to be checked for entries
,dInfos ##<< checked list of dInfo. See the argument in \code{\link{twDEMCBlockInt}}.
,parNames ##<< parameter dimension names
){
##details<<
## Similar to dInfos (see \code{\link{.completeDInfo}}), there is a list with ample information on each block.
## As the user is required to specify only the minimal information or information in more convenient format
## The block list is checked and information is completed.
#
if( 0 == length(block$compPos) ) stop(".completeBlockInfo: entry compPos must be specified, either as indices or as variable names")
##details<<
## If the user specified \code{compPos} as a character string,
## then convert it to positions within vector of all parameter
compPosOrig <- block$compPos
if( is.character(compPosOrig) ){
block$compPos <- match( compPosOrig, parNames)
}
if( any(is.na(block$compPos)) ){
i <- which(!(names(compPosOrig) %in% parNames))
stop(paste(".completeBlockInfo: unknown variable names in entry compPos in block:",paste(parNames[i],collapse=",") ))
}
##details<<
## If no other block update function is given, assume Metropolis update
if( 0 == length(block$useMetropolisUpdate) )
block$useMetropolisUpdate = (0==length(block$fUpdateBlock))
if( block$useMetropolisUpdate ){
block$fUpdateBlock <- updateBlockTwDEMC
dInfoPosOrig <- block$dInfoPos
if( 1 != length(block$dInfoPos) ) stop(".checkBlock: dInfoPos must be of length 1.")
if( is.character(dInfoPosOrig) ) block$dInfoPos <- match(dInfoPosOrig,names(dInfos))
if( is.na(block$dInfoPos) ){
stop(paste0(".checkBlock: specified density '",dInfoPosOrig,
"' not found in names of list of density information (argument dInfos)."))
}
dInfo <- dInfos[[block$dInfoPos]]
#
# positions of parameters to update in list of parameters to logDensity function
block$compPosInDen <- match(block$compPos,dInfo$compPosDen )
if( any(is.na(block$compPosInDen)) ){
stop(".checkBlock: some of compPos specified in block are not in compPosDen of fLogDenInfo.")
}
} # useMetropolisUpdate
#
if( 1!=length(block$fUpdateBlock) || !is.function(block$fUpdateBlock) ) stop(".checkBlock: fUpdateBlock must be a single function.")
#
if( 0 == length(block$nStates) ){
# number of states to choose from for proposal generations on conditional sampling
block$nStates <- calcM0twDEMC( length(block$compPos), 1)
}
#
##value<< argument block with entries compPos, dInfoPos, fUpdateBlock, and useMetropolisUpdate definded
if( length(i <- which(!(names( block) %in% c("compPos","dInfoPos","fUpdateBlock","useMetropolisUpdate","argsFUpdate"
,"nStates","compPosInDen","intermediateId")))) )
stop(paste(".checkBlock: unknown entries in block:",paste(names(block)[i],collapse=",") ))
block
}
.getResFLogDenNames <- function(
### Extracting/Creating names for result components of logDenComp
logDenComp ##<< vector or array (vars in rows) of results of logDensity function
,logDenName ##<< name of the density
){
##details<<
## Names of the result components of fLogDen are used to distinguish internal components.
## If the density function does not provide names, they are created.
## If the density function has only one result component, the name is lost during parallel execution.
## Hence, a single name is also re-created, to avoid errors on checking.
## The created names are of the form <logDenName>_<iComp>
if( is.array(logDenComp)){
res <- rownames(logDenComp)
return( if( !is.null(res) ) res else paste(logDenName,1:nrow(logDenComp),sep="_") )
}else{
res <- names(logDenComp)
#return( if( !is.null(res) && length(res) > 1) res else paste(logDenName,1:length(logDenComp),sep="") )
return( if( !is.null(res) ) res else paste(logDenName,1:length(logDenComp),sep="_") )
}
}
.sampleStatesConditional <- function(
### Sample records and initial state from past records.
Z ##<< population (steps x parms x chains)
,mZ ##<< number of steps in Z that have been already sampled (Z is bigger)
,nGenBack ##<< number of samples in recent past, used to generate the proposals
,nSample ##<< number of proposals to generate
,blocks ##<< information on blocks with item \code{compPos} specifiy the parameter positions
##<< and item \code{nStates} specifying the number of closest states to choose from, usually calculated by CalcM0 with one chain
#,rLogDen
){
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.generateXPropThin}}
nChainPop = dim(Z)[3]
nParm = dim(Z)[2]
# blocks <- lapply(blocks, function(block){ within(block, nStates <- 8) })
# calculate the number of states
# blocks <- lapply(blocks, function(block){ within(block, nStates <- m0b*nChainPop) })
Zt <- t(apply(Z[1:(mZ),,],2,function(Zs){Zs})) # matrix (nParm x nState): all states except the last (current) one
#iChain <- 1
zxbL <- lapply( 1:nChainPop, function(iChain){
Xc <- Z[mZ,,iChain] # current state
d <- abs(Xc - Zt) # absolute deviation from current state (as parameters are the first dimension, we can apply array diff)
#block <- blocks[[1]]
ZtsubL <- lapply( blocks, function(block){
db <- d[-block$compPos, ,drop=FALSE] # subset of parameters not in block
rbp <- apply( db, 1, rank) # rank of difference for each parameter
rb <- apply( rbp,1, max) # maximum over parameters
iStates <- order(rb)[1:block$nStates] # the index of the states to choose from
iS <- matrix( iStates[sample.int(block$nStates, 3*nSample, replace=TRUE)], ncol=3 ) #TODO save invocations to sample.int by sampling before
Ztsub <- Zt[block$compPos, iS,drop=FALSE]
Xcb <- Xc[block$compPos]
zx <- array(c(Ztsub, rep(Xcb,nSample)), dim=c(parms=nrow(Ztsub), step=nSample, iz=4), dimnames=list(parms=rownames(Ztsub),step=NULL,iz=NULL) )
# check that z1 and z2 are different
# check that z3 is different from x
##details<<
## States z2 is attempted to be distinct.
## And z3 is attempted to be distinct from x (assuming steps in z are stacked chains collectives of X)
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,2] <- zx[,sample.int(dim(zx)[2],length(iSame)),2] # sample from z2 of other steps
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i<-i+1
}
#(tmp <- which( zx[,,2] == zx[,,1], arr.ind=TRUE))
#when z3 is the same as x, i.e. zx[,,4]
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
#iSame <- unique(which( (zx[-nrow(zx),,3]==Xs), arr.ind=TRUE )[,2])
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,3] <- zx[,sample.int(dim(zx)[2],length(iSame)),3] # sample from z3 of other steps
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
i<-i+1
}
zx
}) # updates for each block
zxb <- abind(ZtsubL, along=1)
}) # updates for each chain
zxbc <- abind(zxbL, along=2) # stack chains across samples
### numeric array: (nParm x (nSample*nChainPop) x 4) of randomly choosen states
### states of different blocks are choosen independently of each other, i.e. are combinations from different states
### 4th state is the currently accepted state of the chain (for Snooker update)
zxc <- zxbc[ colnames(Z),, ] # order parameters
}
.sampleStates <- function(
### Sample records and initial state from past records.
Z ##<< population (steps x parms x chains)
,mZ ##<< number of steps in Z that have been sampled
,nGenBack ##<< number of samples in recent past, used to generate the proposals
,nSample ##<< number of proposals to generate
#,rLogDen
){
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.generateXPropThin}}
nChainPop = dim(Z)[3]
nParm = dim(Z)[2]
##details<<
## Random states for chains for difference vector generation are within subsets of populations.
## This allows simulating independent population of chains.
## The acceptance rate may differ among populations. Hence, the set of previous generations to
## randomly select from also differs between poplations.
# integer array (thinSteps*nChain*4) sample of chains, within populations
X <- adrop(Z[mZ,,,drop=FALSE],1) # assume current state X as beginning of the interval
nStates <- nSample*nChainPop*3 # need to sample three states for snooker update
sChains <- 1:nChainPop
sGens <- (mZ-nGenBack+1):mZ
rrGenPop <- sample(sGens, nStates, replace = TRUE)
rrChainsPop <- sample(sChains, nStates, replace = TRUE)
# in order to constrain two dimensions at the same time use the [] subset with an array see ?"["
rr <- cbind(rrGenPop,rrChainsPop)
#zLogDen <- array(rLogDen[rr], dim=c(1,nSample*nChainPop,3), dimnames=list(parms="logDen", steps=NULL, zi=NULL) )
rrParms <- cbind( rep(rrGenPop,each=nParm), rep(1:nParm, nStates), rep(rrChainsPop,each=nParm) )
zParms <- array(Z[rrParms], dim=c(nParm,nSample*nChainPop,3), dimnames=list( parms=rownames(Z), steps=NULL, zi=NULL) )
# note that parameters are the first dimension
#rrParms <- cbind( rep(rrGenPop,nParm), rep(1:nParm, each=nStates), rep(rrChainsPop,nParm) )
#zParms <- matrix( Z[rrParms], ncol=nParm)
#z0 <- abind( zParms, zLogDen, along=1)
#
# append as forth initial state x vector to z
# assume that chains is the last dimension, for each step assume the same initial state x
chainsZ <- rep(sChains,each=nSample)
Xs <- array(X[,chainsZ], dim=c(nParm,nSample*nChainPop,1), dimnames=list(parms=rownames(Z), steps=NULL, zi=NULL) )
#XsLogDen <- array( rLogDen[mZ,chainsZ],dim=c(1,nSample*nChainPop,1), dimnames=list(parms="logDen", steps=NULL, zi=NULL) )
#Xs <- abind( Xs, XsLogDen, along=1) #logDen row never used for X
zx <- abind( zParms, Xs, along=3)
#
##details<<
## States z2 is attempted to be distinct.
## And z3 is attempted to be distinct from x (assuming steps in z are stacked chains collectives of X)
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,2] <- zx[,sample.int(dim(zx)[2],length(iSame)),2]
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i<-i+1
}
#(tmp <- which( zx[,,2] == zx[,,1], arr.ind=TRUE))
#when z3 is the same as x, i.e. zx[,,4]
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
#iSame <- unique(which( (zx[-nrow(zx),,3]==Xs), arr.ind=TRUE )[,2])
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,3] <- zx[,sample.int(dim(zx)[2],length(iSame)),3]
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
i<-i+1
}
#(tmp <- which( zx[,,3] == zx[,,4], arr.ind=TRUE))
### random states (Nparms+1,(steps*nChainPop), 4)
### first dimension is the state vector
### random states for each step and chains (chains stacked to be vectorized)
### chain is last dimensionion in stack (consequtive steps for one chain) in order to abind across populations
### zx dim: three random vectors, forth dimension is the initial state x of the chain
rownames(zx) <- colnames(Z) # preserve parameter names
zx
}
.generateXPropThin <- function(
### Generate Proposals from given random states.
zx ##<< output of \code{\link{.sampleStates}} (maybe stackes states)
,ctrl ##<< list with tuning parameters, see \code{twDEMCBlockInt}
,nChain ##<< number of chains
){
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.sampleStates}}
## \code{\link{.xStepSnooker}}
## \code{\link{.xStepParallel}}
#
nParm <- nrow(zx)
nState <- ncol(zx)
nSample <- nState/nChain
#zxl <- lapply( 1:nPop, fzPop )
#zxAndLogDen <- structure( abind(zxl, along=2) , dimnames=list(parms=c(rownames(Z),"logDen"),steps=NULL,zi=NULL))
#zx <- zxAndLogDen[1:nParm,,,drop=FALSE]
z <- zx[,,1:3,drop=FALSE] #three random state vectors per step
X <- adrop(zx[,,4,drop=FALSE],3) #initial state vector for step
#zLogDen <- adrop(zxAndLogDen[,,1:3,drop=FALSE][nParm+1,,,drop=FALSE],1)
dz <- as.list( structure(dim(z), names=c("parms","steps","zi")) )
#
res <- matrix( NA_real_, nrow=dz$parms+1, ncol=dz$steps, dimnames=c(list(c(rownames(z),"rExtra")),list(NULL)) )
boSnooker <- runif(dz$steps)< ctrl$pSnooker
if( 0 < sum(boSnooker) ){
res[,boSnooker] <- .xStepSnooker(z[,boSnooker,,drop=FALSE],X[,boSnooker,drop=FALSE],ctrl=ctrl)
}
if( 0 < sum(!boSnooker) )
#res[,!boSnooker] <- .xStepParallel(z[,!boSnooker,,drop=FALSE],zLogDen=zLogDen[!boSnooker,,drop=FALSE],ctrl=ctrl)
res[,!boSnooker] <- .xStepParallel(z[,!boSnooker,,drop=FALSE],ctrl=ctrl)
#
#second dimension is Nsteps*nChain (we can set nChains as last dimension
#array(chainsZ,dim=c(nSample,d$chains))
xStepAndExtra <- resArraySteps <- array(res, dim=c(nParm+1,nSample,nChain), dimnames=list(parms=rownames(res),steps=NULL,chains=NULL) )
#deprecated expected steps as last dimension
#xStepAndExtra <- aperm(resArraySteps,c(1,3,2))
# numeric array (Npar+1,Nchains,Nsteps): difference vectors in parameter space for steps and chains
# last row is the extra LogDen associated with snooker update
#
xStep <- xStepAndExtra[-nrow(xStepAndExtra),,,drop=FALSE] #dito
rExtra <- adrop(xStepAndExtra[nrow(xStepAndExtra),,,drop=FALSE],1) #second dim (columns step within Thinning interval)
list(xStep=xStep, rExtra=rExtra)
### List with components \describe{
### \item{xStep}{numeric array (Npar x NSteps x Nchain): difference vectors in parameter space}
### \item{rExtra}{numeric matrix (Nsteps x NChain): some extra LogDensity from snooker update}}
### Nsteps=ctrl$thin
}
.xStepSnooker <- function(
### Generates Snooker updates based on given random numbers.
z, ##<< numeric array (Nparms,(nsteps), 3) of random states, dimnames parms,steps,zi
X, ##<< current state (Nparms,(nsteps)) corresponding to chain of second dimension in z
ctrl
){
# DE-Snooker update
##seealso<<
## \code{\link{.generateXPropThin}}
nParm <- nrow(z)
nState <- ncol(z)
gamma_snooker = runif(nState, min=1.2,max=2.2)
res <- matrix( as.numeric(NA), nrow=nParm+1, ncol=nState, dimnames=c(list(c(rownames(z),"rExtra")),list(NULL)) )
#need loop because of inner products
for( i in 1:nState){
x_z = X[,i] - z[,i,3]
D2 = max(1.0e-300, x_z %*% x_z)
#gamma_snooker =1.7
projdiff = ((z[,i,1] -z[,i,2]) %*% x_z)/D2 # inner_product of difference with x_z / squared norm x_z
res[-(nParm+1),i] <- xStepChain <- (gamma_snooker[i] * projdiff) * x_z
xPropChain = X[,i] + xStepChain
x_z = xPropChain - z[,i,3]
D2prop = max((x_z%*%x_z), 1.0e-30)
res[(nParm+1),i] <- rExtra <- ctrl$Npar12 * (log(D2prop) - log(D2)) # Npar12 =(nParm - 1)/2 # extra term in logr for accept - reject ratio
}
res
} # DE-Snooker update
.xStepParallel <- function(
### DE-parallel direction update based on given random numbers.
z, ##<< numeric array (Nparms,(nsteps), 3) of random states, dimnames parms,steps,zi
#zLogDen, ##<< numeric matrix (nsteps,3): logDen corresponding to the random states z
ctrl
){
##seealso<<
## \code{\link{.generateXPropThin}}
nParm <- nrow(z)
nState <- ncol(z)
dz <- adrop((z[,,1,drop=FALSE]-z[,,2,drop=FALSE]),3) #jump vector as the difference between two random states (from z2 towards z1)
boGammasF1 <- (runif(nState) < ctrl$pGamma1)
gamma_par <- matrix( ctrl$F1, nrow=nParm, ncol=nState) # to be able to jump between modes
gamma_par[,!boGammasF1] <- ctrl$F2 * runif( nParm*sum(!boGammasF1), min=1-ctrl$epsMult, max=1+ctrl$epsMult) # multiplicative error to be applied to the difference
xStepChain <- if (ctrl$epsAdd ==0) { # avoid generating normal random variates if possible
gamma_par * dz
} else {
gamma_par * dz + rnorm(length(dz),0,ctrl$epsAdd)
}
xStepChainAndRExtra <- rbind(xStepChain,rExtra=0)
### Numeric matrix (Nparms, nsteps): difference vectors in parameter space.
}
.tmp.f <- function(){
if( !is.null(ctrl$probUpDir) && (ctrl$probUpDir != 1/2) ){
iFinites <- which(is.finite(zLogDen[,1]) & is.finite(zLogDen[,2]))
}
}
.updateIntervalTwDEMCParChains <- function(
### Perform update steps within next thinning interval parellelizing by chains.
# first three arguments relate to current state (stored in chainState in twDEMCBlockInt)
X ##<< numeric matrix: current location of chains rows: parameters, columns: chains
,logDenCompX ##<< numeric matrix: result of calls to density functions from previous step
,intermediatesChains = lapply(1:nsChain, function(iChain){ list()} ) ##<< list(nChain) intermediate states
,parUpdateDen ##<< for each parameter/density combination: is the density up to date from previous step
,xStep ##<< array rows: difference vectors in parameter space, cols: steps, zdim: chains
,rExtra ##<< numeric matrix: row: step within thinning interval, col: chain
,tempDenCompSteps ##<< numeric array: (nGenThin x nResComp x nsPop)
,nsPop ##<< the number of populations that take part in this step
,pAcceptPar ##<< numeric matrix: current acceptance rate for each (block x population)
#argsUpdateBlocksTwDEMC, ##<< list of arguments to updateBlocksTwDEMC
#argsDEMCStep, ##<< see \code{\link{.doDEMCStep}}
,remoteFunArgs ##<< arguments that do not change across calls
## may be specified as a name instead of complete object (see code{\link[twSnowfall]{sfRemoteWrapper}})
## Must include entries \itemize{
## \item{remoteFun=.updateBlocksTwDEMCChains}
## \item{argsUpdateBlocksTwDEMC
## \item{remoteDumpfileBasename}
## }
## Can be a name of a previously exported variable.
,debugSequential=FALSE ##<< see \code{\link[twSnowfall]{sfFArgsApplyDep}}
,isRecordProposalsPop=rep(FALSE,nsPop) ##<< logical vector: for each population: if TRUE then proposals and results of rLogDen are recorded in result$Y.
, chainsPop=matrix(1:ncol(X),ncol=nsPop, dimnames=list(iChainInPop=NULL, iPop=NULL)) ##<< list of integer vectors:
, isPops ##<< integer vector of populations that have not yet dropped out
){
# .updateIntervalTwDEMCParChains
dimXStep <- dim(xStep)
nStep <- dimXStep[2]
nsChain <- dimXStep[3]
nChainPop = nsChain / nsPop
#
# if( length(intermediatesChains[[1]]) ) recover()
#
#iChain <- 1
argsListsChain <- lapply( 1:nsChain, function(iChain){
iiPop <- ((iChain-1) %/% nChainPop) +1
iPop <- isPops[ iiPop ] # refer to population in non-dropped array
list( X = X[,iChain]
,logDenCompX = logDenCompX[,iChain]
, intermediates = intermediatesChains[[iiPop]]
, parUpdateDen = adrop(parUpdateDen[,,iChain,drop=FALSE],3)
, xStep = adrop(xStep[,,iChain,drop=FALSE],3)
, rExtra = rExtra[,iChain]
, tempDenCompSteps = adrop(tempDenCompSteps[,,iiPop,drop=FALSE],3)
, pAcceptPar = pAcceptPar[,iiPop]
, isRecordProposalsPop = isRecordProposalsPop[iiPop]
, iPop = iPop
)
})
#
#remoteFunArgs$remoteFun <- .updateIntervalTwDEMCChain # on changed function
res <- if( isTRUE(debugSequential) ){
lapply( argsListsChain, remoteFunArgs$remoteFun, argsUpdateBlocksTwDEMC=remoteFunArgs$argsUpdateBlocksTwDEMC )
}else{
# if( sfParallel() ) ( list=c("nChainPop", "X", "logDenCompX", "parUpdateDen", "xStep"
# , "rExtra", "tempDenCompSteps", "pAcceptPar"
# , "remoteFunArgsChain", "isRecordProposalsPop","returnIntermediateChains"))
sfLapply( argsListsChain, sfRemoteWrapper, remoteFunArgs = remoteFunArgs )
}
##value<< list with components
res1 <- res[[1]]
X <- array( NA_real_, c(length(res1$X), nsChain), dimnames=list(parms=names(res1$X), iChain=NULL))
logDenCompX <- array( NA_real_, c(length(res1$logDenCompX), nsChain), dimnames=list(parms=names(res1$logDenCompX), iChain=NULL))
parUpdateDen <- array( NA, c(dim(res1$parUpdateDen), nsChain), dimnames=c(dimnames(res1$parUpdateDen), iChain=NULL))
accepted <- array( NA_real_, c(length(res1$accepted), nsChain), dimnames=list(iBlock=names(res1$accepted), iChain=NULL))
Y <- array( NA_real_, c(dim(res1$Y), nsChain), dimnames=c(dimnames(res1$Y), iChain=NULL))
for( iChain in 1:nsChain ){
resi <- res[[iChain]]
X[,iChain] <- resi$X
logDenCompX[,iChain] <- resi$logDenCompX
parUpdateDen[,,iChain] <- resi$parUpdateDen
accepted[,iChain] <- resi$accepted
Y[,,iChain] <- resi$Y
}
resDo <- list( ##describe<<
X=X ##<< matrix current position, column for each chain
, logDenCompX=logDenCompX ##<< matrix: result components of fLogDen current position, column for each chain
#, logDenX=logDenX ##<< vector current logDen of chains
, parUpdateDen=parUpdateDen ##<< numeric matrix (par x den x chain): for each parameter/density combination: is the density up to date
, accepted=accepted ##<< numeric matrix (blocks x chains): number of accepted steps for each chain
, Y=if( !is.null(res1$Y) ) Y else NULL ##<< numeric matrix (steps x 2+params+result x chains): accepted, rLogDen, parms, and all fLogDen result components for each proposal
, intermediatesChains = lapply(res,"[[","intermediatesX")
) ##end<<
}
.updateIntervalTwDEMCChain <- function(
### Perform update steps within next thinning interval for a single chain (called remotely from .updateIntervalTwDEMCParChains)
# first three arguments relate to current state (stored in chainState in twDEMCBlockInt)
argsUpdateIntervalTwDEMCChain = list( ##<< argument that change across invocations of .updateIntervalTwDEMCChain
## supplied as list, so that sfLapply can be used
##describe<<
X=c() ##<< numeric vector: current location
,logDenCompX=c()##<< numeric vector: result of calls to density functions for current location
,intermediates=list() ##<< list: intermediate results for current location
,parUpdateDen=c() ##<< for each parameter/density combination: is the density up to date from previous step
,xStep=c() ##<< matrix rows: difference vectors in parameter space, cols: steps
,rExtra=c() ##<< numeric vector: row: step within thinning interval
,tempDenCompSteps=c() ##<< numeric matrix: (nGenThin x nResComp)
#,nsPop ##<< the number of populations that take part in this step
,pAcceptPar=c() ##<< numeric vector: current acceptance rate for each (block)
,isRecordProposalsPop=FALSE ##<< logical
, iPop=1 ##<< the population that the current chain is in
)
##end<<
,argsUpdateBlocksTwDEMC ##<< list: further arguments to \code{\link{.updateBlocksTwDEMC}}
){
a <- argsUpdateIntervalTwDEMCChain # shortcut
# .updateIntervalTwDEMCPar
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.doDEMCStep}}
dimXStep <- dim(a$xStep)
nStep <- dimXStep[2]
#d <- as.list(structure(dimXStep,names=c("parms","chains","steps")))
iGenT <- (1:nStep)
nCases = nStep
if( !(nCases == length(a$rExtra)) )
stop("number of cases in steps must correspond to length of rExtra")
res0 <- list(
xC=a$X
,logDenCompC=a$logDenCompX
,intermediatesC=a$intermediates
,parUpdateDenC=a$parUpdateDen
)
#
# record the proposals, their density results and their acceptance for each generation in thinning interval
parNames <- names(a$X)
resCompNames <- names(a$logDenCompX)
if( is.null(resCompNames) ) resCompNames <- paste("lDen",seq_along(a$logDenCompX),sep="")
tmp <- c(parNames,paste("accepted",seq_along(a$pAcceptPar),sep=""),resCompNames)
Y <- array( NA_real_, dim=c(nStep=nStep,comp=length(tmp)), dimnames=list(steps=NULL,comp=tmp) )
#
# recored acceptance state
acceptedSteps <- array( NA_real_, dim=c(nStep=nStep,block=length(a$pAcceptPar)) )
#
#i <- 1
prevRes <- res0
for( iStep in 1:nStep){
args <-c( prevRes[1:4]
,list( step=a$xStep[,iStep]
, rExtra=a$rExtra[iStep ,drop=TRUE]
#, tempGlobalC=tempGlobalSteps[iStep0+1,isPop ,drop=TRUE]
, tempDenCompC=a$tempDenCompSteps[iStep,, drop=TRUE]
, pAccept=a$pAcceptPar
#, isPop=isPop
, iPop=a$iPop
, argsUpdateBlocksTwDEMC=argsUpdateBlocksTwDEMC
)
)
prevRes <- do.call( .updateBlocksTwDEMC, args )
Y[iStep,] <- c( prevRes$xP, prevRes$accepted, prevRes$logDenCompP )
acceptedSteps[iStep,] <- prevRes$accepted
prevRes$accepted
}
#
#if( (iPop==2) && (X["a"] < 10.8) ) recover()
##value<< list with components
ctrl <- argsUpdateBlocksTwDEMC$argsFUpdateBlocks[[1]]$ctrl
resDo <- list( ##describe<<
X=prevRes$xC ##<< vector (nParm) current position
, logDenCompX=prevRes$logDenCompC ##<< vector (nResComp): result components of fLogDen current position
, intermediatesX=if( ctrl$returnIntermediate ) ##<< list: intermediate results for current position
prevRes$intermediatesC else list()
#, logDenX=logDenX ##<< vector current logDen of chains
, parUpdateDen=prevRes$parUpdateDenC ##<< numeric matrix (par x den): for each parameter/density combination: is the density up to date
, accepted=colSums(acceptedSteps) ##<< numeric vector (blocks): number of accepted steps for each chain
, Y=Y ##<< numeric matrix (steps x 2+params+result): accepted, rLogDen, parms, and all fLogDen result components for each proposal
) ##end<<
}
.updateBlocksTwDEMC <- function(
### One step of updating all blocks for one chain
xC ##<< current state
,logDenCompC ##<< result of calls to density functions of current state
,intermediatesC=list() ##<< intermediate results of current state
,parUpdateDenC ##<< matrix (nPar x nDen): FALSE: logDenCompC is not uptodate and needs recalculation
## (because depends on parameters that was updated against other density)
,step ##<< proposed step
,rExtra ##<< adjustment of Metropolis state
,tempDenCompC ##<< numeric vector of length(logDenCompC): current temperature for each result component
,pAcceptC ##<< numeric vector (nBlock) current acceptance rates of the blocks
,iPop ##<< population in total array of populations
,argsUpdateBlocksTwDEMC = list( ##<< further arguments that do not change between calls.
## These can be transferred by sfExport once and are handled by sfRemoteWrapper.
##describe<<
argsFUpdateBlocks=list() ##<< list(nBlock) arguments for each block: dInfo, block
,nBlock=1 ##<< scalar integer: the number of blocks
,upperParBoundsPop=c() ##<< upper parameter bounds
,lowerParBoundsPop=c() ##<< lower parameter bounds
##end<<
)
){
#if( !all(is.finite(pAcceptC)))
# stop(".updateBlocksTwDEMC: encountered non-finited pAcceptC")
a <- argsUpdateBlocksTwDEMC
#with( argsUpdateBlocksTwDEMC,{
acceptedBlock <- numeric(a$nBlock)
xP <- xC # place for recording proposal
logDenCompP <- logDenCompC # place for recording density results for proposal
#
#recover()
#iBlock<-1
#iBlock<-2
#dInfo1 <- a$argsFUpdateBlocks[[1]]
for( iBlock in seq_along(a$argsFUpdateBlocks)){
#block <- blocks[[iBlock]]
blockArgs <- a$argsFUpdateBlocks[[iBlock]] #merged the information in argsFUpdateBlocks
usesIntermediate <- length(blockArgs$intermediateId) != 0
#
argsFUpdateBlock <- c( blockArgs, list(
iPop=iPop
, upperParBounds=a$upperParBoundsPop[[iPop]]
, lowerParBounds=a$lowerParBoundsPop[[iPop]]
)) # will be extendend if Metropolis is used and by intermediate further down
#
##details<< \describe{}{\item{Recalculating logDensity of accepted state.}{
## If one of the components used by the current density function has been
## updated against another density function, then the current density of
## the accepted state is out of date, and needs to be recalculated.
##}}
if( blockArgs$useMetropolisUpdate ){
dInfo <- blockArgs$dInfo
cd <- dInfo$compPosDen # components used by the current density function
if( !all( parUpdateDenC[blockArgs$dInfoPos,cd] ) ){
# intermediate state may have been calculated, no need of recalculation
argsFLogDen <- dInfo$argsFLogDen
argsFLogDen$intermediate <- if( usesIntermediate ) intermediatesC[[blockArgs$intermediateId]] else NULL
# here do not allow for internal rejection
logDenCompC[dInfo$resCompPos] <- tmp <- .calcFLogDen( dInfo$fLogDen, xP, compInt=dInfo$logDenCompC[dInfo$intResCompPos], TiInt = dInfo$tempDenCompC[dInfo$intResCompPos]
,argsFLogDen = argsFLogDen
,maxLogDen = dInfo$maxLogDen
)
#do.call( dInfo$fLogDen, c(list(xC[cd]), dInfo$argsFLogDen) ) # evaluate logDen
#undebug(ofDalecPhen)
#tmp2 <- do.call(ofDalecPhen, c(list(xC[cd]), dInfo$argsFLogDen) ) # evaluate logDen
if( any(!is.finite(tmp)) ){
dump.frames("updateBlocksTwDEMC",TRUE)
stop(".updateBlocksTwDEMC (25): encountered non-finite resLogDen of in chainState. Dump.frames to file updateBlocksTwDEMC.rda. Implement better consistency checks across all fLogDen.")
}
parUpdateDenC[blockArgs$dInfoPos,cd] <- TRUE
if( usesIntermediate ){
intermediatesC[[ blockArgs$intermediateId ]] <- attr(tmp,"intermediate")
}
} # end if intermediate state changed
argsFUpdateBlock$metInfoStep <- list(
step=step
, rExtra=rExtra
, logDenCompC=logDenCompC[ dInfo$resCompPos ]
#, tempC=tempGlobalC
, tempDenCompC=tempDenCompC[ dInfo$resCompPos ]
, pAccept=pAcceptC
)
} # end if( blockArgs$useMetropolisUpdate )
#
# intermediate may have changed in Metropolis-Block
argsFUpdateBlock$intermediate <- if( usesIntermediate ) intermediatesC[[blockArgs$intermediateId]] else NULL
argsFUpdateBlock$intermediateId <- NULL # else this will be uses for intermediate: R bug
#
# make sure to use names different from already existing in argsFUpdateBlocks[[iBlock]]
# in c latter entries overwrite previous entries, so use argsFUpdateBlocks first
# updateBlockTwDEMC will set it to NULL for new parameter proposal
#resUpdate <- updateBlockTwDEMC( xC[cd], argsFUpdateBlock=argsFUpdateBlock )
# call the update
resUpdateB <- do.call( blockArgs$fUpdateBlock, c(list(xC[cd], argsFUpdateBlock=argsFUpdateBlock), argsFUpdateBlock$argsFUpdate) )
#resUpdate <- blockArgs$fUpdateBlock( xC[cd], argsFUpdateBlock=argsFUpdateBlock )
acceptedBlock[iBlock] <- resUpdateB$accepted # record acceptance
if( !all(is.finite(resUpdateB$accepted)))
stop(".updateBlocksTwDEMC: encountered non-finited acceptance rate")
if( any(!is.finite(resUpdateB$logDenCompC)) ){
stop(".updateBlocksTwDEMC (35): encountered non-finite resLogDen in accepted state")
}
cu <- blockArgs$compPos # components to be updated in this block
if( resUpdateB$accepted != 0){
#if( resUpdate$xC["a"] > 10.8) recover()
xC[cu] <- xP[cu] <- resUpdateB$xC # update the components of current state
if( usesIntermediate )
intermediatesC[[ blockArgs$intermediateId]] <- resUpdateB$intermediate # update the intermediate
#undebug(setIntermediate)
parUpdateDenC[,cu] <- FALSE # indicate all density results out of date for the updated parameters
if( 0 != length(resUpdateB$logDenCompC) ){
#if( iBlock==1 && resUpdate$logDenCompC[2] > -5) recover()
logDenCompC[ dInfo$resCompPos ] <- logDenCompP[ dInfo$resCompPos ] <- resUpdateB$logDenCompC # update the result of the used density function
parUpdateDenC[blockArgs$dInfoPos,cu] <- TRUE # if density has been calculated, indicate that it is up to date
}
}else{ # not accepted
if( 0 != length(resUpdateB$xP) ) xP[ cu ] <- resUpdateB$xP # proposal
if( 0 != length(resUpdateB$logDenCompP) ) logDenCompP[ dInfo$resCompPos ] <- resUpdateB$logDenCompP # density results for proposal
# keep with previous intermediate
}
} #iBlock
##value<< list with components
resUpdateStep <- list( ##describe <<
xC=xC ##<< numeric vecotr: current accepted state
,logDenCompC=logDenCompC ##<< numeric vector: result components of fLogDen for current position
,intermediatesC= intermediatesC ##<< the list of intermediate results of current state
,parUpdateDenC=parUpdateDenC ##<< integer vector: logDensity that recently updated parameter at given index
,accepted=acceptedBlock ##<< numeric vector: acceptance of each block (0: not, 1: did, (0..1): DRstep)
,xP=xP ##<< numeric vector: result components of fLogDen for proposal
#,intermediate
,logDenCompP=logDenCompP
) ##enf<<
#}) # with
}
### \item{tempC}{global temperature}
### \item{TFix}{named numeric vector (comp): components with fixed Temperature}
### \item{posTFix}{integer vector (comp): =match(TFix, compNames): positions of TFix within comp provided for performance reasons}
### \item{ctrl$useMultiT}{boolean wheter to scale temperature for different data streams}
### \item{TProp}{numeric vector of length(resCompPos): proportions of temperature for different data streams}
.calcFLogDen <- function(
## calculating the logDensity and adjusting results
fLogDen ##<< logDensity function
, xP ##<< parameter vector (first argument to fLogDen)
, compInt ##<< internal components
, TiInt ##<< numeric vector (nCompInt): Temperature for internal components
, argsFLogDen ##<< further named arguments to fLogDen
, maxLogDen ##<< integer vector (nResComp): maximum logDensity (overfitting control, usually -1/2 nObs)
, fLogDenScale=1 ##<< multiplicator to convert logDen result (maybe a cost) to true logDensity
){
l <- try( argsFLogDen ); if( inherits(l,"try-error")) recover()
LpOverfit <- fLogDenScale * if(length(compInt)){
do.call( fLogDen, c(list(xP, compInt/fLogDenScale, TiInt), argsFLogDen) ) # evaluate logDen
}else{
do.call( fLogDen, c(list(xP), argsFLogDen) ) # evaluate logDen
#a$fLogDen( xProp, blockIndices=a$argsFLogDen$blockIndices, fModel=a$argsFLogDen$fModel, obs=a$argsFLogDen$obs, invCovar=a$argsFLogDen$invCovar, xval=a$argsFLogDen$xval )
}
Lp <- if( length(maxLogDen) ) pmin( LpOverfit, maxLogDen ) else LpOverfit
##value<< numeric vector of logDensity, corrected for logDenScale and maximum LogDensity, i.e. overfitting control
}
#undebug(updateBlockTwDEMC)
updateBlockTwDEMC <- function(
### Block update function by a Metropolis decision
xC ##<< numeric vector: current state that is used in density function of the block
,argsFUpdateBlock
### further argument provided for generating the update \describe{
### \item{upperParBounds}{ named numeric vector, upper limits of the parameters to generate, see \code{\link{twDEMCBlockInt}} }
### \item{lowerParBounds}{ named numeric vector, lower limits of the parameters to generate, see \code{\link{twDEMCBlockInt}} }
### \item{iPop}{ the population for which the sample is generated (just for debugging) }
### \item{intermediate}{ intermediate result for current state xC, see end of vignette on using two Densities }
### }
# # # \item{fCalcComponentTemp}{ functiont to calculate temperature of result components, (way of transporting calcComponentTemp to remote process) }
){
#updateBlockTwDEMC
#
##details<<
## Function arguments above are provided to all block update functions.
## If Metropolis update is used, argument \code{argsFUpdateBlock} constains the following additional entries \describe{
## \item{compPosInDen}{ positions of the dimensions in x that are updated in this block}
## \item{ctrl$DRgamma}{ if given and >0 delayed Rejection (DR) (Haario06) is applied:
## If proposal is rejected, then a second proposal is tried jumping only DRgamma distance along the proposal }
## \item{dInfo}{ list on used logDensity function static troughout inversion with items \describe{
## \item{argsFLogDen}{additional arguments to fLogDen and scalar factor applied to result of fLogDen}
## \item{posLogDenInt}{the matching positions of intResCompNames within the the results components that are handled internally}
## \item{maxLogDen}{ integer vector (nResComp): maximum logDensity (overfitting control, usually -1/2 nObs) }
## \item{fDiscrProp,argsFDiscrProp}{function and additional arguments applied to xProp, e.g. to round it to discrete values}
## }}
## \item{metInfoStep}{ list on current Metropolis update variable troughout inversion with items \describe{
## \item{step}{proposed jump}
## \item{rExtra}{extra portion in Metropolis decision because of selecting the jump}
## \item{logDenCompC}{numeric vector: former result of call to same fLogDen}
## \item{tempDenCompC}{numeric vector of length(logDenCompC): temperature for each density result component}
## }}
## }
#
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{.updateBlocksTwDEMC}
# # \code{\link{.updateBlocksTwDEMC}}
#attach( argsDEMCStep )
#stop(".doDEMCStep: stop to trace error in remote R invocation.")
# collect information into a single big shorthand list for shorter code
a <- c(argsFUpdateBlock, argsFUpdateBlock$dInfo, argsFUpdateBlock$metInfoStep )
a$argsFLogDen$intermediate <- NULL # intermediate refers to current state, not to new proposal
#recover()
#with( argsFUpdateBlock, { # with produces significant performance loss
boResFLogDenX <- (length(a$intResCompPos) > 0)
# LogDensity of accepted state
#La <- logDenCompC #logDensity components of accepted state
if( 0 == length(a$logDenCompC)) stop("updateBlockTwDEMC: encountered empty logDenCompC")
#assume that all is.finite(logDenCompC), make sure in twDEMCBlockInt
LaExt <- La <- a$logDenCompC
#logDenC <- sum(La)
#a$tempDenCompC <- a$fCalcComponentTemp( temp=a$tempC, TFix=a$TFix, TProp=a$TProp, useMultiT=a$ctrl$useMultiT,posTFix=a$posTFix)
#
accepted<-FALSE
xP <- xC
xP[a$compPosInDen] <- xP[a$compPosInDen] + a$step[a$compPosInDen]
#
#if( a$iPop==2 && (a$lowerParBounds != 10.8) ) recover()
boOutside <-
any( sapply( names(a$upperParBounds), function(pname){ xP[pname] >= a$upperParBounds[pname] })) ||
any( sapply( names(a$lowerParBounds), function(pname){ xP[pname] < a$lowerParBounds[pname] }))
#if(xProp[1] < 10.8) recover()
#sapply( names(a$lowerParBounds), function(pname){ pname })
#sapply( names(a$lowerParBounds), function(pname){ xP[pname] })
#sapply( names(a$lowerParBounds), function(pname){ a$lowerParBounds[pname] })
if( boOutside ){
# bp(11)
# if it is still outside (maybe opposite border) reject step and give -Inf as logDenResult
#logDenP=logAlpha10=-Inf
logAlpha10=-Inf #logAlpha10 is log of the initial acceptance ratio for DR step (-Inf no chance of acceptance)
Lp <- a$logDenCompC #results for the proposal
Lp[] <- -Inf
}else{
# discrtize proposal
if( is.function(a$fDiscrProp)) xP = do.call(a$fDiscrProp,xP,a$argsFDiscrProp, quote=TRUE)
Lp <- .calcFLogDen( a$fLogDen, xP, compInt=a$logDenCompC[a$intResCompPos], TiInt = a$tempDenCompC[a$intResCompPos]
,argsFLogDen = a$argsFLogDen
,maxLogDen = a$maxLogDen
)
#take care that the result has always the same sames, even when if fails
#if( 0==length(names(res)))
# stop("encountered result of fLogDen without names")
#if( !identical(names(logDenCompC),names(res)))
# stop("encountered result with different names")
#strip attributes other than names, else twDynamicClusterApplyDep fails with big data chunks
#attributes(Lp) <- list(names=names(Lp))
#logDenP=-Inf
#make sure Lp, La have the same order and legnth
#if( !identical( names(Lp), names(La)) ) stop(".doDEMCStep: logDenCompC must contain the same components and the order of result of fLogDen." )
if( all(is.finite(Lp))){
#logDenP <- sum(Lp)
##details<< \describe{\item{internal Metropolis step}{
## if posLogDenInt is given, then these components of result of fLogDen are handled
## internally. Hence, for Metropolis step here operates only on other remaining components.
## }}
posTExt <- setdiff( seq_along(a$resCompPos), a$intResCompPos ) #externally handled components
nExt <- length(posTExt)
#posTFixExt <- setdiff(posTFix,posLogDenInt) #externally handled components with fixed temperature
#posTVarExt <- setdiff(seq_along(Lp), c(posTFix,posLogDenInt)) #externally handled componetns with variable temperature
#nFixExt <- length(posTFixExt)
#nVarExt <- length(posTVarExt)
#nExt <- nFixExt + nVarExt
#
#Metropolis step
#logr = (logDenPropExt+rExtra - logDenXExt) / tempC
logrDS10 <- (Lp[posTExt]-La[posTExt])/a$tempDenCompC[posTExt]
logAlpha10 <- a$rExtra + sum(logrDS10)
accepted <- if( is.finite(logAlpha10) && (logAlpha10 > log(runif(1)) ) ) 1 else 0
if(accepted != 0){
xC <- xP
#logDenC <- logDenP
a$logDenCompC <- Lp
}
}else logAlpha10 <- -Inf
} # end check outside parBounds
# DR-step
if(!accepted && !is.null(a$ctrl$DRgamma) && (a$ctrl$DRgamma > 0) &&
( boOutside || (!is.null(a$ctrl$minPCompAcceptTempDecr) && (a$pAccept < 1.2*a$ctrl$minPCompAcceptTempDecr)))
) {
#----- delayed rejection (DR) step
# only if across parBoundEdge or acceptance rate drops below 1.2*minAcceptrate
# repeat all above with delayed rejection (DR) step, only adjust DRfac after calculating Lp
Lp1 <- Lp
xProp1 <- xP # former proposal
xP <- xC
xP[a$compPosInDen] <- xP[a$compPosInDen] + a$ctrl$DRgamma*a$step[a$compPosInDen]
boOutside <-
any( sapply( names(a$upperParBounds), function(pname){ xP[pname] > a$upperParBounds[pname] })) ||
any( sapply( names(a$lowerParBounds), function(pname){ xP[pname] < a$lowerParBounds[pname] }))
if( !boOutside ){
if( is.function(a$fDiscrProp)) xP = do.call(a$fDiscrProp,xP,a$argsFDiscrProp, quote=TRUE)
Lp <- .calcFLogDen( a$fLogDen, xP, compInt=a$logDenCompC[a$intResCompPos], TiInt = a$tempDenCompC[a$intResCompPos]
,argsFLogDen = a$argsFLogDen
,maxLogDen = a$maxLogDen
)
#take care that the result has always the same sames, even when if fails
#if( 0==length(names(res)))
# stop("encountered result of fLogDen without names")
#if( !identical(names(logDenCompC),names(res)))
# stop("encountered result with different names")
#strip attributes other than names, else twDynamicClusterApplyDep fails with big data chunks
attributes(Lp) <- list(names=names(Lp))
#logDenP=-Inf
#make sure Lp, La have the same order and legnth
#if( !identical( names(Lp), names(La)) ) stop(".doDEMCStep: logDenCompC must contain the same components and the order of result of fLogDen." )
if( all(is.finite(Lp))){
#logDenP <- sum(Lp)
posTExt <- setdiff( seq_along(a$tempDenCompC), a$intResCompPos ) #externally handled components
nExt <- length(posTExt)
#posTFixExt <- setdiff(posTFix,posLogDenInt) #externally handled components with fixed temperature
#posTVarExt <- setdiff(seq_along(Lp), c(posTFix,posLogDenInt)) #externally handled componetns with variable temperature
#nFixExt <- length(posTFixExt)
#nVarExt <- length(posTVarExt)
#nExt <- nFixExt + nVarExt
#Metropolis step
#logr = (logDenPropExt+rExtra - logDenXExt) / tempC
logrDS20 <- (Lp[posTExt]-La[posTExt])/a$tempDenCompC[posTExt]
logAlpha20 <- a$rExtra + sum(logrDS20)
#--- here correct with first stage DR factor (1-alpha21)/(1-alpha10) with meaning 0:accepted 1:first proposal 2:second proposal
logrDS21 <- (Lp1[posTExt]-Lp[posTExt])/a$tempDenCompC[posTExt]
logAlpha21 <- sum(logrDS21)
logAlpha2 <- suppressWarnings( logAlpha20 +log(1-exp(logAlpha21)) -log(1-exp(logAlpha10)) ) # log and exp may produce NaNs
accepted <- if( is.finite(logAlpha2) && ( logAlpha2 > log(runif(1)) ) ) a$ctrl$DRgamma else 0
if(accepted != 0){
xC <- xP
#logDenC <- logDenP
a$logDenCompC <- Lp # also transferres the "intermediate" attribute
}
}
} # end !boOutside in DR step
} # end DR step
#will invoke prevRes[c("x", "logDenCompC", "logDenAcc")]
if(!is.numeric(xC) | !is.numeric(a$logDenCompC) ) #| !is.numeric(logDenC) )
stop(".twDEMCBlocks: x, logDenAcc and logDenCompC must be numeric")
##value<< list with components
list( ##describe<<
accepted=accepted ##<< boolean scalar: if step was accepted
, xC=xC[a$compPosInDen] ##<< numeric vector: components of position in parameter space that are being updated (if accepted then the same as xP)
, logDenCompC=a$logDenCompC ##<< numeric vector: result components of fLogDen for current position (if accepted then the same as logDenCompP)
, intermediate=attr(a$logDenCompC,"intermediate") ##<< intermediate result of current state (if provided by logDensity)
#, logDenC =logDenC ##<< numeric vector: summed fLogDen for current accepted position
, xP=xP[a$compPosInDen] ##<< numeric vector: components of proposal that are being
, logDenCompP=Lp ##<< numeric vector: result components of fLogDen for proposal, may inlcude attr "intermediate"
#, logDenP=logDenP ##<< numeric vector: summed fLogDen for proposal
) ##end<<
#})
}
R.methodsS3::setMethodS3("twDEMCBlock","array", function(
### Initialize \code{\link{twDEMCBlockInt}} by array of initial population and remove those generations from results afterwards
x ##<< initial population: a numeric array (M0 x d x nChain) see details in \code{\link{twDEMCBlockInt}}
,... ##<< further arguments to \code{\link{twDEMCBlockInt}}
,nPop = 1 ##<< number of populations in x
#,TProp = NULL ##<< numeric vector (nResComp): proportions of temperature for result components: equal for all populations
, X=NULL ##<< initial state (nParm x nChain)
, logDenCompX = NULL ##<< numeric matrix (nResComp x nChains) initial state
, upperParBounds = vector("list",nPop)
### list of named numeric vectors: giving upper parameter bounds for each population
### for exploring subspaces of the limiting distribution, see details
### , Alternatively a single numeric vector can be supplied, which is replicated for each population.
, lowerParBounds = vector("list",nPop)
### similar to upperParBounds
){
#twDEMCBlock.array
#if( 1 == length(T0) ) T0 <- rep(T0,nPop)
#if( 1 == length(TEnd) ) TEnd <- rep(TEnd,nPop)
if( is.numeric(upperParBounds) ) upperParBounds <- lapply(1:nPop, function(iPop) upperParBounds )
if( is.numeric(lowerParBounds) ) lowerParBounds <- lapply(1:nPop, function(iPop) lowerParBounds )
nChains <- dim(x)[3]
nChainPop <- nChains %/% nPop
isX <- (0 != length(X))
isLogDenCompX <- (0 != length(logDenCompX))
#isTProp <- (0 != length(TProp))
#if( isTProp && !is.matrix(TProp) ) TProp <- matrix( TProp, nrow=length(TProp), ncol=nPop, dimnames=list(resComp=names(TProp),pop=NULL) )
pops <- lapply( 1:nPop, function(iPop){
chainsPopI <- (iPop-1)*nChainPop + 1:nChainPop
pop <- list(
parms=x[,,chainsPopI, drop=FALSE]
#,T0=T0[iPop]
#,TEnd=TEnd[iPop]
,upperParBounds=upperParBounds[[iPop]]
,lowerParBounds=lowerParBounds[[iPop]]
)
if( isX ) pop$X <- X[,chainsPopI , drop=FALSE]
if( isLogDenCompX ) pop$logDenCompX <- logDenCompX[,chainsPopI ,drop=FALSE]
#if( isTProp ) pop$TProp <- TProp[,iPop ,drop=TRUE]
#pop$spaceInd <- spacePop[iPop]
pop
})
#res <- twDEMCBlockInt(pops)
res <- twDEMCBlockInt(pops=pops,...)
#.dots <- list(...)
#do.call( twDEMCBlockInt, c(list(pops=pops),.dots))
#resc <- concatPops.twDEMCPops(res) # all have the same length, so allow concatenate population results
})
attr(twDEMCBlock.array,"ex") <- function(){
data(twLinreg1)
attach( twLinreg1 )
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nPop = 2
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
dim(Zinit)
.nGen=100
#nGen=3
#mtrace(twDEMC.array)
#mtrace(.updateIntervalTwDEMCPar)
#mtrace(twDEMCBlockInt)
tmp1 <- tmp <- twDEMCBlock( Zinit, nPop=.nPop
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
,nGen=.nGen
)
plot( as.mcmc.list(tmp), smooth=FALSE )
tmp2 <- tmp <- twDEMCBlock( tmp1, nGen=200 )
str(tmp)
}
.concatTwDEMCRuns <- function(
x ##<< initial run (of class twDEMCPops)
,res ##<< extended run
, doRecordProposals=FALSE ##<< if TRUE then an array of each proposal together with the results of fLogDen are recorded and returned in component Y
){
##details<<
## initial populations are appended by results
# iPop = length(x$pops)
for( iPop in 1:length(x$pops) ){
xPop <- x$pops[[iPop]]
resPop <- res$pops[[iPop]]
if( nrow(resPop$parms) < 2){
# no samples to append
res$pops[[iPop]] <- xPop
}else{
res$pops[[iPop]]$parms <- abind(xPop$parms, resPop$parms[-1,, ,drop=FALSE], along=1)
res$pops[[iPop]]$temp <- abind(xPop$temp, resPop$temp[-1, ,drop=FALSE], along=1 )
res$pops[[iPop]]$pAccept <- abind(xPop$pAccept, resPop$pAccept[-1,, ,drop=FALSE], along=1)
res$pops[[iPop]]$resLogDen <- abind( xPop$resLogDen, resPop$resLogDen[-1,, ,drop=FALSE], along=1)
res$pops[[iPop]]$logDen <- abind( xPop$logDen, resPop$logDen[-1,, ,drop=FALSE], along=1)
nrY <- nrow(resPop$Y)
if( doRecordProposals || (nrY < 128) ){
# if new Y has less than 128 rows append previous Y
res$pops[[iPop]]$Y <- abind( xPop$Y, resPop$Y, along=1)
nrY <- nrow(res$pops[[iPop]]$Y) # rows did change
if( !doRecordProposals && (nrY > 1) )
res$pops[[iPop]]$Y <- res$pops[[iPop]]$Y[ 1:min(128,nrY),, ,drop=FALSE] # cut to 128 cases
}
} # if res has samples
}
# copy other items present in x to result
namesAdd <- names(x)[ !(names(x) %in% names(res)) ]
res[ namesAdd ] <- x[ namesAdd ]
res
}
R.methodsS3::setMethodS3("twDEMCBlock","twDEMCPops", function(
### initialize \code{\link{twDEMCBlockInt}} by former run and append results to former run
x ##<< list of class twDEMCPops, result of \code{\link{twDEMCBlockInt}}
,... ##<< further arguments to \code{\link{twDEMCBlockInt}}
, TEnd=numeric(0) ##<< numeric vector (nResComp) of end temperatures for each result component, if not given stays at current temperature, if in addition TSpec is specified, it has no effect
, TEnd0 ##<< numeric scalar: alternative way of specifying end tempertaure: by baseTemperature for all data streams
, TStart=numeric(0) ##<< numeric vector (nResComp): specifing starting temperature, if neigher TSpec, TStart, or TStart0 are given stays at the current temperature
, TStart0 ##<< numeric scalar: specifing starting temperature by baseTemperature for all data streams
, TFix ##<< see \code{\link{twDEMCSACont}}, argument \code{ctrlT$TFix}. Must be given if specifying TEnd0 or TStart0
, nObs ##<< see \code{\link{twDEMCSACont}}, argument must be given if specifying TEnd0 or TStart0
, doRecordProposals=FALSE ##<< if TRUE then an array of each proposal together with the results of fLogDen are recorded and returned in component Y
, extendRun=TRUE ##<< if set to FALSE then only the new samples are returned
){
#twDEMCBlock.twDEMCPops
##details<<
## pops, dInfos, and blocks are reused from x or overwritten by arguments
nSamplePop <- getNSamples(x) # used in several cases
#
.dots <- argsList <- list(...)
argsList$pops <- x$pops
argsList$doRecordProposals <- doRecordProposals
if( 0 == length(.dots$dInfos) )
argsList$dInfos <- x$dInfos
if( 0 == length(.dots$blocks) )
argsList$blocks <- x$blocks
#for( iPop in seq_along(x$pops) ){
# if( 0 == length(argsList$pops[[iPop]]$T0) )
# argsList$pops[[iPop]]$T0 <-
# x$pops[[iPop]]$temp[ nSamplePop[iPop] ]
#}
if( 0 == length(argsList$TSpec) ){
##details<< \describe{\item{TSpec}{
## If TSpec is not explicitly given as an array for \code{\link{twDEMCBlockInt}} it is constructed as follows:
## , If TEnd is given use it to explicitely define end Temperature for all data streams.
## , If TEnd0 is given calculate TEnd for all data streams by scaling for number of observations (see \code{\link{calcStreamTemp}}).
## , If TStart0 is given calcualte TCurr for all data streams
## , else use current Temperature of streams in x as starting temperature
## }}
if( (!missing(TStart0)) || (!missing(TEnd0)) ){
if( missing(nObs) || missing(TFix) ) stop("twDEMCBlock.twDEMCPops: if TStart0 or TENd0 are specified, one needs to specify nObs and TFix too.")
resCompNames <- colnames(x$pops[[1]]$resLogDen)
nObsA <- nObs[ resCompNames ]
TFix <- completeResCompVec( TFix, resCompNames )
iFixTemp <- which( is.finite(TFix) )
}
TCurr <- if( (!missing(TStart))){
TStart
}else if( (!missing(TStart0))){
calcStreamTemp(TStart0, nObsA , TFix, iFixTemp )
}else{
x$pops[[1]]$temp[ nrow(x$pops[[1]]$temp), ]
}
if( 0 == length(TEnd) )
TEnd <- if( !missing(TEnd0)){
calcStreamTemp(TEnd0, nObsA , TFix, iFixTemp )
} else TCurr
if( 1 == length(TEnd) ) TEnd <- rep(TEnd, length(TCurr) )
argsList$TSpec <- cbind( T0=TCurr, TEnd=TEnd )
}
#
if( 0 == length(.dots$controlTwDEMC) )
argsList$controlTwDEMC <- list()
if( 0 == length(argsList$controlTwDEMC$initialAcceptanceRate) ){
ar <- abind( lapply( seq_along(x$pops), function(iPop){
adrop( x$pops[[iPop]]$pAccept[ nSamplePop[iPop],, ,drop=FALSE ],1)
}),along=2)
argsList$controlTwDEMC$initialAcceptanceRate <- ar
}
res <- res0 <- do.call( twDEMCBlockInt, argsList )
if( extendRun ) res <- .concatTwDEMCRuns(x,res0,doRecordProposals=doRecordProposals)
res
})
attr(twDEMCBlock.twDEMCPops,"ex") <- function(){
data(twdemcEx1) # previous run of twDEMCBlock
class(twdemcEx1)
twdemcEx1$thin # thinning interval
(nGen0 <- getNGen(twdemcEx1)) # number of generations
# extend by 16 generations
nGen <- 16
#mtrace(twDEMCBlock.twDEMCPops)
res <- twDEMCBlock( twdemcEx1, nGen=nGen )
identical( nGen0+nGen, getNGen(res) )
plot( as.mcmc.list(res), smooth=FALSE )
}
R.methodsS3::setMethodS3("twDEMCBlock","twDEMC", function(
### initialize \code{\link{twDEMCBlockInt}} by former run and append results to former run
x, ##<< list of class twDEMC, result of \code{\link{twDEMCBlockInt}}
... ##<< further arguments to \code{\link{twDEMCBlockInt}}
){
#twDEMC.twDEMC
.dots <- list(...)
argsList <- list(x=x$parms)
M0 <- nrow(x$rLogDen)
#extract X, logDenX, and logDenCompX from parms, but only if not given with \dots
if( is.null(.dots$X) ) #use is.null, because if provided zero length vector, we want to use it
argsList$X <- adrop(x$parms[M0,,,drop=FALSE],2)
if( is.null(.dots$logDenX) )
argsList$logDenX <- x$rLogDen[M0,,drop=TRUE]
if( is.null(.dots$logDenCompX) )
argsList$logDenCompX <- adrop(x$logDenComp[M0,,,drop=FALSE],2)
if( is.null(.dots$nPop) )
argsList$nPop <- getNPops.twDEMC(x)
res <- do.call( twDEMCBlock.array, c(argsList,.dots))
res$rLogDen[1:M0,] <- x$rLogDen
res$logDenComp[,1:M0,] <- x$logDenComp
res$pAccept[1:M0,] <- x$pAccept
res$temp[1:M0,] <- x$temp
res
### parms appended with the further generations of \code{\link{twDEMCBlockInt}}
})
#mtrace(twDEMC.twDEMC)
#mtrace(twDEMCBlockInt)
#twUtest(twDEMC,"test.ZinittwDEMC")
Try the twDEMC package in your browser
Any scripts or data that you put into this service are public.
twDEMC documentation built on May 2, 2019, 5:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
twDEMCBlockInt <- function(
### Parallelized differential Evolution Markov Chain with blocked parameter update
pops = list( list( ##<< list of population infos for each population, each info is a list with components
##describe<<
parms=list() ##<< list of matrices (nState x nParm x nChain) initial states for each population see details and \code{\link{initZtwDEMCNormal}}.
, nGen = 12 ##<< number of generations, i.e. steps to proceed
, TProp=NULL ##<< numeric vector (nResComp) temperature proportions of result components.
## It can also be given as character vector with names of result components, however, be aware that htis fails if several logDen may return components of the same name
, X=NULL ##<< numeric matrix (nParm x nChainPop) initial state
, logDenCompX=NULL ##<< numeric matrix (nComp x nChain): logDen components of initial state X, see details
, spaceInd = 1 ##<< the space replicate that this population belongs to
, upperParBounds = numeric(0) ##<< named numeric vectors: giving upper parameter bounds: lowerBound < par <= upperBound.
## For exploring subspaces of the limiting distribution, see details
, lowerParBounds = numeric(0) ##<< similar to upperParBounds: sample > bound
, splits=numeric(0) ##<< named numeric vector of splitting points, used to remerge divided subspaces
))
##end<<
,dInfos = list( list( ##<< named list of used density functions. Each entry is a list with components
##describe<<
fLogDen=NULL ##<< \code{function(theta, ...)} calculates a vector of logDensities
## corresponding to different data streams of parameter vector theta
## , or \code{function(theta, logDenCompX, metropolisStepTemp, ...)}
, argsFLogDen=list() ##<< further arguments passed to fLogDen
, compPosDen #=dInfos[[1]]$compPosDen ##<< names or index of the parameters used in with this density, defaults to all parameters
, intResComp=vector("integer",0) ##<< integer or character vector: indices or names of results components of fLogDen
## that are used for internal Metropolis decisions
#, TFix = vector("numeric",0) ##<< named numeric vector with Temperature for result components that have fixed temperature
, maxLogDen=c() ##<< integer vector (nRespComp): maximum logDensity (=usually -1/2* number of observations) for each result component, used in control of overfitting
, fDiscrProp=NULL ##<< function applied to proposal, e.g. to round proposals to to discrete possible values function(theta,...)
, argsFDiscrProp=list() ##<< further arguments to fDiscrProp
))
##end<<
, blocks = list( list( ##<< list of parameter blocks, each block is a list with entries
##describe<<
compPos = dInfos[[1]]$compPosDen ##<< names or index of the parameter components to be updated
, fUpdateBlock = NULL ##<< function to update the parameters.
## It must return a list with first named components \code{accepted}, \code{xC}, and optionally \code{intermediate}
## as described in \code{\link{updateBlockTwDEMC}}. If not specified (default), then \code{\link{updateBlockTwDEMC}} is used and \code{dInfoPos} must be specified
, argsFUpdate=list() ##<< further arguments passed to fUpdate
, intermediateId=character(0) ##<< string: identifier in list of intermediate results that are shared between block udpate functions. See details.
, dInfoPos=1 ##<< name or position to \code{fLogDenInfo} if using Metropolis update.
## Several blocks may share the same density but update different parameters
))
##end<<
, stateCurrent = list(
,xC=c() ##<< numeric vector: current accepted state
,logDenCompC=c() ##<< numeric vector: result components of all fLogDens for current position
,parUpdateDenC=c() ##<< integer vector: logDensity that recently updated parameter at given index
## to handle first steps in Multi-step Metropolis decision internally.
## See details.
)
, TSpec = matrix(1, ncol=2, nrow=1, dimnames=list(NULL, c("T0","TEnd")) )
### numeric matrix (nResComp x 2): specifying Initial and End Temperature of each fLogDen result component.
### If only one row is specified, the Temperature is taken for all result components
, m0 = numeric(0) ##<< scalar integer: number of samples in initial population, if length==0 will be calculated from number of chains and number of parameters
, nGen = integer(0) ##<< scalar integer: number of generations, if given overwrites \code{pops[[i]]$nGen}
, spacePop = 1:nPop ##<< the space replicate that each population belongs to. By default assume only one population per space, overridden by entry in pops
, controlTwDEMC = list() ##<< control parameters influencing the update and the convergens of the chains (see details)
, debugSequential=FALSE ##<< if TRUE apply is used instead of sfApply, for easier debugging
, remoteDumpfileBasename=NULL ##<< the basename of a dumpfile that is created on error on remote process (see example section)
, doRecordProposals=FALSE ##<< if TRUE then an array of each proposal together with the results of fLogDen are recorded and returned in component Y
, progressOutput="." ##<< output after each thinning interval
){
##details<<
## This is the central method for applying a Differential Evolution Markov Chain given a function of
## an unnormalized probability density.
## It is invoked usually by (\code{\link{twDEMCBlock.array}} or \code{\link{twDEMCBlock.twDEMCPops}})
##seealso<<
## \code{\link{calcDEMCTemp}}
## \code{\link{logDenGaussian}}
# # \code{\link{.generateXPropThin}}
# # \code{\link{.doDEMCSteps}}
# # \code{\link{.doDEMCStep}}
## #details<<
## # further elements in pop that are calculated are \itemize{
## # \item resCompPos: position of the result components in the concatenation across all density functions
## # }
##details<< \describe{ \item{recognized control parameters in \code{controlTwDEMC}: }{
##describe<<
ctrl = list(
thin = 4 ##<< thinning interval, Default: 4
,F = 2.38 ##<< related to multiplicative error (F2=F/sqrt(2*Npar), see eps.mult, Default: 2.38
,pSnooker= 0.1 ##<< probability of a snooker update (others parallel updates), Default 0.1
,pGamma1 = 0.1 ##<< probability of jumping to state of another chain (different modes), Default 0.1
,epsMult =0.2 ##<< >0 gives d-dimensional variation around gamma.
## It adds scaled uncorrelated noise to the proposal. Default: 0.2 ,
## Its advantage over eps.add is
## that its effect scales with the differences of vectors in the population whereas eps.add does not.
## if the variance of a dimensions is close to 0, eps.mult gives smaller changes. ,
## A uniformly distributed error, i.e. F2*runif(1+-epsMult*prop) multiplied to difference vector
## from parallel update
,epsAdd = 0 ##<< >0 is needed to ensure that all positions in the space can be reached.
## For targets without gaps, it can set small or even to 0. Default 0 ,
## sd of normally distributed error added to proposal by parallel or snooker update.
,pAcceptWindowWidth = 8 ##<< number of thinning intervals back over which the acceptance rate is calculated Default 8
,probUpDir=0.5 ##<< probability of direction between two states of increasing Density, Default 0.5
## , Increasing this during burin may accelerate convergence
,initialAcceptanceRate=0.25 ##<< numeric matrix (nBlock x nChains) initially assumed acceptance rate. Default 0.25
## , Used to calculate the number of generations backwards to sample from
,DRgamma=0 ##<< factor for reducing step length [0..1) in delayed rejection step, 0 means no DR step, Default 0
,minPCompAcceptTempDecr=0.15 ##<< if acceptance rate drops below minPCompAcceptTempDecr+0.02 times this level, employ delayed rejection (DR), Default 0.15
,pIndStep = 1.5 ##<< assume state to be independent, after on average about those number of accepted steps, Default 1.5
,nPastGen = 10 ##<< factor for determining the number of recent past states to sample during burnin. Default 10
## It is multiplied by the number of parameters. Past generations are calculated
## by deviding by the number of chains per population
#,useLoadBalancing=FALSE ##<< if set to TRUE, submits each step separaetely to a free node, else each node gets an entire chain process
#,freeMasterNode=FALSE ##<< if set to TRUE, no job is submitted to first node, so that this node can dispatch jobs without waiting,, see \code{\link[twSnowfall]{sfFArgsApplyDep}}
,returnIntermediates=TRUE ##<< set to FALSE if intermediate result is large and transferring it between slaves slows down calculation
#,useMultiT = FALSE ##<< whether to downscale Temperature of result components during burnin
#moved to block,TFix = vector("numeric",0) ##<< named numeric vector: result components for which temperature shoudl not change
,useConditionalProposal = FALSE ##<< if set to TRUE, poposal steps are generated conditional on the state of the other blocks (experimental)
#,controlOverfittingMinNObs = Inf ##<< scalar positive integer: minimum number of observations where bias is controlled. Need at least about 20 observations to work. Default Inf says no bias correction.
)
##details<< }}
ctrl[names(controlTwDEMC)] <- controlTwDEMC
#
#-- dimensions of pops
nPop <- length(pops)
iPops <- 1:nPop
pops <- if( 0 == length(nGen) ){
lapply( iPops, function(iPop){ .completePopDescription( pops[[iPop]], spaceInd=spacePop[iPop] )}) # fill in default values for missing entries
}else{
if( length(nGen)==1) nGen=rep(nGen[1],nPop)
if( length(nGen) != nPop ) stop("twDEMCBlockInt: lenght of nGen must equal the number of populations.")
lapply( iPops, function(iPop){ .completePopDescription( pops[[iPop]], nGen=nGen[iPop], spaceInd=spacePop[iPop])}) # fill in default values for missing entries
}
ZinitPops <- lapply(pops,"[[","parms")
parNames <- colnames(ZinitPops[[1]])
nParm <- ncol(ZinitPops[[1]])
nChainPop <- dim(ZinitPops[[1]])[3]
if( nChainPop < 2) stop(paste("Initialized with only ",nChainPop," per population, but must have at least two chains per population (recommened is at least 4)"))
nChain <- nChainPop*nPop
popChain <- rep(1:nPop,each=nChainPop) # population for chain at given index
chainsPop <- lapply( iPops, function(iPop){ (iPop-1)*nChainPop+(1:nChainPop)}) # chains for given population
#
# check that have enough values for initial populations of each population
if( 0 == length(m0) ) m0 = calcM0twDEMC(nParm,nChainPop) # number of samples in initial population
M0Pops <- sapply( ZinitPops, nrow )
if( any(M0Pops < 0.8*m0) ) warning(paste("twDEMCBlockInt: too few initial cases for populations",paste(which(M0Pops < m0),collapse=",")) )
#
# determine number of thinning intervals and create result arrays
nGenPops <- sapply( pops, "[[", "nGen")
nThinnedGenPops = sapply( nGenPops %/% ctrl$thin, max, 0 ) #number of thinning intervals
nGenPops = nThinnedGenPops * ctrl$thin #number of generations in this batch (do not need to move behind last thinning interval)
Nz <- M0Pops +(nThinnedGenPops) #number of rows in Z after batch
XPops <- lapply( pops, "[[", "X" )
XChains <- do.call(cbind,XPops)
logDenCompXPops <- lapply( pops, "[[", "logDenCompX" ) # see initial states for initializing missing entries
upperParBoundsPop <- lapply( pops, "[[", "upperParBounds" )
lowerParBoundsPop <- lapply( pops, "[[", "lowerParBounds" )
#
#-- dimensions of dInfo
nDen <- length(dInfos)
iDens <- 1:nDen
if( is.null(names(dInfos))) names(dInfos) <- paste("dInfo",iDens,sep="") # make up names if not given
posInfosEmptyString <- which(names(dInfos) == "")
names(dInfos)[posInfosEmptyString] <- paste("dInfo",posInfosEmptyString,sep="")
#compPosDens <- lapply(dInfos, "[[", "compPosDen")
#calculate each logDen once to obtain result component names required for proper dInfos initialization
#dInfo <- dInfos[[1]]
resCompNamesList <- lapply(iDens, function(iInfo){
dInfo <- dInfos[[iInfo]]
.resFLogDen <- .calcFLogDen(dInfo$fLogDen, xP=XChains[,1], compInt=c(), TiInt=c()
, argsFLogDen=dInfo$argsFLogDen, maxLogDen=c() )
resCompNames <- .getResFLogDenNames( .resFLogDen, logDenName=names(dInfos)[iInfo] )
})
for( i in seq_along(dInfos)){
dInfos[[i]] <- .completeDInfo( dInfos[[i]], parNames=parNames, resCompNames = resCompNamesList[[i]] )
}
#
#-- calculate positions in resCompNames of results of different blocks
resCompNamesU <- lapply( dInfos, "[[", "resCompNames" )
resCompNamesUFlat <- as.vector(do.call(c, resCompNamesU)) # concatenate result names
nResComp <- length(resCompNamesUFlat)
resCompPosPops <- { # list: for each dInfo: position of resCompNames in resCompNamesUFlat
tmp.length <- sapply(resCompNamesU, length)
tmp.end <- cumsum(tmp.length)
tmp.i <- vector("list",nDen)
for( iDen in iDens ){
tmp.i[[iDen]] <- dInfos[[iDen]]$resCompPos <- tmp.end[iDen]+1-(tmp.length[iDen]:1)
}
tmp.i
}
#
#-- dimensions of blocks
#mtrace(.checkBlock)
#recover()
blocks <- lapply( blocks, .completeBlockInfo, dInfos=dInfos, parNames=parNames)
nBlock <- length(blocks)
iBlocks <- 1:length(blocks)
#block <- blocks[[1]]
dInfoPosBlock <- sapply(blocks, "[[", "dInfoPos")
# following refers to position in parNames
compPosBlock <- lapply( blocks, "[[", "compPos" ) # already transformed to positions in .checkBlock
# following refers to position in dInfo$compPosDen
compPosInDenBlock <- lapply( blocks, "[[", "compPosInDen" )
# make sure that all parameters are updated at least once
if( 0==length(compPosBlock) || (length(sBlock <- sort(unique(do.call(c, compPosBlock)))) != nParm) || !all.equal( sBlock, 1:nParm, check.attributes = FALSE) )
stop(paste("each parameter (columns of parms) must be updated in at least one block. nParm=",nParm,"updatedVars=",paste(sBlock,collapse=",")) )
#
#-- initialize further parameters to parallel and snooker steps
# that depend on pops (nParm, nChainPop)
ctrl$Npar12 =(nParm - 1)/2 # factor for Metropolis ratio DE Snooker update
ctrl$F2 = ctrl$F/sqrt(2*nParm)
ctrl$F1 = 1
# terBraak report that may not sample the distribution, if not using the full past
# but together with decreasing temperature acceptance rate drops very low if using the full past
# hence constrain to recent past during burnin
ctrl$gInd <- ctrl$nPastGen*nParm/nChainPop #number of independent generations to sample from, similarly to M0 (usually ctrl$nPastGen=10)
#
#-- evaluate logDen of initial states and get names of the result components (dInfo$resCompNames)
##details<< \describe{ \item{Initial state: \code{X}}{
## If initial state X is not specified, the last column (generation) of Z is utilized.
## If in addition to X, logDenX is specified as a numeric vector, the fLogDen will not be avaluated for the initial state of the chains.
## All the results of fLogDen for the initial state must be finite.
## }}
# find those chains where no logDenCompX is given, or not all components are finite
popsMissingResX <- ( 0 == sapply( logDenCompXPops, length) )
iMissingPops <- which(popsMissingResX)
iNonMissingPops <- which(!popsMissingResX)
missingPopInfo <- do.call( rbind, lapply( iMissingPops, function(iPop){
cbind(iPop=iPop,iChainInPop=1:nChainPop)
} ))
missingPopInfo <- rbind( missingPopInfo, do.call( rbind, lapply( iNonMissingPops, function(iPop){
#all( is.finite(logDenCompXPops[[iPop]][,1]))
#chainIsAllFinite <- apply( logDenCompXPops[[iPop]], 2, all, is.finite )
chainIsAllFinite <- apply( logDenCompXPops[[iPop]], 2, function(cChain){ all(is.finite(cChain))} )
iChainInPop=which(!chainIsAllFinite)
if( 0 != length(iChainInPop) )
cbind(iPop=iPop,iChainInPop=iChainInPop)
else
cbind(iPop=iPop, iChainInPop=NA_integer_ )[FALSE,]
})))
if( 0 < nrow(missingPopInfo) ){
# evaluate logDen for those
missingPopInfo <- cbind(missingPopInfo, iChain=(missingPopInfo[,"iPop"]-1)*nChainPop + missingPopInfo[,"iChainInPop"])
XChainsMissing <- XChains[,missingPopInfo[,"iChain"] ,drop=FALSE] # initial states of the missings
logDenCompDen <- vector("list", nDen )
#XXTODO calculate using intermediate
#XXTODO on restart provide existing resLogDen
logDenRes <- twCalcLogDensPar(
### Invokes all fLogDens with proposal in a parallel way, taking care of intermediates between densities
dInfos ##<< list describing the logDensities, see \code{\link{twDEMCBlockInt}}
,xProp=t(XChainsMissing[,,drop=FALSE])
, debugSequential=debugSequential
, remoteDumpfileBasename=remoteDumpfileBasename
)
if( any(!is.finite(logDenRes$logDen)) )
stop(paste("twDEMCBlockInt: non-finite logDensity of starting value ",sep=""))
#iDen <- 1
for( iDen in seq_along(dInfos) ){
#calculated further down, for all dInfos, not just the missing ones: dInfos[[iDen]]$resCompPos <-
.resCompPosI <- which(logDenRes$logDenCompPos==iDen)
logDenCompDen[[iDen]] <- .logDenCompI <- logDenRes$logDenComp[,.resCompPosI]
}
logDenCompXMiss <- logDenRes$logDenComp
# replace statesin initial history for proposal generation with infinite logDen by
# ones with finite states from other chains/populations
#jPop <- 1
for( jPop in seq_along(iMissingPops) ){
iiPop <- iMissingPops[jPop]
iInfo <- which( missingPopInfo[,"iPop"] == iiPop )
logDenCompXPops[[ iMissingPops[jPop] ]] <-
t(logDenCompXMiss[iInfo,])
}
for( jPop in seq_along(iNonMissingPops) ){
iiPop <- iNonMissingPops[jPop]
iInfo <- which( missingPopInfo[,"iPop"] == iiPop )
logDenCompXPops[[iiPop]][,missingPopInfo$iChainInPop[iInfo] ] <-
t(logDenCompXMiss[iInfo,])
}
}
#-- do overfitting control on the intial logDenComp
maxLogDen <- structure( numeric(length(resCompNamesUFlat)), names=resCompNamesUFlat)
#dInfo <- dInfos[[1]]
for( dInfo in dInfos){
maxLogDen[dInfo$resCompNames] <- dInfo$maxLogDen
}
#logDenCompXPop <- logDenCompXPops[[2]]
tmp <- logDenCompXPops
logDenCompXPops <- lapply(logDenCompXPops, function(logDenCompXPop){
logDenCompXPop[] <- apply( logDenCompXPop, 2, function(logDenComp){ pmin( maxLogDen, logDenComp) })
logDenCompXPop # keep structure
})
#
logDenCompXChains <- do.call( cbind, logDenCompXPops )
# update the names in intial conditons to hold unique names
rownames(logDenCompXChains) <- resCompNamesUFlat
#
#-- get indices of internal components and TFix (dInfo$ intResCompPosWithin, intResCompPos, posTFix)
# depends on proper resCompNames (inferred from invoking fLogDen with initial states)
intResCompDen <- vector("list",nDen)
#.TFixL <- vector("list", nDen )
#.posTFixL <- vector("list", nDen )
for( iDen in iDens ){
dInfo = dInfos[[iDen]]
#.TFixL[[iDen]] <- dInfo$TFix
#.posTFixL[[iDen]] <- dInfos[[iDen]]$posTFix <- dInfo$resCompPos[ match( names(dInfo$TFix), dInfo$resCompNames ) ] # positions in concatenated result components
intResCompI <- dInfo$intResComp
iPosW <- if( is.character(intResCompI) ){
iPosW <- match(intResCompI, dInfo$resCompNames )
if( any(is.na(iPosW)))
stop(paste("not all names of intResCompNames (",paste(intResCompI,collapse=","),") in return of fLogDen: ",paste(dInfo$resCompNames,collapse=",")))
iPosW
} else intResCompI
dInfos[[iDen]]$intResCompPosWithin <- iPosW
intResCompDen[[iDen]] <- dInfos[[iDen]]$intResCompPos <- dInfo$resCompPos[iPosW] # position in concatenation of all result components
} #iDen
#TFixAll <- do.call(c, .TFixL )
#posTFixAll <- do.call(c, .posTFixL ) # position in concatenation of all resComp
#-- Temperature of result components, i.e. data streams
if( nrow(TSpec) == 1 ){
TSpec <- do.call( rbind, lapply(1:nResComp,function(i){TSpec[1,]}) )
rownames(TSpec) <- resCompNamesUFlat
}
if( nrow(TSpec) != nResComp ){
.missingResCompNames <- unique(resCompNamesUFlat[ !(resCompNamesUFlat %in% rownames(TSpec)) ])
if( length(.missingResCompNames) ){
stop(paste("twDEMCBlockInt: TSpec missing for result components",paste(.missingResCompNames,collapse=",")))
}
TSpec <- TSpec[ resCompNamesUFlat, ]
}
tempResCompPops <- vector("list", nPop )
for( iPop in iPops){
.tempResL <- lapply( 1:nResComp, function(iResComp){
.tempG <- pmax(1,c( TSpec[iResComp,"T0"],
calcDEMCTemp( TSpec[iResComp,"T0"], TSpec[iResComp,"TEnd"], nGenPops[iPop] )))
})
tempResCompPops[[iPop]] <- do.call( cbind, .tempResL)
colnames(tempResCompPops[[iPop]]) <- resCompNamesUFlat
}
#-- setup acceptance rate recording
# number of requried independent generations to choose from (10*d independent states: TerBraak06 TerBraak08) devided by number of chains
if( length(ctrl$initialAcceptanceRate) == 1)
ctrl$initialAcceptanceRate <- matrix(ctrl$initialAcceptanceRate, nrow=nBlock, ncol=nChain, dimnames=list(blocks=NULL, chains=NULL) )
if( !is.matrix(ctrl$initialAcceptanceRate) || dim(ctrl$initialAcceptanceRate) != c(nBlock,nChain))
stop("ctrl$initialAcceptanceRate must be matrix of dim nBlock x nChain")
pAcceptPops <- popMeansTwDEMC(ctrl$initialAcceptanceRate, nPop) # current acceptance rate of population (block x pop)
pAcceptPops[] <- pmax(1/ctrl$thin,pAcceptPops ) # lower bound
#nGenBack <- pmin(mZ,ceiling(ctrl$rIndToAcc * ctrl$gInd * pmax(1,ctrl$pIndStep/(ar * ctrl$thin)))) #number of genrations to select states from for each population, ctrl$gInd is multiplied by number of rows for one step depending on acceptance rate and thinning but at least one
pAcceptMinPops <- apply(pAcceptPops,2,min) # minimum of acceptance rate across blocks for each population
nGenBackPops <- pmin(M0Pops,ceiling(ctrl$gInd * pmax(1,ctrl$pIndStep/(pAcceptMinPops * ctrl$thin)))) #number of genrations to select states from for each population, ctrl$gInd is multiplied by number of rows for one step depending on acceptance rate and thinning but at least one
##details<< \describe{ \item{Acceptance rate}{
## The acceptance rate is tracked for each chain across ctrl$pAcceptWindowWidth thinning intervals. ,
## If acceptance rates drops to low values, this might be because of bad localization
## ,i.e the requirement to use states from more distant past.
## In order to improve localization, less parameters or more chains per population are required.
## }}
# acceptWindow holds the number of accepted steps for each thinning interval for each chain
# if its boundaries are filled up, the last nGenBack states are copied to the first interval
aThinWinW <- ctrl$pAcceptWindowWidth #%/% ctrl$thin # number of thinning intervals over which to calculate average acceptance rate
#ctrl$pAcceptWindowWidth <- aThinWinW * ctrl$thin
# acceptWindow records number of accepted steps in thinnging interval
# (stepslots x blocks x chains)
#acceptWindow <- matrix( rep(ctrl$thin*ar, each=2*aThinWinW*nChainPop), nrow=2*aThinWinW, ncol=d$chains ) #record number of accepted steps in thinning interval
# (width x block x chain )
acceptWindow <- array( rep(ctrl$thin*ctrl$initialAcceptanceRate, each=2*aThinWinW)
, dim=c(2*aThinWinW, dim(ctrl$initialAcceptanceRate) )
, dimnames=c(list(steps=NULL),dimnames(ctrl$initialAcceptanceRate)) ) #record number of accepted steps in thinning interval
#-- preallocate output of parameters, calculated LogDen, indices of accepted steps, and next proposal
# need lists for populations because Nz (number of samples including intial population ) and nGen (samlples without initial) and nThinnedGen is differ across populations in general
ZPops <- lapply( iPops, function(iPop){
dNames <- if( !is.null(dimnames(ZinitPops[[iPop]])) ) dimnames(ZinitPops[[iPop]]) else list( samples=NULL, parms=paste("p",1:nParm,sep="_"), chains=NULL )
Z <- array(NA_real_, dim=c(Nz[iPop],nParm,nChainPop), dimnames=dNames)
Z[1:M0Pops[iPop],,] <- ZinitPops[[iPop]]
Z
})
# several resLogDen per chain (dInfos are concatenated)
resLogDen <- lapply(iPops,function(iPop){ res <- array( NA_real_
, dim=c(1+nThinnedGenPops[iPop], nResComp, nChainPop)
, dimnames=list(steps=NULL, resComp=resCompNamesUFlat, chains=NULL))
res[1,,] <- logDenCompXPops[[iPop]]
res
})
# one rLogDen per chain and per dInfo
logDen <- lapply(iPops, function(iPop){ res <- array(NA_real_
, dim=c(1+nThinnedGenPops[iPop], nDen, nChainPop)
, dimnames=list(steps=NULL, den=names(dInfos), chains=NULL) )
for( iChainPop in 1:nChainPop)
for( iDen in iDens )
res[ 1,iDen,iChainPop] <- sum(logDenCompXPops[[iPop]][resCompPosPops[[iDen]],iChainPop ])
res
})
# one pAcceptance indication per chain and per block
pAccept <- lapply(iPops, function(iPop){ res <- array(NA_real_
, dim=c(1+nThinnedGenPops[iPop], nBlock, nChainPop)
, dimnames=list(steps=NULL, block=names(blocks), chains=NULL) )
})
for( iPop in iPops ){
pAccept[[iPop]][1,,] <- ctrl$initialAcceptanceRate[,chainsPop[[iPop]] ]
}
# record of proposals and fLogDen results, rows c(accepted, parNames, resCompNames, rLogDen)
# if doRecordProposals is FALSE record only the thinning intervals for the last 128 generations
# +1 for the initial state
nThinLastPops <- if(doRecordProposals) nThinnedGenPops else pmin(nThinnedGenPops, ceiling(128/ctrl$thin)) #number of thinning steps to record
nThinOmitRecordPops = nThinnedGenPops-nThinLastPops #the Thinning intervals with no recording of outputs
nGenOmitRecordPops = nThinOmitRecordPops*ctrl$thin
nGenYPops <- nThinLastPops*ctrl$thin #+1?
YPops <- lapply(iPops, function(iPop){ res <- array( NA_real_
, dim=c( nGenYPops[iPop], nParm+nBlock+nResComp, nChainPop)
, dimnames=list(steps=NULL, vars=c(rownames(XPops[[iPop]]),paste("accepted",iBlocks,sep=""),colnames(resLogDen[[iPop]])), chains=NULL ) )
# already initialized first row in nGen cycle
#if( nGenYPops[iPop] == nGenPops[iPop] ) # record first state as accepted
# res[1,,] <- rbind( XPops[[iPop]], matrix(TRUE, nrow=nBlock, ncol=nChainPop), logDenCompXPops[[iPop]] )
res
}) # records proposals, indication of acceptance per Block for each Chain
# here accepted is 0: not accepted, 1: full step accepted, between DR step accepted at DRGamma
#
##details<< \describe{ \item{Overfitting control}{
## Cost, i.e. -2*logDensity, below the number of observations indicates an overfitting to the data stream.
##
## By supplying a maximum logDensity for one result component with \code{dInfo$maxLogDen = -1/2*nObs} one can
## activate the control for this overfitting. This works only, however, if there are enough observations.
## When a result component is based on less than about 30 observartions, one should give a maximum Densitiy of 0.
##
## With overfitting control, there is an upper bound bound for the logDensity of one result compoent-
## Calculated higher logDensities indicate overfitting and are decreases to the maxLogDen value.
##}}
#
#
#-- arguments to fUpdate argument argsFUpdate that do not change within thinning interval
argsFUpdateDefault = list(
ctrl=ctrl
#,fCalcComponentTemp=calcComponentTemp
#,fDiscrProp=fDiscrProp
#,argsFDiscrProp=argsFDiscrProp
)
# construct an list for each block that includes compPos, and intResCompPos
#block <- blocks[[1]]
argsFUpdateBlocks <- lapply( blocks, function(block){
dInfo <- if( block$useMetropolisUpdate) dInfos[[block$dInfoPos]] else NULL
c( argsFUpdateDefault, dInfo=list(dInfo), block
#, list(maxLogDen=maxLogDenL[[block$dInfoPos]]) # defined in dInfo on .checkDInfoResCompNames
#, list(TPropPops=TPropPops[dInfo$compPosDen,])
)
})
remoteArgsFUpdateBlocksTwDEMC <- list(
remoteFun = .updateIntervalTwDEMCChain # this function is called remotely and invokes the blockupdate for all steps of a chain within a thinning interval.
,remoteDumpfileBasename = remoteDumpfileBasename #if null delete
,argsUpdateBlocksTwDEMC=list(
argsFUpdateBlocks = argsFUpdateBlocks
,nBlock=nBlock #, iBlocks=iBlocks
,upperParBoundsPop=upperParBoundsPop
,lowerParBoundsPop=lowerParBoundsPop
,returnIntermediate=TRUE # here stays on slave, so no overhead of returning it between block updates
#,temp=temp
)
)
#
if( !debugSequential && sfParallel() ){
sfExport("remoteArgsFUpdateBlocksTwDEMC") #evaluated in remote process, only passed one time
sfExport(".updateIntervalTwDEMCChain",".updateBlocksTwDEMC",".calcFLogDen", namespace="twDEMC") # used in remote process
suppressWarnings(sfLibrary(abind))
}
#else( assign("remoteArgsFUpdateBlocksTwDEMC", remoteArgsFUpdateBlocksTwDEMC, pos=1)) #export to current
# pass remoteArgsFUpdateBlocksTwDEMC as a name instead of each time complete object in a call
# by using sfRemoteWrapper it is obtained from exported object on the remote node
nameRemoteArgs <- if( !debugSequential && sfParallel() ) as.name("remoteArgsFUpdateBlocksTwDEMC") else remoteArgsFUpdateBlocksTwDEMC # eval.parent fails for remoteArgsFUpdateBlocksTwDEMC within sequential function
#
# arguments to demc that change between thinning intervals
# substract logDen components of logDenCompX from logDenXExt
chainState <- list(
X = XChains
,logDenCompX = logDenCompXChains
,intermediatesChains=lapply( 1:nChain, function(iChain){list()} ) # initially empty
,parUpdateDen = array(TRUE, dim=c(nDen,nParm,nChain), dimnames=list(dens=NULL,parms=parNames,chains=NULL)) ##<< for each parameter/density combination: is the density up to date
#,logDenX = logDenX
)
parUpdateDenLast = chainState$parUpdateDen # to store the last parUpdateDen
#
maxNThinned <- max(nThinnedGenPops)
isSamePopLength = all( nThinnedGenPops == maxNThinned)
isPops <- 1:nPop # populations that take part in this step (some may need less steps)
isChains <- 1:nChain # chains that take part in this step
#iThin0 <- maxNThinned-1
#
#--- main cycle of updates (each ctrl$thin generations)
for( iThin0 in (0:(maxNThinned-1)) ){
if( length(progressOutput)!=0 && nchar(progressOutput)!=0 ){ cat(progressOutput) }
iGen = iThin0*ctrl$thin+(1:ctrl$thin) # the generations within this thinning step
mZPops <- M0Pops + iThin0
iPopsOut <- which(nThinnedGenPops == iThin0) # those pops drop out
#if( length(iPopsOut)) recover()
if( 0 != length(iPopsOut) ){
iChainsOut <- do.call( c, lapply( iPopsOut, function(iPop){ chainsPop[[iPop]] }) )
iiPopsOut <- which(isPops %in% iPopsOut ) # index in set of current populations
iiChainsOut <- rep((iiPopsOut-1)*nChainPop, each=nChainPop) + (1:nChainPop)
isPops <- isPops[-iiPopsOut]
isChains <- isChains[-iiChainsOut]
# recored the last parUpdateDen for dropouts
parUpdateDenLast[,,iChainsOut] <- chainState$parUpdateDen[,,iiChainsOut]
# adapt chainState that refer to current populations only
chainState <- within(chainState,{
X <- X[,-iiChainsOut ,drop=FALSE]
intermediatesChains <- intermediatesChains[-iiChainsOut]
logDenCompX <- logDenCompX[,-iiChainsOut ,drop=FALSE]
parUpdateDen <- parUpdateDen[,,-iiChainsOut ,drop=FALSE]
})
#acceptance rates and nGenBack refer to all chains
}
#
# sample random vectors (note here parameters in rows)
# calculate proposed steps (differences not destinations) within next thinning interval
zxl <- if( (length(blocks)!=1) && ctrl$useConditionalProposal ){
lapply( iPops[isPops], function(iPop){
.sampleStatesConditional(Z=ZPops[[iPop]],mZ=mZPops[iPop],nGenBack=nGenBackPops[[iPop]], nSample=ctrl$thin, blocks=blocks )
})
} else {
lapply( iPops[isPops], function(iPop){
.sampleStates(Z=ZPops[[iPop]],mZ=mZPops[iPop],nGenBack=nGenBackPops[[iPop]], nSample=ctrl$thin )
})
}
zx <- abind( zxl, along= 2 ) # combine chains (and steps within each chain) of all populations
# second dim: nSample x nChainPop x nPop
genPropRes <- .generateXPropThin(zx,ctrl=ctrl, nChain=length(isPops)*nChainPop)
#
# check if proposals should be recorded
# for the first time
iPopsRecord <- which( (1:nPop %in% isPops) & (iThin0 == nThinOmitRecordPops))
for( iiPop in iPopsRecord ){
# record the initial state of the proposals record
iisPop <- match(iiPop,isPops) # index in iPop
YPops[[iiPop]][1,seq_along(parNames),] <- chainState$X[ ,chainsPop[[ iisPop ]] ]
YPops[[iiPop]][1,nParm+iBlocks,] <- 1 #TRUE
YPops[[iiPop]][1,nParm+nBlock+seq_along(resCompNamesUFlat),] <- chainState$logDenCompX[ ,chainsPop[[ iisPop ]] ]
}
# after the proposal has been made
isRecordProposalsPop = ((iThin0 >= nThinOmitRecordPops))[isPops] # only record the last proposals after nThinOmitRecord
#
# numeric matrix (nGenThin x nPop)
#tempGlobalThinStepsL <- lapply( iPops[isPops], function(iPop){ tempGlobalPops[[iPop]][ 1+iGen ] })
#tempGlobalThinSteps <- structure( do.call( cbind,tempGlobalThinStepsL ), dimnames=list(steps=NULL,pops=NULL) )
# numeric array (nGenThin x nResComp x nPop)
tempDenCompThinStepsL <- lapply( iPops[isPops], function(iPop){ tempResCompPops[[iPop]][1+iGen, ,drop=FALSE] })
tempDenCompThinSteps <- structure( abind( tempDenCompThinStepsL, along=3 ), dimnames=list(steps=NULL,resComp=resCompNamesUFlat, pops=NULL) )
#
# here may use code in .tmp.f.testStep
#
#-- do the steps of next thinning interval
resUpdate <- .updateIntervalTwDEMCParChains( X=chainState$X, logDenCompX=chainState$logDenCompX
, intermediatesChains=chainState$intermediatesChains
, parUpdateDen=chainState$parUpdateDen
, xStep=genPropRes$xStep
, rExtra=genPropRes$rExtra
#, tempGlobalSteps=tempGlobalThinSteps
, tempDenCompSteps=tempDenCompThinSteps
, nsPop=length(isPops)
, pAcceptPar=pAcceptPops[,isPops ,drop=FALSE]
,remoteFunArgs=nameRemoteArgs
,debugSequential=debugSequential
,isRecordProposalsPop=isRecordProposalsPop
,isPops=isPops
#,freeMasterNode=ctrl$freeMasterNode
)
# if( any(!is.finite(resUpdate$logDenCompX)) ){
# if( length(remoteDumpfileBasename) ) dump.frames(remoteDumpfileBasename, TRUE)
# stop("twDEMCBlocks (10): encountered non-finite resLogDen of accepted states")
# }
chainState <- resUpdate[ names(chainState) ]
#
#-- calculate accepteance rate
#row in acceptance Window to record acceptance, if exceeds window, copy second part to first (rewind)
acceptPos0 <- aThinWinW + (iThin0 %% aThinWinW)
if( acceptPos0 == aThinWinW ){ #interval exeeded or new, copy second half to first half
acceptWindow[ 1:aThinWinW,, ] <- acceptWindow[ aThinWinW+(1:aThinWinW),, ]
acceptWindow[ aThinWinW+(1:aThinWinW),, ] <- NA
}
acceptWindow[ acceptPos0+1,, isChains] <- resUpdate$accepted
curAcceptRows <- (acceptPos0+1)-(min(aThinWinW-1,max(iThin0,8)):0) # at the start average over at least 4 slots,
#acceptWindow[curAcceptRows,]
pAcceptChains <- apply(acceptWindow[curAcceptRows,, ,drop=FALSE],c(2,3),sum) / (length(curAcceptRows)*ctrl$thin)
#if( !all(is.finite(pAcceptChains[,isChains ,drop=FALSE])))
# stop("twDEMCBlocks: encountered non-finite pAcceptChains")
pAcceptPops <- popMeansTwDEMC(pAcceptChains, nPop) # current acceptance rate of population (block x pop)
pAcceptPops[] <- pmax(1/ctrl$thin,pAcceptPops ) # lower bound
pAcceptMinPops <- apply(pAcceptPops,2,min) # minimum of acceptance rate across blocks for each population
nGenBackPops <- pmin(M0Pops,ceiling(ctrl$gInd * pmax(1,ctrl$pIndStep/(pAcceptMinPops * ctrl$thin)))) #number of genrations to select states from for each population, ctrl$gInd is multiplied by number of rows for one step depending on acceptance rate and thinning but at least one
#
#-- record current thinning step for each population
for( iiPop in seq_along(isPops) ){
iPop <- isPops[iiPop] # translate index in current populations to all populations
mZ <- mZPops[ iPop ]
ipChains <- chainsPop[[iiPop]] # chains within current chains
igChains <- chainsPop[[iPop]] # chains within all chains
ZPops[[iPop]][mZ+1,,] <- chainState$X[,ipChains]
resLogDen[[iPop]][iThin0+2,,] <- chainState$logDenCompX[,ipChains] # first state is initial state
pAccept[[iPop]][iThin0+2,,] <- pAcceptChains[,igChains]
#if( iPop==13 && any(chainState$X["a",ipChains] < -0.6) ) recover()
}
#
#-- record poprosals
iYSteps0 <- (iThin0 - nThinOmitRecordPops)*ctrl$thin
for( iiPop in which(isRecordProposalsPop) ){
# recordProposals refers to current populations (some may have dropped out already)
iPop <- isPops[iiPop]
iiChains = chainsPop[[iiPop]]
YPops[[iPop]][iYSteps0[iPop]+(1:ctrl$thin),,] <- resUpdate$Y[,,iiChains]
}
}# for iThin0
if( length(progressOutput)!=0 && nchar(progressOutput)!=0 ){ cat("\n") }
# recored the last parUpdateDen for dropouts
iChainsOut <- do.call( c, lapply( isPops, function(iPop){ chainsPop[[iPop]] }) )
parUpdateDenLast[,,iChainsOut] <- chainState$parUpdateDen
#-- calculate the log-Density of last states that are not up to date, because they might be used for reinitialization
for( iDen in iDens ){
dInfo <- dInfos[[iDen]]
cd <- dInfo$compPosDen # components used by the current density function
udi <- adrop(parUpdateDenLast[iDen,, ,drop=FALSE],1) # (iDen= x par x chain
iChains <- which( apply(udi,2,function(ud){ !all(ud) }) )
# TODO parallelize those calculations
for( iChain in iChains ){
iPop <- popChain[iChain]
iChainInPop <- ((iChain-1) %% nChainPop) +1
Xc <- ZPops[[iPop]][(M0Pops+nThinnedGenPops)[iPop],cd,iChainInPop ]
#resCompDenC <- do.call( dInfo$fLogDen, c(list(Xc), dInfo$argsFLogDen) ) # evaluate logDen
resCompDenC <- .calcFLogDen( dInfo$fLogDen, xP=Xc, argsFLogDen=dInfo$argsFLogDen, compInt=c(), TiInt=c(), maxLogDen = dInfo$maxLogDen ) # evaluate logDen
resLogDen[[iPop]][ 1+nThinnedGenPops[iPop],dInfo$resCompPos,iChainInPop] <- resCompDenC
} # end need update
} # for iDen
#-- sum over resLogDen components to yield logDens for each given density
# (might be slightly off because not up to date, but work for inspecting the trend)
for( iDen in iDens ){
dInfo <- dInfos[[iDen]]
rcp <- dInfo$resCompPos # components used by the current density function
if( length(rcp) == 1){
for( iPop in iPops ){
logDen[[iPop]][,iDen,] <- resLogDen[[iPop]][,rcp,]
}
}else{
for( iPop in iPops ){
resLogDenI <- resLogDen[[iPop]][,rcp, ,drop=FALSE]
logDen[[iPop]][,iDen,] <- apply(resLogDenI,c(1,3),sum)
}
} # length(resComp) == 1
} # iDen
##details<< \describe{ \item{subspace localized sampling}{
## The covariance matix for new proposals might look very different in different locations of the parameter space.
## For example, overall CO2 efflux might be composed of two parts microbial respiration and inorganic CO2 production.
## that are relevant in different parts of the parameter space (governed by another parameter). The mirobial growth
## rate will be very sensitive in the respiration part and is allowed to to only small jumps, while it might vary large
## in the inorganic part. Hence, it might be beneficial to use different jumping proposal and different population
## samples in different subspaces.
##
## In order to support such supspace-localized sampling XX
##
## When using parallel execution on multiple cores, the transfer of very big intermediates can cause
## performance issues. Use argument \code{controlTwDEMC$returnIntermediate=FALSE} to ensure that
## intemediates are stored only within process boundaries. When executing the next interval or the next round of
## updates of all blocks, the intermediate will then be the empty list initially.
##}}
##details<< \describe{ \item{intermediate results}{
## Different density functions can share intermediate results, so that these do not need to be recomputed.
## For example two densities can be based on the same model predictions for given parameters, but compare different
## subsets of the predictions against different observations. If for a given parameter vector, the model output has been
## computed in one density, it does not need to be recomputed in the second density function.
##
## This is an advanced option. Care must be taken that the intermediate is really the same between
## densities. And care must be taken that intermediate is recalculated if the parameters change on
## which the intermediate is based on.
##
## Description of both block update in argument \code{blocks} (or argument \code{dInfos} if default Metropolis is used)
## should specify the same identifier string in argument \code{dInfo$intermediateId}.
## The densities or block update functions should return the model output in attribute "intermediate"
## with the result vector of logDensity
## components. \code{twDEMC} ensures that this result is provided with argument \code{intermediate}
## in another call to the density or block update function.
## With a new parameter proposal in Metropolis update, an empty list will be given with this argument.
##
## For an example see test case \code{ofMultiIntermediate} in file unitTests/runittwDEMC.R.
## Using densities \code{\link{denSparsePrior}} and \code{\link{denRichPrior}}.
##
## When using parallel execution on multiple cores, the transfer of very big intermediates can cause
## performance issues. Use argument \code{controlTwDEMC$returnIntermediate=FALSE} to ensure that
## intemediates are stored only within process boundaries. When executing the next interval or the next round of
## updates of all blocks, the intermediate will then be the empty list initially.
##}}
#-- return
for( iPop in iPops ){
if( any(!is.finite(resLogDen[[iPop]])) ){
if( length(remoteDumpfileBasename) ) dump.frames(remoteDumpfileBasename, TRUE)
stop("twDEMCBlocks (20): encountered non-finite resLogDen of accepted states")
}
}
spacesPop <- sapply(pops,"[[","spaceInd")
##value<< a list with entries of populations, each entry is a list
res <- list(
thin=ctrl$thin ##<< thinning interval that has been used
,dInfos=dInfos ##<< list of information on densities (argument \code{dInfos})
,blocks=blocks ##<< list of information on blocks (argument \code{blocks})
,YPos = list( ##<< list of column positions in pops[i]$Y, a list with entries \describe{
##describe<<
accepted = nParm+iBlocks ##<< integer vector of positions of acceptance indication of block at given index
,resLogDen0 = nParm+nBlock ##<< integer scalar: postion before first column of results of fLogDen
)
##end<<
,pops = lapply( iPops, function(iPop){ list( ##<< info on each population. A list with entries: \describe{
##describe<<
upperParBounds = upperParBoundsPop[[iPop]] ##<< upper parameter bounds for sampling
,lowerParBounds = lowerParBoundsPop[[iPop]] ##<< lower parameter bounds for sampling
,splits=pops[[iPop]]$splits ##<< named numeric vector: splitting history
,spaceInd = spacesPop[iPop] ##<< the space replicate that the population belongs to
,parms = ZPops[[iPop]][ M0Pops[iPop]:nrow(ZPops[[iPop]]),, ,drop=FALSE] ##<< numeric array (steps x parms x chains): collected states, including the initial states
,temp = tempResCompPops[[iPop]][seq(1,nGenPops[iPop]+1,by=ctrl$thin), ,drop=FALSE] ##<< numeric array (nSample+1, nResComp): temperature, i.e. cost reduction factor in each step for each datastream
,pAccept= pAccept[[iPop]] ##<< acceptance rate of chains (nStep x nChainPop)
,resLogDen = resLogDen[[iPop]] ##<< numeric array (steps x resComps x chains): results components of fLogDen of blocks
,logDen = logDen[[iPop]] ##<< numberic array (steps x iDen x chains): results summed over blocks
,Y = YPops[[iPop]] ##<< numeric matrix (steps x nParam+nBlock+nResComp x chains): parms, isAccepted for each block, and all fLogDen result components for each proposal
)})
##end<<
)
##end<<
#names(res$pops) <- names(pops)
class(res) <- c( class(res), "twDEMCPops" ) #monte carlo proposal list
res
}
attr(twDEMCBlockInt,"ex") <- function(){
# see unit test
# twUtestF("twDEMCPops","distinctLogDen")
data(twLinreg1); attach( twLinreg1 )
# tracing error in remote session:
# Provide argument \code{remoteDumpfileBasename="dumpFile"}
# Then a dumpfile is created on remote error by \code{\link{sfRemoteWrapper}}.
# In order to trace the density function, the following can be done
if( FALSE ){
.remoteDumpfileBasename="dumpfile"
.remoteDumpfile <- paste(.remoteDumpfileBasename,".rda",sep="")
load(.remoteDumpfile)
debugger(get(.remoteDumpfileBasename))
# choose last step (18)
#require(debug)
fDen <- remoteFunArgs$argsUpdateBlocksTwDEMC$argsFUpdateBlocks[[1]]$fLogDen
mtrace(fDen); remoteFunArgs$argsUpdateBlocksTwDEMC$argsFUpdateBlocks[[1]]$fLogDen <- fDen
do.call(remoteFun, c(remoteFunArgs, list(...)))
# go()
# qqq()
}
}
.tmp.f <- function(){ # testStep
# test one step for one chain
iStep <- 1
iChain <- 1
# those are arguments, generated updateBlocksTwDEMCPar
iPop <- popChain[iChain]
step <- genPropRes$xStep[,iChain,iStep]
rExtra <- genPropRes$rExtra[iChain,iStep]
Xc <- chainState$X[,iChain]
logDenCompC <- chainState$logDenCompX[,iChain]
parUpdateDenC <- chainState$parUpdateDen[,,iChain]
tempC <- tempThinSteps[iChain,iStep]
.updateBlocksTwDEMC(iPop = iPop,step=step,rExtra=rExtra,Xc=Xc,logDenCompC=logDenCompC
,parUpdateDenC=parUpdateDenC,pAccept=pAccept, temp=tempC
, argsUpdateBlocksTwDEMC=remoteArgsFUpdateBlocksTwDEMC$argsUpdateBlocksTwDEMC )
}
.tmp.f <- function(){
# before (code here) calculated temperature for each chain, but actually need it only by population
tempThinStepsL <- lapply( iPops, function(iPop){
matrix( temp[[iPop]][ 1+iGen ], nrow=length(iGen), ncol=nChainPop )
})
tempThinSteps <- t(do.call( cbind, tempThinStepsL)) # steps need to be in last dimension (so chains in first)
}
.getRescompNameBlock <- function(
### appending iBlock to the resCompNames to make the names unique among blocks
resCompNames
,iBlock
){
paste(resCompNames,iBlock,sep="_")
}
.completePopDescription <- function(
### filling default values in pop description
pop
,nGen=numeric(0) ##<< number of generations to take: overwrites settings given in pop
,spaceInd=numeric(0) ##<< the space replicate that this population belongs to, overwritten by setting in pop
){
##details<<
## If several entries are null, they are set to default values
if( 0==length(pop$parms) || !is.numeric(pop$parms) || length(dim(pop$parms)) != 3 )
stop(".checkPop: entry parms must be a numeric 3D array")
if( 0 != length(nGen) ) pop$nGen <- nGen
if( (1 != length(pop$nGen)) )
stop(".checkPop: scalar integer nGen must be given as entry nGen in population or in invocation of twDEMCBlockInt.")
if( 0 == length(pop$spaceInd) ) pop$spaceInd = spaceInd
if( (1 != length(pop$spaceInd)) )
stop(".checkPop: scalar integer spaceInd must be given as entry nGen in population or in invocation of twDEMCBlockInt.")
# if( 0==length(pop$T0) ) pop$T0 <- 1
# if( 0==length(pop$TEnd) ) pop$TEnd <- 1
if( 0==length(pop$X) ){
pop$X <- adrop( pop$parms[nrow(pop$parms),, ,drop=FALSE],1)
}else{
if( !all.equal( rownames(pop$X), colnames(pop$parms), check.attributes = FALSE) )
stop(".checkPop: rownames of X and colnames of parms must match")
}
if( 0 != length(pop$upperParBounds) && any(is.na(pop$upperParBounds)))
stop(".checkPop: encountered NA in upperParBounds")
if( 0 != length(pop$lowerParBounds) && any(is.na(pop$lowerParBounds)))
stop(".checkPop: encountered NA in lowerParBounds")
### pop with entries nGen, and X defined
pop
}
.completeDInfo <- function(
### filling default values in information on logDensities (dInfo)
dInfo ##<< list of information on logDensities to be checked for entries
,parNames ##<< string vector: names of all parameters. A subset (compPosDen) of these is required by the logDensity
,resCompNames ##<< string vector: list of names of result vector of logDensity function.
){
##details<<
## \code{dInfo} stores ample information on a log-Density function in order to
## help fast execution of calling the log-Density function and the metropolis update. However, the user
## is required to specify only minimal detail (see \code{\link{twDEMCBlockInt}}) and most of the information is
## filled in this function.
## Some of the information can only be completed after the names of logDensity results are
## know after calling the function. Therefore there is a second function \code{\link{.completeDInfoResCompNames}}.
#
if( 0 == length(dInfo$fLogDen) || !is.function(dInfo$fLogDen))
stop(".checkLogDenInfo: entry fLogDen (function of log of unnomalized density) must be specified")
#
##details<<
## \code{compPosDen} specifies the parameters provided to the given logDensity.
## If entry compPosDen is not given, then assume all parameter dimensions are required.
## If compPosDen is given as a string vector, then it is translated to positions within
## vector of all parameter names.
if( 0 == length(dInfo$compPosDen) ) dInfo$compPosDen <- seq_along(parNames)
compPosDenOrig <- dInfo$compPosDen
if( is.character(compPosDenOrig) ) dInfo$compPosDen <- match(compPosDenOrig,parNames)
if( any(is.na(dInfo$compPosDen)) ){
stop(".completeDInfo: some of compPosDen specified in dInfo are among parameter names (of initial values).")
}
#
if( 0==length(dInfo$intResComp) ) dInfo$intResComp <- vector("character",0)
#if( 0==length(dInfo$TFix) ) dInfo$TFix <- vector("numeric",0)
#
if( 0!=length(dInfo$TFix) ) warning("Usage of TFix is deprecated. Use arguemt TSpec of twDEMCBlockInt with T0 and TEnd of one for given components.")
if( 0!=length(dInfo$intermediate) || 0!=length(dInfo$intermediateId) ) warning("intermediateId now must be specified with blocks.")
if( length(i <- which(!(names( dInfo) %in% c("fLogDen","xC","logDenCompC","parUpdateDenC","compPosDen","intResComp","argsFLogDen"
,"fDiscrProp","argsFDiscrProp"
,"resCompNames","resCompPos","intResCompPos","intResCompPosWithin","nObs","maxLogDen")))) )
stop(paste(".completeDInfo: unknown entries in dInfo:",paste(names(dInfo)[i],collapse=",") ))
if( 0 != length(dInfo$nObs))
warning(".checkDInfoResCompNames: usage of nObs in dInfos is deprecated and not considered. Instead use: maxLogDen=-1/2*nObs")
#
dInfo$resCompNames <- resCompNames
#
if( 0 == length(dInfo$maxLogDen) ) dInfo$maxLogDen <- rep( 0, length(resCompNames) )
if( length(names(dInfo$maxLogDen)) ){
# check that all names matches the resCompNames, and sort along resCompNames
iM <- match( resCompNames, names(dInfo$maxLogDen) )
if( any(is.na(iM)) )
stop(paste(".checkLogDenInfo: maximum logDensitiy, maxLogDen, of observations missing for result components ",paste(resCompNames[which(is.na(iM))], collapse=","),sep="") )
dInfo$maxLogDen <- dInfo$maxLogDen[iM]
}
if( length(dInfo$maxLogDen) != length(resCompNames) )
stop(paste("dInfo maxLogDen (length=",length(dInfo$maxLogDen),") must have same lenght as number of result components=",length(resCompNames),sep=""))
### argument \code{dfo} with entries fLogDen, compPosDen, intResComp, resCompNames, and maxLogDen properly specified
dInfo
}
.completeBlockInfo <- function(
### filling default values in block description
block ##<< list to be checked for entries
,dInfos ##<< checked list of dInfo. See the argument in \code{\link{twDEMCBlockInt}}.
,parNames ##<< parameter dimension names
){
##details<<
## Similar to dInfos (see \code{\link{.completeDInfo}}), there is a list with ample information on each block.
## As the user is required to specify only the minimal information or information in more convenient format
## The block list is checked and information is completed.
#
if( 0 == length(block$compPos) ) stop(".completeBlockInfo: entry compPos must be specified, either as indices or as variable names")
##details<<
## If the user specified \code{compPos} as a character string,
## then convert it to positions within vector of all parameter
compPosOrig <- block$compPos
if( is.character(compPosOrig) ){
block$compPos <- match( compPosOrig, parNames)
}
if( any(is.na(block$compPos)) ){
i <- which(!(names(compPosOrig) %in% parNames))
stop(paste(".completeBlockInfo: unknown variable names in entry compPos in block:",paste(parNames[i],collapse=",") ))
}
##details<<
## If no other block update function is given, assume Metropolis update
if( 0 == length(block$useMetropolisUpdate) )
block$useMetropolisUpdate = (0==length(block$fUpdateBlock))
if( block$useMetropolisUpdate ){
block$fUpdateBlock <- updateBlockTwDEMC
dInfoPosOrig <- block$dInfoPos
if( 1 != length(block$dInfoPos) ) stop(".checkBlock: dInfoPos must be of length 1.")
if( is.character(dInfoPosOrig) ) block$dInfoPos <- match(dInfoPosOrig,names(dInfos))
if( is.na(block$dInfoPos) ){
stop(paste0(".checkBlock: specified density '",dInfoPosOrig,
"' not found in names of list of density information (argument dInfos)."))
}
dInfo <- dInfos[[block$dInfoPos]]
#
# positions of parameters to update in list of parameters to logDensity function
block$compPosInDen <- match(block$compPos,dInfo$compPosDen )
if( any(is.na(block$compPosInDen)) ){
stop(".checkBlock: some of compPos specified in block are not in compPosDen of fLogDenInfo.")
}
} # useMetropolisUpdate
#
if( 1!=length(block$fUpdateBlock) || !is.function(block$fUpdateBlock) ) stop(".checkBlock: fUpdateBlock must be a single function.")
#
if( 0 == length(block$nStates) ){
# number of states to choose from for proposal generations on conditional sampling
block$nStates <- calcM0twDEMC( length(block$compPos), 1)
}
#
##value<< argument block with entries compPos, dInfoPos, fUpdateBlock, and useMetropolisUpdate definded
if( length(i <- which(!(names( block) %in% c("compPos","dInfoPos","fUpdateBlock","useMetropolisUpdate","argsFUpdate"
,"nStates","compPosInDen","intermediateId")))) )
stop(paste(".checkBlock: unknown entries in block:",paste(names(block)[i],collapse=",") ))
block
}
.getResFLogDenNames <- function(
### Extracting/Creating names for result components of logDenComp
logDenComp ##<< vector or array (vars in rows) of results of logDensity function
,logDenName ##<< name of the density
){
##details<<
## Names of the result components of fLogDen are used to distinguish internal components.
## If the density function does not provide names, they are created.
## If the density function has only one result component, the name is lost during parallel execution.
## Hence, a single name is also re-created, to avoid errors on checking.
## The created names are of the form <logDenName>_<iComp>
if( is.array(logDenComp)){
res <- rownames(logDenComp)
return( if( !is.null(res) ) res else paste(logDenName,1:nrow(logDenComp),sep="_") )
}else{
res <- names(logDenComp)
#return( if( !is.null(res) && length(res) > 1) res else paste(logDenName,1:length(logDenComp),sep="") )
return( if( !is.null(res) ) res else paste(logDenName,1:length(logDenComp),sep="_") )
}
}
.sampleStatesConditional <- function(
### Sample records and initial state from past records.
Z ##<< population (steps x parms x chains)
,mZ ##<< number of steps in Z that have been already sampled (Z is bigger)
,nGenBack ##<< number of samples in recent past, used to generate the proposals
,nSample ##<< number of proposals to generate
,blocks ##<< information on blocks with item \code{compPos} specifiy the parameter positions
##<< and item \code{nStates} specifying the number of closest states to choose from, usually calculated by CalcM0 with one chain
#,rLogDen
){
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.generateXPropThin}}
nChainPop = dim(Z)[3]
nParm = dim(Z)[2]
# blocks <- lapply(blocks, function(block){ within(block, nStates <- 8) })
# calculate the number of states
# blocks <- lapply(blocks, function(block){ within(block, nStates <- m0b*nChainPop) })
Zt <- t(apply(Z[1:(mZ),,],2,function(Zs){Zs})) # matrix (nParm x nState): all states except the last (current) one
#iChain <- 1
zxbL <- lapply( 1:nChainPop, function(iChain){
Xc <- Z[mZ,,iChain] # current state
d <- abs(Xc - Zt) # absolute deviation from current state (as parameters are the first dimension, we can apply array diff)
#block <- blocks[[1]]
ZtsubL <- lapply( blocks, function(block){
db <- d[-block$compPos, ,drop=FALSE] # subset of parameters not in block
rbp <- apply( db, 1, rank) # rank of difference for each parameter
rb <- apply( rbp,1, max) # maximum over parameters
iStates <- order(rb)[1:block$nStates] # the index of the states to choose from
iS <- matrix( iStates[sample.int(block$nStates, 3*nSample, replace=TRUE)], ncol=3 ) #TODO save invocations to sample.int by sampling before
Ztsub <- Zt[block$compPos, iS,drop=FALSE]
Xcb <- Xc[block$compPos]
zx <- array(c(Ztsub, rep(Xcb,nSample)), dim=c(parms=nrow(Ztsub), step=nSample, iz=4), dimnames=list(parms=rownames(Ztsub),step=NULL,iz=NULL) )
# check that z1 and z2 are different
# check that z3 is different from x
##details<<
## States z2 is attempted to be distinct.
## And z3 is attempted to be distinct from x (assuming steps in z are stacked chains collectives of X)
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,2] <- zx[,sample.int(dim(zx)[2],length(iSame)),2] # sample from z2 of other steps
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i<-i+1
}
#(tmp <- which( zx[,,2] == zx[,,1], arr.ind=TRUE))
#when z3 is the same as x, i.e. zx[,,4]
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
#iSame <- unique(which( (zx[-nrow(zx),,3]==Xs), arr.ind=TRUE )[,2])
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,3] <- zx[,sample.int(dim(zx)[2],length(iSame)),3] # sample from z3 of other steps
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
i<-i+1
}
zx
}) # updates for each block
zxb <- abind(ZtsubL, along=1)
}) # updates for each chain
zxbc <- abind(zxbL, along=2) # stack chains across samples
### numeric array: (nParm x (nSample*nChainPop) x 4) of randomly choosen states
### states of different blocks are choosen independently of each other, i.e. are combinations from different states
### 4th state is the currently accepted state of the chain (for Snooker update)
zxc <- zxbc[ colnames(Z),, ] # order parameters
}
.sampleStates <- function(
### Sample records and initial state from past records.
Z ##<< population (steps x parms x chains)
,mZ ##<< number of steps in Z that have been sampled
,nGenBack ##<< number of samples in recent past, used to generate the proposals
,nSample ##<< number of proposals to generate
#,rLogDen
){
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.generateXPropThin}}
nChainPop = dim(Z)[3]
nParm = dim(Z)[2]
##details<<
## Random states for chains for difference vector generation are within subsets of populations.
## This allows simulating independent population of chains.
## The acceptance rate may differ among populations. Hence, the set of previous generations to
## randomly select from also differs between poplations.
# integer array (thinSteps*nChain*4) sample of chains, within populations
X <- adrop(Z[mZ,,,drop=FALSE],1) # assume current state X as beginning of the interval
nStates <- nSample*nChainPop*3 # need to sample three states for snooker update
sChains <- 1:nChainPop
sGens <- (mZ-nGenBack+1):mZ
rrGenPop <- sample(sGens, nStates, replace = TRUE)
rrChainsPop <- sample(sChains, nStates, replace = TRUE)
# in order to constrain two dimensions at the same time use the [] subset with an array see ?"["
rr <- cbind(rrGenPop,rrChainsPop)
#zLogDen <- array(rLogDen[rr], dim=c(1,nSample*nChainPop,3), dimnames=list(parms="logDen", steps=NULL, zi=NULL) )
rrParms <- cbind( rep(rrGenPop,each=nParm), rep(1:nParm, nStates), rep(rrChainsPop,each=nParm) )
zParms <- array(Z[rrParms], dim=c(nParm,nSample*nChainPop,3), dimnames=list( parms=rownames(Z), steps=NULL, zi=NULL) )
# note that parameters are the first dimension
#rrParms <- cbind( rep(rrGenPop,nParm), rep(1:nParm, each=nStates), rep(rrChainsPop,nParm) )
#zParms <- matrix( Z[rrParms], ncol=nParm)
#z0 <- abind( zParms, zLogDen, along=1)
#
# append as forth initial state x vector to z
# assume that chains is the last dimension, for each step assume the same initial state x
chainsZ <- rep(sChains,each=nSample)
Xs <- array(X[,chainsZ], dim=c(nParm,nSample*nChainPop,1), dimnames=list(parms=rownames(Z), steps=NULL, zi=NULL) )
#XsLogDen <- array( rLogDen[mZ,chainsZ],dim=c(1,nSample*nChainPop,1), dimnames=list(parms="logDen", steps=NULL, zi=NULL) )
#Xs <- abind( Xs, XsLogDen, along=1) #logDen row never used for X
zx <- abind( zParms, Xs, along=3)
#
##details<<
## States z2 is attempted to be distinct.
## And z3 is attempted to be distinct from x (assuming steps in z are stacked chains collectives of X)
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,2] <- zx[,sample.int(dim(zx)[2],length(iSame)),2]
iSame <- twWhichColsEqual( adrop(zx[,,2 ,drop=FALSE],3), adrop(zx[,,1 ,drop=FALSE],3) )
i<-i+1
}
#(tmp <- which( zx[,,2] == zx[,,1], arr.ind=TRUE))
#when z3 is the same as x, i.e. zx[,,4]
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
#iSame <- unique(which( (zx[-nrow(zx),,3]==Xs), arr.ind=TRUE )[,2])
i <- 0
while( 0 < length(iSame) && i<5){
zx[,iSame,3] <- zx[,sample.int(dim(zx)[2],length(iSame)),3]
iSame <- twWhichColsEqual( adrop(zx[,,3 ,drop=FALSE],3), adrop(zx[,,4 ,drop=FALSE],3) )
i<-i+1
}
#(tmp <- which( zx[,,3] == zx[,,4], arr.ind=TRUE))
### random states (Nparms+1,(steps*nChainPop), 4)
### first dimension is the state vector
### random states for each step and chains (chains stacked to be vectorized)
### chain is last dimensionion in stack (consequtive steps for one chain) in order to abind across populations
### zx dim: three random vectors, forth dimension is the initial state x of the chain
rownames(zx) <- colnames(Z) # preserve parameter names
zx
}
.generateXPropThin <- function(
### Generate Proposals from given random states.
zx ##<< output of \code{\link{.sampleStates}} (maybe stackes states)
,ctrl ##<< list with tuning parameters, see \code{twDEMCBlockInt}
,nChain ##<< number of chains
){
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.sampleStates}}
## \code{\link{.xStepSnooker}}
## \code{\link{.xStepParallel}}
#
nParm <- nrow(zx)
nState <- ncol(zx)
nSample <- nState/nChain
#zxl <- lapply( 1:nPop, fzPop )
#zxAndLogDen <- structure( abind(zxl, along=2) , dimnames=list(parms=c(rownames(Z),"logDen"),steps=NULL,zi=NULL))
#zx <- zxAndLogDen[1:nParm,,,drop=FALSE]
z <- zx[,,1:3,drop=FALSE] #three random state vectors per step
X <- adrop(zx[,,4,drop=FALSE],3) #initial state vector for step
#zLogDen <- adrop(zxAndLogDen[,,1:3,drop=FALSE][nParm+1,,,drop=FALSE],1)
dz <- as.list( structure(dim(z), names=c("parms","steps","zi")) )
#
res <- matrix( NA_real_, nrow=dz$parms+1, ncol=dz$steps, dimnames=c(list(c(rownames(z),"rExtra")),list(NULL)) )
boSnooker <- runif(dz$steps)< ctrl$pSnooker
if( 0 < sum(boSnooker) ){
res[,boSnooker] <- .xStepSnooker(z[,boSnooker,,drop=FALSE],X[,boSnooker,drop=FALSE],ctrl=ctrl)
}
if( 0 < sum(!boSnooker) )
#res[,!boSnooker] <- .xStepParallel(z[,!boSnooker,,drop=FALSE],zLogDen=zLogDen[!boSnooker,,drop=FALSE],ctrl=ctrl)
res[,!boSnooker] <- .xStepParallel(z[,!boSnooker,,drop=FALSE],ctrl=ctrl)
#
#second dimension is Nsteps*nChain (we can set nChains as last dimension
#array(chainsZ,dim=c(nSample,d$chains))
xStepAndExtra <- resArraySteps <- array(res, dim=c(nParm+1,nSample,nChain), dimnames=list(parms=rownames(res),steps=NULL,chains=NULL) )
#deprecated expected steps as last dimension
#xStepAndExtra <- aperm(resArraySteps,c(1,3,2))
# numeric array (Npar+1,Nchains,Nsteps): difference vectors in parameter space for steps and chains
# last row is the extra LogDen associated with snooker update
#
xStep <- xStepAndExtra[-nrow(xStepAndExtra),,,drop=FALSE] #dito
rExtra <- adrop(xStepAndExtra[nrow(xStepAndExtra),,,drop=FALSE],1) #second dim (columns step within Thinning interval)
list(xStep=xStep, rExtra=rExtra)
### List with components \describe{
### \item{xStep}{numeric array (Npar x NSteps x Nchain): difference vectors in parameter space}
### \item{rExtra}{numeric matrix (Nsteps x NChain): some extra LogDensity from snooker update}}
### Nsteps=ctrl$thin
}
.xStepSnooker <- function(
### Generates Snooker updates based on given random numbers.
z, ##<< numeric array (Nparms,(nsteps), 3) of random states, dimnames parms,steps,zi
X, ##<< current state (Nparms,(nsteps)) corresponding to chain of second dimension in z
ctrl
){
# DE-Snooker update
##seealso<<
## \code{\link{.generateXPropThin}}
nParm <- nrow(z)
nState <- ncol(z)
gamma_snooker = runif(nState, min=1.2,max=2.2)
res <- matrix( as.numeric(NA), nrow=nParm+1, ncol=nState, dimnames=c(list(c(rownames(z),"rExtra")),list(NULL)) )
#need loop because of inner products
for( i in 1:nState){
x_z = X[,i] - z[,i,3]
D2 = max(1.0e-300, x_z %*% x_z)
#gamma_snooker =1.7
projdiff = ((z[,i,1] -z[,i,2]) %*% x_z)/D2 # inner_product of difference with x_z / squared norm x_z
res[-(nParm+1),i] <- xStepChain <- (gamma_snooker[i] * projdiff) * x_z
xPropChain = X[,i] + xStepChain
x_z = xPropChain - z[,i,3]
D2prop = max((x_z%*%x_z), 1.0e-30)
res[(nParm+1),i] <- rExtra <- ctrl$Npar12 * (log(D2prop) - log(D2)) # Npar12 =(nParm - 1)/2 # extra term in logr for accept - reject ratio
}
res
} # DE-Snooker update
.xStepParallel <- function(
### DE-parallel direction update based on given random numbers.
z, ##<< numeric array (Nparms,(nsteps), 3) of random states, dimnames parms,steps,zi
#zLogDen, ##<< numeric matrix (nsteps,3): logDen corresponding to the random states z
ctrl
){
##seealso<<
## \code{\link{.generateXPropThin}}
nParm <- nrow(z)
nState <- ncol(z)
dz <- adrop((z[,,1,drop=FALSE]-z[,,2,drop=FALSE]),3) #jump vector as the difference between two random states (from z2 towards z1)
boGammasF1 <- (runif(nState) < ctrl$pGamma1)
gamma_par <- matrix( ctrl$F1, nrow=nParm, ncol=nState) # to be able to jump between modes
gamma_par[,!boGammasF1] <- ctrl$F2 * runif( nParm*sum(!boGammasF1), min=1-ctrl$epsMult, max=1+ctrl$epsMult) # multiplicative error to be applied to the difference
xStepChain <- if (ctrl$epsAdd ==0) { # avoid generating normal random variates if possible
gamma_par * dz
} else {
gamma_par * dz + rnorm(length(dz),0,ctrl$epsAdd)
}
xStepChainAndRExtra <- rbind(xStepChain,rExtra=0)
### Numeric matrix (Nparms, nsteps): difference vectors in parameter space.
}
.tmp.f <- function(){
if( !is.null(ctrl$probUpDir) && (ctrl$probUpDir != 1/2) ){
iFinites <- which(is.finite(zLogDen[,1]) & is.finite(zLogDen[,2]))
}
}
.updateIntervalTwDEMCParChains <- function(
### Perform update steps within next thinning interval parellelizing by chains.
# first three arguments relate to current state (stored in chainState in twDEMCBlockInt)
X ##<< numeric matrix: current location of chains rows: parameters, columns: chains
,logDenCompX ##<< numeric matrix: result of calls to density functions from previous step
,intermediatesChains = lapply(1:nsChain, function(iChain){ list()} ) ##<< list(nChain) intermediate states
,parUpdateDen ##<< for each parameter/density combination: is the density up to date from previous step
,xStep ##<< array rows: difference vectors in parameter space, cols: steps, zdim: chains
,rExtra ##<< numeric matrix: row: step within thinning interval, col: chain
,tempDenCompSteps ##<< numeric array: (nGenThin x nResComp x nsPop)
,nsPop ##<< the number of populations that take part in this step
,pAcceptPar ##<< numeric matrix: current acceptance rate for each (block x population)
#argsUpdateBlocksTwDEMC, ##<< list of arguments to updateBlocksTwDEMC
#argsDEMCStep, ##<< see \code{\link{.doDEMCStep}}
,remoteFunArgs ##<< arguments that do not change across calls
## may be specified as a name instead of complete object (see code{\link[twSnowfall]{sfRemoteWrapper}})
## Must include entries \itemize{
## \item{remoteFun=.updateBlocksTwDEMCChains}
## \item{argsUpdateBlocksTwDEMC
## \item{remoteDumpfileBasename}
## }
## Can be a name of a previously exported variable.
,debugSequential=FALSE ##<< see \code{\link[twSnowfall]{sfFArgsApplyDep}}
,isRecordProposalsPop=rep(FALSE,nsPop) ##<< logical vector: for each population: if TRUE then proposals and results of rLogDen are recorded in result$Y.
, chainsPop=matrix(1:ncol(X),ncol=nsPop, dimnames=list(iChainInPop=NULL, iPop=NULL)) ##<< list of integer vectors:
, isPops ##<< integer vector of populations that have not yet dropped out
){
# .updateIntervalTwDEMCParChains
dimXStep <- dim(xStep)
nStep <- dimXStep[2]
nsChain <- dimXStep[3]
nChainPop = nsChain / nsPop
#
# if( length(intermediatesChains[[1]]) ) recover()
#
#iChain <- 1
argsListsChain <- lapply( 1:nsChain, function(iChain){
iiPop <- ((iChain-1) %/% nChainPop) +1
iPop <- isPops[ iiPop ] # refer to population in non-dropped array
list( X = X[,iChain]
,logDenCompX = logDenCompX[,iChain]
, intermediates = intermediatesChains[[iiPop]]
, parUpdateDen = adrop(parUpdateDen[,,iChain,drop=FALSE],3)
, xStep = adrop(xStep[,,iChain,drop=FALSE],3)
, rExtra = rExtra[,iChain]
, tempDenCompSteps = adrop(tempDenCompSteps[,,iiPop,drop=FALSE],3)
, pAcceptPar = pAcceptPar[,iiPop]
, isRecordProposalsPop = isRecordProposalsPop[iiPop]
, iPop = iPop
)
})
#
#remoteFunArgs$remoteFun <- .updateIntervalTwDEMCChain # on changed function
res <- if( isTRUE(debugSequential) ){
lapply( argsListsChain, remoteFunArgs$remoteFun, argsUpdateBlocksTwDEMC=remoteFunArgs$argsUpdateBlocksTwDEMC )
}else{
# if( sfParallel() ) ( list=c("nChainPop", "X", "logDenCompX", "parUpdateDen", "xStep"
# , "rExtra", "tempDenCompSteps", "pAcceptPar"
# , "remoteFunArgsChain", "isRecordProposalsPop","returnIntermediateChains"))
sfLapply( argsListsChain, sfRemoteWrapper, remoteFunArgs = remoteFunArgs )
}
##value<< list with components
res1 <- res[[1]]
X <- array( NA_real_, c(length(res1$X), nsChain), dimnames=list(parms=names(res1$X), iChain=NULL))
logDenCompX <- array( NA_real_, c(length(res1$logDenCompX), nsChain), dimnames=list(parms=names(res1$logDenCompX), iChain=NULL))
parUpdateDen <- array( NA, c(dim(res1$parUpdateDen), nsChain), dimnames=c(dimnames(res1$parUpdateDen), iChain=NULL))
accepted <- array( NA_real_, c(length(res1$accepted), nsChain), dimnames=list(iBlock=names(res1$accepted), iChain=NULL))
Y <- array( NA_real_, c(dim(res1$Y), nsChain), dimnames=c(dimnames(res1$Y), iChain=NULL))
for( iChain in 1:nsChain ){
resi <- res[[iChain]]
X[,iChain] <- resi$X
logDenCompX[,iChain] <- resi$logDenCompX
parUpdateDen[,,iChain] <- resi$parUpdateDen
accepted[,iChain] <- resi$accepted
Y[,,iChain] <- resi$Y
}
resDo <- list( ##describe<<
X=X ##<< matrix current position, column for each chain
, logDenCompX=logDenCompX ##<< matrix: result components of fLogDen current position, column for each chain
#, logDenX=logDenX ##<< vector current logDen of chains
, parUpdateDen=parUpdateDen ##<< numeric matrix (par x den x chain): for each parameter/density combination: is the density up to date
, accepted=accepted ##<< numeric matrix (blocks x chains): number of accepted steps for each chain
, Y=if( !is.null(res1$Y) ) Y else NULL ##<< numeric matrix (steps x 2+params+result x chains): accepted, rLogDen, parms, and all fLogDen result components for each proposal
, intermediatesChains = lapply(res,"[[","intermediatesX")
) ##end<<
}
.updateIntervalTwDEMCChain <- function(
### Perform update steps within next thinning interval for a single chain (called remotely from .updateIntervalTwDEMCParChains)
# first three arguments relate to current state (stored in chainState in twDEMCBlockInt)
argsUpdateIntervalTwDEMCChain = list( ##<< argument that change across invocations of .updateIntervalTwDEMCChain
## supplied as list, so that sfLapply can be used
##describe<<
X=c() ##<< numeric vector: current location
,logDenCompX=c()##<< numeric vector: result of calls to density functions for current location
,intermediates=list() ##<< list: intermediate results for current location
,parUpdateDen=c() ##<< for each parameter/density combination: is the density up to date from previous step
,xStep=c() ##<< matrix rows: difference vectors in parameter space, cols: steps
,rExtra=c() ##<< numeric vector: row: step within thinning interval
,tempDenCompSteps=c() ##<< numeric matrix: (nGenThin x nResComp)
#,nsPop ##<< the number of populations that take part in this step
,pAcceptPar=c() ##<< numeric vector: current acceptance rate for each (block)
,isRecordProposalsPop=FALSE ##<< logical
, iPop=1 ##<< the population that the current chain is in
)
##end<<
,argsUpdateBlocksTwDEMC ##<< list: further arguments to \code{\link{.updateBlocksTwDEMC}}
){
a <- argsUpdateIntervalTwDEMCChain # shortcut
# .updateIntervalTwDEMCPar
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{.doDEMCStep}}
dimXStep <- dim(a$xStep)
nStep <- dimXStep[2]
#d <- as.list(structure(dimXStep,names=c("parms","chains","steps")))
iGenT <- (1:nStep)
nCases = nStep
if( !(nCases == length(a$rExtra)) )
stop("number of cases in steps must correspond to length of rExtra")
res0 <- list(
xC=a$X
,logDenCompC=a$logDenCompX
,intermediatesC=a$intermediates
,parUpdateDenC=a$parUpdateDen
)
#
# record the proposals, their density results and their acceptance for each generation in thinning interval
parNames <- names(a$X)
resCompNames <- names(a$logDenCompX)
if( is.null(resCompNames) ) resCompNames <- paste("lDen",seq_along(a$logDenCompX),sep="")
tmp <- c(parNames,paste("accepted",seq_along(a$pAcceptPar),sep=""),resCompNames)
Y <- array( NA_real_, dim=c(nStep=nStep,comp=length(tmp)), dimnames=list(steps=NULL,comp=tmp) )
#
# recored acceptance state
acceptedSteps <- array( NA_real_, dim=c(nStep=nStep,block=length(a$pAcceptPar)) )
#
#i <- 1
prevRes <- res0
for( iStep in 1:nStep){
args <-c( prevRes[1:4]
,list( step=a$xStep[,iStep]
, rExtra=a$rExtra[iStep ,drop=TRUE]
#, tempGlobalC=tempGlobalSteps[iStep0+1,isPop ,drop=TRUE]
, tempDenCompC=a$tempDenCompSteps[iStep,, drop=TRUE]
, pAccept=a$pAcceptPar
#, isPop=isPop
, iPop=a$iPop
, argsUpdateBlocksTwDEMC=argsUpdateBlocksTwDEMC
)
)
prevRes <- do.call( .updateBlocksTwDEMC, args )
Y[iStep,] <- c( prevRes$xP, prevRes$accepted, prevRes$logDenCompP )
acceptedSteps[iStep,] <- prevRes$accepted
prevRes$accepted
}
#
#if( (iPop==2) && (X["a"] < 10.8) ) recover()
##value<< list with components
ctrl <- argsUpdateBlocksTwDEMC$argsFUpdateBlocks[[1]]$ctrl
resDo <- list( ##describe<<
X=prevRes$xC ##<< vector (nParm) current position
, logDenCompX=prevRes$logDenCompC ##<< vector (nResComp): result components of fLogDen current position
, intermediatesX=if( ctrl$returnIntermediate ) ##<< list: intermediate results for current position
prevRes$intermediatesC else list()
#, logDenX=logDenX ##<< vector current logDen of chains
, parUpdateDen=prevRes$parUpdateDenC ##<< numeric matrix (par x den): for each parameter/density combination: is the density up to date
, accepted=colSums(acceptedSteps) ##<< numeric vector (blocks): number of accepted steps for each chain
, Y=Y ##<< numeric matrix (steps x 2+params+result): accepted, rLogDen, parms, and all fLogDen result components for each proposal
) ##end<<
}
.updateBlocksTwDEMC <- function(
### One step of updating all blocks for one chain
xC ##<< current state
,logDenCompC ##<< result of calls to density functions of current state
,intermediatesC=list() ##<< intermediate results of current state
,parUpdateDenC ##<< matrix (nPar x nDen): FALSE: logDenCompC is not uptodate and needs recalculation
## (because depends on parameters that was updated against other density)
,step ##<< proposed step
,rExtra ##<< adjustment of Metropolis state
,tempDenCompC ##<< numeric vector of length(logDenCompC): current temperature for each result component
,pAcceptC ##<< numeric vector (nBlock) current acceptance rates of the blocks
,iPop ##<< population in total array of populations
,argsUpdateBlocksTwDEMC = list( ##<< further arguments that do not change between calls.
## These can be transferred by sfExport once and are handled by sfRemoteWrapper.
##describe<<
argsFUpdateBlocks=list() ##<< list(nBlock) arguments for each block: dInfo, block
,nBlock=1 ##<< scalar integer: the number of blocks
,upperParBoundsPop=c() ##<< upper parameter bounds
,lowerParBoundsPop=c() ##<< lower parameter bounds
##end<<
)
){
#if( !all(is.finite(pAcceptC)))
# stop(".updateBlocksTwDEMC: encountered non-finited pAcceptC")
a <- argsUpdateBlocksTwDEMC
#with( argsUpdateBlocksTwDEMC,{
acceptedBlock <- numeric(a$nBlock)
xP <- xC # place for recording proposal
logDenCompP <- logDenCompC # place for recording density results for proposal
#
#recover()
#iBlock<-1
#iBlock<-2
#dInfo1 <- a$argsFUpdateBlocks[[1]]
for( iBlock in seq_along(a$argsFUpdateBlocks)){
#block <- blocks[[iBlock]]
blockArgs <- a$argsFUpdateBlocks[[iBlock]] #merged the information in argsFUpdateBlocks
usesIntermediate <- length(blockArgs$intermediateId) != 0
#
argsFUpdateBlock <- c( blockArgs, list(
iPop=iPop
, upperParBounds=a$upperParBoundsPop[[iPop]]
, lowerParBounds=a$lowerParBoundsPop[[iPop]]
)) # will be extendend if Metropolis is used and by intermediate further down
#
##details<< \describe{}{\item{Recalculating logDensity of accepted state.}{
## If one of the components used by the current density function has been
## updated against another density function, then the current density of
## the accepted state is out of date, and needs to be recalculated.
##}}
if( blockArgs$useMetropolisUpdate ){
dInfo <- blockArgs$dInfo
cd <- dInfo$compPosDen # components used by the current density function
if( !all( parUpdateDenC[blockArgs$dInfoPos,cd] ) ){
# intermediate state may have been calculated, no need of recalculation
argsFLogDen <- dInfo$argsFLogDen
argsFLogDen$intermediate <- if( usesIntermediate ) intermediatesC[[blockArgs$intermediateId]] else NULL
# here do not allow for internal rejection
logDenCompC[dInfo$resCompPos] <- tmp <- .calcFLogDen( dInfo$fLogDen, xP, compInt=dInfo$logDenCompC[dInfo$intResCompPos], TiInt = dInfo$tempDenCompC[dInfo$intResCompPos]
,argsFLogDen = argsFLogDen
,maxLogDen = dInfo$maxLogDen
)
#do.call( dInfo$fLogDen, c(list(xC[cd]), dInfo$argsFLogDen) ) # evaluate logDen
#undebug(ofDalecPhen)
#tmp2 <- do.call(ofDalecPhen, c(list(xC[cd]), dInfo$argsFLogDen) ) # evaluate logDen
if( any(!is.finite(tmp)) ){
dump.frames("updateBlocksTwDEMC",TRUE)
stop(".updateBlocksTwDEMC (25): encountered non-finite resLogDen of in chainState. Dump.frames to file updateBlocksTwDEMC.rda. Implement better consistency checks across all fLogDen.")
}
parUpdateDenC[blockArgs$dInfoPos,cd] <- TRUE
if( usesIntermediate ){
intermediatesC[[ blockArgs$intermediateId ]] <- attr(tmp,"intermediate")
}
} # end if intermediate state changed
argsFUpdateBlock$metInfoStep <- list(
step=step
, rExtra=rExtra
, logDenCompC=logDenCompC[ dInfo$resCompPos ]
#, tempC=tempGlobalC
, tempDenCompC=tempDenCompC[ dInfo$resCompPos ]
, pAccept=pAcceptC
)
} # end if( blockArgs$useMetropolisUpdate )
#
# intermediate may have changed in Metropolis-Block
argsFUpdateBlock$intermediate <- if( usesIntermediate ) intermediatesC[[blockArgs$intermediateId]] else NULL
argsFUpdateBlock$intermediateId <- NULL # else this will be uses for intermediate: R bug
#
# make sure to use names different from already existing in argsFUpdateBlocks[[iBlock]]
# in c latter entries overwrite previous entries, so use argsFUpdateBlocks first
# updateBlockTwDEMC will set it to NULL for new parameter proposal
#resUpdate <- updateBlockTwDEMC( xC[cd], argsFUpdateBlock=argsFUpdateBlock )
# call the update
resUpdateB <- do.call( blockArgs$fUpdateBlock, c(list(xC[cd], argsFUpdateBlock=argsFUpdateBlock), argsFUpdateBlock$argsFUpdate) )
#resUpdate <- blockArgs$fUpdateBlock( xC[cd], argsFUpdateBlock=argsFUpdateBlock )
acceptedBlock[iBlock] <- resUpdateB$accepted # record acceptance
if( !all(is.finite(resUpdateB$accepted)))
stop(".updateBlocksTwDEMC: encountered non-finited acceptance rate")
if( any(!is.finite(resUpdateB$logDenCompC)) ){
stop(".updateBlocksTwDEMC (35): encountered non-finite resLogDen in accepted state")
}
cu <- blockArgs$compPos # components to be updated in this block
if( resUpdateB$accepted != 0){
#if( resUpdate$xC["a"] > 10.8) recover()
xC[cu] <- xP[cu] <- resUpdateB$xC # update the components of current state
if( usesIntermediate )
intermediatesC[[ blockArgs$intermediateId]] <- resUpdateB$intermediate # update the intermediate
#undebug(setIntermediate)
parUpdateDenC[,cu] <- FALSE # indicate all density results out of date for the updated parameters
if( 0 != length(resUpdateB$logDenCompC) ){
#if( iBlock==1 && resUpdate$logDenCompC[2] > -5) recover()
logDenCompC[ dInfo$resCompPos ] <- logDenCompP[ dInfo$resCompPos ] <- resUpdateB$logDenCompC # update the result of the used density function
parUpdateDenC[blockArgs$dInfoPos,cu] <- TRUE # if density has been calculated, indicate that it is up to date
}
}else{ # not accepted
if( 0 != length(resUpdateB$xP) ) xP[ cu ] <- resUpdateB$xP # proposal
if( 0 != length(resUpdateB$logDenCompP) ) logDenCompP[ dInfo$resCompPos ] <- resUpdateB$logDenCompP # density results for proposal
# keep with previous intermediate
}
} #iBlock
##value<< list with components
resUpdateStep <- list( ##describe <<
xC=xC ##<< numeric vecotr: current accepted state
,logDenCompC=logDenCompC ##<< numeric vector: result components of fLogDen for current position
,intermediatesC= intermediatesC ##<< the list of intermediate results of current state
,parUpdateDenC=parUpdateDenC ##<< integer vector: logDensity that recently updated parameter at given index
,accepted=acceptedBlock ##<< numeric vector: acceptance of each block (0: not, 1: did, (0..1): DRstep)
,xP=xP ##<< numeric vector: result components of fLogDen for proposal
#,intermediate
,logDenCompP=logDenCompP
) ##enf<<
#}) # with
}
### \item{tempC}{global temperature}
### \item{TFix}{named numeric vector (comp): components with fixed Temperature}
### \item{posTFix}{integer vector (comp): =match(TFix, compNames): positions of TFix within comp provided for performance reasons}
### \item{ctrl$useMultiT}{boolean wheter to scale temperature for different data streams}
### \item{TProp}{numeric vector of length(resCompPos): proportions of temperature for different data streams}
.calcFLogDen <- function(
## calculating the logDensity and adjusting results
fLogDen ##<< logDensity function
, xP ##<< parameter vector (first argument to fLogDen)
, compInt ##<< internal components
, TiInt ##<< numeric vector (nCompInt): Temperature for internal components
, argsFLogDen ##<< further named arguments to fLogDen
, maxLogDen ##<< integer vector (nResComp): maximum logDensity (overfitting control, usually -1/2 nObs)
, fLogDenScale=1 ##<< multiplicator to convert logDen result (maybe a cost) to true logDensity
){
l <- try( argsFLogDen ); if( inherits(l,"try-error")) recover()
LpOverfit <- fLogDenScale * if(length(compInt)){
do.call( fLogDen, c(list(xP, compInt/fLogDenScale, TiInt), argsFLogDen) ) # evaluate logDen
}else{
do.call( fLogDen, c(list(xP), argsFLogDen) ) # evaluate logDen
#a$fLogDen( xProp, blockIndices=a$argsFLogDen$blockIndices, fModel=a$argsFLogDen$fModel, obs=a$argsFLogDen$obs, invCovar=a$argsFLogDen$invCovar, xval=a$argsFLogDen$xval )
}
Lp <- if( length(maxLogDen) ) pmin( LpOverfit, maxLogDen ) else LpOverfit
##value<< numeric vector of logDensity, corrected for logDenScale and maximum LogDensity, i.e. overfitting control
}
#undebug(updateBlockTwDEMC)
updateBlockTwDEMC <- function(
### Block update function by a Metropolis decision
xC ##<< numeric vector: current state that is used in density function of the block
,argsFUpdateBlock
### further argument provided for generating the update \describe{
### \item{upperParBounds}{ named numeric vector, upper limits of the parameters to generate, see \code{\link{twDEMCBlockInt}} }
### \item{lowerParBounds}{ named numeric vector, lower limits of the parameters to generate, see \code{\link{twDEMCBlockInt}} }
### \item{iPop}{ the population for which the sample is generated (just for debugging) }
### \item{intermediate}{ intermediate result for current state xC, see end of vignette on using two Densities }
### }
# # # \item{fCalcComponentTemp}{ functiont to calculate temperature of result components, (way of transporting calcComponentTemp to remote process) }
){
#updateBlockTwDEMC
#
##details<<
## Function arguments above are provided to all block update functions.
## If Metropolis update is used, argument \code{argsFUpdateBlock} constains the following additional entries \describe{
## \item{compPosInDen}{ positions of the dimensions in x that are updated in this block}
## \item{ctrl$DRgamma}{ if given and >0 delayed Rejection (DR) (Haario06) is applied:
## If proposal is rejected, then a second proposal is tried jumping only DRgamma distance along the proposal }
## \item{dInfo}{ list on used logDensity function static troughout inversion with items \describe{
## \item{argsFLogDen}{additional arguments to fLogDen and scalar factor applied to result of fLogDen}
## \item{posLogDenInt}{the matching positions of intResCompNames within the the results components that are handled internally}
## \item{maxLogDen}{ integer vector (nResComp): maximum logDensity (overfitting control, usually -1/2 nObs) }
## \item{fDiscrProp,argsFDiscrProp}{function and additional arguments applied to xProp, e.g. to round it to discrete values}
## }}
## \item{metInfoStep}{ list on current Metropolis update variable troughout inversion with items \describe{
## \item{step}{proposed jump}
## \item{rExtra}{extra portion in Metropolis decision because of selecting the jump}
## \item{logDenCompC}{numeric vector: former result of call to same fLogDen}
## \item{tempDenCompC}{numeric vector of length(logDenCompC): temperature for each density result component}
## }}
## }
#
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{.updateBlocksTwDEMC}
# # \code{\link{.updateBlocksTwDEMC}}
#attach( argsDEMCStep )
#stop(".doDEMCStep: stop to trace error in remote R invocation.")
# collect information into a single big shorthand list for shorter code
a <- c(argsFUpdateBlock, argsFUpdateBlock$dInfo, argsFUpdateBlock$metInfoStep )
a$argsFLogDen$intermediate <- NULL # intermediate refers to current state, not to new proposal
#recover()
#with( argsFUpdateBlock, { # with produces significant performance loss
boResFLogDenX <- (length(a$intResCompPos) > 0)
# LogDensity of accepted state
#La <- logDenCompC #logDensity components of accepted state
if( 0 == length(a$logDenCompC)) stop("updateBlockTwDEMC: encountered empty logDenCompC")
#assume that all is.finite(logDenCompC), make sure in twDEMCBlockInt
LaExt <- La <- a$logDenCompC
#logDenC <- sum(La)
#a$tempDenCompC <- a$fCalcComponentTemp( temp=a$tempC, TFix=a$TFix, TProp=a$TProp, useMultiT=a$ctrl$useMultiT,posTFix=a$posTFix)
#
accepted<-FALSE
xP <- xC
xP[a$compPosInDen] <- xP[a$compPosInDen] + a$step[a$compPosInDen]
#
#if( a$iPop==2 && (a$lowerParBounds != 10.8) ) recover()
boOutside <-
any( sapply( names(a$upperParBounds), function(pname){ xP[pname] >= a$upperParBounds[pname] })) ||
any( sapply( names(a$lowerParBounds), function(pname){ xP[pname] < a$lowerParBounds[pname] }))
#if(xProp[1] < 10.8) recover()
#sapply( names(a$lowerParBounds), function(pname){ pname })
#sapply( names(a$lowerParBounds), function(pname){ xP[pname] })
#sapply( names(a$lowerParBounds), function(pname){ a$lowerParBounds[pname] })
if( boOutside ){
# bp(11)
# if it is still outside (maybe opposite border) reject step and give -Inf as logDenResult
#logDenP=logAlpha10=-Inf
logAlpha10=-Inf #logAlpha10 is log of the initial acceptance ratio for DR step (-Inf no chance of acceptance)
Lp <- a$logDenCompC #results for the proposal
Lp[] <- -Inf
}else{
# discrtize proposal
if( is.function(a$fDiscrProp)) xP = do.call(a$fDiscrProp,xP,a$argsFDiscrProp, quote=TRUE)
Lp <- .calcFLogDen( a$fLogDen, xP, compInt=a$logDenCompC[a$intResCompPos], TiInt = a$tempDenCompC[a$intResCompPos]
,argsFLogDen = a$argsFLogDen
,maxLogDen = a$maxLogDen
)
#take care that the result has always the same sames, even when if fails
#if( 0==length(names(res)))
# stop("encountered result of fLogDen without names")
#if( !identical(names(logDenCompC),names(res)))
# stop("encountered result with different names")
#strip attributes other than names, else twDynamicClusterApplyDep fails with big data chunks
#attributes(Lp) <- list(names=names(Lp))
#logDenP=-Inf
#make sure Lp, La have the same order and legnth
#if( !identical( names(Lp), names(La)) ) stop(".doDEMCStep: logDenCompC must contain the same components and the order of result of fLogDen." )
if( all(is.finite(Lp))){
#logDenP <- sum(Lp)
##details<< \describe{\item{internal Metropolis step}{
## if posLogDenInt is given, then these components of result of fLogDen are handled
## internally. Hence, for Metropolis step here operates only on other remaining components.
## }}
posTExt <- setdiff( seq_along(a$resCompPos), a$intResCompPos ) #externally handled components
nExt <- length(posTExt)
#posTFixExt <- setdiff(posTFix,posLogDenInt) #externally handled components with fixed temperature
#posTVarExt <- setdiff(seq_along(Lp), c(posTFix,posLogDenInt)) #externally handled componetns with variable temperature
#nFixExt <- length(posTFixExt)
#nVarExt <- length(posTVarExt)
#nExt <- nFixExt + nVarExt
#
#Metropolis step
#logr = (logDenPropExt+rExtra - logDenXExt) / tempC
logrDS10 <- (Lp[posTExt]-La[posTExt])/a$tempDenCompC[posTExt]
logAlpha10 <- a$rExtra + sum(logrDS10)
accepted <- if( is.finite(logAlpha10) && (logAlpha10 > log(runif(1)) ) ) 1 else 0
if(accepted != 0){
xC <- xP
#logDenC <- logDenP
a$logDenCompC <- Lp
}
}else logAlpha10 <- -Inf
} # end check outside parBounds
# DR-step
if(!accepted && !is.null(a$ctrl$DRgamma) && (a$ctrl$DRgamma > 0) &&
( boOutside || (!is.null(a$ctrl$minPCompAcceptTempDecr) && (a$pAccept < 1.2*a$ctrl$minPCompAcceptTempDecr)))
) {
#----- delayed rejection (DR) step
# only if across parBoundEdge or acceptance rate drops below 1.2*minAcceptrate
# repeat all above with delayed rejection (DR) step, only adjust DRfac after calculating Lp
Lp1 <- Lp
xProp1 <- xP # former proposal
xP <- xC
xP[a$compPosInDen] <- xP[a$compPosInDen] + a$ctrl$DRgamma*a$step[a$compPosInDen]
boOutside <-
any( sapply( names(a$upperParBounds), function(pname){ xP[pname] > a$upperParBounds[pname] })) ||
any( sapply( names(a$lowerParBounds), function(pname){ xP[pname] < a$lowerParBounds[pname] }))
if( !boOutside ){
if( is.function(a$fDiscrProp)) xP = do.call(a$fDiscrProp,xP,a$argsFDiscrProp, quote=TRUE)
Lp <- .calcFLogDen( a$fLogDen, xP, compInt=a$logDenCompC[a$intResCompPos], TiInt = a$tempDenCompC[a$intResCompPos]
,argsFLogDen = a$argsFLogDen
,maxLogDen = a$maxLogDen
)
#take care that the result has always the same sames, even when if fails
#if( 0==length(names(res)))
# stop("encountered result of fLogDen without names")
#if( !identical(names(logDenCompC),names(res)))
# stop("encountered result with different names")
#strip attributes other than names, else twDynamicClusterApplyDep fails with big data chunks
attributes(Lp) <- list(names=names(Lp))
#logDenP=-Inf
#make sure Lp, La have the same order and legnth
#if( !identical( names(Lp), names(La)) ) stop(".doDEMCStep: logDenCompC must contain the same components and the order of result of fLogDen." )
if( all(is.finite(Lp))){
#logDenP <- sum(Lp)
posTExt <- setdiff( seq_along(a$tempDenCompC), a$intResCompPos ) #externally handled components
nExt <- length(posTExt)
#posTFixExt <- setdiff(posTFix,posLogDenInt) #externally handled components with fixed temperature
#posTVarExt <- setdiff(seq_along(Lp), c(posTFix,posLogDenInt)) #externally handled componetns with variable temperature
#nFixExt <- length(posTFixExt)
#nVarExt <- length(posTVarExt)
#nExt <- nFixExt + nVarExt
#Metropolis step
#logr = (logDenPropExt+rExtra - logDenXExt) / tempC
logrDS20 <- (Lp[posTExt]-La[posTExt])/a$tempDenCompC[posTExt]
logAlpha20 <- a$rExtra + sum(logrDS20)
#--- here correct with first stage DR factor (1-alpha21)/(1-alpha10) with meaning 0:accepted 1:first proposal 2:second proposal
logrDS21 <- (Lp1[posTExt]-Lp[posTExt])/a$tempDenCompC[posTExt]
logAlpha21 <- sum(logrDS21)
logAlpha2 <- suppressWarnings( logAlpha20 +log(1-exp(logAlpha21)) -log(1-exp(logAlpha10)) ) # log and exp may produce NaNs
accepted <- if( is.finite(logAlpha2) && ( logAlpha2 > log(runif(1)) ) ) a$ctrl$DRgamma else 0
if(accepted != 0){
xC <- xP
#logDenC <- logDenP
a$logDenCompC <- Lp # also transferres the "intermediate" attribute
}
}
} # end !boOutside in DR step
} # end DR step
#will invoke prevRes[c("x", "logDenCompC", "logDenAcc")]
if(!is.numeric(xC) | !is.numeric(a$logDenCompC) ) #| !is.numeric(logDenC) )
stop(".twDEMCBlocks: x, logDenAcc and logDenCompC must be numeric")
##value<< list with components
list( ##describe<<
accepted=accepted ##<< boolean scalar: if step was accepted
, xC=xC[a$compPosInDen] ##<< numeric vector: components of position in parameter space that are being updated (if accepted then the same as xP)
, logDenCompC=a$logDenCompC ##<< numeric vector: result components of fLogDen for current position (if accepted then the same as logDenCompP)
, intermediate=attr(a$logDenCompC,"intermediate") ##<< intermediate result of current state (if provided by logDensity)
#, logDenC =logDenC ##<< numeric vector: summed fLogDen for current accepted position
, xP=xP[a$compPosInDen] ##<< numeric vector: components of proposal that are being
, logDenCompP=Lp ##<< numeric vector: result components of fLogDen for proposal, may inlcude attr "intermediate"
#, logDenP=logDenP ##<< numeric vector: summed fLogDen for proposal
) ##end<<
#})
}
R.methodsS3::setMethodS3("twDEMCBlock","array", function(
### Initialize \code{\link{twDEMCBlockInt}} by array of initial population and remove those generations from results afterwards
x ##<< initial population: a numeric array (M0 x d x nChain) see details in \code{\link{twDEMCBlockInt}}
,... ##<< further arguments to \code{\link{twDEMCBlockInt}}
,nPop = 1 ##<< number of populations in x
#,TProp = NULL ##<< numeric vector (nResComp): proportions of temperature for result components: equal for all populations
, X=NULL ##<< initial state (nParm x nChain)
, logDenCompX = NULL ##<< numeric matrix (nResComp x nChains) initial state
, upperParBounds = vector("list",nPop)
### list of named numeric vectors: giving upper parameter bounds for each population
### for exploring subspaces of the limiting distribution, see details
### , Alternatively a single numeric vector can be supplied, which is replicated for each population.
, lowerParBounds = vector("list",nPop)
### similar to upperParBounds
){
#twDEMCBlock.array
#if( 1 == length(T0) ) T0 <- rep(T0,nPop)
#if( 1 == length(TEnd) ) TEnd <- rep(TEnd,nPop)
if( is.numeric(upperParBounds) ) upperParBounds <- lapply(1:nPop, function(iPop) upperParBounds )
if( is.numeric(lowerParBounds) ) lowerParBounds <- lapply(1:nPop, function(iPop) lowerParBounds )
nChains <- dim(x)[3]
nChainPop <- nChains %/% nPop
isX <- (0 != length(X))
isLogDenCompX <- (0 != length(logDenCompX))
#isTProp <- (0 != length(TProp))
#if( isTProp && !is.matrix(TProp) ) TProp <- matrix( TProp, nrow=length(TProp), ncol=nPop, dimnames=list(resComp=names(TProp),pop=NULL) )
pops <- lapply( 1:nPop, function(iPop){
chainsPopI <- (iPop-1)*nChainPop + 1:nChainPop
pop <- list(
parms=x[,,chainsPopI, drop=FALSE]
#,T0=T0[iPop]
#,TEnd=TEnd[iPop]
,upperParBounds=upperParBounds[[iPop]]
,lowerParBounds=lowerParBounds[[iPop]]
)
if( isX ) pop$X <- X[,chainsPopI , drop=FALSE]
if( isLogDenCompX ) pop$logDenCompX <- logDenCompX[,chainsPopI ,drop=FALSE]
#if( isTProp ) pop$TProp <- TProp[,iPop ,drop=TRUE]
#pop$spaceInd <- spacePop[iPop]
pop
})
#res <- twDEMCBlockInt(pops)
res <- twDEMCBlockInt(pops=pops,...)
#.dots <- list(...)
#do.call( twDEMCBlockInt, c(list(pops=pops),.dots))
#resc <- concatPops.twDEMCPops(res) # all have the same length, so allow concatenate population results
})
attr(twDEMCBlock.array,"ex") <- function(){
data(twLinreg1)
attach( twLinreg1 )
argsFLogDen <- list(
fModel=dummyTwDEMCModel, ### the model function, which predicts the output based on theta
obs=obs, ### vector of data to compare with
invCovar=invCovar, ### the inverse of the Covariance of obs (its uncertainty)
thetaPrior = thetaTrue, ### the prior estimate of the parameters
invCovarTheta = invCovarTheta, ### the inverse of the Covariance of the prior parameter estimates
xval=xval
)
do.call( logDenGaussian, c(list(theta=theta0),argsFLogDen))
.nPop = 2
Zinit <- initZtwDEMCNormal( theta0, diag(sdTheta^2), nChainPop=4, nPop=.nPop)
dim(Zinit)
.nGen=100
#nGen=3
#mtrace(twDEMC.array)
#mtrace(.updateIntervalTwDEMCPar)
#mtrace(twDEMCBlockInt)
tmp1 <- tmp <- twDEMCBlock( Zinit, nPop=.nPop
,dInfos=list(list(fLogDen=logDenGaussian, argsFLogDen=argsFLogDen))
,nGen=.nGen
)
plot( as.mcmc.list(tmp), smooth=FALSE )
tmp2 <- tmp <- twDEMCBlock( tmp1, nGen=200 )
str(tmp)
}
.concatTwDEMCRuns <- function(
x ##<< initial run (of class twDEMCPops)
,res ##<< extended run
, doRecordProposals=FALSE ##<< if TRUE then an array of each proposal together with the results of fLogDen are recorded and returned in component Y
){
##details<<
## initial populations are appended by results
# iPop = length(x$pops)
for( iPop in 1:length(x$pops) ){
xPop <- x$pops[[iPop]]
resPop <- res$pops[[iPop]]
if( nrow(resPop$parms) < 2){
# no samples to append
res$pops[[iPop]] <- xPop
}else{
res$pops[[iPop]]$parms <- abind(xPop$parms, resPop$parms[-1,, ,drop=FALSE], along=1)
res$pops[[iPop]]$temp <- abind(xPop$temp, resPop$temp[-1, ,drop=FALSE], along=1 )
res$pops[[iPop]]$pAccept <- abind(xPop$pAccept, resPop$pAccept[-1,, ,drop=FALSE], along=1)
res$pops[[iPop]]$resLogDen <- abind( xPop$resLogDen, resPop$resLogDen[-1,, ,drop=FALSE], along=1)
res$pops[[iPop]]$logDen <- abind( xPop$logDen, resPop$logDen[-1,, ,drop=FALSE], along=1)
nrY <- nrow(resPop$Y)
if( doRecordProposals || (nrY < 128) ){
# if new Y has less than 128 rows append previous Y
res$pops[[iPop]]$Y <- abind( xPop$Y, resPop$Y, along=1)
nrY <- nrow(res$pops[[iPop]]$Y) # rows did change
if( !doRecordProposals && (nrY > 1) )
res$pops[[iPop]]$Y <- res$pops[[iPop]]$Y[ 1:min(128,nrY),, ,drop=FALSE] # cut to 128 cases
}
} # if res has samples
}
# copy other items present in x to result
namesAdd <- names(x)[ !(names(x) %in% names(res)) ]
res[ namesAdd ] <- x[ namesAdd ]
res
}
R.methodsS3::setMethodS3("twDEMCBlock","twDEMCPops", function(
### initialize \code{\link{twDEMCBlockInt}} by former run and append results to former run
x ##<< list of class twDEMCPops, result of \code{\link{twDEMCBlockInt}}
,... ##<< further arguments to \code{\link{twDEMCBlockInt}}
, TEnd=numeric(0) ##<< numeric vector (nResComp) of end temperatures for each result component, if not given stays at current temperature, if in addition TSpec is specified, it has no effect
, TEnd0 ##<< numeric scalar: alternative way of specifying end tempertaure: by baseTemperature for all data streams
, TStart=numeric(0) ##<< numeric vector (nResComp): specifing starting temperature, if neigher TSpec, TStart, or TStart0 are given stays at the current temperature
, TStart0 ##<< numeric scalar: specifing starting temperature by baseTemperature for all data streams
, TFix ##<< see \code{\link{twDEMCSACont}}, argument \code{ctrlT$TFix}. Must be given if specifying TEnd0 or TStart0
, nObs ##<< see \code{\link{twDEMCSACont}}, argument must be given if specifying TEnd0 or TStart0
, doRecordProposals=FALSE ##<< if TRUE then an array of each proposal together with the results of fLogDen are recorded and returned in component Y
, extendRun=TRUE ##<< if set to FALSE then only the new samples are returned
){
#twDEMCBlock.twDEMCPops
##details<<
## pops, dInfos, and blocks are reused from x or overwritten by arguments
nSamplePop <- getNSamples(x) # used in several cases
#
.dots <- argsList <- list(...)
argsList$pops <- x$pops
argsList$doRecordProposals <- doRecordProposals
if( 0 == length(.dots$dInfos) )
argsList$dInfos <- x$dInfos
if( 0 == length(.dots$blocks) )
argsList$blocks <- x$blocks
#for( iPop in seq_along(x$pops) ){
# if( 0 == length(argsList$pops[[iPop]]$T0) )
# argsList$pops[[iPop]]$T0 <-
# x$pops[[iPop]]$temp[ nSamplePop[iPop] ]
#}
if( 0 == length(argsList$TSpec) ){
##details<< \describe{\item{TSpec}{
## If TSpec is not explicitly given as an array for \code{\link{twDEMCBlockInt}} it is constructed as follows:
## , If TEnd is given use it to explicitely define end Temperature for all data streams.
## , If TEnd0 is given calculate TEnd for all data streams by scaling for number of observations (see \code{\link{calcStreamTemp}}).
## , If TStart0 is given calcualte TCurr for all data streams
## , else use current Temperature of streams in x as starting temperature
## }}
if( (!missing(TStart0)) || (!missing(TEnd0)) ){
if( missing(nObs) || missing(TFix) ) stop("twDEMCBlock.twDEMCPops: if TStart0 or TENd0 are specified, one needs to specify nObs and TFix too.")
resCompNames <- colnames(x$pops[[1]]$resLogDen)
nObsA <- nObs[ resCompNames ]
TFix <- completeResCompVec( TFix, resCompNames )
iFixTemp <- which( is.finite(TFix) )
}
TCurr <- if( (!missing(TStart))){
TStart
}else if( (!missing(TStart0))){
calcStreamTemp(TStart0, nObsA , TFix, iFixTemp )
}else{
x$pops[[1]]$temp[ nrow(x$pops[[1]]$temp), ]
}
if( 0 == length(TEnd) )
TEnd <- if( !missing(TEnd0)){
calcStreamTemp(TEnd0, nObsA , TFix, iFixTemp )
} else TCurr
if( 1 == length(TEnd) ) TEnd <- rep(TEnd, length(TCurr) )
argsList$TSpec <- cbind( T0=TCurr, TEnd=TEnd )
}
#
if( 0 == length(.dots$controlTwDEMC) )
argsList$controlTwDEMC <- list()
if( 0 == length(argsList$controlTwDEMC$initialAcceptanceRate) ){
ar <- abind( lapply( seq_along(x$pops), function(iPop){
adrop( x$pops[[iPop]]$pAccept[ nSamplePop[iPop],, ,drop=FALSE ],1)
}),along=2)
argsList$controlTwDEMC$initialAcceptanceRate <- ar
}
res <- res0 <- do.call( twDEMCBlockInt, argsList )
if( extendRun ) res <- .concatTwDEMCRuns(x,res0,doRecordProposals=doRecordProposals)
res
})
attr(twDEMCBlock.twDEMCPops,"ex") <- function(){
data(twdemcEx1) # previous run of twDEMCBlock
class(twdemcEx1)
twdemcEx1$thin # thinning interval
(nGen0 <- getNGen(twdemcEx1)) # number of generations
# extend by 16 generations
nGen <- 16
#mtrace(twDEMCBlock.twDEMCPops)
res <- twDEMCBlock( twdemcEx1, nGen=nGen )
identical( nGen0+nGen, getNGen(res) )
plot( as.mcmc.list(res), smooth=FALSE )
}
R.methodsS3::setMethodS3("twDEMCBlock","twDEMC", function(
### initialize \code{\link{twDEMCBlockInt}} by former run and append results to former run
x, ##<< list of class twDEMC, result of \code{\link{twDEMCBlockInt}}
... ##<< further arguments to \code{\link{twDEMCBlockInt}}
){
#twDEMC.twDEMC
.dots <- list(...)
argsList <- list(x=x$parms)
M0 <- nrow(x$rLogDen)
#extract X, logDenX, and logDenCompX from parms, but only if not given with \dots
if( is.null(.dots$X) ) #use is.null, because if provided zero length vector, we want to use it
argsList$X <- adrop(x$parms[M0,,,drop=FALSE],2)
if( is.null(.dots$logDenX) )
argsList$logDenX <- x$rLogDen[M0,,drop=TRUE]
if( is.null(.dots$logDenCompX) )
argsList$logDenCompX <- adrop(x$logDenComp[M0,,,drop=FALSE],2)
if( is.null(.dots$nPop) )
argsList$nPop <- getNPops.twDEMC(x)
res <- do.call( twDEMCBlock.array, c(argsList,.dots))
res$rLogDen[1:M0,] <- x$rLogDen
res$logDenComp[,1:M0,] <- x$logDenComp
res$pAccept[1:M0,] <- x$pAccept
res$temp[1:M0,] <- x$temp
res
### parms appended with the further generations of \code{\link{twDEMCBlockInt}}
})
#mtrace(twDEMC.twDEMC)
#mtrace(twDEMCBlockInt)
#twUtest(twDEMC,"test.ZinittwDEMC")
Try the twDEMC package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.