Nothing
R/MCMC_conjugacy.R
In nimble: MCMC, Particle Filtering, and Programmable Hierarchical Modeling
Defines functions
dynamicConjugateSamplerWrite
dynamicConjugateSamplerGet
dynamicConjugateSamplerAdd
dynamicConjugateSamplerExists
makeIndexedVariable
makeDeclareSizeField
createDynamicConjugateSamplerName
compareConjugacyLists
cc_checkStickbreaking
cc_combineExprsDivision
cc_combineExprsMultiplication
cc_combineExprsSubtraction
cc_combineExprsAddition
cc_negateExpr
cc_checkLinearity
cc_replace01
cc_vectorizedComponentCheck
cc_getNodesInExpr
cc_nodeInExpr
cc_otherParamsCheck
cc_stripExpr
cc_checkScalar
cc_linkCheck
cc_createStructureExpr
cc_expandDetermNodesInExpr
cc_makeRDistributionName
cc_makeConjugateSamplerName
cc_makeSamplerTypeName
Documented in
cc_getNodesInExpr
conjugacyRelationshipsInputList <- list(
## beta
list(prior = 'dbeta',
link = 'identity',
dependents = list(
dbern = list(param = 'prob', contribution_shape1 = 'value', contribution_shape2 = '1 - value' ),
dbin = list(param = 'prob', contribution_shape1 = 'value', contribution_shape2 = 'size - value'),
dnegbin = list(param = 'prob', contribution_shape1 = 'size', contribution_shape2 = 'value' ),
## using 'stickbreaking' not 'stick_breaking' as link is now appended to dist with a '_' in realized conjugacy system
## Note that 'stick_breaking' is still the user-facing **function** name.
dcat = list(param = 'prob', link = 'stickbreaking', contribution_shape1 = 'value == offset',
contribution_shape2 = 'value > offset')), ## offset is set to be where in the broken stick the target is
posterior = 'dbeta(shape1 = prior_shape1 + contribution_shape1,
shape2 = prior_shape2 + contribution_shape2)'),
## dirichlet
list(prior = 'ddirch',
link = 'identity',
dependents = list(
dmulti = list(param = 'prob', contribution_alpha = 'value'),
dcat = list(param = 'prob', contribution_alpha = 'calc_dcatConjugacyContributions(k, value)')),
posterior = 'ddirch(alpha = prior_alpha + contribution_alpha)'),
## flat
list(prior = 'dflat',
link = 'linear',
dependents = list(
dnorm = list(param = 'mean', contribution_mean = 'coeff * (value-offset) * tau', contribution_tau = 'coeff^2 * tau'),
dlnorm = list(param = 'meanlog', contribution_mean = 'coeff * (log(value)-offset) * taulog', contribution_tau = 'coeff^2 * taulog')),
posterior = 'dnorm(mean = contribution_mean / contribution_tau,
sd = contribution_tau^(-0.5))'),
## halfflat - first possible conjugacy; user can achieve this with gamma(1, 0) as we handle that
## list(prior = 'dhalfflat',
## link = 'multiplicative',
## dependents = list(
## dpois = list(param = 'lambda', contribution_shape = 'value', contribution_rate = 'coeff' ),
## dnorm = list(param = 'tau', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (value-mean)^2' ),
## dlnorm = list(param = 'taulog', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (log(value)-meanlog)^2'),
## dgamma = list(param = 'rate', contribution_shape = 'shape', contribution_rate = 'coeff * value' ),
## dinvgamma = list(param = 'scale', contribution_shape = 'shape', contribution_rate = 'coeff / value' ),
## dexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * value' ),
## dweib = list(param = 'lambda', contribution_shape = '1', contribution_rate = 'coeff * value^shape' )),
## ## ddexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * abs(value-location)' )
## ## dpar = list(...) ## contribution_shape=1; contribution_rate=coeff*log(value/c) 'c is 2nd param of pareto'
## posterior = 'dgamma(shape = 1 + contribution_shape,
## scale = 1 / contribution_rate)'),
## halfflat - second possible conjugacy; note that user will be able to do invgamma(-1,0) and
## get conjugacy once we fix issue #314
## list(prior = 'dinvgamma',
## link = 'multiplicative',
## dependents = list(
## dnorm = list(param = 'var', contribution_shape = '1/2', contribution_scale = '(value-mean)^2 / (coeff * 2)' ),
## dlnorm = list(param = 'varlog', contribution_shape = '1/2', contribution_scale = '(log(value)-meanlog)^2 / (coeff*2)'),
## dgamma = list(param = 'scale', contribution_shape = 'shape', contribution_scale = 'value / coeff' ),
## dinvgamma = list(param = 'rate', contribution_shape = 'shape', contribution_scale = '1 / (coeff * value)' ),
## dexp = list(param = 'scale', contribution_shape = '1', contribution_scale = 'value / coeff' )),
## ## add ddexp
## posterior = 'dinvgamma(shape = -1 + contribution_shape,
## rate = 1 / contribution_scale)'),
## halfflat - third possible conjugacy - we use this because it corresponds to the
## Gelman (2006) recommended uniform on sd scale prior for variance components
## and current NIMBLE conjugacy system only allows one possible form of conjugacy.
## Note that sd ~ U(0,Inf) equivalent to var ~ IG(-1/2, 0).
## Also note that if conj system could detect 'squared' dependency, then
## we could allow dnorm with param = 'var'.
list(prior = 'dhalfflat',
link = 'multiplicative',
dependents = list(
dnorm = list(param = 'sd', contribution_shape = '1/2', contribution_scale = '(value-mean)^2 / (coeff^2 * 2)' ),
dlnorm = list(param = 'sdlog', contribution_shape = '1/2', contribution_scale = '(log(value)-meanlog)^2 / (coeff^2 * 2)')),
posterior = 'dsqrtinvgamma(shape = -1/2 + contribution_shape,
rate = 1 / contribution_scale)'),
## gamma
list(prior = 'dgamma',
link = 'multiplicative',
dependents = list(
dpois = list(param = 'lambda', contribution_shape = 'value', contribution_rate = 'coeff' ),
dnorm = list(param = 'tau', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (value-mean)^2' ),
dlnorm = list(param = 'taulog', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (log(value)-meanlog)^2'),
dgamma = list(param = 'rate', contribution_shape = 'shape', contribution_rate = 'coeff * value' ),
dinvgamma = list(param = 'scale', contribution_shape = 'shape', contribution_rate = 'coeff / value' ),
dexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * value' ),
dweib = list(param = 'lambda', contribution_shape = '1', contribution_rate = 'coeff * value^shape' ),
ddexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * abs(value-location)' )),
## dpar = list(...) ## contribution_shape=1; contribution_rate=coeff*log(value/c) 'c is 2nd param of pareto'
posterior = 'dgamma(shape = prior_shape + contribution_shape,
scale = 1 / (prior_rate + contribution_rate))'),
## invgamma
list(prior = 'dinvgamma',
link = 'multiplicative',
dependents = list(
dnorm = list(param = 'var', contribution_shape = '1/2', contribution_scale = '(value-mean)^2 / (coeff * 2)' ),
dlnorm = list(param = 'varlog', contribution_shape = '1/2', contribution_scale = '(log(value)-meanlog)^2 / (coeff*2)'),
dgamma = list(param = 'scale', contribution_shape = 'shape', contribution_scale = 'value / coeff' ),
dinvgamma = list(param = 'rate', contribution_shape = 'shape', contribution_scale = '1 / (coeff * value)' ),
dexp = list(param = 'scale', contribution_shape = '1', contribution_scale = 'value / coeff' ),
ddexp = list(param = 'scale', contribution_shape = '1', contribution_scale = 'abs(value-location) / coeff' )),
posterior = 'dinvgamma(shape = prior_shape + contribution_shape,
rate = 1 / (prior_scale + contribution_scale))'),
## normal
list(prior = 'dnorm',
link = 'linear',
dependents = list(
dnorm = list(param = 'mean', contribution_mean = 'coeff * (value-offset) * tau', contribution_tau = 'coeff^2 * tau' ),
dlnorm = list(param = 'meanlog', contribution_mean = 'coeff * (log(value)-offset) * taulog', contribution_tau = 'coeff^2 * taulog')),
posterior = 'dnorm(mean = (prior_mean*prior_tau + contribution_mean) / (prior_tau + contribution_tau),
sd = (prior_tau + contribution_tau)^(-0.5))'),
#####
## pareto
## these are idiosyncratic enough that we probably want to skip them
## list(prior = 'dpar', ##### waiting for dpar() distribution
## link = 'multiplicative',
## dependents = list(
## dunif = list(param = 'max', contribution_alpha = '1', contribution_c = 'value/coeff'), # only works if 0 is min of the unif so this will be hard to do
## careful with next one - data seem to impose upper bound on posterior, so not clear this goes through
## dpar = list(param = 'c', contribution_alpha = '-alpha'), contribution_c = '0'),
## posterior = 'dpar(alpha = prior_alpha + contribution_alpha,
## c = max(prior_c, contribution_c)'),
#####
## multivariate-normal; note that to avoid unnecessary computations, there is special processing of the contribution_mean and contribution_prec
## in the generation of the conjugate sampler run method, so the specification here should not be changed without looking at that processing.
list(prior = 'dmnorm',
link = 'linear',
dependents = list(
##dmnorm = list(param = 'mean', contribution_mean = '(t(coeff) %*% prec %*% asCol(value-offset))[,1]', contribution_prec = 't(coeff) %*% prec %*% coeff')),
dmnorm = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)')),
## LINK will be replaced with appropriate link via code processing
## original less efficient posterior definition:
## posterior = 'dmnorm_chol(mean = (inverse(prior_prec + contribution_prec) %*% (prior_prec %*% asCol(prior_mean) + asCol(contribution_mean)))[,1],
## cholesky = chol(prior_prec + contribution_prec),
## prec_param = 1)'),
posterior = '{ R <- chol(prior_prec + contribution_prec)
A <- prior_prec %*% asCol(prior_mean) + asCol(contribution_mean)
mu <- backsolve(R, forwardsolve(t(R), A))[,1]
dmnorm_chol(mean = mu, cholesky = R, prec_param = 1) }'),
## wishart
list(prior = 'dwish',
## changing to only use link='identity' case, since the link='linear' case was not correct.
## -DT March 2017
## link = 'linear',
link = 'multiplicativeScalar', # we only handle scalar 'coeff'; this naming is slightly awkward since for univar dists, link is of course scalar
dependents = list(
## parentheses added to the contribution_R calculation:
## colVec * (rowVec * matrix)
## Chris is checking to see whether this makes a difference for Eigen
## -DT April 2016
## changing to only use link='identity' case, since the link='linear' case was not correct
## -DT March 2017
## dmnorm = list(param = 'prec', contribution_R = 'asCol(value-mean) %*% (asRow(value-mean) %*% coeff)', contribution_df = '1')),
dmnorm = list(param = 'prec', contribution_R = 'coeff * asCol(value-mean) %*% asRow(value-mean)', contribution_df = '1')),
posterior = 'dwish_chol(cholesky = chol(prior_R + contribution_R),
df = prior_df + contribution_df,
scale_param = 0)'),
## inverse wishart
list(prior = 'dinvwish',
link = 'multiplicativeScalar', # we only handle scalar 'coeff'; this naming is slightly awkward since for univar dists, link is of course scalar
dependents = list(
dmnorm = list(param = 'cov', contribution_S = 'asCol(value-mean) %*% asRow(value-mean) / coeff', contribution_df = '1')),
posterior = 'dinvwish_chol(cholesky = chol(prior_S + contribution_S),
df = prior_df + contribution_df,
scale_param = 1)')
)
##############################################################################################
##############################################################################################
## reference class definitions:
## conjugacyRelationshipsClass
## conjugacyClass
## dependentClass
## posteriorClass
##############################################################################################
##############################################################################################
conjugacyRelationshipsClass <- setRefClass(
Class = 'conjugacyRelationshipsClass',
fields = list(
conjugacys = 'ANY' ## a (named) list of conjugacyClass objects, each describes the conjugacies for a particular prior distribution (name is prior distribution name)
),
methods = list(
initialize = function(crl) {
conjugacys <<- list()
for(i in seq_along(crl)) {
conjugacys[[i]] <<- conjugacyClass(crl[[i]])
}
names(conjugacys) <<- unlist(lapply(conjugacys, function(cr) cr$prior))
},
checkConjugacy = function(model, nodeIDs, restrictLink = NULL) {
if(isTRUE(getNimbleOption("oldConjugacyChecking")))
checkConjugacy_original(model, nodeIDs, restrictLink)
else
checkConjugacy_new(model, nodeIDs, restrictLink)
},
checkConjugacy_new = function(model, nodeIDs, restrictLink = NULL) {
maps <- model$modelDef$maps
nodeDeclIDs <- maps$graphID_2_declID[nodeIDs] ## declaration IDs of the nodeIDs
declID2nodeIDs <- split(nodeIDs, nodeDeclIDs) ## nodeIDs grouped by declarationID
ansList <- ansList2 <- list()
for(i in seq_along(declID2nodeIDs)) { ## For each group of nodeIDs from the same declarationID
nodeIDsFromOneDecl <- declID2nodeIDs[[i]]
firstNodeName <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[1]]
if(model$isTruncated(firstNodeName)) next ## we say non-conjugate if the targetNode is truncated
dist <- model$getDistribution(firstNodeName)
conjugacyObj <- conjugacys[[dist]]
if(is.null(conjugacyObj)) next
# NO: insert logic here to check a single dependency and do next if can't be conjugate
#model$getDependencies('mu[1]',self=F,stochOnly=T)
#for( loop through deps )
# conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode)
# if(not conj) next
# now try to guess if finding paths will be more intensive than simply looking at target-dependent pairs, to avoid path finding when there is nested structure such as stickbreaking
deps <- lapply(nodeIDsFromOneDecl, function(x) model$getDependencies(x, stochOnly = TRUE, self = FALSE))
numDeps <- sapply(deps, length)
sumNumDeps <- sum(numDeps)
maxNumPaths <- 0
## We only need check until we get to as many paths as `sumNumDeps`.
## This also avoids integer overflow that can occur with recursive indexing.
for(nodeID in nodeIDsFromOneDecl) {
numPaths <- model$getDependencyPathCountOneNode(nodeID, max = sumNumDeps + 1)
if(numPaths > maxNumPaths)
maxNumPaths <- numPaths
if(maxNumPaths > sumNumDeps)
break
}
if(maxNumPaths > sumNumDeps) {
# max(numPaths) is reasonable guess at number of unique (by node) paths (though it overestimates number of unique (by declaration ID) paths; if we have to evaluate conjugacy for more paths than we would by simply looking at all pairs of target-dependent nodes, then just use node pairs
# note that it's not clear what criterion to use here since computational time is combination of time for finding all paths and then for evaluating conjugacy for unique (by declaration ID) paths, but the hope is to make a crude cut here that avoids path calculations when there would be a lot of them
ansList[[length(ansList)+1]] <- lapply(seq_along(nodeIDsFromOneDecl),
function(index) {
targetNode <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[index]]
depEnds <- deps[[index]]
depTypes <- sapply(depEnds, function(x) conjugacyObj$checkConjugacyOneDep(model, targetNode, x, restrictLink))
if(!length(depTypes)) return(NULL)
if(!any(sapply(depTypes, is.null))) {
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
depEnds[boolMatch]
})
names(control) <- uniqueDepTypes
return(list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control))
} else return(NULL)
})
names(ansList[[length(ansList)]]) <- maps$graphID_2_nodeName[nodeIDsFromOneDecl]
} else {
# determine conjugacy based on unique (by declaration ID) paths
depPathsByNode <- lapply(nodeIDsFromOneDecl, model$getDependencyPaths) ## make list (by nodeID) of lists of paths through graph
depPathsByNode <- depPathsByNode[!unlist(lapply(depPathsByNode, function(x) is.null(x) || (length(x)==0)))]
depPathsByNodeLabels <- lapply(depPathsByNode, function(z) ## make character labels that match for same path through graph
unlist(lapply(z,
function(x)
paste(maps$graphID_2_declID[x[,1]], x[,2], collapse = '\r', sep='\r'))))
depPathsByNodeUnlisted <- unlist(depPathsByNode, recursive = FALSE)
depPathsByNodeLabelsUnlisted <- unlist(depPathsByNodeLabels)
## uniquePaths <- unique(depPathsByNodeLabelsUnlisted)
uniquePathsUnlistedIndices <- split(seq_along(depPathsByNodeLabelsUnlisted), depPathsByNodeLabelsUnlisted)
conjDepTypes <- character(length(uniquePathsUnlistedIndices))
for(j in seq_along(uniquePathsUnlistedIndices)) {
firstDepPath <- depPathsByNodeUnlisted[[ uniquePathsUnlistedIndices[[j]][1] ]]
targetNode <- maps$graphID_2_nodeName[firstDepPath[1,1]]
depNode <- maps$graphID_2_nodeName[firstDepPath[nrow(firstDepPath), 1]]
oneDepType <- conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode, restrictLink)
conjDepTypes[j] <- if(is.null(oneDepType)) "" else oneDepType
}
conjBool <- conjDepTypes != ""
names(conjDepTypes) <- names(conjBool) <- names(uniquePathsUnlistedIndices)
if(any(conjBool)) {
targetNodes <- unlist(lapply(depPathsByNode, function(x) if(is.null(x)) '_NO_DEPS_' else maps$graphID_2_nodeName[x[[1]][1,1]]))
ansList[[length(ansList)+1]] <- mapply(
function(targetNode, depPathsOneNode, depPathsLabelsOneNode) {
if(targetNode == '_NO_DEPS_') return(NULL) ## these should have already been weeded out
if(all(conjBool[depPathsLabelsOneNode])) {
depTypes <- conjDepTypes[depPathsLabelsOneNode]
depEnds <- maps$graphID_2_nodeName[ unlist(lapply(depPathsOneNode, function(x) x[nrow(x)])) ]
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
## prevent multiple instances of same dependent node name
## (via different graph dependency paths -DT Oct 2016
##depEnds[boolMatch]
unique(depEnds[boolMatch])
})
names(control) <- uniqueDepTypes
list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control)
}
},
targetNodes, depPathsByNode, depPathsByNodeLabels, USE.NAMES = TRUE, SIMPLIFY = FALSE)
}
}
}
if(length(ansList) > 0) ansList <- do.call('c', ansList)
ansList <- ansList[!sapply(ansList, is.null)] # strips out any NULL values
if(!length(ansList)) return(list()) else return(ansList) # replaces empty named list with empty list
},
checkConjugacy_original = function(model, nodeIDs, restrictLink = NULL) {
maps <- model$modelDef$maps
nodeDeclIDs <- maps$graphID_2_declID[nodeIDs] ## declaration IDs of the nodeIDs
declID2nodeIDs <- split(nodeIDs, nodeDeclIDs) ## nodeIDs grouped by declarationID
ansList <- ansList2 <- list()
for(i in seq_along(declID2nodeIDs)) { ## For each group of nodeIDs from the same declarationID
nodeIDsFromOneDecl <- declID2nodeIDs[[i]]
firstNodeName <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[1]]
if(model$isTruncated(firstNodeName)) next ## we say non-conjugate if the targetNode is truncated
dist <- model$getDistribution(firstNodeName)
conjugacyObj <- conjugacys[[dist]]
if(is.null(conjugacyObj)) next
# NO: insert logic here to check a single dependency and do next if can't be conjugate
#model$getDependencies('mu[1]',self=F,stochOnly=T)
#for( loop through deps )
# conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode)
# if(not conj) next
# now try to guess if finding paths will be more intensive than simply looking at target-dependent pairs, to avoid path finding when there is nested structure such as stickbreaking
numPaths <- sapply(nodeIDsFromOneDecl, model$getDependencyPathCountOneNode)
deps <- lapply(nodeIDsFromOneDecl, function(x) model$getDependencies(x, stochOnly = TRUE, self = FALSE))
numDeps <- sapply(deps, length)
if(max(numPaths) > sum(numDeps)) {
# max(numPaths) is reasonable guess at number of unique (by node) paths (though it overestimates number of unique (by declaration ID) paths; if we have to evaluate conjugacy for more paths than we would by simply looking at all pairs of target-dependent nodes, then just use node pairs
# note that it's not clear what criterion to use here since computational time is combination of time for finding all paths and then for evaluating conjugacy for unique (by declaration ID) paths, but the hope is to make a crude cut here that avoids path calculations when there would be a lot of them
ansList[[length(ansList)+1]] <- lapply(seq_along(nodeIDsFromOneDecl),
function(index) {
targetNode <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[index]]
depEnds <- deps[[index]]
depTypes <- sapply(depEnds, function(x) conjugacyObj$checkConjugacyOneDep(model, targetNode, x, restrictLink))
if(!length(depTypes)) return(NULL)
if(!any(sapply(depTypes, is.null))) {
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
depEnds[boolMatch]
})
names(control) <- uniqueDepTypes
return(list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control))
} else return(NULL)
})
names(ansList[[length(ansList)]]) <- maps$graphID_2_nodeName[nodeIDsFromOneDecl]
} else {
# determine conjugacy based on unique (by declaration ID) paths
depPathsByNode <- lapply(nodeIDsFromOneDecl, getDependencyPaths, maps = maps) ## make list (by nodeID) of lists of paths through graph
depPathsByNode <- depPathsByNode[!unlist(lapply(depPathsByNode, function(x) is.null(x) || (length(x)==0)))]
depPathsByNodeLabels <- lapply(depPathsByNode, function(z) ## make character labels that match for same path through graph
unlist(lapply(z,
function(x)
paste(maps$graphID_2_declID[x[,1]], x[,2], collapse = '\r', sep='\r'))))
depPathsByNodeUnlisted <- unlist(depPathsByNode, recursive = FALSE)
depPathsByNodeLabelsUnlisted <- unlist(depPathsByNodeLabels)
## uniquePaths <- unique(depPathsByNodeLabelsUnlisted)
uniquePathsUnlistedIndices <- split(seq_along(depPathsByNodeLabelsUnlisted), depPathsByNodeLabelsUnlisted)
conjDepTypes <- character(length(uniquePathsUnlistedIndices))
for(j in seq_along(uniquePathsUnlistedIndices)) {
firstDepPath <- depPathsByNodeUnlisted[[ uniquePathsUnlistedIndices[[j]][1] ]]
targetNode <- maps$graphID_2_nodeName[firstDepPath[1,1]]
depNode <- maps$graphID_2_nodeName[firstDepPath[nrow(firstDepPath), 1]]
oneDepType <- conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode, restrictLink)
conjDepTypes[j] <- if(is.null(oneDepType)) "" else oneDepType
}
conjBool <- conjDepTypes != ""
names(conjDepTypes) <- names(conjBool) <- names(uniquePathsUnlistedIndices)
if(any(conjBool)) {
targetNodes <- unlist(lapply(depPathsByNode, function(x) if(is.null(x)) '_NO_DEPS_' else maps$graphID_2_nodeName[x[[1]][1,1]]))
ansList[[length(ansList)+1]] <- mapply(
function(targetNode, depPathsOneNode, depPathsLabelsOneNode) {
if(targetNode == '_NO_DEPS_') return(NULL) ## these should have already been weeded out
if(all(conjBool[depPathsLabelsOneNode])) {
depTypes <- conjDepTypes[depPathsLabelsOneNode]
depEnds <- maps$graphID_2_nodeName[ unlist(lapply(depPathsOneNode, function(x) x[nrow(x)])) ]
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
## prevent multiple instances of same dependent node name
## (via different graph dependency paths -DT Oct 2016
##depEnds[boolMatch]
unique(depEnds[boolMatch])
})
names(control) <- uniqueDepTypes
list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control)
}
},
targetNodes, depPathsByNode, depPathsByNodeLabels, USE.NAMES = TRUE, SIMPLIFY = FALSE)
}
}
}
if(length(ansList) > 0) ansList <- do.call('c', ansList)
ansList <- ansList[!sapply(ansList, is.null)] # strips out any NULL values
if(!length(ansList)) return(list()) else return(ansList) # replaces empty named list with empty list
},
generateDynamicConjugateSamplerDefinition = function(prior, dependentCounts, doDependentScreen = FALSE) {
## conjugateSamplerDefinitions[[paste0('sampler_conjugate_', conjugacyResult$prior)]] ## using original (non-dynamic) conjugate sampler functions
conjugacys[[prior]]$generateConjugateSamplerDef(dynamic = TRUE, dependentCounts = dependentCounts,
doDependentScreen = doDependentScreen)
}
)
)
setMethod(
'[[',
'conjugacyRelationshipsClass',
function(x, i) return(x$conjugacys[[i]])
)
conjugacyClass <- setRefClass(
Class = 'conjugacyClass',
fields = list(
samplerType = 'ANY', ## name of the sampler for this conjugacy class, e.g. 'conjugate_dnorm'
prior = 'ANY', ## name of the prior distribution, e.g. 'dnorm'
link = 'ANY', ## the link ('linear', 'additive', 'multiplicative', 'multiplicativeScalar', or 'identity')
dependents = 'ANY', ## (named) list of dependentClass objects, each contains conjugacy information specific to a particular sampling distribution (name is sampling distribution name)
dependentDistNames = 'ANY', ## character vector of the names of all allowable dependent sampling distributions. same as: names(dependents)
posteriorObject = 'ANY', ## an object of posteriorClass
## needsLinearityCheck = 'ANY', ## logical specifying whether we need to do the linearity check; if the link is 'multiplicative' or 'linear'
## needsStickbreakingCheck = 'ANY', ## logical specifying whether we need to do the stickbreaking check
model = 'ANY', ## these fields ONLY EXIST TO PREVENT A WARNING for '<<-',
DEP_VALUES_VAR_INDEXED = 'ANY', ## in the code for generating the conjugate sampler functions
DEP_PARAM_VAR_INDEXED = 'ANY', ##
DEP_OFFSET_VAR = 'ANY', ##
DEP_COEFF_VAR = 'ANY', ##
CONTRIB_NAME = 'ANY' ##
),
methods = list(
initialize = function(cr) {
dependents <<- list()
samplerType <<- cc_makeSamplerTypeName(cr$prior)
prior <<- cr$prior
link <<- cr$link
initialize_addDependents(cr$dependents)
## removed because link info in specific conjugacy can now override default link
## needsLinearityCheck <<- link %in% c('multiplicative', 'linear')
## needsStickbreakingCheck <<- link %in% c('stickbreaking')
posteriorObject <<- posteriorClass(cr$posterior, prior)
},
initialize_addDependents = function(depList) {
for(i in seq_along(depList)) {
dependents[[i]] <<- dependentClass(depList[[i]], names(depList)[i])
}
names(dependents) <<- names(depList)
dependentDistNames <<- names(dependents)
},
## used by new checkConjugacy() system
## see checkConjugacy for more explanation of each step
checkConjugacyOneDep = function(model, targetNode, depNode, restrictLink = NULL) {
if(model$getDistribution(targetNode) != prior) return(NULL) # check prior distribution of targetNode
if(model$isTruncated(depNode)) return(NULL) # if depNode is truncated, then not conjugate
depNodeDist <- model$getDistribution(depNode)
if(!(depNodeDist %in% dependentDistNames)) return(NULL) # check sampling distribution of depNode
dependentObj <- dependents[[depNodeDist]]
if(is.null(restrictLink)) {
if(!is.null(dependentObj$link)) currentLink <- dependentObj$link else currentLink <- link # handle multiple link case introduced for beta stickbreaking
} else currentLink <- restrictLink
if(currentLink != 'stickbreaking') {
depNodeParamName <- dependentObj$param
linearityCheckExprRaw <- model$getParamExpr(depNode, depNodeParamName) # extracts the expression for 'param' from 'depNode'
linearityCheckExpr <- cc_expandDetermNodesInExpr(model, linearityCheckExprRaw, targetNode = targetNode)
if(!cc_nodeInExpr(targetNode, linearityCheckExpr)) return(NULL)
if(cc_vectorizedComponentCheck(targetNode, linearityCheckExpr)) return(NULL) # if targetNode is vectorized, make sure none of its components appear in expr
linearityCheck <- cc_checkLinearity(linearityCheckExpr, targetNode) # determines whether paramExpr is linear in targetNode
realizedLink <- cc_linkCheck(linearityCheck, currentLink)
if(is.null(realizedLink)) return(NULL)
## ensure targetNode appears in only *one* depNode parameter expression
if(!cc_otherParamsCheck(model, depNode, targetNode, depNodeExprExpanded = linearityCheckExpr, depParamNodeName = depNodeParamName)) return(NULL)
} else {
stickbreakingCheckExpr <- model$getParamExpr(depNode, dependentObj$param) # extracts the expression for 'param' from 'depNode'
stickbreakingCheckExpr <- cc_expandDetermNodesInExpr(model, stickbreakingCheckExpr, targetNode = targetNode)
if(is.null(cc_checkStickbreaking(stickbreakingCheckExpr, targetNode))) return(NULL)
realizedLink <- 'stickbreaking'
}
return(paste0('dep_', depNodeDist, '_', realizedLink))
},
addDependentNodeToControl = function(control, depNodeDist, depNode) {
listName <- paste0('dep_', depNodeDist)
control[[listName]] <- c(control[[listName]], depNode)
control
},
## workhorse for creating conjugate sampler nimble functions
generateConjugateSamplerDef = function(dynamic = FALSE, dependentCounts, doDependentScreen = FALSE) {
if(!dynamic) stop('something went wrong, should never have dynamic = FALSE here')
substitute(
nimbleFunction(contains = sampler_BASE,
setup = SETUPFUNCTION,
run = RUNFUNCTION,
methods = list(getPosteriorLogDensity = GETPOSTERIORLOGDENSITYFUNCTION,
reset = function() {})
),
list(SETUPFUNCTION = genSetupFunction(dependentCounts = dependentCounts, doDependentScreen = doDependentScreen),
RUNFUNCTION = genRunFunction(dependentCounts = dependentCounts, doDependentScreen = doDependentScreen),
GETPOSTERIORLOGDENSITYFUNCTION = genGetPosteriorLogDensityFunction(dependentCounts = dependentCounts, doDependentScreen = doDependentScreen)
)
)
},
genSetupFunction = function(dependentCounts, doDependentScreen = FALSE) {
functionBody <- codeBlockClass()
functionBody$addCode({
calcNodes <- model$getDependencies(target)
calcNodesDeterm <- model$getDependencies(target, determOnly = TRUE)
})
## if this conjugate sampler is for a multivariate node (i.e., nDim > 0), then we need to determine the size (d)
if(getDimension(prior) > 0) {
functionBody$addCode(d <- max(determineNodeIndexSizes(target)))
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], "_")[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
functionBody$addCode({
DEP_NODENAMES <- control$DEP_CONTROL_NAME
N_DEP <- length(control$DEP_CONTROL_NAME)
}, list(DEP_NODENAMES = as.name(paste0( 'dep_', distLinkName, '_nodeNames')),
N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_CONTROL_NAME = as.name(paste0( 'dep_', distLinkName))))
if(any(getDimension(distName, includeParams = TRUE) > 0)) {
if(getDimension(distName) > 0) { ## usual case: dependent node is multivariate -> get sizes from dependent node names
functionBody$addCode({
DEP_NODESIZES <- sapply(DEP_NODENAMES, function(node) max(determineNodeIndexSizes(node)), USE.NAMES = FALSE)
}, list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames'))))
} else { ## unusual case: dependent is univariate -> get sizes from dependent node's param (e.g., ddirch-dcat conjugacy)
functionBody$addCode({
DEP_NODESIZES <- sapply(DEP_NODENAMES, function(node) max(nimDim(model$getParam(node, MULTIVARIATE_PARAM_NAME))), USE.NAMES = FALSE)
}, list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
MULTIVARIATE_PARAM_NAME = names(which.max(getDimension(distName, includeParams = TRUE)))))
}
functionBody$addCode({
if(length(DEP_NODESIZES) == 1) DEP_NODESIZES <- c(DEP_NODESIZES, -1) ## guarantee to be a vector, for indexing and size processing
DEP_NODESIZEMAX <- max(DEP_NODESIZES)
}, list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_NODESIZEMAX = as.name(paste0('dep_', distLinkName, '_nodeSizeMax'))))
}
## declare() statements are removed from run() code,
## and were replaced with setup output array() calls below.
## July 2017
depNodeValueNdim <- getDimension(distName)
functionBody$addCode(DEP_VALUES_VAR <- array(0, dim = DECLARE_SIZE),
list(DEP_VALUES_VAR = as.name(paste0('dep_', distLinkName, '_values')),
DECLARE_SIZE = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), depNodeValueNdim)))
neededParams <- dependents[[distName]]$neededParamsForPosterior
for(param in neededParams) {
depNodeParamNdim <- getDimension(distName, param)
functionBody$addCode(DEP_PARAM_VAR <- array(0, dim = DECLARE_SIZE),
list(DEP_PARAM_VAR = as.name(paste0('dep_', distLinkName, '_', param)),
DECLARE_SIZE = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), depNodeParamNdim)))
}
}
## more new array() setup outputs, instead of declare() statements, for offset and coeff variables
## July 2017
targetNdim <- getDimension(prior)
targetCoeffNdim <- switch(as.character(targetNdim), `0`=0, `1`=2, `2`=2, stop())
if(targetCoeffNdim == 2 && link == 'multiplicativeScalar') ## Handles wish/invwish. There are no cases where we allow non-scalar 'coeff'.
targetCoeffNdim <- 0
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], "_")[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('additive', 'multiplicative', 'multiplicativeScalar', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
## the 2's here are *only* to prevent warnings about assigning into member variable names using
inputList <- list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')), ## local assignment '<-', so changed the names to "...2"
DEP_COEFF_VAR2 = as.name(paste0('dep_', distLinkName, '_coeff')), ## so it doesn't recognize the ref class field name
DECLARE_SIZE_OFFSET = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), targetNdim),
DECLARE_SIZE_COEFF = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), quote(d), targetCoeffNdim))
if(currentLink == 'additive')
functionBody$addCode(
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET),
inputList)
if(currentLink %in% c('multiplicative', 'multiplicativeScalar'))
functionBody$addCode(
DEP_COEFF_VAR2 <- array(0, dim = DECLARE_SIZE_COEFF),
inputList)
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET)
DEP_COEFF_VAR2 <- array(0, dim = DECLARE_SIZE_COEFF)
}, inputList)
}
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], "_")[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(!is.null(dependents[[distName]]$link)) currentLink <- dependents[[distName]]$link
if(currentLink == 'stickbreaking') {
functionBody$addCode({
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET) ## the 2's here are *only* to prevent warnings about
}, list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')), ## local assignment '<-', so changed the names to "...2"
DECLARE_SIZE_OFFSET = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), 0))
)
functionBody$addCode({
stickbreakingCheckExpr <- model$getValueExpr(calcNodesDeterm)
stickbreakingCheckExpr <- cc_expandDetermNodesInExpr(model, stickbreakingCheckExpr, targetNode = target)
})
functionBody$addCode(DEP_OFFSET_VAR2 <- rep(cc_checkStickbreaking(stickbreakingCheckExpr, target)$offset, DEP_OFFSET_SIZE),
list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')),
DEP_OFFSET_SIZE = as.name(paste0('N_dep_', distLinkName))))
}
}
## adding declarations for the contribution terms, to remove Windows compiler warnings, DT August 2015
## moved these numeric() and array() declarations for contribution terms to setup outputs, July 2017
for(contributionName in posteriorObject$neededContributionNames) {
contribNdim <- posteriorObject$neededContributionDims[[contributionName]]
functionBody$addCode(CONTRIB_NAME2 <- CONTRIB_INITIAL_DECLARATION, ## the 2's here are *only* to prevent warnings about
list(CONTRIB_NAME2 = as.name(contributionName), ## local assignment '<-' versus '<<-'
CONTRIB_INITIAL_DECLARATION = switch(as.character(contribNdim),
`0` = 0, `1` = quote(rep(0, length = d)), `2` = quote(array(0, dim = c(d, d))), stop())))
}
functionDef <- quote(function(model, mvSaved, target, control) {})
functionDef[[3]] <- functionBody$getCode()
functionDef[[4]] <- NULL ## removes the 'scrref' attribute
return(functionDef)
},
genRunFunction = function(dependentCounts, doDependentScreen = FALSE) {
functionBody <- codeBlockClass()
## only if we're verifying conjugate posterior distributions: get initial targetValue, and modelLogProb -- model$getLogProb(calcNodes)
if(getNimbleOption('verifyConjugatePosteriors')) {
functionBody$addCode({
modelLogProb0 <- model$getLogProb(calcNodes)
origTargetValue <- model[[target]]
})
}
## adds code to generate the quantities prior_xxx, and contribution_xxx
addPosteriorQuantitiesGenerationCode(functionBody = functionBody, dependentCounts = dependentCounts, doDependentScreen = doDependentScreen)
## generate new value, store, calculate, copy, etc...
functionBody$addCode(posteriorObject$prePosteriorCodeBlock, quote = FALSE)
functionBody$addCode({
newTargetValue <- RPOSTERIORCALL
model[[target]] <<- newTargetValue
model$calculate(calcNodes)
nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
}, list(RPOSTERIORCALL = posteriorObject$rCallExpr))
## only if we're verifying conjugate posterior distributions: figure out if conjugate posterior distribution is correct
if(getNimbleOption('verifyConjugatePosteriors')) {
functionBody$addCode({
modelLogProb1 <- model$getLogProb(calcNodes)
posteriorLogDensity0 <- DPOSTERIORCALL_ORIG
posteriorLogDensity1 <- DPOSTERIORCALL_NEW
posteriorVerification <- modelLogProb0 - posteriorLogDensity0 - modelLogProb1 + posteriorLogDensity1
if(abs(posteriorVerification) > 1e-8) {
nimPrint('conjugate posterior density appears to be wrong, off by ', posteriorVerification)
}
}, list(DPOSTERIORCALL_ORIG = eval(substitute(substitute(expr, list(VALUE=quote(origTargetValue))), list(expr=posteriorObject$dCallExpr))),
DPOSTERIORCALL_NEW = eval(substitute(substitute(expr, list(VALUE=quote(newTargetValue))), list(expr=posteriorObject$dCallExpr)))))
}
functionDef <- quote(function() {})
functionDef[[3]] <- functionBody$getCode()
functionDef[[4]] <- NULL ## removes the 'scrref' attribute
return(functionDef)
},
genGetPosteriorLogDensityFunction = function(dependentCounts, doDependentScreen = FALSE) {
functionBody <- codeBlockClass()
functionBody$addCode(origTargetValue <- model[[target]])
## adds code to generate the quantities prior_xxx, and contribution_xxx
addPosteriorQuantitiesGenerationCode(functionBody = functionBody, dependentCounts = dependentCounts, doDependentScreen = doDependentScreen)
## calculate and return the (log)density for the current value of target
functionBody$addCode(posteriorObject$prePosteriorCodeBlock, quote = FALSE)
if(link != 'identity') {
functionBody$addCode({
model[[target]] <<- origTargetValue
model$calculate(calcNodesDeterm)})
}
functionBody$addCode({
posteriorLogDensity <- DPOSTERIORCALL
returnType(double())
return(posteriorLogDensity)
}, list(DPOSTERIORCALL = eval(substitute(substitute(expr, list(VALUE=quote(origTargetValue))), list(expr=posteriorObject$dCallExpr)))))
functionDef <- quote(function() {})
functionDef[[3]] <- functionBody$getCode()
functionDef[[4]] <- NULL ## removes the 'scrref' attribute
return(functionDef)
},
addPosteriorQuantitiesGenerationCode = function(functionBody = functionBody, dependentCounts = dependentCounts, doDependentScreen = FALSE) {
## get current value of prior parameters which appear in the posterior expression
for(priorParam in posteriorObject$neededPriorParams) {
functionBody$addCode(PRIOR_PARAM_VAR <- model$getParam(target[1], PARAM_NAME),
list(PRIOR_PARAM_VAR = as.name(paste0('prior_', priorParam)),
PARAM_NAME = priorParam))
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
neededParams <- dependents[[distName]]$neededParamsForPosterior
depNodeValueNdim <- getDimension(distName)
forLoopBody <- codeBlockClass()
## get *value* of each dependent node
if(any(getDimension(distName, includeParams = TRUE) > 0)) {
forLoopBody$addCode(thisNodeSize <- DEP_NODESIZES[iDep],
list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes'))))
}
forLoopBody$addCode(DEP_VALUES_VAR_INDEXED <<- model$getParam(DEP_NODENAMES[iDep], 'value'),
list(DEP_VALUES_VAR_INDEXED = makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_values')), depNodeValueNdim, indexExpr = quote(iDep), secondSize = quote(thisNodeSize), thirdSize = quote(thisNodeSize)),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames'))))
for(param in neededParams) {
depNodeParamNdim <- getDimension(distName, param)
## get *parameter values* for each dependent node
forLoopBody$addCode(DEP_PARAM_VAR_INDEXED <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME),
list(DEP_PARAM_VAR_INDEXED = makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_', param)), depNodeParamNdim, indexExpr = quote(iDep), secondSize = quote(thisNodeSize), thirdSize = quote(thisNodeSize)),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames')),
PARAM_NAME = param))
}
functionBody$addCode(for(iDep in 1:N_DEP) FORLOOPBODY,
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
FORLOOPBODY = forLoopBody$getCode()))
}
targetNdim <- getDimension(prior)
allCurrentLinks <- sapply(names(dependentCounts), function(x) strsplit(x, '_')[[1]][[2]])
switch(as.character(targetNdim),
`0` = {
if(any(allCurrentLinks %in% c('additive', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- 0
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('additive', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode(
for(iDep in 1:N_DEP)
DEP_OFFSET_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME),
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
}
if(any(allCurrentLinks == 'multiplicativeScalar'))
stop("Found 'multiplicativeScalar' link for 0-dimensional case.")
if(any(allCurrentLinks %in% c('multiplicative', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- 1
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('multiplicative', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
inputList <- list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param,
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')))
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode(
for(iDep in 1:N_DEP)
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME) - DEP_OFFSET_VAR[iDep],
inputList)
} else {
functionBody$addCode(
for(iDep in 1:N_DEP)
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME),
inputList)
}
}
}
}
},
`1` = {
if(any(allCurrentLinks == 'multiplicativeScalar'))
stop("Found 'multiplicativeScalar' link for 1-dimensional case.")
if(any(allCurrentLinks %in% c('additive', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- rep(0, d)
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('additive', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
for(iDep in 1:N_DEP) {
thisNodeSize <- DEP_NODESIZES[iDep]
DEP_OFFSET_VAR[iDep, 1:thisNodeSize] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)
}
},
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
}
if(any(allCurrentLinks %in% c('multiplicative', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode(unitVector <- rep(0, d))
forLoopBody <- codeBlockClass()
forLoopBody$addCode({
unitVector[sizeIndex] <- 1
model[[target]] <<- unitVector
unitVector[sizeIndex] <- 0
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('multiplicative', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
inputList <- list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param,
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')))
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
forLoopBody$addCode(
for(iDep in 1:N_DEP) {
thisNodeSize <- DEP_NODESIZES[iDep]
DEP_COEFF_VAR[iDep, 1:thisNodeSize, sizeIndex] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME) - DEP_OFFSET_VAR[iDep, 1:thisNodeSize]
}, inputList)
} else {
forLoopBody$addCode(
for(iDep in 1:N_DEP) {
thisNodeSize <- DEP_NODESIZES[iDep]
DEP_COEFF_VAR[iDep, 1:thisNodeSize, sizeIndex] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)
}, inputList)
}
}
}
functionBody$addCode(for(sizeIndex in 1:d) FORLOOPBODY,
list(FORLOOPBODY = forLoopBody$getCode()))
}
},
`2` = {
if(!all(allCurrentLinks %in% c('identity', 'multiplicativeScalar')))
stop("Found non-multiplicative link for 2-d variable.")
if(any(allCurrentLinks == 'multiplicativeScalar') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- diag(d)
model$calculate(calcNodesDeterm)
})
## Use _COEFF_VAR to store value; we need this for stoch indexing case
## where we determine that coeff = 0 because the potential dependency is not a dependency
## given current index values.
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink == 'multiplicativeScalar' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
for(iDep in 1:N_DEP) {
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)[1, 1] ## DEP_COEFF_VAR = (A+2B)-(A+B) = B
}
},
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
functionBody$addCode({
## Can't use zeros matrix as Cholesky fails; this is solely to determine if
## potential dependency is not a dependency of target (due to stochastic indexing).
model[[target]] <<- 2*diag(d)
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink == 'multiplicativeScalar'|| (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
for(iDep in 1:N_DEP) {
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)[1, 1] - DEP_COEFF_VAR[iDep]
}
},
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
}
},
stop()
)
targetNdim <- getDimension(prior)
## contribution terms have been moved to be setup function outputs,
## but we still need to *zero these variables out* before adding contribution terms into them
for(contributionName in posteriorObject$neededContributionNames) {
contribNdim <- posteriorObject$neededContributionDims[[contributionName]]
functionBody$addCode(CONTRIB_NAME <<- CONTRIB_ZERO_OUT,
list(CONTRIB_NAME = as.name(contributionName),
CONTRIB_ZERO_OUT = switch(as.character(contribNdim),
`0` = 0, `1` = quote(rep(0, length = d)), `2` = quote(array(0, dim = c(d, d))), stop())))
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
targetCoeffNdim <- switch(as.character(targetNdim), `0`=0, `1`=2, `2`=2, stop())
if(targetCoeffNdim == 2 && link == 'multiplicativeScalar') ## There are no cases where we allow non-scalar 'coeff'.
targetCoeffNdim <- 0
if(!any(posteriorObject$neededContributionNames %in% dependents[[distName]]$contributionNames)) next
depParamsAvailable <- dependents[[distName]]$neededParamsForPosterior
## don't allow ragged dependencies for 2D conjugate case.
## no such cases exist, and it causes a runtime size check compiler warning.
## nonRaggedSizeExpr used to replace quote(thisNodeSize) below.
## August 2016
nonRaggedSizeExpr <- if(targetNdim < 2) quote(thisNodeSize) else quote(d)
subList <- lapply(depParamsAvailable, function(param)
makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_', param)), getDimension(distName, param), indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = nonRaggedSizeExpr))
names(subList) <- depParamsAvailable
subList$value <- makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_values')), getDimension(distName), indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = nonRaggedSizeExpr)
if(currentLink %in% c('additive', 'linear', 'stickbreaking'))
subList$offset <- makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_offset')), targetNdim, indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = nonRaggedSizeExpr)
if(currentLink %in% c('multiplicative', 'multiplicativeScalar', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
subList$coeff <- makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_coeff')), targetCoeffNdim, indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = quote(d))
forLoopBody <- codeBlockClass()
if(any(getDimension(distName, includeParams = TRUE) > 0)) {
if(targetNdim == 1) ## 1D
forLoopBody$addCode(thisNodeSize <- DEP_NODESIZES[iDep],
list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes'))))
if(targetNdim == 2) ## 2D ## formerly this was 'else', but for 'dcat' we have targetNdim=0 while max(getDimension(distName, includeParams = TRUE)) is 1 so need explicit check for 2D
forLoopBody$addCode(if(DEP_NODESIZES[iDep] != d) print('runtime error with sizes of 2D conjugate sampler'),
list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes'))))
}
for(contributionName in posteriorObject$neededContributionNames) {
if(!(contributionName %in% dependents[[distName]]$contributionNames)) next
contributionExpr <- dependents[[distName]]$contributionExprs[[contributionName]]
if(distName == 'dmnorm' && prior == 'dmnorm') {
## need to deal with [,1] in contribution_mean
if(contributionName == 'contribution_mean') tmpExpr <- contributionExpr[[2]][[2]] else tmpExpr <- contributionExpr
if(getNimbleOption('allowDynamicIndexing') && doDependentScreen) {
## Will always have 'coeff' but may not need to use it
if(currentLink %in% c('identity', 'additive'))
tmpExpr[[6]] <- quote(0) else tmpExpr[[6]] <- quote(1)
} else {
if(currentLink %in% c('identity', 'additive')) {
tmpExpr[[2]] <- quote(prec) ## hack, since 'coeff' won't exist
tmpExpr[[6]] <- quote(0)
} else tmpExpr[[6]] <- quote(1)
}
tmpExpr[[4]] <- cc_stripExpr(tmpExpr[[4]], offset = currentLink %in% c('identity','multiplicative'), coeff = FALSE) # strip 'offset'
if(contributionName == 'contribution_mean') contributionExpr[[2]][[2]] <- tmpExpr else contributionExpr <- tmpExpr
} else contributionExpr <- cc_stripExpr(contributionExpr, offset = currentLink %in% c('identity','multiplicative','multiplicativeScalar'),
coeff = currentLink %in% c('identity','additive'))
contributionExpr <- eval(substitute(substitute(EXPR, subList), list(EXPR=contributionExpr)))
if(getNimbleOption('allowDynamicIndexing') && doDependentScreen) { ## FIXME: would be nice to only have one if() here when we loop through multiple parameters
if(targetCoeffNdim == 0)
forLoopBody$addCode(if(COEFF_EXPR != 0) CONTRIB_NAME <<- CONTRIB_NAME + CONTRIB_EXPR,
list(COEFF_EXPR = subList$coeff, CONTRIB_NAME = as.name(contributionName), CONTRIB_EXPR = contributionExpr))
else forLoopBody$addCode(if(min(COEFF_EXPR) != 0 | max(COEFF_EXPR) != 0) CONTRIB_NAME <<- CONTRIB_NAME + CONTRIB_EXPR,
list(COEFF_EXPR = subList$coeff, CONTRIB_NAME = as.name(contributionName), CONTRIB_EXPR = contributionExpr))
} else forLoopBody$addCode(CONTRIB_NAME <<- CONTRIB_NAME + CONTRIB_EXPR,
list(CONTRIB_NAME = as.name(contributionName), CONTRIB_EXPR = contributionExpr))
}
functionBody$addCode(for(iDep in 1:N_DEP) FORLOOPBODY,
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
FORLOOPBODY = forLoopBody$getCode()))
}
}
)
)
dependentClass <- setRefClass(
Class = 'dependentClass',
fields = list(
distribution = 'ANY', ## the name of the (dependent) sampling distribution, e.g. 'dnorm'
param = 'ANY', ## the name of the sampling distribution parameter in which target must appear
contributionExprs = 'ANY', ## a (named) list of expressions, giving the (additive) contribution to any parameters of the posterior. names correspond to variables in the posterior expressions
contributionNames = 'ANY', ## names of the contributions to the parameters of the posterior distribution. same as names(posteriorExprs)
neededParamsForPosterior = 'ANY', ## names of all parameters appearing in the posteriorExprs
link = 'ANY' ## optional specific link for a given conjugacy; currently only used for stickbreaking case
),
methods = list(
initialize = function(depInfoList, depDistName) {
contributionExprs <<- list()
distribution <<- depDistName
param <<- depInfoList$param
link <<- depInfoList$link ## will be NULL unless specific conjugacy overrides default link
initialize_contributionExprs(depInfoList)
initialize_neededParamsForPosterior()
},
initialize_contributionExprs = function(depInfoList) {
depInfoList['param'] <- NULL
depInfoList <- lapply(depInfoList, function(di) parse(text=di)[[1]])
contributionExprs <<- depInfoList
contributionExprs[['link']] <<- NULL ## so link not included in neededParamsForPosterior
contributionNames <<- names(depInfoList)
},
initialize_neededParamsForPosterior = function() {
posteriorVars <- unique(unlist(lapply(contributionExprs, all.vars)))
## Formerly, we always included the target 'param' in 'neededParamsForPosterior'.
## Not certain why, it wasn't necessary in many cases, e.g., dbeta-dbinom conjugacy. Excluding it now.
## -DT, August 2017
##neededParamsForPosterior <<- unique(c(param, posteriorVars[!(posteriorVars %in% c('value', 'coeff', 'offset'))]))
neededParamsForPosterior <<- setdiff(posteriorVars, c('value', 'coeff', 'offset'))
}
)
)
posteriorClass <- setRefClass(
Class = 'posteriorClass',
fields = list(
prePosteriorCodeBlock = 'ANY', ## a quoted {...} code block containing DSL code to execute before making posterior call, possibly empty
posteriorExpr = 'ANY', ## the full, parsed, posterior distribution expression, e.g. dnorm(mean = prior_mean + ..., sd = ...)
rDistribution = 'ANY', ## the *R* name of the posterior distribution, e.g. 'rnorm'
dDistribution = 'ANY', ## the *R* name of the posterior density distribution, e.g. 'dnorm'
argumentExprs = 'ANY', ## (named) list of expressions for each argument to the posterior distribution. names are the posterior distribution argument names
argumentNames = 'ANY', ## character vector of the argument names to the posterior distribution. same as: names(argumentExprs)
rCallExpr = 'ANY', ## the actual 'rnorm(1, ...)' call, which will be substituted into the conjugate sampler function
dCallExpr = 'ANY', ## the 'dnorm(value, ...)' call, which can be used to get values of the posterior density
neededPriorParams = 'ANY', ## the names of any prior parameters (e.g., 'mean') which appear in the posterior expression as 'prior_mean'
neededContributionNames = 'ANY', ## the names of contributions from dependent nodes, such as 'contribution_scale'
neededContributionDims = 'ANY' ## a named list of contribution dimensions (0, 1, or 2). List element names are, e.g., 'contribution_scale'
),
methods = list(
initialize = function(posteriorText, prior) {
parsedTotalPosterior <- parse(text = posteriorText)[[1]]
if(parsedTotalPosterior[[1]] != '{') parsedTotalPosterior <- substitute({POST}, list(POST = parsedTotalPosterior))
prePosteriorCodeBlock <<- parsedTotalPosterior[-length(parsedTotalPosterior)]
posteriorExpr <<- parsedTotalPosterior[[length(parsedTotalPosterior)]]
rDistribution <<- cc_makeRDistributionName(as.character(posteriorExpr[[1]]))
dDistribution <<- as.character(posteriorExpr[[1]])
argumentExprs <<- as.list(posteriorExpr)[-1]
argumentNames <<- names(argumentExprs)
rCallExpr <<- as.call(c(as.name(rDistribution), 1, argumentExprs))
dCallExpr <<- as.call(c(as.name(dDistribution), quote(VALUE), argumentExprs, log = 1))
posteriorVars <- all.vars(parsedTotalPosterior)
neededPriorParams <<- gsub('^prior_', '', posteriorVars[grepl('^prior_', posteriorVars)])
neededContributionNames <<- posteriorVars[grepl('^contribution_', posteriorVars)]
neededContributionDims <<- inferContributionTermDimensions(prior)
},
inferContributionTermDimensions = function(prior) {
distToLookup <- if(dDistribution %in% distributions$namesVector) dDistribution else if(prior %in% distributions$namesVector) prior else stop('cannot locate prior or posterior distribution in conjugacy processing')
targetNdim <- getDimension(distToLookup)
## if posterior distribution is univariate, assume all contributions are scalar
if(targetNdim == 0) {
theDims <- lapply(neededContributionNames, function(x) 0)
names(theDims) <- neededContributionNames
return(theDims)
}
## if posterior distribution is multivariate, attempt to infer contribution dimensionality from the *name* of each contribution term
theDims <- list()
typeNamesAvailable <- getParamNames(distToLookup)
for(contribName in neededContributionNames) {
contribNameBase <- gsub('^contribution_', '', contribName)
if(contribNameBase %in% typeNamesAvailable) {
## contribution base name matches a parameter name of the posterior
theDims[[contribName]] <- getDimension(distToLookup, contribNameBase)
} else {
## contribution base name doesn't match any parameter; can't easily infer the dimensionality
stop(message('The NIMBLE conjugacy system is attempting to infer the dimensionality of the contribution term: ',
contribName, '. However, since the posterior distribution is multivariate, and the contribution name doesn\'t match any parameter names of the posterior distribution, NIMBLE can\'t infer this one. This means the conjugacy system might need to be extended, to allow users to provide the dimensionality of contribution terms. Or perhaps something more clever. -DT August 2015'), call. = FALSE)
}
}
return(theDims)
}
)
)
##############################################################################################
##############################################################################################
## utility functions
##############################################################################################
##############################################################################################
cc_makeSamplerTypeName <- function(distName) return(paste0('conjugate_', distName)) ## 'dnorm' --> 'conjugate_dnorm'
cc_makeConjugateSamplerName <- function(samplerType) return(paste0('sampler_', samplerType)) ## 'conjugate_dnorm' --> 'sampler_conjugate_dnorm'
cc_makeRDistributionName <- function(distName) return(paste0('r', substring(distName, 2))) ## 'dnorm' --> 'rnorm'
## expands all deterministic nodes in expr, to create a single expression with only stochastic nodes
cc_expandDetermNodesInExpr <- function(model, expr, targetNode = NULL, skipExpansionsNode = NULL, prevExpr = NULL) {
## prevents an infinite-recursion case, where structureExpr(EXPR) is repeatedly processed, indefinitely (DT Aug 2022):
if(!is.null(prevExpr) && identical(expr, prevExpr)) return(expr)
if(is.numeric(expr)) return(expr) # return numeric
if(is.name(expr) || (is.call(expr) && (expr[[1]] == '[') && is.name(expr[[2]]))) { # expr is a name, or an indexed name
if(getNimbleOption('allowDynamicIndexing')) {
## this deals with having mu[k[1]] (which won't pass through expandNodeNames), replacing k[1] with the index from targetNode
if(!is.name(expr)) {
indexExprs <- expr[3:length(expr)]
numericOrVectorIndices <- sapply(indexExprs,
function(x) is.numeric(x) || (length(x) == 3 && x[[1]] == ':'))
if(!all(numericOrVectorIndices)) {
if(model$getVarNames(nodes = targetNode) == model$getVarNames(nodes = safeDeparse(expr, warn = TRUE))) {
## expr var is same as target var, so plug in target indexes for
## non-constant expr indexes to allow for possibility that
## dynamic index will be that of target (in some cases based on allowed
## values of dynamic index, this might actually not be possible, but
## the only disdvantage should be failing to detect conjugacy where it is
## present).
expr[which(!numericOrVectorIndices)+2] <- parse(text = targetNode)[[1]][3:length(expr)][!numericOrVectorIndices]
## now check that target is not actually used in index expressions
newExpr <- as.call(c(cc_structureExprName, sapply(indexExprs[!numericOrVectorIndices], function(x) x)))
for(i in seq_along(newExpr)[-1])
newExpr[[i]] <- cc_expandDetermNodesInExpr(model, newExpr[[i]], targetNode, skipExpansionsNode)
if(cc_nodeInExpr(targetNode, newExpr))
return(as.call(c(cc_structureExprName, expr, sapply(indexExprs[!numericOrVectorIndices], function(x) x)))) # put 'expr' back in though shouldn't be needed downstream
## otherwise continue with processing as in non-dynamic index case
} else {
## in this case the expr var is a different var, so shouldn't really matter what
## indexes get plugged in; plug in mins of indexes as a default
varInfo <- model$modelDef$varInfo[[safeDeparse(expr[[2]], warn = TRUE)]]
expr[which(!numericOrVectorIndices)+2] <- varInfo$mins[!numericOrVectorIndices]
## sapply business gets rid of () at end of index expression
newExpr <- as.call(c(cc_structureExprName, expr, sapply(indexExprs[!numericOrVectorIndices], function(x) x)))
for(i in seq_along(newExpr)[-1])
newExpr[[i]] <- cc_expandDetermNodesInExpr(model, newExpr[[i]], targetNode, skipExpansionsNode, prevExpr = expr)
return(newExpr)
}
} ## else continue with processing as in non-dynamic index case
}
}
exprText <- safeDeparse(expr, warn = TRUE)
expandedNodeNamesRaw <- model$expandNodeNames(exprText)
## skipExpansionsNode was added specifically for CAR model target nodes:
if(!is.null(skipExpansionsNode) && (exprText %in% model$expandNodeNames(skipExpansionsNode, returnScalarComponents=TRUE))) return(expr)
## if exprText is a node itself (and also part of a larger node), then we only want the expansion to be the exprText node:
if(exprText %in% expandedNodeNamesRaw) {
expandedNodeNames <- exprText
} else {
## if exprText is itself *not* a formal node name, but rather a subset of a larger formal node,
## then we don't want it to be expanded to include the subsubming full node
if(length(setdiff(model$expandNodeNames(expandedNodeNamesRaw, returnScalarComponents = TRUE), model$expandNodeNames(exprText, returnScalarComponents = TRUE)))) {
expandedNodeNames <- exprText
} else {
expandedNodeNames <- expandedNodeNamesRaw
}
}
if(length(expandedNodeNames) == 1 && (expandedNodeNames == exprText) && (exprText %in% expandedNodeNamesRaw)) {
## expr is a single node in the model
type <- model$getNodeType(exprText)
## this next block covers a **real corner case** (from a legacy BUGS model: biops) where dimensions are inferred
## (since not provided), thus causing a valid stochastic node to rather have dimensions larger than it acually is
## (being defined as part of a ragged declaration of multivariate stochastic nodes),
## thus the "type" of this expanded node includes both 'stoch' and also 'RHSonly', leading to a "node"
## with two types... which triggers the check below. So, if a "node" we're processing has multiple types, then
## we remove the extraneous 'RHSonly' types here:
if((length(type) > 1) && ('RHSonly' %in% type) && !all(type == 'RHSonly')) type <- setdiff(type, 'RHSonly')
if(length(type) > 1) {
## if exprText is a node itself (and also part of a larger node), then we only want the expansion to be the exprText node:
type <- type[which(exprText == expandedNodeNamesRaw)]
}
if(type == 'stoch') return(expr)
if(type == 'determ') {
if(!(exprText %in% model$getNodeNames(determOnly = TRUE))) return(expr) ## exprText is a single element of a multivariate deterministic node
newExpr <- model$getValueExpr(exprText)
return(cc_expandDetermNodesInExpr(model, newExpr, targetNode, skipExpansionsNode))
}
if(type == 'RHSonly') return(expr)
stop(paste0('Unexpected model structure for expression: ', exprText, '. Please report this to the NIMBLE development team.'))
}
newExpr <- cc_createStructureExpr(model, exprText)
if(is.call(newExpr) && newExpr[[1]] == 'structureExpr') { ## recurse, if there's a newly created structureExpr()
for(i in seq_along(newExpr)[-1])
newExpr[[i]] <- cc_expandDetermNodesInExpr(model, newExpr[[i]], targetNode, skipExpansionsNode, prevExpr = expr)
}
return(newExpr)
}
if(is.call(expr)) {
if(any(sapply(expr[-1], function(x) x=='')))
stop('Found missing indexing in ', safeDeparse(expr), ' that prevents conjugacy processing for this particular model structure.')
for(i in seq_along(expr)[-1])
expr[[i]] <- cc_expandDetermNodesInExpr(model, expr[[i]], targetNode, skipExpansionsNode)
return(expr)
}
stop(paste0('something went wrong processing: ', safeDeparse(expr)))
}
## special name used to represent vectors / arrays defined in terms of other stoch/determ nodes
cc_structureExprName <- quote(structureExpr)
## creates an expression of the form [cc_structureExprName](element11, element12, etc...) to represent vectors / arrays defined in terms of other stoch/determ nodes,
cc_createStructureExpr <- function(model, exprText) {
expandedNodeNamesVector <- model$expandNodeNames(exprText)
## remove expanded node names which are not fully represented in the original scalar components of exprText:
exprTextScalarComponents <- model$expandNodeNames(exprText, returnScalarComponents=TRUE)
expandedNodeNamesToKeepBool <- sapply(expandedNodeNamesVector, function(n) all(model$expandNodeNames(n, returnScalarComponents=TRUE) %in% exprTextScalarComponents))
expandedNodeNamesVector <- expandedNodeNamesVector[expandedNodeNamesToKeepBool]
## now, add in any scalar components of original exprText which are not appearing in the expanded node names:
scalarComponentsToAdd <- setdiff(exprTextScalarComponents, model$expandNodeNames(expandedNodeNamesVector, returnScalarComponents=TRUE))
nodesForStructureExpr <- c(expandedNodeNamesVector, scalarComponentsToAdd)
## if exprText is a non-node scalar component, then just return that parse(exprText); (rather than structureExpr(exprText))
if((length(nodesForStructureExpr) == 1) && identical(nodesForStructureExpr, scalarComponentsToAdd)) return(parse(text=exprText)[[1]])
## otherwise, create and return a structureExpr
structureExprExprList <- lapply(nodesForStructureExpr, function(x) parse(text=x)[[1]])
structureExpr <- c(cc_structureExprName, structureExprExprList)
structureExprCall <- as.call(structureExpr)
return(structureExprCall)
}
## verifies that 'link' is satisfied by the results of linearityCheck
cc_linkCheck <- function(linearityCheck, link) {
linearityLinks <- c('identity', 'additive', 'multiplicative', 'multiplicativeScalar', 'linear')
if(!(link %in% c(linearityLinks, 'stickbreaking')))
stop(paste0('unknown link: \'', link, '\''))
if(is.null(linearityCheck)) return(NULL)
offset <- linearityCheck$offset
scale <- linearityCheck$scale
if(link %in% linearityLinks && offset == 0 && scale == 1)
return('identity')
if(link %in% c('additive', 'linear') && scale == 1)
return('additive')
## For Wishart/inverse-Wishart, only handle scalar 'scale'.
if(link == 'multiplicativeScalar' && offset == 0 && cc_checkScalar(scale))
return('multiplicativeScalar')
if(link %in% c('multiplicative', 'linear') && offset == 0)
return('multiplicative')
if(link == 'linear')
return('linear')
return(NULL)
}
cc_checkScalar <- function(expr) {
if(length(expr) == 1) return(TRUE)
if(expr[[1]] == ":") return(FALSE)
## We don't know the output dims of all functions,
## as there can be user-defined functions that produce non-scalar output from scalar inputs,
## so only say it could be scalar (based on input args) in specific known cases we enumerate.
if(safeDeparse(expr[[1]], warn = TRUE) %in% c('(','[','*','+','-','/','exp','log','^','pow','pow_int','sqrt')) {
return(all(sapply(expr[2:length(expr)], cc_checkScalar)))
} else return(FALSE)
}
## Removes 'coeff' and/or 'offset' from expression for use in simplifying conjugacy contribution
## expressions when realized link is 'simpler' than possible links.
cc_stripExpr <- function(expr, offset = TRUE, coeff = TRUE) {
if(length(expr) == 1) {
if(coeff && expr == 'coeff') return(1)
if(offset && expr == 'offset') return(0)
return(expr)
}
if(offset && (expr[[1]] == '+' || expr[[1]] == '-')) {
if(expr[[2]] == 'offset' || expr[[2]] == 0)
return(cc_stripExpr(expr[[3]], offset, coeff))
if(expr[[3]] == 'offset' || expr[[3]] == 0)
return(cc_stripExpr(expr[[2]], offset, coeff))
}
if(coeff && expr[[1]] == '*') {
if(expr[[2]] == 'coeff' || expr[[2]] == 1)
return(cc_stripExpr(expr[[3]], offset, coeff))
if(expr[[3]] == 'coeff' || expr[[3]] == 1)
return(cc_stripExpr(expr[[2]], offset, coeff))
}
if(coeff && expr[[1]] == '/' && (expr[[3]] == 'coeff' || expr[[3]] == 1))
return(cc_stripExpr(expr[[2]], offset, coeff))
if(coeff && expr[[1]] == '/' && expr[[2]] == 'coeff') {
expr[[2]] <- 1
expr[[3]] <- cc_stripExpr(expr[[3]], offset, coeff)
return(expr)
}
if(coeff && (expr[[1]] == '^' || expr[[1]] == 'sqrt') && expr[[2]] == 'coeff') {
return(1)
}
if(is.call(expr))
for(i in seq_along(expr)[-1])
expr[[i]] <- cc_stripExpr(expr[[i]], offset, coeff)
return(expr)
}
## checks the parameter expressions in the stochastic distribution of depNode
## returns FALSE if we find 'targetNode' in ***more than one*** of these expressions
cc_otherParamsCheck <- function(model, depNode, targetNode, skipExpansionsNode = NULL, depNodeExprExpanded, depParamNodeName) {
paramsList <- as.list(model$getValueExpr(depNode)[-1]) # extracts the list of all parameters, for the distribution of depNode
timesFound <- 0 ## for success, we'll find targetNode in only *one* parameter expression
for(i in seq_along(paramsList)) {
if(!missing(depParamNodeName) && (names(paramsList)[i] == depParamNodeName)) {
expr <- depNodeExprExpanded
} else { expr <- cc_expandDetermNodesInExpr(model, paramsList[[i]], targetNode, skipExpansionsNode) }
if(cc_vectorizedComponentCheck(targetNode, expr)) return(FALSE)
if(cc_nodeInExpr(targetNode, expr)) { timesFound <- timesFound + 1 } ## we found 'targetNode'
}
if(timesFound == 0) stop('something went wrong; targetNode not found in any parameter expressions')
if(timesFound == 1) return(TRUE)
if(timesFound > 1) return(FALSE)
}
## determines whether node appears anywhere in expr
cc_nodeInExpr <- function(node, expr) { return(node %in% cc_getNodesInExpr(expr)) }
## determines which nodes apppear in an expression
## exporting as used in setup code of a nf called directly by user
#' @export
cc_getNodesInExpr <- function(expr) {
if(is.numeric(expr)) return(character(0)) ## expr is numeric
if(is.logical(expr)) return(character(0)) ## expr is logical
if(is.name(expr) || (is.call(expr) && (expr[[1]] == '[') && is.name(expr[[2]]))) return(safeDeparse(expr, warn = TRUE)) ## expr is a node name
if(is.call(expr)) return(unlist(lapply(expr[-1], cc_getNodesInExpr))) ## expr is some general call
stop(paste0('something went wrong processing: ', safeDeparse(expr)))
}
## if targetNode is vectorized: determines if any components of targetNode appear in expr
cc_vectorizedComponentCheck <- function(targetNode, expr) {
if(!is.vectorized(targetNode)) return(FALSE)
targetNodesExpanded <- nl_expandNodeIndex(targetNode)
if(any(unlist(lapply(targetNodesExpanded, function(node) cc_nodeInExpr(node, expr))))) return(TRUE) ## any components of *vectorized* targetNode appear in expr
return(FALSE)
}
##############################################################################################
##############################################################################################
## general, flexible, check for linear relationship
## returns list(offset, scale), or NULL
##############################################################################################
##############################################################################################
cc_replace01 <- function(expr) {
## replace {0,1} with {0,1}+1i so distinguished from offset,scale of 0,1 when match target
if(is.numeric(expr) && (expr %in% c(0, 1)))
return(expr + 1i)
if(is.call(expr))
for(i in 2:length(expr))
expr[[i]] <- cc_replace01(expr[[i]])
return(expr)
}
cc_checkLinearity <- function(expr, targetNode) {
## targetNode doesn't appear in expr
if(!cc_nodeInExpr(targetNode, expr)) {
if(is.call(expr) && expr[[1]] == '(') return(cc_checkLinearity(expr[[2]], targetNode))
## add +1i to tags 0s and 1s as not being from exact match to target
return(list(offset = cc_replace01(expr), scale = 0))
}
## Check if expr is exactly the targetNode.
## `deparse()` may return multiple lines, which will cause `FALSE` here.
## That should always be correct as multiple lines would only occur
## when `expr` is not a node name. It's possible one could have a node name
## with many indices and large integers as indices (e.g., x[13434, 1342352,...])
## but to get to a length that would get deparsed to multiple lines
## would entail a gigantic object.
if(identical(targetNode, deparse(expr)))
return(list(offset = 0, scale = 1))
if(!is.call(expr)) stop('cc_checkLinearity: expression is not a call object')
## process the expression contents of the parentheses
if(expr[[1]] == '(')
return(cc_checkLinearity(expr[[2]], targetNode))
if(expr[[1]] == '[')
return(cc_checkLinearity(expr[[2]], targetNode))
## Look for individual nodes in vectorized use or other strange cases.
if(expr[[1]] == 'structureExpr') {
## Can't have target appear multiple times.
if(sum(sapply(expr[2:length(expr)], function(x)
cc_nodeInExpr(targetNode, x))) != 1)
return(NULL)
wh <- which(sapply(expr[2:length(expr)], function(x)
cc_nodeInExpr(targetNode, x)))
checkLinearityStrucExpr <- cc_checkLinearity(expr[[wh+1]], targetNode)
if(is.null(checkLinearityStrucExpr)) return(NULL)
if(is.indexed(targetNode) && grepl(':', targetNode)) {
## if this is a structureExpression, and also
## if targetNode is multivariate (is indexed and contains ":"),
## and is also part of a larger structureExpr,
## then we do not handle conjugacy for those cases
return(NULL)
}
return(list(offset = cc_combineExprsAddition(expr, checkLinearityStrucExpr$offset), # was expr?,
scale = checkLinearityStrucExpr$scale))
}
## we'll just have to skip over asRow() and asCol(), so they don't mess up the linearity check
if(expr[[1]] == 'asRow' || expr[[1]] == 'asCol') {
return(cc_checkLinearity(expr[[2]], targetNode))
}
## minus sign: change to a plus sign, and invert the sign of the RHS
if(expr[[1]] == '-') {
if(length(expr) == 3) {
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
return(list(offset = cc_combineExprsSubtraction(checkLinearityLHS$offset, checkLinearityRHS$offset),
scale = cc_combineExprsSubtraction(checkLinearityLHS$scale, checkLinearityRHS$scale)))
}
if(length(expr) == 2) {
checkLin <- cc_checkLinearity(expr[[2]], targetNode)
if(is.null(checkLin)) return(NULL)
return(list(offset = cc_negateExpr(checkLin$offset),
scale = cc_negateExpr(checkLin$scale)))
}
stop('cc_checkLinearity: problem with negation expression')
}
if(expr[[1]] == '+') {
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
return(list(offset = cc_combineExprsAddition(checkLinearityLHS$offset, checkLinearityRHS$offset),
scale = cc_combineExprsAddition(checkLinearityLHS$scale, checkLinearityRHS$scale)))
}
if(expr[[1]] == '*' || expr[[1]] == '%*%' || expr[[1]] == 'inprod' || expr[[1]] == 'sum') {
## X[,] %*% beta[] where beta[i] are nodes is not standard matrix multiplication, so there is offset and scale
isMatrixMult <- ifelse(expr[[1]] == '%*%' &&
length(expr[[2]]) > 1 && expr[[2]][[1]] != 'structureExpr' &&
length(expr[[3]]) > 1 && expr[[3]][[1]] != 'structureExpr',
TRUE, FALSE)
if(expr[[1]] == 'sum')
if(length(expr[[2]]) == 3 && expr[[2]][[1]] == '*') {
tmpExpr <- quote(inprod(a, b))
tmpExpr[[2]] <- expr[[2]][[2]]
tmpExpr[[3]] <- expr[[2]][[3]]
expr <- tmpExpr
} else {
return(NULL) # cases such as sum(p[1:5]); there may be some unusual conjugacy cases here that we don't detect
}
if(length(expr) != 3) return(NULL) # avoid error at next step for unexpected cases
if(cc_nodeInExpr(targetNode, expr[[2]]) && cc_nodeInExpr(targetNode, expr[[3]])) return(NULL)
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
if((checkLinearityLHS$scale != 0) && (checkLinearityRHS$scale != 0)) stop('cc_checkLinearity: incompatible scales in * operation')
if((checkLinearityLHS$scale == 0) && (checkLinearityRHS$scale == 0)) {
return(list(offset = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$offset, isMatrixMult),
scale = 0)) }
if(checkLinearityLHS$scale != 0) {
return(list(offset = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$offset, isMatrixMult),
scale = cc_combineExprsMultiplication(checkLinearityLHS$scale, checkLinearityRHS$offset, isMatrixMult))) }
if(checkLinearityRHS$scale != 0) {
return(list(offset = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$offset, isMatrixMult),
scale = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$scale , isMatrixMult))) }
stop('cc_checkLinearity: something went wrong')
}
if(expr[[1]] == '/') {
if(cc_nodeInExpr(targetNode, expr[[3]])) return(NULL)
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
if(checkLinearityLHS$scale == 0) stop('left hand side scale == 0')
if(checkLinearityRHS$scale != 0) stop('right hand side scale is non-zero')
return(list(offset = cc_combineExprsDivision(checkLinearityLHS$offset, checkLinearityRHS$offset),
scale = cc_combineExprsDivision(checkLinearityLHS$scale, checkLinearityRHS$offset)))
}
## returns not conjugate (NULL) if we don't recognize the call (expr[[1]])
return(NULL)
}
cc_negateExpr <- function(expr) {
if(expr == 0) return(0)
if(is.numeric(expr)) return(-1*expr)
return(substitute(-EXPR, list(EXPR = expr)))
}
cc_combineExprsAddition <- function(expr1, expr2) {
if((expr1 == 0) && (expr2 == 0)) return(0)
if((expr1 != 0) && (expr2 == 0)) return(expr1)
if((expr1 == 0) && (expr2 != 0)) return(expr2)
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 + expr2)
return(substitute(EXPR1 + EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_combineExprsSubtraction <- function(expr1, expr2) {
if((expr1 == 0) && (expr2 == 0)) return(0)
if((expr1 != 0) && (expr2 == 0)) return(expr1)
if((expr1 == 0) && (expr2 != 0)) return(cc_negateExpr(expr2))
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 - expr2)
return(substitute(EXPR1 - EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_combineExprsMultiplication <- function(expr1, expr2, isMatrixMult) {
if((expr1 == 0) || (expr2 == 0)) return(0)
if((expr1 == 1) && (expr2 == 1)) return(1)
if((expr1 != 1) && (expr2 == 1)) return(expr1)
if((expr1 == 1) && (expr2 != 1)) return(expr2)
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 * expr2)
if(!isMatrixMult) return(substitute(EXPR1 * EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
if( isMatrixMult) return(substitute(EXPR1 %*% EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_combineExprsDivision <- function(expr1, expr2) {
if( (expr2 == 0)) stop('error, division by 0')
if( (expr2 == 1)) return(expr1)
if((expr1 == 0) ) return(0)
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 / expr2)
return(substitute(EXPR1 / EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_checkStickbreaking <- function(expr, targetNode) {
if(!is.call(expr) || expr[[1]] != 'stick_breaking' || !cc_nodeInExpr(targetNode, expr))
return(NULL)
expr <- expr[[2]]
if(!is.call(expr) || expr[[1]] != 'structureExpr')
return(NULL)
## offset <- which(targetNode == cc_getNodesInExpr(expr)) # Issue 1409.
## Need that target is the the entirety of one element of the structureExpr.
offset <- which(sapply(expr[2:length(expr)], function(x) x == targetNode))
if(length(offset) != 1) return(NULL) else return(list(offset = offset, scale = 1))
}
compareConjugacyLists <- function(C1, C2) {
if(identical(C1, C2)) return(TRUE)
if(!identical(names(C1), names(C2))) {cat('Names do not match\n'); return(FALSE)}
for(i in seq_along(C1)) {
if(!identical(C1[[i]]$type, C2[[i]]$type)) cat(paste0('type mismatch for i =',i))
if(!identical(C1[[i]]$target, C2[[i]]$target)) cat(paste0('target mismatch for i =',i))
if(!identical(C1[[i]]$target, C2[[i]]$target)) cat(paste0('target mismatch for i =',i))
if(!identical(names(C1[[i]]$control), names(C2[[i]]$control))) cat(paste0('control names mismatch for i =',i,'. Skipping node comparison'))
else {
for(j in seq_along(C1[[i]]$control)) {
if(!identical(sort(C1[[i]]$control[[j]]), sort(C2[[i]]$control[[j]]))) cat(paste0('target mismatch for i =',i, 'j =', j))
}
}
}
}
createDynamicConjugateSamplerName <- function(prior, dependentCounts, dynamicallyIndexed = FALSE) {
##depString <- paste0(dependentCounts, names(dependentCounts), collapse='_') ## including the numbers of dependents
depString <- paste0(sort(names(dependentCounts)), collapse='_') ## without the numbers of each type of dependent node
paste0('sampler_conjugate_', prior, '_', depString, ifelse(dynamicallyIndexed, '_dynamicDeps', ''))
}
makeDeclareSizeField <- function(firstSize, secondSize, thirdSize, nDim) {
eval(substitute(switch(as.character(nDim),
`0` = quote(FIRSTSIZE),
`1` = quote(c(FIRSTSIZE, SECONDSIZE)),
`2` = quote(c(FIRSTSIZE, SECONDSIZE, THIRDSIZE)),
stop()),
list(FIRSTSIZE = firstSize,
SECONDSIZE = secondSize,
THIRDSIZE = thirdSize)))
}
makeIndexedVariable <- function(varName, nDim, indexExpr, secondSize, thirdSize) {
eval(substitute(switch(as.character(nDim),
`0` = quote(VARNAME[INDEXEXPR]),
`1` = quote(VARNAME[INDEXEXPR, 1:SECONDSIZE]),
`2` = quote(VARNAME[INDEXEXPR, 1:SECONDSIZE, 1:THIRDSIZE]),
stop()),
list(VARNAME = varName,
INDEXEXPR = indexExpr,
SECONDSIZE = secondSize,
THIRDSIZE = thirdSize)))
}
##############################################################################################
##############################################################################################
## create object: conjugacyRelationshipsObject
## also, generate all conjugate sampler nimbleFunctions
## and a function to rebuild conjugate sampler functions
##############################################################################################
##############################################################################################
## this is still *necessary* (and exported):
conjugacyRelationshipsObject <- conjugacyRelationshipsClass(conjugacyRelationshipsInputList)
## here after is for handling of dynamic conjugate sampler function
dynamicConjugateSamplerDefinitionsEnv <- new.env()
dynamicConjugateSamplerFunctionsEnv <- new.env()
dynamicConjugateSamplerExists <- function(name) {
return(name %in% ls(dynamicConjugateSamplerDefinitionsEnv))
}
dynamicConjugateSamplerAdd <- function(name, def) {
dynamicConjugateSamplerDefinitionsEnv[[name]] <- def
dynamicConjugateSamplerFunctionsEnv[[name]] <- eval(def)
}
dynamicConjugateSamplerGet <- function(name) {
return(dynamicConjugateSamplerFunctionsEnv[[name]])
}
dynamicConjugateSamplerWrite <- function(file = 'TEMP_dynamicConjugateSamplerDefinitions.R') {
## environment to create functions in is throw-away, since they've all already been created in dynamicConjugateSamplerFunctionsEnv
createNamedObjectsFromList(as.list(dynamicConjugateSamplerDefinitionsEnv), writeToFile = file, envir = new.env())
}
Try the nimble package in your browser
Any scripts or data that you put into this service are public.
nimble documentation built on Sept. 11, 2024, 7:10 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.
Defines functions dynamicConjugateSamplerWrite dynamicConjugateSamplerGet dynamicConjugateSamplerAdd dynamicConjugateSamplerExists makeIndexedVariable makeDeclareSizeField createDynamicConjugateSamplerName compareConjugacyLists cc_checkStickbreaking cc_combineExprsDivision cc_combineExprsMultiplication cc_combineExprsSubtraction cc_combineExprsAddition cc_negateExpr cc_checkLinearity cc_replace01 cc_vectorizedComponentCheck cc_getNodesInExpr cc_nodeInExpr cc_otherParamsCheck cc_stripExpr cc_checkScalar cc_linkCheck cc_createStructureExpr cc_expandDetermNodesInExpr cc_makeRDistributionName cc_makeConjugateSamplerName cc_makeSamplerTypeName
Documented in cc_getNodesInExpr
conjugacyRelationshipsInputList <- list(
## beta
list(prior = 'dbeta',
link = 'identity',
dependents = list(
dbern = list(param = 'prob', contribution_shape1 = 'value', contribution_shape2 = '1 - value' ),
dbin = list(param = 'prob', contribution_shape1 = 'value', contribution_shape2 = 'size - value'),
dnegbin = list(param = 'prob', contribution_shape1 = 'size', contribution_shape2 = 'value' ),
## using 'stickbreaking' not 'stick_breaking' as link is now appended to dist with a '_' in realized conjugacy system
## Note that 'stick_breaking' is still the user-facing **function** name.
dcat = list(param = 'prob', link = 'stickbreaking', contribution_shape1 = 'value == offset',
contribution_shape2 = 'value > offset')), ## offset is set to be where in the broken stick the target is
posterior = 'dbeta(shape1 = prior_shape1 + contribution_shape1,
shape2 = prior_shape2 + contribution_shape2)'),
## dirichlet
list(prior = 'ddirch',
link = 'identity',
dependents = list(
dmulti = list(param = 'prob', contribution_alpha = 'value'),
dcat = list(param = 'prob', contribution_alpha = 'calc_dcatConjugacyContributions(k, value)')),
posterior = 'ddirch(alpha = prior_alpha + contribution_alpha)'),
## flat
list(prior = 'dflat',
link = 'linear',
dependents = list(
dnorm = list(param = 'mean', contribution_mean = 'coeff * (value-offset) * tau', contribution_tau = 'coeff^2 * tau'),
dlnorm = list(param = 'meanlog', contribution_mean = 'coeff * (log(value)-offset) * taulog', contribution_tau = 'coeff^2 * taulog')),
posterior = 'dnorm(mean = contribution_mean / contribution_tau,
sd = contribution_tau^(-0.5))'),
## halfflat - first possible conjugacy; user can achieve this with gamma(1, 0) as we handle that
## list(prior = 'dhalfflat',
## link = 'multiplicative',
## dependents = list(
## dpois = list(param = 'lambda', contribution_shape = 'value', contribution_rate = 'coeff' ),
## dnorm = list(param = 'tau', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (value-mean)^2' ),
## dlnorm = list(param = 'taulog', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (log(value)-meanlog)^2'),
## dgamma = list(param = 'rate', contribution_shape = 'shape', contribution_rate = 'coeff * value' ),
## dinvgamma = list(param = 'scale', contribution_shape = 'shape', contribution_rate = 'coeff / value' ),
## dexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * value' ),
## dweib = list(param = 'lambda', contribution_shape = '1', contribution_rate = 'coeff * value^shape' )),
## ## ddexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * abs(value-location)' )
## ## dpar = list(...) ## contribution_shape=1; contribution_rate=coeff*log(value/c) 'c is 2nd param of pareto'
## posterior = 'dgamma(shape = 1 + contribution_shape,
## scale = 1 / contribution_rate)'),
## halfflat - second possible conjugacy; note that user will be able to do invgamma(-1,0) and
## get conjugacy once we fix issue #314
## list(prior = 'dinvgamma',
## link = 'multiplicative',
## dependents = list(
## dnorm = list(param = 'var', contribution_shape = '1/2', contribution_scale = '(value-mean)^2 / (coeff * 2)' ),
## dlnorm = list(param = 'varlog', contribution_shape = '1/2', contribution_scale = '(log(value)-meanlog)^2 / (coeff*2)'),
## dgamma = list(param = 'scale', contribution_shape = 'shape', contribution_scale = 'value / coeff' ),
## dinvgamma = list(param = 'rate', contribution_shape = 'shape', contribution_scale = '1 / (coeff * value)' ),
## dexp = list(param = 'scale', contribution_shape = '1', contribution_scale = 'value / coeff' )),
## ## add ddexp
## posterior = 'dinvgamma(shape = -1 + contribution_shape,
## rate = 1 / contribution_scale)'),
## halfflat - third possible conjugacy - we use this because it corresponds to the
## Gelman (2006) recommended uniform on sd scale prior for variance components
## and current NIMBLE conjugacy system only allows one possible form of conjugacy.
## Note that sd ~ U(0,Inf) equivalent to var ~ IG(-1/2, 0).
## Also note that if conj system could detect 'squared' dependency, then
## we could allow dnorm with param = 'var'.
list(prior = 'dhalfflat',
link = 'multiplicative',
dependents = list(
dnorm = list(param = 'sd', contribution_shape = '1/2', contribution_scale = '(value-mean)^2 / (coeff^2 * 2)' ),
dlnorm = list(param = 'sdlog', contribution_shape = '1/2', contribution_scale = '(log(value)-meanlog)^2 / (coeff^2 * 2)')),
posterior = 'dsqrtinvgamma(shape = -1/2 + contribution_shape,
rate = 1 / contribution_scale)'),
## gamma
list(prior = 'dgamma',
link = 'multiplicative',
dependents = list(
dpois = list(param = 'lambda', contribution_shape = 'value', contribution_rate = 'coeff' ),
dnorm = list(param = 'tau', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (value-mean)^2' ),
dlnorm = list(param = 'taulog', contribution_shape = '1/2', contribution_rate = 'coeff/2 * (log(value)-meanlog)^2'),
dgamma = list(param = 'rate', contribution_shape = 'shape', contribution_rate = 'coeff * value' ),
dinvgamma = list(param = 'scale', contribution_shape = 'shape', contribution_rate = 'coeff / value' ),
dexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * value' ),
dweib = list(param = 'lambda', contribution_shape = '1', contribution_rate = 'coeff * value^shape' ),
ddexp = list(param = 'rate', contribution_shape = '1', contribution_rate = 'coeff * abs(value-location)' )),
## dpar = list(...) ## contribution_shape=1; contribution_rate=coeff*log(value/c) 'c is 2nd param of pareto'
posterior = 'dgamma(shape = prior_shape + contribution_shape,
scale = 1 / (prior_rate + contribution_rate))'),
## invgamma
list(prior = 'dinvgamma',
link = 'multiplicative',
dependents = list(
dnorm = list(param = 'var', contribution_shape = '1/2', contribution_scale = '(value-mean)^2 / (coeff * 2)' ),
dlnorm = list(param = 'varlog', contribution_shape = '1/2', contribution_scale = '(log(value)-meanlog)^2 / (coeff*2)'),
dgamma = list(param = 'scale', contribution_shape = 'shape', contribution_scale = 'value / coeff' ),
dinvgamma = list(param = 'rate', contribution_shape = 'shape', contribution_scale = '1 / (coeff * value)' ),
dexp = list(param = 'scale', contribution_shape = '1', contribution_scale = 'value / coeff' ),
ddexp = list(param = 'scale', contribution_shape = '1', contribution_scale = 'abs(value-location) / coeff' )),
posterior = 'dinvgamma(shape = prior_shape + contribution_shape,
rate = 1 / (prior_scale + contribution_scale))'),
## normal
list(prior = 'dnorm',
link = 'linear',
dependents = list(
dnorm = list(param = 'mean', contribution_mean = 'coeff * (value-offset) * tau', contribution_tau = 'coeff^2 * tau' ),
dlnorm = list(param = 'meanlog', contribution_mean = 'coeff * (log(value)-offset) * taulog', contribution_tau = 'coeff^2 * taulog')),
posterior = 'dnorm(mean = (prior_mean*prior_tau + contribution_mean) / (prior_tau + contribution_tau),
sd = (prior_tau + contribution_tau)^(-0.5))'),
#####
## pareto
## these are idiosyncratic enough that we probably want to skip them
## list(prior = 'dpar', ##### waiting for dpar() distribution
## link = 'multiplicative',
## dependents = list(
## dunif = list(param = 'max', contribution_alpha = '1', contribution_c = 'value/coeff'), # only works if 0 is min of the unif so this will be hard to do
## careful with next one - data seem to impose upper bound on posterior, so not clear this goes through
## dpar = list(param = 'c', contribution_alpha = '-alpha'), contribution_c = '0'),
## posterior = 'dpar(alpha = prior_alpha + contribution_alpha,
## c = max(prior_c, contribution_c)'),
#####
## multivariate-normal; note that to avoid unnecessary computations, there is special processing of the contribution_mean and contribution_prec
## in the generation of the conjugate sampler run method, so the specification here should not be changed without looking at that processing.
list(prior = 'dmnorm',
link = 'linear',
dependents = list(
##dmnorm = list(param = 'mean', contribution_mean = '(t(coeff) %*% prec %*% asCol(value-offset))[,1]', contribution_prec = 't(coeff) %*% prec %*% coeff')),
dmnorm = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)')),
## LINK will be replaced with appropriate link via code processing
## original less efficient posterior definition:
## posterior = 'dmnorm_chol(mean = (inverse(prior_prec + contribution_prec) %*% (prior_prec %*% asCol(prior_mean) + asCol(contribution_mean)))[,1],
## cholesky = chol(prior_prec + contribution_prec),
## prec_param = 1)'),
posterior = '{ R <- chol(prior_prec + contribution_prec)
A <- prior_prec %*% asCol(prior_mean) + asCol(contribution_mean)
mu <- backsolve(R, forwardsolve(t(R), A))[,1]
dmnorm_chol(mean = mu, cholesky = R, prec_param = 1) }'),
## wishart
list(prior = 'dwish',
## changing to only use link='identity' case, since the link='linear' case was not correct.
## -DT March 2017
## link = 'linear',
link = 'multiplicativeScalar', # we only handle scalar 'coeff'; this naming is slightly awkward since for univar dists, link is of course scalar
dependents = list(
## parentheses added to the contribution_R calculation:
## colVec * (rowVec * matrix)
## Chris is checking to see whether this makes a difference for Eigen
## -DT April 2016
## changing to only use link='identity' case, since the link='linear' case was not correct
## -DT March 2017
## dmnorm = list(param = 'prec', contribution_R = 'asCol(value-mean) %*% (asRow(value-mean) %*% coeff)', contribution_df = '1')),
dmnorm = list(param = 'prec', contribution_R = 'coeff * asCol(value-mean) %*% asRow(value-mean)', contribution_df = '1')),
posterior = 'dwish_chol(cholesky = chol(prior_R + contribution_R),
df = prior_df + contribution_df,
scale_param = 0)'),
## inverse wishart
list(prior = 'dinvwish',
link = 'multiplicativeScalar', # we only handle scalar 'coeff'; this naming is slightly awkward since for univar dists, link is of course scalar
dependents = list(
dmnorm = list(param = 'cov', contribution_S = 'asCol(value-mean) %*% asRow(value-mean) / coeff', contribution_df = '1')),
posterior = 'dinvwish_chol(cholesky = chol(prior_S + contribution_S),
df = prior_df + contribution_df,
scale_param = 1)')
)
##############################################################################################
##############################################################################################
## reference class definitions:
## conjugacyRelationshipsClass
## conjugacyClass
## dependentClass
## posteriorClass
##############################################################################################
##############################################################################################
conjugacyRelationshipsClass <- setRefClass(
Class = 'conjugacyRelationshipsClass',
fields = list(
conjugacys = 'ANY' ## a (named) list of conjugacyClass objects, each describes the conjugacies for a particular prior distribution (name is prior distribution name)
),
methods = list(
initialize = function(crl) {
conjugacys <<- list()
for(i in seq_along(crl)) {
conjugacys[[i]] <<- conjugacyClass(crl[[i]])
}
names(conjugacys) <<- unlist(lapply(conjugacys, function(cr) cr$prior))
},
checkConjugacy = function(model, nodeIDs, restrictLink = NULL) {
if(isTRUE(getNimbleOption("oldConjugacyChecking")))
checkConjugacy_original(model, nodeIDs, restrictLink)
else
checkConjugacy_new(model, nodeIDs, restrictLink)
},
checkConjugacy_new = function(model, nodeIDs, restrictLink = NULL) {
maps <- model$modelDef$maps
nodeDeclIDs <- maps$graphID_2_declID[nodeIDs] ## declaration IDs of the nodeIDs
declID2nodeIDs <- split(nodeIDs, nodeDeclIDs) ## nodeIDs grouped by declarationID
ansList <- ansList2 <- list()
for(i in seq_along(declID2nodeIDs)) { ## For each group of nodeIDs from the same declarationID
nodeIDsFromOneDecl <- declID2nodeIDs[[i]]
firstNodeName <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[1]]
if(model$isTruncated(firstNodeName)) next ## we say non-conjugate if the targetNode is truncated
dist <- model$getDistribution(firstNodeName)
conjugacyObj <- conjugacys[[dist]]
if(is.null(conjugacyObj)) next
# NO: insert logic here to check a single dependency and do next if can't be conjugate
#model$getDependencies('mu[1]',self=F,stochOnly=T)
#for( loop through deps )
# conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode)
# if(not conj) next
# now try to guess if finding paths will be more intensive than simply looking at target-dependent pairs, to avoid path finding when there is nested structure such as stickbreaking
deps <- lapply(nodeIDsFromOneDecl, function(x) model$getDependencies(x, stochOnly = TRUE, self = FALSE))
numDeps <- sapply(deps, length)
sumNumDeps <- sum(numDeps)
maxNumPaths <- 0
## We only need check until we get to as many paths as `sumNumDeps`.
## This also avoids integer overflow that can occur with recursive indexing.
for(nodeID in nodeIDsFromOneDecl) {
numPaths <- model$getDependencyPathCountOneNode(nodeID, max = sumNumDeps + 1)
if(numPaths > maxNumPaths)
maxNumPaths <- numPaths
if(maxNumPaths > sumNumDeps)
break
}
if(maxNumPaths > sumNumDeps) {
# max(numPaths) is reasonable guess at number of unique (by node) paths (though it overestimates number of unique (by declaration ID) paths; if we have to evaluate conjugacy for more paths than we would by simply looking at all pairs of target-dependent nodes, then just use node pairs
# note that it's not clear what criterion to use here since computational time is combination of time for finding all paths and then for evaluating conjugacy for unique (by declaration ID) paths, but the hope is to make a crude cut here that avoids path calculations when there would be a lot of them
ansList[[length(ansList)+1]] <- lapply(seq_along(nodeIDsFromOneDecl),
function(index) {
targetNode <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[index]]
depEnds <- deps[[index]]
depTypes <- sapply(depEnds, function(x) conjugacyObj$checkConjugacyOneDep(model, targetNode, x, restrictLink))
if(!length(depTypes)) return(NULL)
if(!any(sapply(depTypes, is.null))) {
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
depEnds[boolMatch]
})
names(control) <- uniqueDepTypes
return(list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control))
} else return(NULL)
})
names(ansList[[length(ansList)]]) <- maps$graphID_2_nodeName[nodeIDsFromOneDecl]
} else {
# determine conjugacy based on unique (by declaration ID) paths
depPathsByNode <- lapply(nodeIDsFromOneDecl, model$getDependencyPaths) ## make list (by nodeID) of lists of paths through graph
depPathsByNode <- depPathsByNode[!unlist(lapply(depPathsByNode, function(x) is.null(x) || (length(x)==0)))]
depPathsByNodeLabels <- lapply(depPathsByNode, function(z) ## make character labels that match for same path through graph
unlist(lapply(z,
function(x)
paste(maps$graphID_2_declID[x[,1]], x[,2], collapse = '\r', sep='\r'))))
depPathsByNodeUnlisted <- unlist(depPathsByNode, recursive = FALSE)
depPathsByNodeLabelsUnlisted <- unlist(depPathsByNodeLabels)
## uniquePaths <- unique(depPathsByNodeLabelsUnlisted)
uniquePathsUnlistedIndices <- split(seq_along(depPathsByNodeLabelsUnlisted), depPathsByNodeLabelsUnlisted)
conjDepTypes <- character(length(uniquePathsUnlistedIndices))
for(j in seq_along(uniquePathsUnlistedIndices)) {
firstDepPath <- depPathsByNodeUnlisted[[ uniquePathsUnlistedIndices[[j]][1] ]]
targetNode <- maps$graphID_2_nodeName[firstDepPath[1,1]]
depNode <- maps$graphID_2_nodeName[firstDepPath[nrow(firstDepPath), 1]]
oneDepType <- conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode, restrictLink)
conjDepTypes[j] <- if(is.null(oneDepType)) "" else oneDepType
}
conjBool <- conjDepTypes != ""
names(conjDepTypes) <- names(conjBool) <- names(uniquePathsUnlistedIndices)
if(any(conjBool)) {
targetNodes <- unlist(lapply(depPathsByNode, function(x) if(is.null(x)) '_NO_DEPS_' else maps$graphID_2_nodeName[x[[1]][1,1]]))
ansList[[length(ansList)+1]] <- mapply(
function(targetNode, depPathsOneNode, depPathsLabelsOneNode) {
if(targetNode == '_NO_DEPS_') return(NULL) ## these should have already been weeded out
if(all(conjBool[depPathsLabelsOneNode])) {
depTypes <- conjDepTypes[depPathsLabelsOneNode]
depEnds <- maps$graphID_2_nodeName[ unlist(lapply(depPathsOneNode, function(x) x[nrow(x)])) ]
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
## prevent multiple instances of same dependent node name
## (via different graph dependency paths -DT Oct 2016
##depEnds[boolMatch]
unique(depEnds[boolMatch])
})
names(control) <- uniqueDepTypes
list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control)
}
},
targetNodes, depPathsByNode, depPathsByNodeLabels, USE.NAMES = TRUE, SIMPLIFY = FALSE)
}
}
}
if(length(ansList) > 0) ansList <- do.call('c', ansList)
ansList <- ansList[!sapply(ansList, is.null)] # strips out any NULL values
if(!length(ansList)) return(list()) else return(ansList) # replaces empty named list with empty list
},
checkConjugacy_original = function(model, nodeIDs, restrictLink = NULL) {
maps <- model$modelDef$maps
nodeDeclIDs <- maps$graphID_2_declID[nodeIDs] ## declaration IDs of the nodeIDs
declID2nodeIDs <- split(nodeIDs, nodeDeclIDs) ## nodeIDs grouped by declarationID
ansList <- ansList2 <- list()
for(i in seq_along(declID2nodeIDs)) { ## For each group of nodeIDs from the same declarationID
nodeIDsFromOneDecl <- declID2nodeIDs[[i]]
firstNodeName <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[1]]
if(model$isTruncated(firstNodeName)) next ## we say non-conjugate if the targetNode is truncated
dist <- model$getDistribution(firstNodeName)
conjugacyObj <- conjugacys[[dist]]
if(is.null(conjugacyObj)) next
# NO: insert logic here to check a single dependency and do next if can't be conjugate
#model$getDependencies('mu[1]',self=F,stochOnly=T)
#for( loop through deps )
# conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode)
# if(not conj) next
# now try to guess if finding paths will be more intensive than simply looking at target-dependent pairs, to avoid path finding when there is nested structure such as stickbreaking
numPaths <- sapply(nodeIDsFromOneDecl, model$getDependencyPathCountOneNode)
deps <- lapply(nodeIDsFromOneDecl, function(x) model$getDependencies(x, stochOnly = TRUE, self = FALSE))
numDeps <- sapply(deps, length)
if(max(numPaths) > sum(numDeps)) {
# max(numPaths) is reasonable guess at number of unique (by node) paths (though it overestimates number of unique (by declaration ID) paths; if we have to evaluate conjugacy for more paths than we would by simply looking at all pairs of target-dependent nodes, then just use node pairs
# note that it's not clear what criterion to use here since computational time is combination of time for finding all paths and then for evaluating conjugacy for unique (by declaration ID) paths, but the hope is to make a crude cut here that avoids path calculations when there would be a lot of them
ansList[[length(ansList)+1]] <- lapply(seq_along(nodeIDsFromOneDecl),
function(index) {
targetNode <- maps$graphID_2_nodeName[nodeIDsFromOneDecl[index]]
depEnds <- deps[[index]]
depTypes <- sapply(depEnds, function(x) conjugacyObj$checkConjugacyOneDep(model, targetNode, x, restrictLink))
if(!length(depTypes)) return(NULL)
if(!any(sapply(depTypes, is.null))) {
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
depEnds[boolMatch]
})
names(control) <- uniqueDepTypes
return(list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control))
} else return(NULL)
})
names(ansList[[length(ansList)]]) <- maps$graphID_2_nodeName[nodeIDsFromOneDecl]
} else {
# determine conjugacy based on unique (by declaration ID) paths
depPathsByNode <- lapply(nodeIDsFromOneDecl, getDependencyPaths, maps = maps) ## make list (by nodeID) of lists of paths through graph
depPathsByNode <- depPathsByNode[!unlist(lapply(depPathsByNode, function(x) is.null(x) || (length(x)==0)))]
depPathsByNodeLabels <- lapply(depPathsByNode, function(z) ## make character labels that match for same path through graph
unlist(lapply(z,
function(x)
paste(maps$graphID_2_declID[x[,1]], x[,2], collapse = '\r', sep='\r'))))
depPathsByNodeUnlisted <- unlist(depPathsByNode, recursive = FALSE)
depPathsByNodeLabelsUnlisted <- unlist(depPathsByNodeLabels)
## uniquePaths <- unique(depPathsByNodeLabelsUnlisted)
uniquePathsUnlistedIndices <- split(seq_along(depPathsByNodeLabelsUnlisted), depPathsByNodeLabelsUnlisted)
conjDepTypes <- character(length(uniquePathsUnlistedIndices))
for(j in seq_along(uniquePathsUnlistedIndices)) {
firstDepPath <- depPathsByNodeUnlisted[[ uniquePathsUnlistedIndices[[j]][1] ]]
targetNode <- maps$graphID_2_nodeName[firstDepPath[1,1]]
depNode <- maps$graphID_2_nodeName[firstDepPath[nrow(firstDepPath), 1]]
oneDepType <- conjugacyObj$checkConjugacyOneDep(model, targetNode, depNode, restrictLink)
conjDepTypes[j] <- if(is.null(oneDepType)) "" else oneDepType
}
conjBool <- conjDepTypes != ""
names(conjDepTypes) <- names(conjBool) <- names(uniquePathsUnlistedIndices)
if(any(conjBool)) {
targetNodes <- unlist(lapply(depPathsByNode, function(x) if(is.null(x)) '_NO_DEPS_' else maps$graphID_2_nodeName[x[[1]][1,1]]))
ansList[[length(ansList)+1]] <- mapply(
function(targetNode, depPathsOneNode, depPathsLabelsOneNode) {
if(targetNode == '_NO_DEPS_') return(NULL) ## these should have already been weeded out
if(all(conjBool[depPathsLabelsOneNode])) {
depTypes <- conjDepTypes[depPathsLabelsOneNode]
depEnds <- maps$graphID_2_nodeName[ unlist(lapply(depPathsOneNode, function(x) x[nrow(x)])) ]
uniqueDepTypes <- unique(depTypes)
control <- lapply(uniqueDepTypes,
function(oneType) {
boolMatch <- depTypes == oneType
## prevent multiple instances of same dependent node name
## (via different graph dependency paths -DT Oct 2016
##depEnds[boolMatch]
unique(depEnds[boolMatch])
})
names(control) <- uniqueDepTypes
list(prior = conjugacyObj$prior, type = conjugacyObj$samplerType, target = targetNode, control = control)
}
},
targetNodes, depPathsByNode, depPathsByNodeLabels, USE.NAMES = TRUE, SIMPLIFY = FALSE)
}
}
}
if(length(ansList) > 0) ansList <- do.call('c', ansList)
ansList <- ansList[!sapply(ansList, is.null)] # strips out any NULL values
if(!length(ansList)) return(list()) else return(ansList) # replaces empty named list with empty list
},
generateDynamicConjugateSamplerDefinition = function(prior, dependentCounts, doDependentScreen = FALSE) {
## conjugateSamplerDefinitions[[paste0('sampler_conjugate_', conjugacyResult$prior)]] ## using original (non-dynamic) conjugate sampler functions
conjugacys[[prior]]$generateConjugateSamplerDef(dynamic = TRUE, dependentCounts = dependentCounts,
doDependentScreen = doDependentScreen)
}
)
)
setMethod(
'[[',
'conjugacyRelationshipsClass',
function(x, i) return(x$conjugacys[[i]])
)
conjugacyClass <- setRefClass(
Class = 'conjugacyClass',
fields = list(
samplerType = 'ANY', ## name of the sampler for this conjugacy class, e.g. 'conjugate_dnorm'
prior = 'ANY', ## name of the prior distribution, e.g. 'dnorm'
link = 'ANY', ## the link ('linear', 'additive', 'multiplicative', 'multiplicativeScalar', or 'identity')
dependents = 'ANY', ## (named) list of dependentClass objects, each contains conjugacy information specific to a particular sampling distribution (name is sampling distribution name)
dependentDistNames = 'ANY', ## character vector of the names of all allowable dependent sampling distributions. same as: names(dependents)
posteriorObject = 'ANY', ## an object of posteriorClass
## needsLinearityCheck = 'ANY', ## logical specifying whether we need to do the linearity check; if the link is 'multiplicative' or 'linear'
## needsStickbreakingCheck = 'ANY', ## logical specifying whether we need to do the stickbreaking check
model = 'ANY', ## these fields ONLY EXIST TO PREVENT A WARNING for '<<-',
DEP_VALUES_VAR_INDEXED = 'ANY', ## in the code for generating the conjugate sampler functions
DEP_PARAM_VAR_INDEXED = 'ANY', ##
DEP_OFFSET_VAR = 'ANY', ##
DEP_COEFF_VAR = 'ANY', ##
CONTRIB_NAME = 'ANY' ##
),
methods = list(
initialize = function(cr) {
dependents <<- list()
samplerType <<- cc_makeSamplerTypeName(cr$prior)
prior <<- cr$prior
link <<- cr$link
initialize_addDependents(cr$dependents)
## removed because link info in specific conjugacy can now override default link
## needsLinearityCheck <<- link %in% c('multiplicative', 'linear')
## needsStickbreakingCheck <<- link %in% c('stickbreaking')
posteriorObject <<- posteriorClass(cr$posterior, prior)
},
initialize_addDependents = function(depList) {
for(i in seq_along(depList)) {
dependents[[i]] <<- dependentClass(depList[[i]], names(depList)[i])
}
names(dependents) <<- names(depList)
dependentDistNames <<- names(dependents)
},
## used by new checkConjugacy() system
## see checkConjugacy for more explanation of each step
checkConjugacyOneDep = function(model, targetNode, depNode, restrictLink = NULL) {
if(model$getDistribution(targetNode) != prior) return(NULL) # check prior distribution of targetNode
if(model$isTruncated(depNode)) return(NULL) # if depNode is truncated, then not conjugate
depNodeDist <- model$getDistribution(depNode)
if(!(depNodeDist %in% dependentDistNames)) return(NULL) # check sampling distribution of depNode
dependentObj <- dependents[[depNodeDist]]
if(is.null(restrictLink)) {
if(!is.null(dependentObj$link)) currentLink <- dependentObj$link else currentLink <- link # handle multiple link case introduced for beta stickbreaking
} else currentLink <- restrictLink
if(currentLink != 'stickbreaking') {
depNodeParamName <- dependentObj$param
linearityCheckExprRaw <- model$getParamExpr(depNode, depNodeParamName) # extracts the expression for 'param' from 'depNode'
linearityCheckExpr <- cc_expandDetermNodesInExpr(model, linearityCheckExprRaw, targetNode = targetNode)
if(!cc_nodeInExpr(targetNode, linearityCheckExpr)) return(NULL)
if(cc_vectorizedComponentCheck(targetNode, linearityCheckExpr)) return(NULL) # if targetNode is vectorized, make sure none of its components appear in expr
linearityCheck <- cc_checkLinearity(linearityCheckExpr, targetNode) # determines whether paramExpr is linear in targetNode
realizedLink <- cc_linkCheck(linearityCheck, currentLink)
if(is.null(realizedLink)) return(NULL)
## ensure targetNode appears in only *one* depNode parameter expression
if(!cc_otherParamsCheck(model, depNode, targetNode, depNodeExprExpanded = linearityCheckExpr, depParamNodeName = depNodeParamName)) return(NULL)
} else {
stickbreakingCheckExpr <- model$getParamExpr(depNode, dependentObj$param) # extracts the expression for 'param' from 'depNode'
stickbreakingCheckExpr <- cc_expandDetermNodesInExpr(model, stickbreakingCheckExpr, targetNode = targetNode)
if(is.null(cc_checkStickbreaking(stickbreakingCheckExpr, targetNode))) return(NULL)
realizedLink <- 'stickbreaking'
}
return(paste0('dep_', depNodeDist, '_', realizedLink))
},
addDependentNodeToControl = function(control, depNodeDist, depNode) {
listName <- paste0('dep_', depNodeDist)
control[[listName]] <- c(control[[listName]], depNode)
control
},
## workhorse for creating conjugate sampler nimble functions
generateConjugateSamplerDef = function(dynamic = FALSE, dependentCounts, doDependentScreen = FALSE) {
if(!dynamic) stop('something went wrong, should never have dynamic = FALSE here')
substitute(
nimbleFunction(contains = sampler_BASE,
setup = SETUPFUNCTION,
run = RUNFUNCTION,
methods = list(getPosteriorLogDensity = GETPOSTERIORLOGDENSITYFUNCTION,
reset = function() {})
),
list(SETUPFUNCTION = genSetupFunction(dependentCounts = dependentCounts, doDependentScreen = doDependentScreen),
RUNFUNCTION = genRunFunction(dependentCounts = dependentCounts, doDependentScreen = doDependentScreen),
GETPOSTERIORLOGDENSITYFUNCTION = genGetPosteriorLogDensityFunction(dependentCounts = dependentCounts, doDependentScreen = doDependentScreen)
)
)
},
genSetupFunction = function(dependentCounts, doDependentScreen = FALSE) {
functionBody <- codeBlockClass()
functionBody$addCode({
calcNodes <- model$getDependencies(target)
calcNodesDeterm <- model$getDependencies(target, determOnly = TRUE)
})
## if this conjugate sampler is for a multivariate node (i.e., nDim > 0), then we need to determine the size (d)
if(getDimension(prior) > 0) {
functionBody$addCode(d <- max(determineNodeIndexSizes(target)))
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], "_")[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
functionBody$addCode({
DEP_NODENAMES <- control$DEP_CONTROL_NAME
N_DEP <- length(control$DEP_CONTROL_NAME)
}, list(DEP_NODENAMES = as.name(paste0( 'dep_', distLinkName, '_nodeNames')),
N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_CONTROL_NAME = as.name(paste0( 'dep_', distLinkName))))
if(any(getDimension(distName, includeParams = TRUE) > 0)) {
if(getDimension(distName) > 0) { ## usual case: dependent node is multivariate -> get sizes from dependent node names
functionBody$addCode({
DEP_NODESIZES <- sapply(DEP_NODENAMES, function(node) max(determineNodeIndexSizes(node)), USE.NAMES = FALSE)
}, list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames'))))
} else { ## unusual case: dependent is univariate -> get sizes from dependent node's param (e.g., ddirch-dcat conjugacy)
functionBody$addCode({
DEP_NODESIZES <- sapply(DEP_NODENAMES, function(node) max(nimDim(model$getParam(node, MULTIVARIATE_PARAM_NAME))), USE.NAMES = FALSE)
}, list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
MULTIVARIATE_PARAM_NAME = names(which.max(getDimension(distName, includeParams = TRUE)))))
}
functionBody$addCode({
if(length(DEP_NODESIZES) == 1) DEP_NODESIZES <- c(DEP_NODESIZES, -1) ## guarantee to be a vector, for indexing and size processing
DEP_NODESIZEMAX <- max(DEP_NODESIZES)
}, list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_NODESIZEMAX = as.name(paste0('dep_', distLinkName, '_nodeSizeMax'))))
}
## declare() statements are removed from run() code,
## and were replaced with setup output array() calls below.
## July 2017
depNodeValueNdim <- getDimension(distName)
functionBody$addCode(DEP_VALUES_VAR <- array(0, dim = DECLARE_SIZE),
list(DEP_VALUES_VAR = as.name(paste0('dep_', distLinkName, '_values')),
DECLARE_SIZE = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), depNodeValueNdim)))
neededParams <- dependents[[distName]]$neededParamsForPosterior
for(param in neededParams) {
depNodeParamNdim <- getDimension(distName, param)
functionBody$addCode(DEP_PARAM_VAR <- array(0, dim = DECLARE_SIZE),
list(DEP_PARAM_VAR = as.name(paste0('dep_', distLinkName, '_', param)),
DECLARE_SIZE = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), depNodeParamNdim)))
}
}
## more new array() setup outputs, instead of declare() statements, for offset and coeff variables
## July 2017
targetNdim <- getDimension(prior)
targetCoeffNdim <- switch(as.character(targetNdim), `0`=0, `1`=2, `2`=2, stop())
if(targetCoeffNdim == 2 && link == 'multiplicativeScalar') ## Handles wish/invwish. There are no cases where we allow non-scalar 'coeff'.
targetCoeffNdim <- 0
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], "_")[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('additive', 'multiplicative', 'multiplicativeScalar', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
## the 2's here are *only* to prevent warnings about assigning into member variable names using
inputList <- list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')), ## local assignment '<-', so changed the names to "...2"
DEP_COEFF_VAR2 = as.name(paste0('dep_', distLinkName, '_coeff')), ## so it doesn't recognize the ref class field name
DECLARE_SIZE_OFFSET = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), targetNdim),
DECLARE_SIZE_COEFF = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), quote(d), targetCoeffNdim))
if(currentLink == 'additive')
functionBody$addCode(
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET),
inputList)
if(currentLink %in% c('multiplicative', 'multiplicativeScalar'))
functionBody$addCode(
DEP_COEFF_VAR2 <- array(0, dim = DECLARE_SIZE_COEFF),
inputList)
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET)
DEP_COEFF_VAR2 <- array(0, dim = DECLARE_SIZE_COEFF)
}, inputList)
}
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], "_")[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(!is.null(dependents[[distName]]$link)) currentLink <- dependents[[distName]]$link
if(currentLink == 'stickbreaking') {
functionBody$addCode({
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET) ## the 2's here are *only* to prevent warnings about
}, list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')), ## local assignment '<-', so changed the names to "...2"
DECLARE_SIZE_OFFSET = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), as.name(paste0('dep_', distLinkName, '_nodeSizeMax')), 0))
)
functionBody$addCode({
stickbreakingCheckExpr <- model$getValueExpr(calcNodesDeterm)
stickbreakingCheckExpr <- cc_expandDetermNodesInExpr(model, stickbreakingCheckExpr, targetNode = target)
})
functionBody$addCode(DEP_OFFSET_VAR2 <- rep(cc_checkStickbreaking(stickbreakingCheckExpr, target)$offset, DEP_OFFSET_SIZE),
list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')),
DEP_OFFSET_SIZE = as.name(paste0('N_dep_', distLinkName))))
}
}
## adding declarations for the contribution terms, to remove Windows compiler warnings, DT August 2015
## moved these numeric() and array() declarations for contribution terms to setup outputs, July 2017
for(contributionName in posteriorObject$neededContributionNames) {
contribNdim <- posteriorObject$neededContributionDims[[contributionName]]
functionBody$addCode(CONTRIB_NAME2 <- CONTRIB_INITIAL_DECLARATION, ## the 2's here are *only* to prevent warnings about
list(CONTRIB_NAME2 = as.name(contributionName), ## local assignment '<-' versus '<<-'
CONTRIB_INITIAL_DECLARATION = switch(as.character(contribNdim),
`0` = 0, `1` = quote(rep(0, length = d)), `2` = quote(array(0, dim = c(d, d))), stop())))
}
functionDef <- quote(function(model, mvSaved, target, control) {})
functionDef[[3]] <- functionBody$getCode()
functionDef[[4]] <- NULL ## removes the 'scrref' attribute
return(functionDef)
},
genRunFunction = function(dependentCounts, doDependentScreen = FALSE) {
functionBody <- codeBlockClass()
## only if we're verifying conjugate posterior distributions: get initial targetValue, and modelLogProb -- model$getLogProb(calcNodes)
if(getNimbleOption('verifyConjugatePosteriors')) {
functionBody$addCode({
modelLogProb0 <- model$getLogProb(calcNodes)
origTargetValue <- model[[target]]
})
}
## adds code to generate the quantities prior_xxx, and contribution_xxx
addPosteriorQuantitiesGenerationCode(functionBody = functionBody, dependentCounts = dependentCounts, doDependentScreen = doDependentScreen)
## generate new value, store, calculate, copy, etc...
functionBody$addCode(posteriorObject$prePosteriorCodeBlock, quote = FALSE)
functionBody$addCode({
newTargetValue <- RPOSTERIORCALL
model[[target]] <<- newTargetValue
model$calculate(calcNodes)
nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
}, list(RPOSTERIORCALL = posteriorObject$rCallExpr))
## only if we're verifying conjugate posterior distributions: figure out if conjugate posterior distribution is correct
if(getNimbleOption('verifyConjugatePosteriors')) {
functionBody$addCode({
modelLogProb1 <- model$getLogProb(calcNodes)
posteriorLogDensity0 <- DPOSTERIORCALL_ORIG
posteriorLogDensity1 <- DPOSTERIORCALL_NEW
posteriorVerification <- modelLogProb0 - posteriorLogDensity0 - modelLogProb1 + posteriorLogDensity1
if(abs(posteriorVerification) > 1e-8) {
nimPrint('conjugate posterior density appears to be wrong, off by ', posteriorVerification)
}
}, list(DPOSTERIORCALL_ORIG = eval(substitute(substitute(expr, list(VALUE=quote(origTargetValue))), list(expr=posteriorObject$dCallExpr))),
DPOSTERIORCALL_NEW = eval(substitute(substitute(expr, list(VALUE=quote(newTargetValue))), list(expr=posteriorObject$dCallExpr)))))
}
functionDef <- quote(function() {})
functionDef[[3]] <- functionBody$getCode()
functionDef[[4]] <- NULL ## removes the 'scrref' attribute
return(functionDef)
},
genGetPosteriorLogDensityFunction = function(dependentCounts, doDependentScreen = FALSE) {
functionBody <- codeBlockClass()
functionBody$addCode(origTargetValue <- model[[target]])
## adds code to generate the quantities prior_xxx, and contribution_xxx
addPosteriorQuantitiesGenerationCode(functionBody = functionBody, dependentCounts = dependentCounts, doDependentScreen = doDependentScreen)
## calculate and return the (log)density for the current value of target
functionBody$addCode(posteriorObject$prePosteriorCodeBlock, quote = FALSE)
if(link != 'identity') {
functionBody$addCode({
model[[target]] <<- origTargetValue
model$calculate(calcNodesDeterm)})
}
functionBody$addCode({
posteriorLogDensity <- DPOSTERIORCALL
returnType(double())
return(posteriorLogDensity)
}, list(DPOSTERIORCALL = eval(substitute(substitute(expr, list(VALUE=quote(origTargetValue))), list(expr=posteriorObject$dCallExpr)))))
functionDef <- quote(function() {})
functionDef[[3]] <- functionBody$getCode()
functionDef[[4]] <- NULL ## removes the 'scrref' attribute
return(functionDef)
},
addPosteriorQuantitiesGenerationCode = function(functionBody = functionBody, dependentCounts = dependentCounts, doDependentScreen = FALSE) {
## get current value of prior parameters which appear in the posterior expression
for(priorParam in posteriorObject$neededPriorParams) {
functionBody$addCode(PRIOR_PARAM_VAR <- model$getParam(target[1], PARAM_NAME),
list(PRIOR_PARAM_VAR = as.name(paste0('prior_', priorParam)),
PARAM_NAME = priorParam))
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
neededParams <- dependents[[distName]]$neededParamsForPosterior
depNodeValueNdim <- getDimension(distName)
forLoopBody <- codeBlockClass()
## get *value* of each dependent node
if(any(getDimension(distName, includeParams = TRUE) > 0)) {
forLoopBody$addCode(thisNodeSize <- DEP_NODESIZES[iDep],
list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes'))))
}
forLoopBody$addCode(DEP_VALUES_VAR_INDEXED <<- model$getParam(DEP_NODENAMES[iDep], 'value'),
list(DEP_VALUES_VAR_INDEXED = makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_values')), depNodeValueNdim, indexExpr = quote(iDep), secondSize = quote(thisNodeSize), thirdSize = quote(thisNodeSize)),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames'))))
for(param in neededParams) {
depNodeParamNdim <- getDimension(distName, param)
## get *parameter values* for each dependent node
forLoopBody$addCode(DEP_PARAM_VAR_INDEXED <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME),
list(DEP_PARAM_VAR_INDEXED = makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_', param)), depNodeParamNdim, indexExpr = quote(iDep), secondSize = quote(thisNodeSize), thirdSize = quote(thisNodeSize)),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames')),
PARAM_NAME = param))
}
functionBody$addCode(for(iDep in 1:N_DEP) FORLOOPBODY,
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
FORLOOPBODY = forLoopBody$getCode()))
}
targetNdim <- getDimension(prior)
allCurrentLinks <- sapply(names(dependentCounts), function(x) strsplit(x, '_')[[1]][[2]])
switch(as.character(targetNdim),
`0` = {
if(any(allCurrentLinks %in% c('additive', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- 0
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('additive', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode(
for(iDep in 1:N_DEP)
DEP_OFFSET_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME),
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
}
if(any(allCurrentLinks == 'multiplicativeScalar'))
stop("Found 'multiplicativeScalar' link for 0-dimensional case.")
if(any(allCurrentLinks %in% c('multiplicative', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- 1
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('multiplicative', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
inputList <- list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param,
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')))
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode(
for(iDep in 1:N_DEP)
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME) - DEP_OFFSET_VAR[iDep],
inputList)
} else {
functionBody$addCode(
for(iDep in 1:N_DEP)
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME),
inputList)
}
}
}
}
},
`1` = {
if(any(allCurrentLinks == 'multiplicativeScalar'))
stop("Found 'multiplicativeScalar' link for 1-dimensional case.")
if(any(allCurrentLinks %in% c('additive', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- rep(0, d)
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('additive', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
for(iDep in 1:N_DEP) {
thisNodeSize <- DEP_NODESIZES[iDep]
DEP_OFFSET_VAR[iDep, 1:thisNodeSize] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)
}
},
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName,'_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
}
if(any(allCurrentLinks %in% c('multiplicative', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode(unitVector <- rep(0, d))
forLoopBody <- codeBlockClass()
forLoopBody$addCode({
unitVector[sizeIndex] <- 1
model[[target]] <<- unitVector
unitVector[sizeIndex] <- 0
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink %in% c('multiplicative', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
inputList <- list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes')),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param,
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')))
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
forLoopBody$addCode(
for(iDep in 1:N_DEP) {
thisNodeSize <- DEP_NODESIZES[iDep]
DEP_COEFF_VAR[iDep, 1:thisNodeSize, sizeIndex] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME) - DEP_OFFSET_VAR[iDep, 1:thisNodeSize]
}, inputList)
} else {
forLoopBody$addCode(
for(iDep in 1:N_DEP) {
thisNodeSize <- DEP_NODESIZES[iDep]
DEP_COEFF_VAR[iDep, 1:thisNodeSize, sizeIndex] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)
}, inputList)
}
}
}
functionBody$addCode(for(sizeIndex in 1:d) FORLOOPBODY,
list(FORLOOPBODY = forLoopBody$getCode()))
}
},
`2` = {
if(!all(allCurrentLinks %in% c('identity', 'multiplicativeScalar')))
stop("Found non-multiplicative link for 2-d variable.")
if(any(allCurrentLinks == 'multiplicativeScalar') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
functionBody$addCode({
model[[target]] <<- diag(d)
model$calculate(calcNodesDeterm)
})
## Use _COEFF_VAR to store value; we need this for stoch indexing case
## where we determine that coeff = 0 because the potential dependency is not a dependency
## given current index values.
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink == 'multiplicativeScalar' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
for(iDep in 1:N_DEP) {
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)[1, 1] ## DEP_COEFF_VAR = (A+2B)-(A+B) = B
}
},
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
functionBody$addCode({
## Can't use zeros matrix as Cholesky fails; this is solely to determine if
## potential dependency is not a dependency of target (due to stochastic indexing).
model[[target]] <<- 2*diag(d)
model$calculate(calcNodesDeterm)
})
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
if(currentLink == 'multiplicativeScalar'|| (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
functionBody$addCode({
for(iDep in 1:N_DEP) {
DEP_COEFF_VAR[iDep] <<- model$getParam(DEP_NODENAMES[iDep], PARAM_NAME)[1, 1] - DEP_COEFF_VAR[iDep]
}
},
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
DEP_NODENAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
PARAM_NAME = dependents[[distName]]$param))
}
}
},
stop()
)
targetNdim <- getDimension(prior)
## contribution terms have been moved to be setup function outputs,
## but we still need to *zero these variables out* before adding contribution terms into them
for(contributionName in posteriorObject$neededContributionNames) {
contribNdim <- posteriorObject$neededContributionDims[[contributionName]]
functionBody$addCode(CONTRIB_NAME <<- CONTRIB_ZERO_OUT,
list(CONTRIB_NAME = as.name(contributionName),
CONTRIB_ZERO_OUT = switch(as.character(contribNdim),
`0` = 0, `1` = quote(rep(0, length = d)), `2` = quote(array(0, dim = c(d, d))), stop())))
}
for(iDepCount in seq_along(dependentCounts)) {
distLinkName <- names(dependentCounts)[iDepCount]
tmp <- strsplit(names(dependentCounts)[iDepCount], '_')[[1]]
distName <- tmp[[1]]
currentLink <- tmp[[2]]
targetCoeffNdim <- switch(as.character(targetNdim), `0`=0, `1`=2, `2`=2, stop())
if(targetCoeffNdim == 2 && link == 'multiplicativeScalar') ## There are no cases where we allow non-scalar 'coeff'.
targetCoeffNdim <- 0
if(!any(posteriorObject$neededContributionNames %in% dependents[[distName]]$contributionNames)) next
depParamsAvailable <- dependents[[distName]]$neededParamsForPosterior
## don't allow ragged dependencies for 2D conjugate case.
## no such cases exist, and it causes a runtime size check compiler warning.
## nonRaggedSizeExpr used to replace quote(thisNodeSize) below.
## August 2016
nonRaggedSizeExpr <- if(targetNdim < 2) quote(thisNodeSize) else quote(d)
subList <- lapply(depParamsAvailable, function(param)
makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_', param)), getDimension(distName, param), indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = nonRaggedSizeExpr))
names(subList) <- depParamsAvailable
subList$value <- makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_values')), getDimension(distName), indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = nonRaggedSizeExpr)
if(currentLink %in% c('additive', 'linear', 'stickbreaking'))
subList$offset <- makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_offset')), targetNdim, indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = nonRaggedSizeExpr)
if(currentLink %in% c('multiplicative', 'multiplicativeScalar', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
subList$coeff <- makeIndexedVariable(as.name(paste0('dep_', distLinkName, '_coeff')), targetCoeffNdim, indexExpr = quote(iDep), secondSize = nonRaggedSizeExpr, thirdSize = quote(d))
forLoopBody <- codeBlockClass()
if(any(getDimension(distName, includeParams = TRUE) > 0)) {
if(targetNdim == 1) ## 1D
forLoopBody$addCode(thisNodeSize <- DEP_NODESIZES[iDep],
list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes'))))
if(targetNdim == 2) ## 2D ## formerly this was 'else', but for 'dcat' we have targetNdim=0 while max(getDimension(distName, includeParams = TRUE)) is 1 so need explicit check for 2D
forLoopBody$addCode(if(DEP_NODESIZES[iDep] != d) print('runtime error with sizes of 2D conjugate sampler'),
list(DEP_NODESIZES = as.name(paste0('dep_', distLinkName, '_nodeSizes'))))
}
for(contributionName in posteriorObject$neededContributionNames) {
if(!(contributionName %in% dependents[[distName]]$contributionNames)) next
contributionExpr <- dependents[[distName]]$contributionExprs[[contributionName]]
if(distName == 'dmnorm' && prior == 'dmnorm') {
## need to deal with [,1] in contribution_mean
if(contributionName == 'contribution_mean') tmpExpr <- contributionExpr[[2]][[2]] else tmpExpr <- contributionExpr
if(getNimbleOption('allowDynamicIndexing') && doDependentScreen) {
## Will always have 'coeff' but may not need to use it
if(currentLink %in% c('identity', 'additive'))
tmpExpr[[6]] <- quote(0) else tmpExpr[[6]] <- quote(1)
} else {
if(currentLink %in% c('identity', 'additive')) {
tmpExpr[[2]] <- quote(prec) ## hack, since 'coeff' won't exist
tmpExpr[[6]] <- quote(0)
} else tmpExpr[[6]] <- quote(1)
}
tmpExpr[[4]] <- cc_stripExpr(tmpExpr[[4]], offset = currentLink %in% c('identity','multiplicative'), coeff = FALSE) # strip 'offset'
if(contributionName == 'contribution_mean') contributionExpr[[2]][[2]] <- tmpExpr else contributionExpr <- tmpExpr
} else contributionExpr <- cc_stripExpr(contributionExpr, offset = currentLink %in% c('identity','multiplicative','multiplicativeScalar'),
coeff = currentLink %in% c('identity','additive'))
contributionExpr <- eval(substitute(substitute(EXPR, subList), list(EXPR=contributionExpr)))
if(getNimbleOption('allowDynamicIndexing') && doDependentScreen) { ## FIXME: would be nice to only have one if() here when we loop through multiple parameters
if(targetCoeffNdim == 0)
forLoopBody$addCode(if(COEFF_EXPR != 0) CONTRIB_NAME <<- CONTRIB_NAME + CONTRIB_EXPR,
list(COEFF_EXPR = subList$coeff, CONTRIB_NAME = as.name(contributionName), CONTRIB_EXPR = contributionExpr))
else forLoopBody$addCode(if(min(COEFF_EXPR) != 0 | max(COEFF_EXPR) != 0) CONTRIB_NAME <<- CONTRIB_NAME + CONTRIB_EXPR,
list(COEFF_EXPR = subList$coeff, CONTRIB_NAME = as.name(contributionName), CONTRIB_EXPR = contributionExpr))
} else forLoopBody$addCode(CONTRIB_NAME <<- CONTRIB_NAME + CONTRIB_EXPR,
list(CONTRIB_NAME = as.name(contributionName), CONTRIB_EXPR = contributionExpr))
}
functionBody$addCode(for(iDep in 1:N_DEP) FORLOOPBODY,
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
FORLOOPBODY = forLoopBody$getCode()))
}
}
)
)
dependentClass <- setRefClass(
Class = 'dependentClass',
fields = list(
distribution = 'ANY', ## the name of the (dependent) sampling distribution, e.g. 'dnorm'
param = 'ANY', ## the name of the sampling distribution parameter in which target must appear
contributionExprs = 'ANY', ## a (named) list of expressions, giving the (additive) contribution to any parameters of the posterior. names correspond to variables in the posterior expressions
contributionNames = 'ANY', ## names of the contributions to the parameters of the posterior distribution. same as names(posteriorExprs)
neededParamsForPosterior = 'ANY', ## names of all parameters appearing in the posteriorExprs
link = 'ANY' ## optional specific link for a given conjugacy; currently only used for stickbreaking case
),
methods = list(
initialize = function(depInfoList, depDistName) {
contributionExprs <<- list()
distribution <<- depDistName
param <<- depInfoList$param
link <<- depInfoList$link ## will be NULL unless specific conjugacy overrides default link
initialize_contributionExprs(depInfoList)
initialize_neededParamsForPosterior()
},
initialize_contributionExprs = function(depInfoList) {
depInfoList['param'] <- NULL
depInfoList <- lapply(depInfoList, function(di) parse(text=di)[[1]])
contributionExprs <<- depInfoList
contributionExprs[['link']] <<- NULL ## so link not included in neededParamsForPosterior
contributionNames <<- names(depInfoList)
},
initialize_neededParamsForPosterior = function() {
posteriorVars <- unique(unlist(lapply(contributionExprs, all.vars)))
## Formerly, we always included the target 'param' in 'neededParamsForPosterior'.
## Not certain why, it wasn't necessary in many cases, e.g., dbeta-dbinom conjugacy. Excluding it now.
## -DT, August 2017
##neededParamsForPosterior <<- unique(c(param, posteriorVars[!(posteriorVars %in% c('value', 'coeff', 'offset'))]))
neededParamsForPosterior <<- setdiff(posteriorVars, c('value', 'coeff', 'offset'))
}
)
)
posteriorClass <- setRefClass(
Class = 'posteriorClass',
fields = list(
prePosteriorCodeBlock = 'ANY', ## a quoted {...} code block containing DSL code to execute before making posterior call, possibly empty
posteriorExpr = 'ANY', ## the full, parsed, posterior distribution expression, e.g. dnorm(mean = prior_mean + ..., sd = ...)
rDistribution = 'ANY', ## the *R* name of the posterior distribution, e.g. 'rnorm'
dDistribution = 'ANY', ## the *R* name of the posterior density distribution, e.g. 'dnorm'
argumentExprs = 'ANY', ## (named) list of expressions for each argument to the posterior distribution. names are the posterior distribution argument names
argumentNames = 'ANY', ## character vector of the argument names to the posterior distribution. same as: names(argumentExprs)
rCallExpr = 'ANY', ## the actual 'rnorm(1, ...)' call, which will be substituted into the conjugate sampler function
dCallExpr = 'ANY', ## the 'dnorm(value, ...)' call, which can be used to get values of the posterior density
neededPriorParams = 'ANY', ## the names of any prior parameters (e.g., 'mean') which appear in the posterior expression as 'prior_mean'
neededContributionNames = 'ANY', ## the names of contributions from dependent nodes, such as 'contribution_scale'
neededContributionDims = 'ANY' ## a named list of contribution dimensions (0, 1, or 2). List element names are, e.g., 'contribution_scale'
),
methods = list(
initialize = function(posteriorText, prior) {
parsedTotalPosterior <- parse(text = posteriorText)[[1]]
if(parsedTotalPosterior[[1]] != '{') parsedTotalPosterior <- substitute({POST}, list(POST = parsedTotalPosterior))
prePosteriorCodeBlock <<- parsedTotalPosterior[-length(parsedTotalPosterior)]
posteriorExpr <<- parsedTotalPosterior[[length(parsedTotalPosterior)]]
rDistribution <<- cc_makeRDistributionName(as.character(posteriorExpr[[1]]))
dDistribution <<- as.character(posteriorExpr[[1]])
argumentExprs <<- as.list(posteriorExpr)[-1]
argumentNames <<- names(argumentExprs)
rCallExpr <<- as.call(c(as.name(rDistribution), 1, argumentExprs))
dCallExpr <<- as.call(c(as.name(dDistribution), quote(VALUE), argumentExprs, log = 1))
posteriorVars <- all.vars(parsedTotalPosterior)
neededPriorParams <<- gsub('^prior_', '', posteriorVars[grepl('^prior_', posteriorVars)])
neededContributionNames <<- posteriorVars[grepl('^contribution_', posteriorVars)]
neededContributionDims <<- inferContributionTermDimensions(prior)
},
inferContributionTermDimensions = function(prior) {
distToLookup <- if(dDistribution %in% distributions$namesVector) dDistribution else if(prior %in% distributions$namesVector) prior else stop('cannot locate prior or posterior distribution in conjugacy processing')
targetNdim <- getDimension(distToLookup)
## if posterior distribution is univariate, assume all contributions are scalar
if(targetNdim == 0) {
theDims <- lapply(neededContributionNames, function(x) 0)
names(theDims) <- neededContributionNames
return(theDims)
}
## if posterior distribution is multivariate, attempt to infer contribution dimensionality from the *name* of each contribution term
theDims <- list()
typeNamesAvailable <- getParamNames(distToLookup)
for(contribName in neededContributionNames) {
contribNameBase <- gsub('^contribution_', '', contribName)
if(contribNameBase %in% typeNamesAvailable) {
## contribution base name matches a parameter name of the posterior
theDims[[contribName]] <- getDimension(distToLookup, contribNameBase)
} else {
## contribution base name doesn't match any parameter; can't easily infer the dimensionality
stop(message('The NIMBLE conjugacy system is attempting to infer the dimensionality of the contribution term: ',
contribName, '. However, since the posterior distribution is multivariate, and the contribution name doesn\'t match any parameter names of the posterior distribution, NIMBLE can\'t infer this one. This means the conjugacy system might need to be extended, to allow users to provide the dimensionality of contribution terms. Or perhaps something more clever. -DT August 2015'), call. = FALSE)
}
}
return(theDims)
}
)
)
##############################################################################################
##############################################################################################
## utility functions
##############################################################################################
##############################################################################################
cc_makeSamplerTypeName <- function(distName) return(paste0('conjugate_', distName)) ## 'dnorm' --> 'conjugate_dnorm'
cc_makeConjugateSamplerName <- function(samplerType) return(paste0('sampler_', samplerType)) ## 'conjugate_dnorm' --> 'sampler_conjugate_dnorm'
cc_makeRDistributionName <- function(distName) return(paste0('r', substring(distName, 2))) ## 'dnorm' --> 'rnorm'
## expands all deterministic nodes in expr, to create a single expression with only stochastic nodes
cc_expandDetermNodesInExpr <- function(model, expr, targetNode = NULL, skipExpansionsNode = NULL, prevExpr = NULL) {
## prevents an infinite-recursion case, where structureExpr(EXPR) is repeatedly processed, indefinitely (DT Aug 2022):
if(!is.null(prevExpr) && identical(expr, prevExpr)) return(expr)
if(is.numeric(expr)) return(expr) # return numeric
if(is.name(expr) || (is.call(expr) && (expr[[1]] == '[') && is.name(expr[[2]]))) { # expr is a name, or an indexed name
if(getNimbleOption('allowDynamicIndexing')) {
## this deals with having mu[k[1]] (which won't pass through expandNodeNames), replacing k[1] with the index from targetNode
if(!is.name(expr)) {
indexExprs <- expr[3:length(expr)]
numericOrVectorIndices <- sapply(indexExprs,
function(x) is.numeric(x) || (length(x) == 3 && x[[1]] == ':'))
if(!all(numericOrVectorIndices)) {
if(model$getVarNames(nodes = targetNode) == model$getVarNames(nodes = safeDeparse(expr, warn = TRUE))) {
## expr var is same as target var, so plug in target indexes for
## non-constant expr indexes to allow for possibility that
## dynamic index will be that of target (in some cases based on allowed
## values of dynamic index, this might actually not be possible, but
## the only disdvantage should be failing to detect conjugacy where it is
## present).
expr[which(!numericOrVectorIndices)+2] <- parse(text = targetNode)[[1]][3:length(expr)][!numericOrVectorIndices]
## now check that target is not actually used in index expressions
newExpr <- as.call(c(cc_structureExprName, sapply(indexExprs[!numericOrVectorIndices], function(x) x)))
for(i in seq_along(newExpr)[-1])
newExpr[[i]] <- cc_expandDetermNodesInExpr(model, newExpr[[i]], targetNode, skipExpansionsNode)
if(cc_nodeInExpr(targetNode, newExpr))
return(as.call(c(cc_structureExprName, expr, sapply(indexExprs[!numericOrVectorIndices], function(x) x)))) # put 'expr' back in though shouldn't be needed downstream
## otherwise continue with processing as in non-dynamic index case
} else {
## in this case the expr var is a different var, so shouldn't really matter what
## indexes get plugged in; plug in mins of indexes as a default
varInfo <- model$modelDef$varInfo[[safeDeparse(expr[[2]], warn = TRUE)]]
expr[which(!numericOrVectorIndices)+2] <- varInfo$mins[!numericOrVectorIndices]
## sapply business gets rid of () at end of index expression
newExpr <- as.call(c(cc_structureExprName, expr, sapply(indexExprs[!numericOrVectorIndices], function(x) x)))
for(i in seq_along(newExpr)[-1])
newExpr[[i]] <- cc_expandDetermNodesInExpr(model, newExpr[[i]], targetNode, skipExpansionsNode, prevExpr = expr)
return(newExpr)
}
} ## else continue with processing as in non-dynamic index case
}
}
exprText <- safeDeparse(expr, warn = TRUE)
expandedNodeNamesRaw <- model$expandNodeNames(exprText)
## skipExpansionsNode was added specifically for CAR model target nodes:
if(!is.null(skipExpansionsNode) && (exprText %in% model$expandNodeNames(skipExpansionsNode, returnScalarComponents=TRUE))) return(expr)
## if exprText is a node itself (and also part of a larger node), then we only want the expansion to be the exprText node:
if(exprText %in% expandedNodeNamesRaw) {
expandedNodeNames <- exprText
} else {
## if exprText is itself *not* a formal node name, but rather a subset of a larger formal node,
## then we don't want it to be expanded to include the subsubming full node
if(length(setdiff(model$expandNodeNames(expandedNodeNamesRaw, returnScalarComponents = TRUE), model$expandNodeNames(exprText, returnScalarComponents = TRUE)))) {
expandedNodeNames <- exprText
} else {
expandedNodeNames <- expandedNodeNamesRaw
}
}
if(length(expandedNodeNames) == 1 && (expandedNodeNames == exprText) && (exprText %in% expandedNodeNamesRaw)) {
## expr is a single node in the model
type <- model$getNodeType(exprText)
## this next block covers a **real corner case** (from a legacy BUGS model: biops) where dimensions are inferred
## (since not provided), thus causing a valid stochastic node to rather have dimensions larger than it acually is
## (being defined as part of a ragged declaration of multivariate stochastic nodes),
## thus the "type" of this expanded node includes both 'stoch' and also 'RHSonly', leading to a "node"
## with two types... which triggers the check below. So, if a "node" we're processing has multiple types, then
## we remove the extraneous 'RHSonly' types here:
if((length(type) > 1) && ('RHSonly' %in% type) && !all(type == 'RHSonly')) type <- setdiff(type, 'RHSonly')
if(length(type) > 1) {
## if exprText is a node itself (and also part of a larger node), then we only want the expansion to be the exprText node:
type <- type[which(exprText == expandedNodeNamesRaw)]
}
if(type == 'stoch') return(expr)
if(type == 'determ') {
if(!(exprText %in% model$getNodeNames(determOnly = TRUE))) return(expr) ## exprText is a single element of a multivariate deterministic node
newExpr <- model$getValueExpr(exprText)
return(cc_expandDetermNodesInExpr(model, newExpr, targetNode, skipExpansionsNode))
}
if(type == 'RHSonly') return(expr)
stop(paste0('Unexpected model structure for expression: ', exprText, '. Please report this to the NIMBLE development team.'))
}
newExpr <- cc_createStructureExpr(model, exprText)
if(is.call(newExpr) && newExpr[[1]] == 'structureExpr') { ## recurse, if there's a newly created structureExpr()
for(i in seq_along(newExpr)[-1])
newExpr[[i]] <- cc_expandDetermNodesInExpr(model, newExpr[[i]], targetNode, skipExpansionsNode, prevExpr = expr)
}
return(newExpr)
}
if(is.call(expr)) {
if(any(sapply(expr[-1], function(x) x=='')))
stop('Found missing indexing in ', safeDeparse(expr), ' that prevents conjugacy processing for this particular model structure.')
for(i in seq_along(expr)[-1])
expr[[i]] <- cc_expandDetermNodesInExpr(model, expr[[i]], targetNode, skipExpansionsNode)
return(expr)
}
stop(paste0('something went wrong processing: ', safeDeparse(expr)))
}
## special name used to represent vectors / arrays defined in terms of other stoch/determ nodes
cc_structureExprName <- quote(structureExpr)
## creates an expression of the form [cc_structureExprName](element11, element12, etc...) to represent vectors / arrays defined in terms of other stoch/determ nodes,
cc_createStructureExpr <- function(model, exprText) {
expandedNodeNamesVector <- model$expandNodeNames(exprText)
## remove expanded node names which are not fully represented in the original scalar components of exprText:
exprTextScalarComponents <- model$expandNodeNames(exprText, returnScalarComponents=TRUE)
expandedNodeNamesToKeepBool <- sapply(expandedNodeNamesVector, function(n) all(model$expandNodeNames(n, returnScalarComponents=TRUE) %in% exprTextScalarComponents))
expandedNodeNamesVector <- expandedNodeNamesVector[expandedNodeNamesToKeepBool]
## now, add in any scalar components of original exprText which are not appearing in the expanded node names:
scalarComponentsToAdd <- setdiff(exprTextScalarComponents, model$expandNodeNames(expandedNodeNamesVector, returnScalarComponents=TRUE))
nodesForStructureExpr <- c(expandedNodeNamesVector, scalarComponentsToAdd)
## if exprText is a non-node scalar component, then just return that parse(exprText); (rather than structureExpr(exprText))
if((length(nodesForStructureExpr) == 1) && identical(nodesForStructureExpr, scalarComponentsToAdd)) return(parse(text=exprText)[[1]])
## otherwise, create and return a structureExpr
structureExprExprList <- lapply(nodesForStructureExpr, function(x) parse(text=x)[[1]])
structureExpr <- c(cc_structureExprName, structureExprExprList)
structureExprCall <- as.call(structureExpr)
return(structureExprCall)
}
## verifies that 'link' is satisfied by the results of linearityCheck
cc_linkCheck <- function(linearityCheck, link) {
linearityLinks <- c('identity', 'additive', 'multiplicative', 'multiplicativeScalar', 'linear')
if(!(link %in% c(linearityLinks, 'stickbreaking')))
stop(paste0('unknown link: \'', link, '\''))
if(is.null(linearityCheck)) return(NULL)
offset <- linearityCheck$offset
scale <- linearityCheck$scale
if(link %in% linearityLinks && offset == 0 && scale == 1)
return('identity')
if(link %in% c('additive', 'linear') && scale == 1)
return('additive')
## For Wishart/inverse-Wishart, only handle scalar 'scale'.
if(link == 'multiplicativeScalar' && offset == 0 && cc_checkScalar(scale))
return('multiplicativeScalar')
if(link %in% c('multiplicative', 'linear') && offset == 0)
return('multiplicative')
if(link == 'linear')
return('linear')
return(NULL)
}
cc_checkScalar <- function(expr) {
if(length(expr) == 1) return(TRUE)
if(expr[[1]] == ":") return(FALSE)
## We don't know the output dims of all functions,
## as there can be user-defined functions that produce non-scalar output from scalar inputs,
## so only say it could be scalar (based on input args) in specific known cases we enumerate.
if(safeDeparse(expr[[1]], warn = TRUE) %in% c('(','[','*','+','-','/','exp','log','^','pow','pow_int','sqrt')) {
return(all(sapply(expr[2:length(expr)], cc_checkScalar)))
} else return(FALSE)
}
## Removes 'coeff' and/or 'offset' from expression for use in simplifying conjugacy contribution
## expressions when realized link is 'simpler' than possible links.
cc_stripExpr <- function(expr, offset = TRUE, coeff = TRUE) {
if(length(expr) == 1) {
if(coeff && expr == 'coeff') return(1)
if(offset && expr == 'offset') return(0)
return(expr)
}
if(offset && (expr[[1]] == '+' || expr[[1]] == '-')) {
if(expr[[2]] == 'offset' || expr[[2]] == 0)
return(cc_stripExpr(expr[[3]], offset, coeff))
if(expr[[3]] == 'offset' || expr[[3]] == 0)
return(cc_stripExpr(expr[[2]], offset, coeff))
}
if(coeff && expr[[1]] == '*') {
if(expr[[2]] == 'coeff' || expr[[2]] == 1)
return(cc_stripExpr(expr[[3]], offset, coeff))
if(expr[[3]] == 'coeff' || expr[[3]] == 1)
return(cc_stripExpr(expr[[2]], offset, coeff))
}
if(coeff && expr[[1]] == '/' && (expr[[3]] == 'coeff' || expr[[3]] == 1))
return(cc_stripExpr(expr[[2]], offset, coeff))
if(coeff && expr[[1]] == '/' && expr[[2]] == 'coeff') {
expr[[2]] <- 1
expr[[3]] <- cc_stripExpr(expr[[3]], offset, coeff)
return(expr)
}
if(coeff && (expr[[1]] == '^' || expr[[1]] == 'sqrt') && expr[[2]] == 'coeff') {
return(1)
}
if(is.call(expr))
for(i in seq_along(expr)[-1])
expr[[i]] <- cc_stripExpr(expr[[i]], offset, coeff)
return(expr)
}
## checks the parameter expressions in the stochastic distribution of depNode
## returns FALSE if we find 'targetNode' in ***more than one*** of these expressions
cc_otherParamsCheck <- function(model, depNode, targetNode, skipExpansionsNode = NULL, depNodeExprExpanded, depParamNodeName) {
paramsList <- as.list(model$getValueExpr(depNode)[-1]) # extracts the list of all parameters, for the distribution of depNode
timesFound <- 0 ## for success, we'll find targetNode in only *one* parameter expression
for(i in seq_along(paramsList)) {
if(!missing(depParamNodeName) && (names(paramsList)[i] == depParamNodeName)) {
expr <- depNodeExprExpanded
} else { expr <- cc_expandDetermNodesInExpr(model, paramsList[[i]], targetNode, skipExpansionsNode) }
if(cc_vectorizedComponentCheck(targetNode, expr)) return(FALSE)
if(cc_nodeInExpr(targetNode, expr)) { timesFound <- timesFound + 1 } ## we found 'targetNode'
}
if(timesFound == 0) stop('something went wrong; targetNode not found in any parameter expressions')
if(timesFound == 1) return(TRUE)
if(timesFound > 1) return(FALSE)
}
## determines whether node appears anywhere in expr
cc_nodeInExpr <- function(node, expr) { return(node %in% cc_getNodesInExpr(expr)) }
## determines which nodes apppear in an expression
## exporting as used in setup code of a nf called directly by user
#' @export
cc_getNodesInExpr <- function(expr) {
if(is.numeric(expr)) return(character(0)) ## expr is numeric
if(is.logical(expr)) return(character(0)) ## expr is logical
if(is.name(expr) || (is.call(expr) && (expr[[1]] == '[') && is.name(expr[[2]]))) return(safeDeparse(expr, warn = TRUE)) ## expr is a node name
if(is.call(expr)) return(unlist(lapply(expr[-1], cc_getNodesInExpr))) ## expr is some general call
stop(paste0('something went wrong processing: ', safeDeparse(expr)))
}
## if targetNode is vectorized: determines if any components of targetNode appear in expr
cc_vectorizedComponentCheck <- function(targetNode, expr) {
if(!is.vectorized(targetNode)) return(FALSE)
targetNodesExpanded <- nl_expandNodeIndex(targetNode)
if(any(unlist(lapply(targetNodesExpanded, function(node) cc_nodeInExpr(node, expr))))) return(TRUE) ## any components of *vectorized* targetNode appear in expr
return(FALSE)
}
##############################################################################################
##############################################################################################
## general, flexible, check for linear relationship
## returns list(offset, scale), or NULL
##############################################################################################
##############################################################################################
cc_replace01 <- function(expr) {
## replace {0,1} with {0,1}+1i so distinguished from offset,scale of 0,1 when match target
if(is.numeric(expr) && (expr %in% c(0, 1)))
return(expr + 1i)
if(is.call(expr))
for(i in 2:length(expr))
expr[[i]] <- cc_replace01(expr[[i]])
return(expr)
}
cc_checkLinearity <- function(expr, targetNode) {
## targetNode doesn't appear in expr
if(!cc_nodeInExpr(targetNode, expr)) {
if(is.call(expr) && expr[[1]] == '(') return(cc_checkLinearity(expr[[2]], targetNode))
## add +1i to tags 0s and 1s as not being from exact match to target
return(list(offset = cc_replace01(expr), scale = 0))
}
## Check if expr is exactly the targetNode.
## `deparse()` may return multiple lines, which will cause `FALSE` here.
## That should always be correct as multiple lines would only occur
## when `expr` is not a node name. It's possible one could have a node name
## with many indices and large integers as indices (e.g., x[13434, 1342352,...])
## but to get to a length that would get deparsed to multiple lines
## would entail a gigantic object.
if(identical(targetNode, deparse(expr)))
return(list(offset = 0, scale = 1))
if(!is.call(expr)) stop('cc_checkLinearity: expression is not a call object')
## process the expression contents of the parentheses
if(expr[[1]] == '(')
return(cc_checkLinearity(expr[[2]], targetNode))
if(expr[[1]] == '[')
return(cc_checkLinearity(expr[[2]], targetNode))
## Look for individual nodes in vectorized use or other strange cases.
if(expr[[1]] == 'structureExpr') {
## Can't have target appear multiple times.
if(sum(sapply(expr[2:length(expr)], function(x)
cc_nodeInExpr(targetNode, x))) != 1)
return(NULL)
wh <- which(sapply(expr[2:length(expr)], function(x)
cc_nodeInExpr(targetNode, x)))
checkLinearityStrucExpr <- cc_checkLinearity(expr[[wh+1]], targetNode)
if(is.null(checkLinearityStrucExpr)) return(NULL)
if(is.indexed(targetNode) && grepl(':', targetNode)) {
## if this is a structureExpression, and also
## if targetNode is multivariate (is indexed and contains ":"),
## and is also part of a larger structureExpr,
## then we do not handle conjugacy for those cases
return(NULL)
}
return(list(offset = cc_combineExprsAddition(expr, checkLinearityStrucExpr$offset), # was expr?,
scale = checkLinearityStrucExpr$scale))
}
## we'll just have to skip over asRow() and asCol(), so they don't mess up the linearity check
if(expr[[1]] == 'asRow' || expr[[1]] == 'asCol') {
return(cc_checkLinearity(expr[[2]], targetNode))
}
## minus sign: change to a plus sign, and invert the sign of the RHS
if(expr[[1]] == '-') {
if(length(expr) == 3) {
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
return(list(offset = cc_combineExprsSubtraction(checkLinearityLHS$offset, checkLinearityRHS$offset),
scale = cc_combineExprsSubtraction(checkLinearityLHS$scale, checkLinearityRHS$scale)))
}
if(length(expr) == 2) {
checkLin <- cc_checkLinearity(expr[[2]], targetNode)
if(is.null(checkLin)) return(NULL)
return(list(offset = cc_negateExpr(checkLin$offset),
scale = cc_negateExpr(checkLin$scale)))
}
stop('cc_checkLinearity: problem with negation expression')
}
if(expr[[1]] == '+') {
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
return(list(offset = cc_combineExprsAddition(checkLinearityLHS$offset, checkLinearityRHS$offset),
scale = cc_combineExprsAddition(checkLinearityLHS$scale, checkLinearityRHS$scale)))
}
if(expr[[1]] == '*' || expr[[1]] == '%*%' || expr[[1]] == 'inprod' || expr[[1]] == 'sum') {
## X[,] %*% beta[] where beta[i] are nodes is not standard matrix multiplication, so there is offset and scale
isMatrixMult <- ifelse(expr[[1]] == '%*%' &&
length(expr[[2]]) > 1 && expr[[2]][[1]] != 'structureExpr' &&
length(expr[[3]]) > 1 && expr[[3]][[1]] != 'structureExpr',
TRUE, FALSE)
if(expr[[1]] == 'sum')
if(length(expr[[2]]) == 3 && expr[[2]][[1]] == '*') {
tmpExpr <- quote(inprod(a, b))
tmpExpr[[2]] <- expr[[2]][[2]]
tmpExpr[[3]] <- expr[[2]][[3]]
expr <- tmpExpr
} else {
return(NULL) # cases such as sum(p[1:5]); there may be some unusual conjugacy cases here that we don't detect
}
if(length(expr) != 3) return(NULL) # avoid error at next step for unexpected cases
if(cc_nodeInExpr(targetNode, expr[[2]]) && cc_nodeInExpr(targetNode, expr[[3]])) return(NULL)
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
if((checkLinearityLHS$scale != 0) && (checkLinearityRHS$scale != 0)) stop('cc_checkLinearity: incompatible scales in * operation')
if((checkLinearityLHS$scale == 0) && (checkLinearityRHS$scale == 0)) {
return(list(offset = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$offset, isMatrixMult),
scale = 0)) }
if(checkLinearityLHS$scale != 0) {
return(list(offset = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$offset, isMatrixMult),
scale = cc_combineExprsMultiplication(checkLinearityLHS$scale, checkLinearityRHS$offset, isMatrixMult))) }
if(checkLinearityRHS$scale != 0) {
return(list(offset = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$offset, isMatrixMult),
scale = cc_combineExprsMultiplication(checkLinearityLHS$offset, checkLinearityRHS$scale , isMatrixMult))) }
stop('cc_checkLinearity: something went wrong')
}
if(expr[[1]] == '/') {
if(cc_nodeInExpr(targetNode, expr[[3]])) return(NULL)
checkLinearityLHS <- cc_checkLinearity(expr[[2]], targetNode)
checkLinearityRHS <- cc_checkLinearity(expr[[3]], targetNode)
if(is.null(checkLinearityLHS) || is.null(checkLinearityRHS)) return(NULL)
if(checkLinearityLHS$scale == 0) stop('left hand side scale == 0')
if(checkLinearityRHS$scale != 0) stop('right hand side scale is non-zero')
return(list(offset = cc_combineExprsDivision(checkLinearityLHS$offset, checkLinearityRHS$offset),
scale = cc_combineExprsDivision(checkLinearityLHS$scale, checkLinearityRHS$offset)))
}
## returns not conjugate (NULL) if we don't recognize the call (expr[[1]])
return(NULL)
}
cc_negateExpr <- function(expr) {
if(expr == 0) return(0)
if(is.numeric(expr)) return(-1*expr)
return(substitute(-EXPR, list(EXPR = expr)))
}
cc_combineExprsAddition <- function(expr1, expr2) {
if((expr1 == 0) && (expr2 == 0)) return(0)
if((expr1 != 0) && (expr2 == 0)) return(expr1)
if((expr1 == 0) && (expr2 != 0)) return(expr2)
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 + expr2)
return(substitute(EXPR1 + EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_combineExprsSubtraction <- function(expr1, expr2) {
if((expr1 == 0) && (expr2 == 0)) return(0)
if((expr1 != 0) && (expr2 == 0)) return(expr1)
if((expr1 == 0) && (expr2 != 0)) return(cc_negateExpr(expr2))
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 - expr2)
return(substitute(EXPR1 - EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_combineExprsMultiplication <- function(expr1, expr2, isMatrixMult) {
if((expr1 == 0) || (expr2 == 0)) return(0)
if((expr1 == 1) && (expr2 == 1)) return(1)
if((expr1 != 1) && (expr2 == 1)) return(expr1)
if((expr1 == 1) && (expr2 != 1)) return(expr2)
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 * expr2)
if(!isMatrixMult) return(substitute(EXPR1 * EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
if( isMatrixMult) return(substitute(EXPR1 %*% EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_combineExprsDivision <- function(expr1, expr2) {
if( (expr2 == 0)) stop('error, division by 0')
if( (expr2 == 1)) return(expr1)
if((expr1 == 0) ) return(0)
if(is.numeric(expr1) && is.numeric(expr2)) return(expr1 / expr2)
return(substitute(EXPR1 / EXPR2, list(EXPR1=expr1, EXPR2=expr2)))
}
cc_checkStickbreaking <- function(expr, targetNode) {
if(!is.call(expr) || expr[[1]] != 'stick_breaking' || !cc_nodeInExpr(targetNode, expr))
return(NULL)
expr <- expr[[2]]
if(!is.call(expr) || expr[[1]] != 'structureExpr')
return(NULL)
## offset <- which(targetNode == cc_getNodesInExpr(expr)) # Issue 1409.
## Need that target is the the entirety of one element of the structureExpr.
offset <- which(sapply(expr[2:length(expr)], function(x) x == targetNode))
if(length(offset) != 1) return(NULL) else return(list(offset = offset, scale = 1))
}
compareConjugacyLists <- function(C1, C2) {
if(identical(C1, C2)) return(TRUE)
if(!identical(names(C1), names(C2))) {cat('Names do not match\n'); return(FALSE)}
for(i in seq_along(C1)) {
if(!identical(C1[[i]]$type, C2[[i]]$type)) cat(paste0('type mismatch for i =',i))
if(!identical(C1[[i]]$target, C2[[i]]$target)) cat(paste0('target mismatch for i =',i))
if(!identical(C1[[i]]$target, C2[[i]]$target)) cat(paste0('target mismatch for i =',i))
if(!identical(names(C1[[i]]$control), names(C2[[i]]$control))) cat(paste0('control names mismatch for i =',i,'. Skipping node comparison'))
else {
for(j in seq_along(C1[[i]]$control)) {
if(!identical(sort(C1[[i]]$control[[j]]), sort(C2[[i]]$control[[j]]))) cat(paste0('target mismatch for i =',i, 'j =', j))
}
}
}
}
createDynamicConjugateSamplerName <- function(prior, dependentCounts, dynamicallyIndexed = FALSE) {
##depString <- paste0(dependentCounts, names(dependentCounts), collapse='_') ## including the numbers of dependents
depString <- paste0(sort(names(dependentCounts)), collapse='_') ## without the numbers of each type of dependent node
paste0('sampler_conjugate_', prior, '_', depString, ifelse(dynamicallyIndexed, '_dynamicDeps', ''))
}
makeDeclareSizeField <- function(firstSize, secondSize, thirdSize, nDim) {
eval(substitute(switch(as.character(nDim),
`0` = quote(FIRSTSIZE),
`1` = quote(c(FIRSTSIZE, SECONDSIZE)),
`2` = quote(c(FIRSTSIZE, SECONDSIZE, THIRDSIZE)),
stop()),
list(FIRSTSIZE = firstSize,
SECONDSIZE = secondSize,
THIRDSIZE = thirdSize)))
}
makeIndexedVariable <- function(varName, nDim, indexExpr, secondSize, thirdSize) {
eval(substitute(switch(as.character(nDim),
`0` = quote(VARNAME[INDEXEXPR]),
`1` = quote(VARNAME[INDEXEXPR, 1:SECONDSIZE]),
`2` = quote(VARNAME[INDEXEXPR, 1:SECONDSIZE, 1:THIRDSIZE]),
stop()),
list(VARNAME = varName,
INDEXEXPR = indexExpr,
SECONDSIZE = secondSize,
THIRDSIZE = thirdSize)))
}
##############################################################################################
##############################################################################################
## create object: conjugacyRelationshipsObject
## also, generate all conjugate sampler nimbleFunctions
## and a function to rebuild conjugate sampler functions
##############################################################################################
##############################################################################################
## this is still *necessary* (and exported):
conjugacyRelationshipsObject <- conjugacyRelationshipsClass(conjugacyRelationshipsInputList)
## here after is for handling of dynamic conjugate sampler function
dynamicConjugateSamplerDefinitionsEnv <- new.env()
dynamicConjugateSamplerFunctionsEnv <- new.env()
dynamicConjugateSamplerExists <- function(name) {
return(name %in% ls(dynamicConjugateSamplerDefinitionsEnv))
}
dynamicConjugateSamplerAdd <- function(name, def) {
dynamicConjugateSamplerDefinitionsEnv[[name]] <- def
dynamicConjugateSamplerFunctionsEnv[[name]] <- eval(def)
}
dynamicConjugateSamplerGet <- function(name) {
return(dynamicConjugateSamplerFunctionsEnv[[name]])
}
dynamicConjugateSamplerWrite <- function(file = 'TEMP_dynamicConjugateSamplerDefinitions.R') {
## environment to create functions in is throw-away, since they've all already been created in dynamicConjugateSamplerFunctionsEnv
createNamedObjectsFromList(as.list(dynamicConjugateSamplerDefinitionsEnv), writeToFile = file, envir = new.env())
}
Try the nimble 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.