avoid sampling parameters of empty clusters in BNP

packages/nimble/R/BNP_samplers.R

@@ -438,6 +438,10 @@ CRP_nonconjugate <- nimbleFunction(
   name = "CRP_nonconjugate",
   contains = CRP_helper,
   setup = function(model, marginalizedNodes, dataNodes, p, nTilde) {
+      savedIdx <- 1
+      saved <- nimNumeric(2) # treated as scalar if length 1
+      saved2 <- saved
+      saved3 <- saved
   },
   methods = list(
     storeParams = function() {},  ## nothing needed for non-conjugate
@@ -446,14 +450,26 @@ CRP_nonconjugate <- nimbleFunction(
       return(model$getLogProb(dataNodes[i]))
     },
     sample = function(i = integer(), j = integer() ) {
-      ## sample from prior
-      if( p == 1 ) {
-        model$simulate(marginalizedNodes[j])
-      } else {
-        for(l in 1:p) {  ## marginalized nodes should be in correct order based on findClusterNodes.
-          model$simulate(marginalizedNodes[(l-1)*nTilde + j])
+      if(j == 0) {  ## reset to stored values (for case of new cluster not opened)
+            values(model, marginalizedNodes[savedIdx]) <<- saved
+            if(p > 1) {
+                values(model, marginalizedNodes[nTilde + savedIdx]) <<- saved2
+                if(p > 2)
+                    values(model, marginalizedNodes[2*nTilde + savedIdx]) <<- saved3
+            }
+        } else {
+            savedIdx <<- j
+            saved <<- values(model, marginalizedNodes[j])
+            model$simulate(marginalizedNodes[j])
+            if(p > 1) {
+                saved2 <<- values(model, marginalizedNodes[nTilde + j])
+                if(p > 2)
+                    saved3 <<- values(model, marginalizedNodes[2*nTilde + j])
+                for(l in 2:p) {  ## marginalized nodes should be in correct order based on findClusterNodes.
+                    model$simulate(marginalizedNodes[(l-1)*nTilde + j])
+                }
+            }
         }
-      }
     }
   )
 )
@@ -984,6 +1000,8 @@ sampler_CRP <- nimbleFunction(
         ## p and nTilde only needed for non-conjugate currently.
         ## Note that the elements of tildeNodes will be in order such that the first element corresponds to the cluster
         ## obtained when xi[i] = 1, the second when xi[i] = 2, etc.
+      if(sampler == 'CRP_nonconjugate' && p > 3)
+          stop("sampler_CRP: CRP_nonconjugate sampler not yet set up to handle clustering of more than three variables.")  ## This is because of how we put old values back into model when proposing a new cluster that is not accepted.
       marginalizedNodes <- unlist(clusterVarInfo$clusterNodes)
       helperFunctions[[1]] <- eval(as.name(sampler))(model, marginalizedNodes, dataNodes, p, min_nTilde)
       calcNodes <- model$getDependencies(c(target, marginalizedNodes))
@@ -1160,6 +1178,8 @@ sampler_CRP <- nimbleFunction(
         }
         xiCounts[model[[target]][i]] <- 1
       } else { # an existing label is sampled
+        if(sampler == 'CRP_nonconjugate')   # reset to previous marginalized node value
+            helperFunctions[[1]]$sample(i, 0)
         if( xiCounts[xi[i]] == 0 ) { # xi_i is a singleton, a component was deleted
           k <- k - 1
           xiUniques <- reorderXiUniques
@@ -1175,7 +1195,7 @@ sampler_CRP <- nimbleFunction(
       }
     }
 
-      ## We have updated cluster variables but not all logProb values are up-to-date.
+    ## We have updated cluster variables but not all logProb values are up-to-date.
     model$calculate(calcNodes)
     copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
   },
@@ -1560,3 +1580,19 @@ checkNormalInvGammaConjugacy <- function(model, clusterVarInfo) {
     }
     return(conjugate)
 }
+
+sampler_CRP_cluster_wrapper <- nimbleFunction(
+    name = "CRP_cluster_wrapper", 
+    contains = sampler_BASE,
+    setup = function(model, mvSaved, target, control) {
+        regular_sampler <- nimbleFunctionList(sampler_BASE)
+        regular_sampler[[1]] <- control$wrapped_conf$buildSampler(model, mvSaved)
+        dcrpNode <- control$dcrpNode
+        clusterID <- control$clusterID
+    },
+    run = function() {
+        if(any(model[[dcrpNode]] == clusterID)) regular_sampler[[1]]$run()
+    },
+    methods = list(
+        reset = function() {regular_sampler[[1]]$reset()}
+    ))

packages/nimble/R/MCMC_configuration.R

