Catching numerical instability

[Cole-Monnahan-NOAA]

How should I deal with R functions that evalulate to Inf or NaN when passing them to stan_sample? For instance I ported the gp_pois_reg example from the posteriordb set (i.e., posteriordb::stan_code('gp_pois_regr') to RTMB and occasionally it crashes out presumably due to numeric issues with the Cholesky decomposition for extreme values during warmup. Reprex below. Is there a way around this? My understand is that it should just be ignored by Stan and counted as a divergence and a new trajectory tried at the next iteration.

remotes::install_github("https://github.com/kaskr/RTMB", subdir="RTMB")
library(RTMB)
dat <- list(x = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10),
     k = c(40, 37, 29, 12, 4, 3, 9, 19, 77, 82, 33))
pars <- list(logrho=1, logalpha=1,  f_tilde=rep(0,11))
func <- function(pars){
  getAll(pars,dat)
  alpha <- exp(logalpha)
  rho <- exp(logrho)
  cov <- matrix(NA, 11,11)
  for(i in 1:11){
    for(j in 1:11){
      # https://mc-stan.org/docs/functions-reference/matrix_operations.html#exponentiated-quadratic-kernel
      cov[i,j] <- alpha^2*exp(-(x[i]-x[j])^2/(2*rho^2))
    }
  }
  cov <- cov+diag(1e-10,11)
  L_cov <- t(chol(cov)) # upper tri chol
  f <- as.numeric(L_cov %*% f_tilde) # log-predicted lambda
  lp <-
    # priors
    dgamma(rho, shape=25,scale=1/4, log=TRUE)+
    dnorm(alpha, 0,2,log=TRUE)+
    sum(dnorm(f_tilde,0,1,log=TRUE)) + # hyperdistribution
    logrho + logalpha + # Jacobians
    sum(dpois(k, exp(f), log=TRUE)) # likelihood
  REPORT(f)
  return(-lp) # note the negative log posterior
}
obj0 <- MakeADFun(func=func, parameters=pars, random=NULL, silent=TRUE)

# will escape warmup with seed=1, fail with seed=2
test <- stan_sample(fn=function(x) -obj0$fn(x), 
                    grad_fun = function(x) -obj0$gr(x), num_chains=1, seed = 2,
                    num_samples=1000, num_warmup=1000, par_inits=obj0$par)