diff --git a/.github/workflows/check.yaml b/.github/workflows/check.yaml
index c1a5be9..9e1fe79 100644
--- a/.github/workflows/check.yaml
+++ b/.github/workflows/check.yaml
@@ -43,12 +43,12 @@ jobs:
     steps:
       - uses: actions/checkout@v2
 
-      - uses: r-lib/actions/setup-r@master
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
           http-user-agent: ${{ matrix.config.http-user-agent }}
 
-      - uses: r-lib/actions/setup-pandoc@master
+      - uses: r-lib/actions/setup-pandoc@v2
 
       - name: Query dependencies
         run: |
diff --git a/.github/workflows/check_windows.yaml b/.github/workflows/check_windows.yaml
index 1b11ca2..134e265 100644
--- a/.github/workflows/check_windows.yaml
+++ b/.github/workflows/check_windows.yaml
@@ -39,12 +39,12 @@ jobs:
     steps:
       - uses: actions/checkout@v2
 
-      - uses: r-lib/actions/setup-r@master
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
           http-user-agent: ${{ matrix.config.http-user-agent }}
 
-      - uses: r-lib/actions/setup-pandoc@master
+      - uses: r-lib/actions/setup-pandoc@v2
 
       - name: Query dependencies
         run: |
diff --git a/DESCRIPTION b/DESCRIPTION
index 88b8c1c..4e43bd9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -28,7 +28,7 @@ Collate:
   dOcc.R
   dNmixture.R
   zzz.R
-RoxygenNote: 7.2.0
+RoxygenNote: 7.3.1
 Suggests: 
     rmarkdown,
     knitr,
diff --git a/R/dBetaBinom.R b/R/dBetaBinom.R
index a65f371..62e4db7 100644
--- a/R/dBetaBinom.R
+++ b/R/dBetaBinom.R
@@ -66,7 +66,9 @@ nimBetaFun <- nimbleFunction(
     if (log) return(lgamma(a) + lgamma(b) - lgamma(a + b))
     else return(exp(lgamma(a) + lgamma(b) - lgamma(a + b)))
     returnType(double(0))
-  })
+  },
+  buildDerivs = list(run = list())
+)
 
 #' @rdname dBetaBinom
 #' @export
@@ -87,7 +89,8 @@ dBetaBinom <- nimbleFunction(
     if (log) return(logprob)
     return(exp(logprob))
     returnType(double(0))
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 #' @rdname dBetaBinom
@@ -107,7 +110,8 @@ dBetaBinom_One <- nimbleFunction(
     if (log) return(logprob)
     return(exp(logprob))
     returnType(double(0))
-  }
+  },
+  buildDerivs = list(run = list())
 )
 
 