@@ -197,6 +197,7 @@ print: A logical argument specifying whether to print the ordered list of defaul
                 isEndNode <- model$isEndNode(nodes)
                 if(useConjugacy) conjugacyResultsAll <- model$checkConjugacy(nodes)
 
+                clusterNodeInfo <- NULL; dcrpNode <- NULL; numCRPnodes <- 0
                 for(i in seq_along(nodes)) {
                     node <- nodes[i]
                     discrete <- model$isDiscrete(node)
@@ -222,7 +223,13 @@ print: A logical argument specifying whether to print the ordered list of defaul
                         if(nodeDist == 'dinvwish')     { addSampler(target = node, type = 'RW_wishart');         next }
                         if(nodeDist == 'dcar_normal')  { addSampler(target = node, type = 'CAR_normal');         next }
                         if(nodeDist == 'dcar_proper')  { addSampler(target = node, type = 'CAR_proper');         next }
-                        if(nodeDist == 'dCRP')         { addSampler(target = node, type = 'CRP', control = list(useConjugacy = useConjugacy));                next }
+                        if(nodeDist == 'dCRP')         {
+                            addSampler(target = node, type = 'CRP', control = list(useConjugacy = useConjugacy))
+                            numCRPnodes <- numCRPnodes + 1
+                            clusterNodeInfo[[numCRPnodes]] <- findClusterNodes(model, node)
+                            dcrpNode[numCRPnodes] <- node
+                            next
+                        }
                         if(multivariateNodesAsScalars) {
                             for(scalarNode in nodeScalarComponents) {
                                 if(onlySlice) addSampler(target = scalarNode, type = 'slice')
@@ -262,12 +269,32 @@ print: A logical argument specifying whether to print the ordered list of defaul
                     ## default: 'RW' sampler
                     addSampler(target = node, type = 'RW');     next
                 }
+
+                ## For CRP-based models, wrap samplers for cluster parameters so not sampled if cluster is unoccupied.
+                if(!is.null(clusterNodeInfo)) {
+                    for(k in seq_along(clusterNodeInfo)) {
+                        for(clusterNodes in clusterNodeInfo[[k]]$clusterNodes) {
+                            samplers <- getSamplers(clusterNodes)
+                            removeSamplers(clusterNodes)
+                            for(i in seq_along(samplers)) {
+                                node <- samplers[[i]]$target
+                                addSampler(target = node, type = 'CRP_cluster_wrapper',
+                                           control = list(wrapped_type = samplers[[i]]$name, wrapped_conf = samplers[[i]],
+                                                          dcrpNode = dcrpNode[[k]], clusterID = i))
+                                ## Note for more general clustering: will probably change to
+                                ## 'clusterID=clusterNodeInfo[[k]]$clusterIDs[[??]][i]'
+                                ## which means we probably need to change to for(clusterNodesIdx in seq_along(clusterNodeInfo[[k]]$clusterNodes))
+                            }
+                        }
+                    }
+                }
+
             }
 
             if(print)   printSamplers()
         },
 
-        addConjugateSampler = function(conjugacyResult, dynamicallyIndexed = FALSE, print = FALSE) {
+        addConjugateSampler = function(conjugacyResult, dynamicallyIndexed = FALSE, dcrpNode = NULL, clusterID = NULL, print = FALSE) {
             ## update May 2016: old (non-dynamic) system is no longer supported -DT
             ##if(!getNimbleOption('useDynamicConjugacy')) {
             ##    addSampler(target = conjugacyResult$target, type = conjugacyResult$type, control = conjugacyResult$control)

packages/nimble/R/MCMC_utils.R

@@ -238,8 +238,11 @@ mcmc_generateControlListArgument <- function(control, controlDefaults) {
 
 mcmc_listContentsToStr <- function(ls, displayControlDefaults=FALSE, displayNonScalars=FALSE, displayConjugateDependencies=FALSE) {
     ##if(any(unlist(lapply(ls, is.function)))) warning('probably provided wrong type of function argument')
-    if(!displayConjugateDependencies)
+    if(!displayConjugateDependencies) {
         if(grepl('^conjugate_d', names(ls)[1])) ls <- ls[1]    ## for conjugate samplers, remove all 'dep_dnorm', etc, control elements (don't print them!)
+        if(grepl('^CRP_cluster_wrapper', names(ls)[1]) && 'wrapped_type' %in% names(ls) &&
+           grepl('conjugate_d', ls$wrapped_type)[1]) ls <- ls[!names(ls) %in% c('wrapped_conf')]
+    }
     ls <- lapply(ls, function(el) if(is.nf(el) || is.function(el)) 'function' else el)   ## functions -> 'function'
     ls2 <- list()
     ## to make displayControlDefaults argument work again, would need to code process