Robustify handling of numerical issues in samplers (#1512)

* Make minor edit to REL_INST.

* Error trap adaptation interval < 2 for various block samplers.

* Robustify handling of NaN/Inf/-Inf in samplers.
Add checkLogProb() that converts NaN to -Inf and warns of Inf values.
Use checkLogProb() in various samplers (primarily slice and RW style).

* Fix indentation.

* Fix syntax error.

* updates

* updated indentation

* Add missing checkLogProb in sampler_RW.

* Fix typo in name.

* Refine use of checkLogProb to use directly on log prob values,
to better handle possible cases like "3 - Inf".
Fix tests in light of use of checkLogProb.

* Clean up testing for use of checkLogProb.
Recreate gold files for trunc and dynamicIndexing testing for mcmc
tests that get different samples because of changes in number of random numbers generated.

* Inline checkLogProb.

---------

Co-authored-by: Daniel Turek <danielturek@gmail.com>

packages/nimble/R/MCMC_samplers.R

@@ -124,16 +124,19 @@ sampler_binary <- nimbleFunction(
         if(!model$isBinary(target))     stop('can only use binary sampler on discrete 0/1 (binary) nodes')
     },
     run = function() {
-        currentLogProb <- model$getLogProb(calcNodes)
+        currentLogProb <- checkLogProb(model$getLogProb(calcNodes))
         model[[target]] <<- 1 - model[[target]]
-        otherLogProbPrior <- model$calculate(target)
+        otherLogProbPrior <- checkLogProb(model$calculate(target))
         if(otherLogProbPrior == -Inf) {
             otherLogProb <- otherLogProbPrior
         } else {
-            otherLogProb <- otherLogProbPrior + model$calculate(calcNodesNoSelf)
+            otherLogProb <- otherLogProbPrior + checkLogProb(model$calculate(calcNodesNoSelf))
         }
-        acceptanceProb <- 1/(exp(currentLogProb - otherLogProb) + 1)
-        jump <- (!is.nan(acceptanceProb)) & (runif(1,0,1) < acceptanceProb)
+        if(currentLogProb == -Inf & otherLogProb == -Inf)
+            stop("in binary sampler, all log probability density values are negative infinity and sampling cannot proceed")
+        logProbDiff <- currentLogProb - otherLogProb
+        acceptanceProb <- 1/(exp(logProbDiff) + 1)
+        jump <- (!is.nan(acceptanceProb)) & (runif(1,0,1) < acceptanceProb)  # `is.nan` probably not needed with use of `checkLogProb`.
         if(jump) {
             ##model$calculate(calcNodesPPomitted)
             nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE)
@@ -184,26 +187,23 @@ sampler_categorical <- nimbleFunction(
     },
     run = function() {
         currentValue <- model[[target]]
-        logProbs[currentValue] <<- model$getLogProb(calcNodes)
+        logProbs[currentValue] <<- checkLogProb(model$getLogProb(calcNodes))
         for(i in 1:k) {
             if(i != currentValue) {
                 model[[target]] <<- i
-                logProbPrior <- model$calculate(target)
+                logProbPrior <- checkLogProb(model$calculate(target))
                 if(logProbPrior == -Inf) {
                     logProbs[i] <<- -Inf
                 } else {
-                    if(is.nan(logProbPrior)) {
-                        logProbs[i] <<- -Inf
-                    } else {
-                        logProbs[i] <<- logProbPrior + model$calculate(calcNodesNoSelf)
-                        if(is.nan(logProbs[i])) logProbs[i] <<- -Inf
-                    }
+                    logProbs[i] <<- logProbPrior + checkLogProb(model$calculate(calcNodesNoSelf))
                 }
             }
         }
         maxLP <- max(logProbs)
-        if(maxLP == Inf | is.nan(maxLP))   cat("Warning: categorical sampler for '", target, "' encountered an invalid model density, and sampling results are likely invalid.\n")
+        if(maxLP ==  -Inf) stop("in categorical sampler, all log probability density values are negative infinity and sampling cannot proceed")
+        infLogProbs <- logProbs == Inf
         logProbs <<- logProbs - maxLP
+        logProbs[infLogProbs] <<- 0   ## Prevent NaN inputs into `rcat`.
         probs <<- exp(logProbs)
         newValue <- rcat(1, probs)   ## rcat normalizes the probabilities internally
         if(!is.na(newValue) & newValue != currentValue) {
@@ -363,12 +363,12 @@ sampler_RW <- nimbleFunction(
             }
         }
         model[[target]] <<- propValue
-        logMHR <- model$calculateDiff(target)
+        logMHR <- checkLogProb(model$calculateDiff(target))
         if(logMHR == -Inf) {
             jump <- FALSE
             nimCopy(from = mvSaved, to = model, row = 1, nodes = target, logProb = TRUE)
         } else {
-            logMHR <- logMHR + model$calculateDiff(calcNodesNoSelf) + propLogScale
+            logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + propLogScale
             jump <- decide(logMHR)
             if(jump) {
                 ##model$calculate(calcNodesPPomitted)
@@ -520,14 +520,14 @@ sampler_RW_noncentered <- nimbleFunction(
         }
         model[[target]] <<- propValue
 
-        logMHR <- model$calculateDiff(target)
+        logMHR <- checkLogProb(model$calculateDiff(target))
         if(logMHR == -Inf) {
             jump <- FALSE
             nimCopy(from = mvSaved, to = model, row = 1, nodes = target, logProb = TRUE)
         } else {
             ## Shift effects and add log-determinant of Jacobian of transformation. 
             logMHR <- logMHR + updateNoncentered(propValue, currentValue)
-            logMHR <- logMHR + model$calculateDiff(calcNodesNoSelf) + propLogScale
+            logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + propLogScale
             jump <- decide(logMHR)
             if(jump) {
                 ##model$calculate(calcNodesPPomitted)
@@ -689,12 +689,12 @@ sampler_RW_block <- nimbleFunction(
         for(i in 1:tries) {
             propValueVector <- generateProposalVector()
             values(model, targetAsScalar) <<- propValueVector
-            lpD <- model$calculateDiff(calcNodesProposalStage)
+            lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage))
             if(lpD == -Inf) {
                 jump <- FALSE
                 nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE)
             } else {
-                lpD <- lpD + model$calculateDiff(calcNodesDepStage)
+                lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage))
                 jump <- decide(lpD)
                 if(jump) {
                     ##model$calculate(calcNodesPPomitted)
@@ -941,9 +941,9 @@ sampler_slice <- nimbleFunction(
         setAndCalculateTarget = function(value = double()) {
             if(discrete)     value <- floor(value)
             model[[target]] <<- value
-            lp <- model$calculate(target)
+            lp <- checkLogProb(model$calculate(target))
             if(lp == -Inf) return(-Inf) 
-            lp <- lp + model$calculate(calcNodesNoSelf)
+            lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf))
             returnType(double())
             return(lp)
         },
@@ -1082,10 +1082,10 @@ sampler_slice_noncentered <- nimbleFunction(
         setAndCalculateTarget = function(value = double()) {
             if(discrete)     value <- floor(value)
             model[[target]] <<- value
-            lp <- model$calculate(target)
+            lp <- checkLogProb(model$calculate(target))
             if(lp == -Inf) return(-Inf)
             lp <- lp + updateNoncentered(value)
-            lp <- lp + model$calculate(calcNodesNoSelf)
+            lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf))
             returnType(double())
             return(lp)
         },
@@ -1215,7 +1215,7 @@ sampler_ess <- nimbleFunction(
         theta_min <- theta - 2*Pi
         theta_max <- theta
         values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d]
-        lp <- model$calculate(calcNodesNoSelf)
+        lp <- checkLogProb(model$calculate(calcNodesNoSelf))
         numContractions <- 0
         while((is.nan(lp) | lp < u) & theta_max - theta_min > eps & numContractions < maxContractions) {   # must be is.nan()
             ## The checks for theta_max - theta_min small and max number of contractions are
@@ -1224,7 +1224,7 @@ sampler_ess <- nimbleFunction(
             if(theta < 0)   theta_min <- theta   else   theta_max <- theta
             theta <- runif(1, theta_min, theta_max)
             values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d]
-            lp <- model$calculate(calcNodesNoSelf)
+            lp <- checkLogProb(model$calculate(calcNodesNoSelf))
             numContractions <- numContractions + 1
         }
         if(theta_max - theta_min <= eps | numContractions == maxContractions) {
@@ -1370,9 +1370,9 @@ sampler_AF_slice <- nimbleFunction(
                 for(i in 1:d)
                     if(discrete[i] == 1)   targetValues[i] <- floor(targetValues[i])            
             values(model, target) <<- targetValues
-            lp <- model$calculate(calcNodesProposalStage)
+            lp <- checkLogProb(model$calculate(calcNodesProposalStage))
             if(lp == -Inf) return(lp)
-            lp <- lp + model$calculate(calcNodesDepStage)
+            lp <- lp + checkLogProb(model$calculate(calcNodesDepStage))
             returnType(double())
             return(lp)
         },
@@ -1499,15 +1499,9 @@ sampler_crossLevel <- nimbleFunction(
         propLP0 <- 0
         for(iSF in seq_along(lowConjugateGetLogDensityFunctions))  { propLP0 <- propLP0 + lowConjugateGetLogDensityFunctions[[iSF]]$run() }
         propValueVector <- topRWblockSamplerFunction$generateProposalVector()
-        topLP <- my_setAndCalculateTop$run(propValueVector)
-        if(is.na(topLP)) {
-            logMHR <- -Inf
-            jump <- decide(logMHR)
-            if(jump) {
-                nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
-            } else {
-                nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE)
-            }
+        topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector))
+        if(topLP == -Inf) {
+            nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE)
         }
         else {
             for(iSF in seq_along(lowConjugateSamplerFunctions))
@@ -1516,7 +1510,7 @@ sampler_crossLevel <- nimbleFunction(
             propLP1 <- 0
             for(iSF in seq_along(lowConjugateGetLogDensityFunctions))
                 propLP1 <- propLP1 + lowConjugateGetLogDensityFunctions[[iSF]]$run()
-            logMHR <- modelLP1 - modelLP0 - propLP1 + propLP0
+            logMHR <- checkLogProb(modelLP1) - checkLogProb(modelLP0) - checkLogProb(propLP1) + checkLogProb(propLP0)
             jump <- decide(logMHR)
             if(jump) {
                 nimCopy(from = model,   to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
@@ -1728,7 +1722,7 @@ sampler_RW_dirichlet <- nimbleFunction(
                 thetaVecProp <- thetaVec
                 thetaVecProp[i] <- propValue
                 values(model, target) <<- thetaVecProp / sum(thetaVecProp)
-                logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + model$calculateDiff(calcNodesNoSelf)
+                logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf))
                 jump <- decide(logMHR)
             } else jump <- FALSE
             if(adaptive & jump)   timesAcceptedVec[i] <<- timesAcceptedVec[i] + 1
@@ -1848,7 +1842,7 @@ sampler_RW_wishart <- nimbleFunction(
         ## matrix multiply to get proposal value (matrix)
         model[[target]] <<- t(propValue_chol) %*% propValue_chol
         ## decide and jump
-        logMHR <- model$calculateDiff(calcNodes)
+        logMHR <- checkLogProb(model$calculateDiff(calcNodes))
         deltaDiag <- thetaVec_prop[1:d]-thetaVec[1:d]
         for(i in 1:d)   logMHR <- logMHR + (d+2-i)*deltaDiag[i]  ## took me quite a while to derive this
         jump <- decide(logMHR)
@@ -1971,7 +1965,8 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction(
                     propValue[jprime] <<- z[jprime, i] * sqrt(partialSumsProp[jprime])
                 }
                 model[[target]][j:i, i] <<- propValue[j:i] 
-                logMHR <- calculateDiff(model, calcNodesNoSelf) + calculateDiff(model, target)
+                logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf)) +
+                            checkLogProb(calculateDiff(model, target))
                 ## Adjust MHR to account for non-symmetric proposal by adjusting prior on U to transformed scale (i.e., y).
                 ## cosh component is for dz/dy and other component is for du/dz  where 'u' is the corr matrix.
                 ## This follows Stan reference manual Section 10.12 (for version 2.27).
@@ -2143,12 +2138,12 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
             ## Adjust for log determinant term from initial values
             logMHR <- logMHR - logDetJac
 
-            lpD <- calculateDiff(model, calcNodesProposalStage)
+            lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage))
             if(lpD == -Inf) {
                 nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE)
                 jump <- FALSE
             } else {
-                logMHR <- logMHR + lpD + calculateDiff(model, calcNodesDepStage)
+                logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage))
                 jump <- decide(logMHR)
                 if(jump) {
                     nimCopy(from = model,   to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
@@ -2423,7 +2418,7 @@ CAR_scalar_RW <- nimbleFunction(
         propValue <- rnorm(1, mean = model[[targetScalar]], sd = scale)
         model[[targetScalar]] <<- propValue
         lp1 <- dcarList[[1]]$run() + model$calculate(depNodes)
-        logMHR <- lp1 - lp0
+        logMHR <- checkLogProb(lp1) - checkLogProb(lp0)
         jump <- decide(logMHR)
         if(jump) {
             model$calculate(targetScalar)

packages/nimble/R/MCMC_utils.R

@@ -11,9 +11,9 @@
 #' @export
 decide <- function(logMetropolisRatio) {
     if(is.na(logMetropolisRatio)) return(FALSE)
-  if(logMetropolisRatio > 0) return(TRUE)
-  if(runif(1,0,1) < exp(logMetropolisRatio)) return(TRUE)
-  return(FALSE)
+    if(logMetropolisRatio > 0) return(TRUE)
+    if(runif(1,0,1) < exp(logMetropolisRatio)) return(TRUE)
+    return(FALSE)
 }
 
 #NOTE: DETAILS(WAS BLANK) REMOVED
@@ -55,7 +55,8 @@ decideAndJump <- nimbleFunction(
         copyNodesDeterm <- ccList$copyNodesDeterm; copyNodesStoch <- ccList$copyNodesStoch  # not used: calcNodes, calcNodesNoSelf
     },
     run = function(modelLP1 = double(), modelLP0 = double(), propLP1 = double(), propLP0 = double()) {
-        logMHR <- modelLP1 - modelLP0 - propLP1 + propLP0
+        ## Check each one individually to catch case like `3 - Inf`.
+        logMHR <- checkLogProb(modelLP1) - checkLogProb(modelLP0) - checkLogProb(propLP1) + checkLogProb(propLP0)
         jump <- decide(logMHR)
         if(jump) {
             nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE)
@@ -71,8 +72,26 @@ decideAndJump <- nimbleFunction(
     }
 )
 
+checkLogProb <- function(logProb) {
+   if(is.na(logProb))
+       return(-Inf)
+   if(logProb == Inf)
+         print("MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.")
+   return(logProb)
+}
 
-
+## checkLogProb <- nimbleFunction(
+##     name = "checkLogProb",
+##     run = function(logProb = double()) {
+##         if(is.na(logProb))
+##             return(-Inf)
+##         if(logProb == Inf)
+##             print("MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.")
+##         return(logProb)
+##         returnType(double())
+
+##     }
+## )
 
 
 #' Creates a nimbleFunction for setting the value of a scalar model node,

packages/nimble/R/RCfunction_core.R

@@ -75,6 +75,8 @@ fxnsNotAllowedInAD <- c(
                        'invwish_chol', 'car_normal','car_proper','multi','dirch') ),
   'getLogProb',
   'decide',
+  'checkLogProb',
+  'checkLogProbWarn',
   'rankSample',
   'any_na', 
   'any_nan',

packages/nimble/R/genCpp_sizeProcessing.R

@@ -139,6 +139,8 @@ sizeCalls <- c(
                    'nimArr_rmulti',
                    'nimArr_rdirch'), 'sizeRmultivarFirstArg'),
     makeCallList(c('decide',
+                   'checkLogProb',
+                   'checkLogProbWarn',
                    'size',
                    'getsize',
                    'getNodeFunctionIndexedInfo',
@@ -155,6 +157,7 @@ sizeCalls <- c(
     ADbreak = 'sizeADbreak')
 
 scalarOutputTypes <- list(decide = 'logical',
+                          checkLogProb = 'double',
                           size = 'integer',
                           nimAnyNA = 'logical',
                           nimAnyNaN = 'logical',

packages/nimble/R/genCpp_toEigenize.R

@@ -3,7 +3,8 @@
 
 toEigenizeNoCalls <- c('dim',
                        'run.time',
-                       'nimOptimDefaultControl')
+                       'nimOptimDefaultControl',
+                       'checkLogProbWarn')
 
 toEigenizeYesCalls <- c(paste0('nimDiagonal', c('D','I','B')),
                         'diagonal',
@@ -23,6 +24,7 @@ toEigenizeYesCalls <- c(paste0('nimDiagonal', c('D','I','B')),
 
 toEigenizeMaybeCalls <- c('map',
                           c('decide',
+                            'checkLogProb',
                             'size',
                             'getsize',
                             'getNodeFunctionIndexedInfo',

packages/nimble/inst/CppCode/Utils.cpp

@@ -21,6 +21,7 @@
 
 #include<limits>
 #include "nimble/Utils.h"
+#include "nimble/RcppNimbleUtils.h"
 
 // A utility function that will return floor(x) unless x is within numerical imprecision of an integer, in which case it will return round(x)
 int floorOrEquivalent(double x) {
@@ -62,6 +63,11 @@ bool decide(double lMHr) { // simple function accept or reject based on log Metr
   return(false);
 }
 
+void checkLogProbWarn() {
+  _nimble_global_output<<"MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.\n";
+  nimble_print_to_R(_nimble_global_output);
+}
+
 void nimStop(string msg) {NIMERROR("%s", msg.c_str());}
 void nimStop() {NIMERROR("Error. Exiting from compiled execution.");}
 

packages/nimble/inst/include/nimble/Utils.h

@@ -162,6 +162,11 @@ class nimbleTimerClass_ {
 bool decide(double lMHr);
 //void allocate(vector< vector <double> > *vv, int sampleSize, int variableSize);
 
+void checkLogProbWarn();
+
+static inline double checkLogProb(double logProb) { if(ISNAN(logProb))  return(-std::numeric_limits<double>::infinity()); if(logProb == std::numeric_limits<double>::infinity()) checkLogProbWarn(); return(logProb);}
+
+
 void nimStop(string msg);
 void nimStop();