diff --git a/R/dNmixture.R b/R/dNmixture.R
index 5905610..74df2f1 100644
--- a/R/dNmixture.R
+++ b/R/dNmixture.R
@@ -220,13 +220,13 @@ dNmixture_v <- nimbleFunction(
     run = function(x = double(1),
                    lambda = double(),
                    prob = double(1),
-                   Nmin = integer(0, default = -1),
-                   Nmax = integer(0, default = -1),
+                   Nmin = double(),
+                   Nmax = double(),
                    len = integer(),
                    log = integer(0, default = 0)) {
   if (length(x) != len) stop("in dNmixture_v, len must equal length(x).")
   if (len != length(prob)) stop("in dNmixture_v, len must equal length(prob).")
-
+  
   # Lambda cannot be negative
   if (lambda < 0) {
     if (log) return(-Inf)
@@ -235,33 +235,36 @@ dNmixture_v <- nimbleFunction(
 
   ## For each x, the conditional distribution of (N - x | x) is pois(lambda * (1-p))
   ## We determine the lowest N and highest N at extreme quantiles and sum over those.
-  if (Nmin == -1) {
-    Nmin <- min(x + qpois(0.00001, lambda * (1 - prob)))
-  }
-  if (Nmax == -1) {
-    Nmax <- max(x + qpois(0.99999, lambda * (1 - prob)))
-  }
-  Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
-
+  # if (Nmin == -1) {
+  #   Nmin <- min(x + qpois(0.00001, lambda * (1 - prob)))
+  # }
+  # if (Nmax == -1) {
+  #   Nmax <- max(x + qpois(0.99999, lambda * (1 - prob)))
+  # }
+  Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
+      
   logProb <- -Inf
-
+      
   if (Nmax > Nmin) {
     numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-    prods <- rep(0, numN)
-    for (i in (Nmin + 1):Nmax) {
+    ##prods <- rep(0, numN) ## This cannot compile with AD
+    prods <- numeric(length = ADbreak(numN), value = 0)
+    for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
       prods[i - Nmin] <- prod(i/(i - x)) / i
     }
-
     ff <- log(lambda) + sum(log(1-prob)) + log(prods)
     log_fac <- nimNmixPois_logFac(numN, ff)
+        
     logProb <- dpois(Nmin, lambda, log = TRUE) + sum(dbinom(x, size = Nmin, prob = prob, log = TRUE)) + log_fac
   }
   if (log) return(logProb)
   else return(exp(logProb))
   returnType(double())
-}
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
+
 NULL
 #' @rdname dNmixture
 #' @export
@@ -269,8 +272,8 @@ dNmixture_s <- nimbleFunction(
     run = function(x = double(1),
                    lambda = double(),
                    prob = double(),
-                   Nmin = double(0, default = -1),
-                   Nmax = double(0, default = -1),
+                   Nmin = double(),
+                   Nmax = double(),
                    len = double(),
                    log = integer(0, default = 0)) {
   if (length(x) != len) stop("in dNmixture_s, len must equal length(x).")
@@ -283,20 +286,21 @@ dNmixture_s <- nimbleFunction(
 
   ## For each x, the conditional distribution of (N - x | x) is pois(lambda * (1-p))
   ## We determine the lowest N and highest N at extreme quantiles and sum over those.
-  if (Nmin == -1) {
-    Nmin <- min(x + qpois(0.00001, lambda * (1 - prob)))
-  }
-  if (Nmax == -1) {
-    Nmax <- max(x + qpois(0.99999, lambda * (1 - prob)))
-  }
-  Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+  # if (Nmin == -1) {
+  #   Nmin <- min(x + qpois(0.00001, lambda * (1 - prob)))
+  # }
+  # if (Nmax == -1) {
+  #   Nmax <- max(x + qpois(0.99999, lambda * (1 - prob)))
+  # }
+  Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
   logProb <- -Inf
 
   if (Nmax > Nmin) {
     numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-    prods <- rep(0, numN)
-    for (i in (Nmin + 1):Nmax) {
+    # prods <- rep(0, numN)
+    prods <- numeric(length = ADbreak(numN), value = 0)
+    for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
       prods[i - Nmin] <- prod(i/(i - x)) / i
     }
 
@@ -307,7 +311,8 @@ dNmixture_s <- nimbleFunction(
   if (log) return(logProb)
   else return(exp(logProb))
   returnType(double())
-}
+ },
+ buildDerivs = list(run = list(ignore = c("i")))
 )
 
 NULL
@@ -317,8 +322,8 @@ rNmixture_v <- nimbleFunction(
   run = function(n = double(),
                  lambda = double(),
                  prob = double(1),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double()) {
     if (n != 1) stop("rNmixture_v only works for n = 1")
     if (length(prob) != len) stop("In rNmixture_v, len must equal length(prob).")
@@ -339,8 +344,8 @@ rNmixture_s <- nimbleFunction(
   run = function(n = double(),
                  lambda = double(),
                  prob = double(),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double()) {
     if (n != 1) stop("rNmixture_v only works for n = 1")
     trueN <- rpois(1, lambda)
@@ -363,8 +368,8 @@ dNmixture_BNB_v <- nimbleFunction(
                  lambda = double(),
                  theta = double(),
                  prob = double(1),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double(),
                  log = integer(0, default = 0)) {
     if (length(x) != len) stop("in dNmixture_BNB_v, len must equal length(x).")
@@ -394,21 +399,22 @@ dNmixture_BNB_v <- nimbleFunction(
     lambda_cond <- omega / (theta_cond * (1 - omega))
     r_cond <- 1 / theta_cond
     pNB_cond <- 1 / (1 + theta_cond * lambda_cond)
-    if (Nmin == -1) {
-      Nmin <- min(x + qnbinom(0.00001, size = r_cond, prob = pNB_cond))
-    }
-    if (Nmax == -1) {
-      Nmax <- max(x + qnbinom(0.99999, size = r_cond, prob = pNB_cond))
-    }
+    # if (Nmin == -1) {
+    #   Nmin <- min(x + qnbinom(0.00001, size = r_cond, prob = pNB_cond))
+    # }
+    # if (Nmax == -1) {
+    #   Nmax <- max(x + qnbinom(0.99999, size = r_cond, prob = pNB_cond))
+    # }
 
-    Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+    Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
-      for (i in (Nmin + 1):Nmax) {
+      # prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         prods[i - Nmin] <- (i + r - 1) * prod(i/(i - x)) / i
       }
 
@@ -421,7 +427,8 @@ dNmixture_BNB_v <- nimbleFunction(
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BNB_s #####
@@ -433,8 +440,8 @@ dNmixture_BNB_s <- nimbleFunction(
                  lambda = double(),
                  theta = double(),
                  prob = double(),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double(),
                  log = integer(0, default = 0)) {
     if (length(x) != len) stop("in dNmixture_BNB_s, len must equal length(x).")
@@ -463,20 +470,21 @@ dNmixture_BNB_s <- nimbleFunction(
     lambda_cond <- omega / (theta_cond * (1 - omega))
     r_cond <- 1 / theta_cond
     pNB_cond <- 1 / (1 + theta_cond * lambda_cond)
-    if (Nmin == -1) {
-      Nmin <- min(x + qnbinom(0.00001, size = r_cond, prob = pNB_cond))
-    }
-    if (Nmax == -1) {
-      Nmax <- max(x + qnbinom(0.99999, size = r_cond, prob = pNB_cond))
-    }
-    Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+    # if (Nmin == -1) {
+    #   Nmin <- min(x + qnbinom(0.00001, size = r_cond, prob = pNB_cond))
+    # }
+    # if (Nmax == -1) {
+    #   Nmax <- max(x + qnbinom(0.99999, size = r_cond, prob = pNB_cond))
+    # }
+    Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
-      for (i in (Nmin + 1):Nmax) {
+      # prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         prods[i - Nmin] <- (i + r - 1) * prod(i/(i - x)) / i
       }
 
@@ -489,7 +497,8 @@ dNmixture_BNB_s <- nimbleFunction(
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BNB_oneObs #####
@@ -502,8 +511,8 @@ dNmixture_BNB_oneObs <- nimbleFunction(
                  lambda = double(),
                  theta = double(),
                  prob = double(),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double(),
                  log = integer(0, default = 0)) {
     if (theta <= 0) {
@@ -530,33 +539,36 @@ dNmixture_BNB_oneObs <- nimbleFunction(
     lambda_cond <- omega / (theta_cond * (1 - omega))
     r_cond <- 1 / theta_cond
     pNB_cond <- 1 / (1 + theta_cond * lambda_cond)
-    if (Nmin == -1) {
-      Nmin <- x + qnbinom(0.00001, size = r_cond, prob = pNB_cond)
-    }
-    if (Nmax == -1) {
-      Nmax <- x + qnbinom(0.99999, size = r_cond, prob = pNB_cond)
-    }
-    Nmin <- max(c(x, Nmin)) ## set Nmin to at least the largest x
+    # if (Nmin == -1) {
+    #   Nmin <- x + qnbinom(0.00001, size = r_cond, prob = pNB_cond)
+    # }
+    # if (Nmax == -1) {
+    #   Nmax <- x + qnbinom(0.99999, size = r_cond, prob = pNB_cond)
+    # }
+    Nmin <- ADbreak(max(x, Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
-      for (i in (Nmin + 1):Nmax) {
+      # prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         prods[i - Nmin] <- (i + r - 1) / (i - x)
       }
 
       ff <- log(1 - pNB) + log(1-prob) + log(prods)
       log_fac <- nimNmixPois_logFac(numN, ff)
+      xx <- nimNumeric(length = 2, value = x)
       logProb <- dnbinom(Nmin, size = r, prob = pNB, log = TRUE) +
-        dbinom(x, size = Nmin, prob = prob, log = TRUE) +
+        dbinom(x, size = Nmin, prob = prob, log = TRUE) + 
         log_fac
     }
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BBP_v #####
@@ -594,20 +606,20 @@ dNmixture_BBP_v <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+    Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
+      #prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
 
-      for (i in (Nmin + 1):Nmax) {
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         # prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) / i
         prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) * (lambda / i)
       }
-
-
+      
       ff <- log(prods)
       log_fac <- nimNmixPois_logFac(numN, ff)
       logProb <- dpois(Nmin, lambda, log = TRUE) +
@@ -617,7 +629,8 @@ dNmixture_BBP_v <- nimbleFunction(
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BBP_s #####
@@ -631,7 +644,7 @@ dNmixture_BBP_s <- nimbleFunction(
                  s = double(),
                  Nmin = double(0, default = -1),
                  Nmax = double(0, default = -1),
-                 len = double(),
+                 len = integer(),
                  log = integer(0, default = 0)) {
     if (length(x) != len) stop("in dNmixture_BBP_s, len must equal length(x).")
 
@@ -654,30 +667,36 @@ dNmixture_BBP_s <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+    Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
+      #prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
 
-      for (i in (Nmin + 1):Nmax) {
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         # prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) / i
         prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) * (lambda / i)
       }
 
-
       ff <- log(prods)
       log_fac <- nimNmixPois_logFac(numN, ff)
+      
+      alphaRep <- numeric(length = len, value = alpha)
+      betaRep <- numeric(length = len, value = beta)
+      
       logProb <- dpois(Nmin, lambda, log = TRUE) +
-        dBetaBinom(x, Nmin, rep(alpha, len), rep(beta, len), log = TRUE) +
+        #dBetaBinom(x, Nmin, rep(alpha, len), rep(beta, len), log = TRUE) +
+        dBetaBinom(x, Nmin, alphaRep, betaRep, log = TRUE) +
         log_fac
     }
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BBP_oneObs #####
@@ -712,22 +731,23 @@ dNmixture_BBP_oneObs <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Nmin <- x
+    Nmin <- ADbreak(x)
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
+      #prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
 
-      for (i in (Nmin + 1):Nmax) {
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         # prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) / i
         prods[i - Nmin] <- i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1)) * (lambda / i)
       }
 
-
       ff <- log(prods)
       log_fac <- nimNmixPois_logFac(numN, ff)
+      
       logProb <- dpois(Nmin, lambda, log = TRUE) +
         dBetaBinom_One(x, Nmin, alpha, beta, log = TRUE) +
         log_fac
@@ -735,7 +755,8 @@ dNmixture_BBP_oneObs <- nimbleFunction(
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 
@@ -782,15 +803,16 @@ dNmixture_BBNB_v <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+    Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
+      #prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
 
-      for (i in (Nmin + 1):Nmax) {
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         # prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) / i
         prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) *
           ((1 - pNB) * (i + r - 1) / i)
@@ -806,7 +828,8 @@ dNmixture_BBNB_v <- nimbleFunction(
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BBNB_s #####
@@ -821,7 +844,7 @@ dNmixture_BBNB_s <- nimbleFunction(
                  s = double(),
                  Nmin = double(0, default = -1),
                  Nmax = double(0, default = -1),
-                 len = double(),
+                 len = integer(),
                  log = integer(0, default = 0)) {
     if (length(x) != len) stop("in dNmixture_BBNB_s, len must equal length(x).")
 
@@ -851,31 +874,37 @@ dNmixture_BBNB_s <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
+    Nmin <- ADbreak(max(max(x), Nmin)) ## set Nmin to at least the largest x
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
+      #prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
 
-      for (i in (Nmin + 1):Nmax) {
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         # prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) / i
         prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) *
           ((1 - pNB) * (i + r - 1) / i)
       }
 
       ff <- log(prods)
-
       log_fac <- nimNmixPois_logFac(numN, ff)
+      
+      alphaRep <- nimNumeric(length = len, value = alpha)
+      betaRep <- nimNumeric(length = len, value = beta)
+      
       logProb <- dnbinom(Nmin, size = r, prob = pNB, log = TRUE) +
-        dBetaBinom(x, Nmin, rep(alpha, len), rep(beta, len), log = TRUE) +
+        #dBetaBinom(x, Nmin, rep(alpha, len), rep(beta, len), log = TRUE) +
+        dBetaBinom(x, Nmin, alphaRep, betaRep, log = TRUE) +
         log_fac
     }
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### dNmixture_BBNB_oneObs #####
@@ -919,15 +948,16 @@ dNmixture_BBNB_oneObs <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    if (Nmin < x) Nmin <- x
+    if (Nmin < x) Nmin <- ADbreak(x)
 
     logProb <- -Inf
 
     if (Nmax > Nmin) {
       numN <- Nmax - Nmin + 1 - 1  ## remember: +1 for the count, but -1 because the summation should run from N = maxN to N = minN + 1
-      prods <- rep(0, numN)
+      #prods <- rep(0, numN)
+      prods <- numeric(length = ADbreak(numN), value = 0)
 
-      for (i in (Nmin + 1):Nmax) {
+      for (i in ADbreak(Nmin + 1):ADbreak(Nmax)) {
         # prods[i - Nmin] <- prod(i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1))) / i
         prods[i - Nmin] <- i * (i - 1 + beta - x) / ((i - x) * (alpha + beta + i - 1)) *
           ((1 - pNB) * (i + r - 1) / i)
@@ -935,7 +965,7 @@ dNmixture_BBNB_oneObs <- nimbleFunction(
       ff <- log(prods)
 
       log_fac <- nimNmixPois_logFac(numN, ff)
-
+      
       logProb <- dnbinom(Nmin, size = r, prob = pNB, log = TRUE) +
         dBetaBinom_One(x, Nmin, alpha, beta, log = TRUE) +
         log_fac
@@ -943,7 +973,8 @@ dNmixture_BBNB_oneObs <- nimbleFunction(
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
-  }
+  },
+  buildDerivs = list(run = list(ignore = c("i")))
 )
 
 ##### rNmixture extensions #####
@@ -956,8 +987,8 @@ rNmixture_BNB_v <- nimbleFunction(
                  lambda = double(),
                  theta = double(),
                  prob = double(1),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double()) {
 
     if (n != 1) stop("rNmixture* only works for n = 1")
@@ -984,8 +1015,8 @@ rNmixture_BNB_s <- nimbleFunction(
                  lambda = double(),
                  theta = double(),
                  prob = double(),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double()) {
 
     if (n != 1) stop("rNmixture* only works for n = 1")
@@ -1012,8 +1043,8 @@ rNmixture_BNB_oneObs <- nimbleFunction(
                  lambda = double(),
                  theta = double(),
                  prob = double(),
-                 Nmin = double(0, default = -1),
-                 Nmax = double(0, default = -1),
+                 Nmin = double(),
+                 Nmax = double(),
                  len = double()) {
     if (n != 1) stop("rNmixture* only works for n = 1")
 
@@ -1062,7 +1093,7 @@ rNmixture_BBP_s <- nimbleFunction(
                  s = double(),
                  Nmin = double(0, default = -1),
                  Nmax = double(0, default = -1),
-                 len = double()) {
+                 len = integer()) {
     if (n != 1) stop("rNmixture* only works for n = 1")
 
     alpha <- prob * s
@@ -1136,7 +1167,7 @@ rNmixture_BBNB_s <- nimbleFunction(
                  s = double(),
                  Nmin = double(0, default = -1),
                  Nmax = double(0, default = -1),
-                 len = double()) {
+                 len = integer()) {
     if (n != 1) stop("rNmixture* only works for n = 1")
 
     alpha <- prob * s
diff --git a/R/utils.R b/R/utils.R
index f85a64a..531a140 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -33,20 +33,22 @@ nimNmixPois_logFac <- nimbleFunction(
       max_index <- 1 # not sure this is relevant. it's defensive.
       sum_ff_g1 <- ff[1]
     }
-
-    terms <- numeric(numN + 1)
+    
+    ## terms <- numeric(numN + 1) ## This cannot compile with AD
+    terms_len <- ADbreak(numN) + 1
+    terms <- numeric(length = ADbreak(terms_len), value = 0)
     terms[max_index + 1] <- 1
-
+    
     sumff <- sum_ff_g1 ## should be the same as sum(ff[1:max_index])
-
+    
     for (i in 1:max_index) {
       # terms[i] <- 1 / exp(sum(ff[i:max_index]))
       terms[i] <- 1 / exp(sumff)
       sumff <- sumff - ff[i]
     }
-
+    
     sumff <- 0
-    for (i in (max_index + 1):numN) {
+    for (i in (max_index + 1):ADbreak(numN)) {
       # terms[i + 1] <- exp(sum(ff[(max_index + 1):i]))
       sumff <- sumff + ff[i]
       terms[i + 1] <- exp(sumff)
@@ -54,4 +56,6 @@ nimNmixPois_logFac <- nimbleFunction(
     log_fac <- sum_ff_g1 + log(sum(terms)) # Final factor is the largest term * (all factors / largest term)    }
     return(log_fac)
     returnType(double())
-  })
+  },
+  buildDerivs = list(run = list(ignore = c("i", "max_index")))
+)
diff --git a/R/zzz.R b/R/zzz.R
index 66e4b46..a2bf8c1 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -296,8 +296,8 @@
       types = c('value = double(1)',
                 'lambda = double()',
                 'prob = double(1)',
-                'Nmin = integer(0, default = -1)',
-                'Nmax = integer(0, default = -1)',
+                'Nmin = double()',
+                'Nmax = double()',
                 'len = integer()'
                 ),
       mixedSizes = FALSE,
@@ -313,8 +313,8 @@
       types = c('value = double(1)',
                 'lambda = double()',
                 'prob = double()',
-                'Nmin = integer(0, default = -1)',
-                'Nmax = integer(0, default = -1)',
+                'Nmin = double()',
+                'Nmax = double()',
                 'len = integer()'
                 ),
       mixedSizes = FALSE,
@@ -331,8 +331,8 @@
                 'lambda = double()',
                 'theta = double()',
                 'prob = double(1)',
-                'Nmin = double(0, default = -1)',
-                'Nmax = double(0, default = -1)',
+                'Nmin = double()',
+                'Nmax = double()',
                 'len = double()'
       ),
       mixedSizes = FALSE,
@@ -348,8 +348,8 @@
                 'lambda = double()',
                 'theta = double()',
                 'prob = double()',
-                'Nmin = double(0, default = -1)',
-                'Nmax = double(0, default = -1)',
+                'Nmin = double()',
+                'Nmax = double()',
                 'len = double()'
       ),
       mixedSizes = FALSE,
@@ -365,8 +365,8 @@
                 'lambda = double()',
                 'theta = double()',
                 'prob = double()',
-                'Nmin = double(0, default = -1)',
-                'Nmax = double(0, default = -1)',
+                'Nmin = double()',
+                'Nmax = double()',
                 'len = double()'
       ),
       mixedSizes = FALSE,
@@ -403,7 +403,7 @@
                 's = double()',
                 'Nmin = double(0, default = -1)',
                 'Nmax = double(0, default = -1)',
-                'len = double()'
+                'len = integer()'
       ),
       mixedSizes = FALSE,
       pqAvail = FALSE
@@ -462,7 +462,7 @@
                 'theta = double()',
                 'Nmin = double(0, default = -1)',
                 'Nmax = double(0, default = -1)',
-                'len = double()'
+                'len = integer()'
       ),
       mixedSizes = FALSE,
       pqAvail = FALSE
diff --git a/inst/AD_test_utils.R b/inst/AD_test_utils.R
index 792bde5..2390c2e 100644
--- a/inst/AD_test_utils.R
+++ b/inst/AD_test_utils.R
@@ -1255,7 +1255,7 @@ setup_update_and_constant_nodes_for_tests <- function(model,
   ##  derivNodes <- model$expandNodeNames(derivNodes) # do not do this because do not want vector node names
   nNodes <- length(derivNodes)
   calcNodes <- model$getDependencies(derivNodes)
-  ucNodes <- makeDerivsInfo(model, derivNodes, calcNodes, dataAsConstantNodes = TRUE)
+  ucNodes <- makeModelDerivsInfo(model, derivNodes, calcNodes, dataAsConstantNodes = TRUE)
   updateNodes <- ucNodes$updateNodes
   constantNodes <- ucNodes$constantNodes
   updateNodes <- setdiff(updateNodes, forceConstantNodes) # remove forceConstants from updates
diff --git a/inst/test_utils.R b/inst/test_utils.R
index bbdaa90..7097c59 100644
--- a/inst/test_utils.R
+++ b/inst/test_utils.R
@@ -1314,7 +1314,7 @@ derivsNimbleFunctionMeta <- nimbleFunction(
   setup = function(model, calcNodes, wrt, reset = FALSE) {
     innerWrtVec <- seq_along(model$expandNodeNames(wrt, returnScalarComponents = TRUE))
     d <- length(innerWrtVec)
-    derivsInfo <- makeDerivsInfo(model, wrt, calcNodes)
+    derivsInfo <- makeModelDerivsInfo(model, wrt, calcNodes)
     updateNodes <- derivsInfo$updateNodes
     constantNodes <- derivsInfo$constantNodes
   },
@@ -1375,7 +1375,7 @@ derivsNimbleFunctionParamTransform <- nimbleFunction(
     my_parameterTransform <- parameterTransform(model, wrtNodesAsScalars)
     d <- my_parameterTransform$getTransformedLength()
     nimDerivs_wrt <- 1:d
-    derivsInfo <- makeDerivsInfo(model, wrt, calcNodes)
+    derivsInfo <- makeModelDerivsInfo(model, wrt, calcNodes)
     updateNodes   <- derivsInfo$updateNodes
     constantNodes <- derivsInfo$constantNodes
   },
@@ -1407,7 +1407,7 @@ derivsNimbleFunctionParamTransformMeta <- nimbleFunction(
     my_parameterTransform <- parameterTransform(model, wrtNodesAsScalars)
     d <- my_parameterTransform$getTransformedLength()
     nimDerivs_wrt <- 1:d
-    derivsInfo <- makeDerivsInfo(model, wrt, calcNodes)
+    derivsInfo <- makeModelDerivsInfo(model, wrt, calcNodes)
     updateNodes   <- derivsInfo$updateNodes
     constantNodes <- derivsInfo$constantNodes
   },
@@ -1811,7 +1811,7 @@ test_ADModelCalculate_internal <- function(model, name = 'unknown', xOrig = NULL
       if(!is.null(d) && length(d) == 2 && d[1] == d[2]) return(TRUE) else return(FALSE)
     }
 
-    derivsInfo <- makeDerivsInfo(model, wrt, calcNodes)
+    derivsInfo <- makeModelDerivsInfo(model, wrt, calcNodes)
     updateNodes <- derivsInfo$updateNodes
     constantNodes <- derivsInfo$constantNodes
 
diff --git a/install_adoak.R b/install_adoak.R
deleted file mode 100644
index 264680a..0000000
--- a/install_adoak.R
+++ /dev/null
@@ -1,4 +0,0 @@
-# For AD testing purposes only.
-# This file should be removed once AD functions are added to a nimble release.
-
-devtools::install_github("nimble-dev/nimble/packages/nimble", ref = "ADoak")
diff --git a/man/dNmixture.Rd b/man/dNmixture.Rd
index 414f14c..c06a090 100644
--- a/man/dNmixture.Rd
+++ b/man/dNmixture.Rd
@@ -26,28 +26,19 @@
 \alias{dNmixture_BBNB_oneObs}
 \title{N-mixture distribution for use in \code{nimble} models}
 \usage{
-dNmixture_v(x, lambda, prob, Nmin = -1, Nmax = -1, len, log = 0)
+dNmixture_v(x, lambda, prob, Nmin, Nmax, len, log = 0)
 
-dNmixture_s(x, lambda, prob, Nmin = -1, Nmax = -1, len, log = 0)
+dNmixture_s(x, lambda, prob, Nmin, Nmax, len, log = 0)
 
-rNmixture_v(n, lambda, prob, Nmin = -1, Nmax = -1, len)
+rNmixture_v(n, lambda, prob, Nmin, Nmax, len)
 
-rNmixture_s(n, lambda, prob, Nmin = -1, Nmax = -1, len)
+rNmixture_s(n, lambda, prob, Nmin, Nmax, len)
 
-dNmixture_BNB_v(x, lambda, theta, prob, Nmin = -1, Nmax = -1, len, log = 0)
+dNmixture_BNB_v(x, lambda, theta, prob, Nmin, Nmax, len, log = 0)
 
-dNmixture_BNB_s(x, lambda, theta, prob, Nmin = -1, Nmax = -1, len, log = 0)
+dNmixture_BNB_s(x, lambda, theta, prob, Nmin, Nmax, len, log = 0)
 
-dNmixture_BNB_oneObs(
-  x,
-  lambda,
-  theta,
-  prob,
-  Nmin = -1,
-  Nmax = -1,
-  len,
-  log = 0
-)
+dNmixture_BNB_oneObs(x, lambda, theta, prob, Nmin, Nmax, len, log = 0)
 
 dNmixture_BBP_v(x, lambda, prob, s, Nmin = -1, Nmax = -1, len, log = 0)
 
@@ -71,11 +62,11 @@ dNmixture_BBNB_oneObs(
   log = 0
 )
 
-rNmixture_BNB_v(n, lambda, theta, prob, Nmin = -1, Nmax = -1, len)
+rNmixture_BNB_v(n, lambda, theta, prob, Nmin, Nmax, len)
 
-rNmixture_BNB_s(n, lambda, theta, prob, Nmin = -1, Nmax = -1, len)
+rNmixture_BNB_s(n, lambda, theta, prob, Nmin, Nmax, len)
 
-rNmixture_BNB_oneObs(n, lambda, theta, prob, Nmin = -1, Nmax = -1, len)
+rNmixture_BNB_oneObs(n, lambda, theta, prob, Nmin, Nmax, len)
 
 rNmixture_BBP_v(n, lambda, prob, s, Nmin = -1, Nmax = -1, len)
 
diff --git a/tests/testthat/test-AD.R b/tests/testthat/test-AD.R
index 51c9d90..063309a 100644
--- a/tests/testthat/test-AD.R
+++ b/tests/testthat/test-AD.R
@@ -1,8 +1,9 @@
 # Testing examples:
 
-# install nimble from branch ADoak:
-# devtools::install_github("nimble-dev/nimble", ref = "ADoak", subdir = "packages/nimble")
-# install nimbleEcology from branch AD_0.3: devtools::install_github("nimble-dev/nimbleEcology", ref = "AD_0.3")
+# install nimble from CRAN:
+devtools::install_cran("nimble", force = TRUE)
+# devtools::install_github("nimble-dev/nimble", ref = "devel", subdir = "packages/nimble")
+# install nimbleEcology from branch Nmixture-AD: devtools::install_github("nimble-dev/nimbleEcology", ref = "Nmixture-AD")
 
 # load nimble's testing tools
 library(nimble)
@@ -33,8 +34,8 @@ nc <- nimbleCode({
 })
 
 Rmodel <- nimbleModel(nc, data = list(x = dat),
-                         inits = list(probOcc = probOcc,
-                                      probDetect = probDetect),
+                      inits = list(probOcc = probOcc,
+                                   probDetect = probDetect),
                       buildDerivs=TRUE)
 
 Cmodel <- compileNimble(Rmodel)
@@ -85,7 +86,8 @@ model_calculate_test_case(Rmodel, Cmodel,
                           v1_case1, v2_case1,
                           0:2)
 })
-test_that ("dNmixture errors on build with AD", {
+
+test_that ("dNmixture works with AD", {
 ##########################
 #### dNmixture_s case ####
 
@@ -104,26 +106,22 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                  inits = list(prob = prob,
                               lambda = lambda),
                  buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda'))
-# v1_case1 <- list(arg1 = c(prob, lambda)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda'))
+v1_case1 <- list(arg1 = c(prob, lambda)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2)) # testing values for prob and lambda
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##########################
 #### dNmixture_BNB_s case ####
@@ -144,26 +142,24 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda,
                                    theta = theta),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-})
-# Cmodel$calculate()
-#
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta'))
-# v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+Cmodel$calculate()
+
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta'))
+v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##############################
 #### dNmixture_BBP_s case ####
@@ -178,31 +174,30 @@ prob2 <- 0.5
 s2 <- 1.111
 
 nc <- nimbleCode({
-  x[1:5] ~ dNmixture_BBNB_s(lambda, prob, s = s,
+  x[1:5] ~ dNmixture_BBP_s(lambda, prob, s = s,
                             Nmin = 0, Nmax = 100, len = 5)
   prob ~ dunif(0, 1)
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda,
                                    s = s),
-                      buildDerivs=TRUE)
-# Rmodel$calculate()
+                      buildDerivs = TRUE)
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's'))
-# v1_case1 <- list(arg1 = c(prob, lambda, s)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+Cmodel$calculate()
+
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's'))
+v1_case1 <- list(arg1 = c(prob, lambda, s)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##############################
 #### dNmixture_BBNB_s case ####
@@ -225,26 +220,25 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda,
                                    theta = theta, s = s),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's'))
-# v1_case1 <- list(arg1 = c(prob, lambda, theta, s)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's'))
+v1_case1 <- list(arg1 = c(prob, lambda, theta, s)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##########################
 #### dNmixture_v case ####
@@ -266,26 +260,23 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda'))
+v1_case1 <- list(arg1 = c(prob, lambda)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2)) # testing values for prob and lambda
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda'))
-# v1_case1 <- list(arg1 = c(prob, lambda)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##############################
 #### dNmixture_BNB_v case ####
@@ -310,26 +301,24 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda,
                                    theta = theta),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta'))
-# v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta'))
+v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##############################
 #### dNmixture_BBP_v case ####
@@ -354,26 +343,24 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda,
                                    s = s),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's'))
-# v1_case1 <- list(arg1 = c(prob, lambda, s)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's'))
+v1_case1 <- list(arg1 = c(prob, lambda, s)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##############################
 #### dNmixture_BBNB_v case ####
@@ -400,26 +387,24 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob,
                                    lambda = lambda,
                                    theta = theta, s = s),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's'))
-# v1_case1 <- list(arg1 = c(prob, lambda, theta, s)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's'))
+v1_case1 <- list(arg1 = c(prob, lambda, theta, s)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##########################
 #### dNmixture_BNB_oneObs case ####
@@ -440,25 +425,23 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob, theta = theta,
                                    lambda = lambda),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta'))
-# v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta'))
+v1_case1 <- list(arg1 = c(prob, lambda, theta)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, theta2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##########################
 #### dNmixture_BBP_oneObs case ####
@@ -479,25 +462,23 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob, s=s,
                                    lambda = lambda),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's'))
-# v1_case1 <- list(arg1 = c(prob, lambda, s2)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 's'))
+v1_case1 <- list(arg1 = c(prob, lambda, s2)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, s2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 
 ##########################
 #### dNmixture_BBNB_oneObs case ####
@@ -520,25 +501,23 @@ nc <- nimbleCode({
   lambda ~ dunif(0, 100)
 })
 
-expect_error({
 Rmodel <- nimbleModel(nc, data = list(x = x),
                       inits = list(prob = prob, theta = theta, s=s,
                                    lambda = lambda),
                       buildDerivs=TRUE)
-# Rmodel$calculate()
+Rmodel$calculate()
 
 Cmodel <- compileNimble(Rmodel)
-# Cmodel$calculate()
-})
+Cmodel$calculate()
 
-# nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's'))
-# v1_case1 <- list(arg1 = c(prob, lambda, theta, s2)) # taping values for prob and lambda
-# v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda
-#
-# model_calculate_test_case(Rmodel, Cmodel,
-#                           model_calculate_test, nodesList_case1,
-#                           v1_case1, v2_case1,
-#                           0:2)
+nodesList_case1 <- setup_update_and_constant_nodes_for_tests(Rmodel, c('prob', 'lambda', 'theta', 's'))
+v1_case1 <- list(arg1 = c(prob, lambda, theta, s2)) # taping values for prob and lambda
+v2_case1 <- list(arg1 = c(prob2, lambda2, theta2, s2)) # testing values for prob and lambda
+
+model_calculate_test_case(Rmodel, Cmodel,
+                          model_calculate_test, nodesList_case1,
+                          v1_case1, v2_case1,
+                          0:2, RCrelTol = c(2e-15, 1e-8, 1e-3, 1e-14))
 })