tests/testthat/test-DynOcc.R

In nimble-dev/nimbleEcology: Distributions for Ecological Models in 'nimble'

context("Testing dDynOcc-related functions.")

test_that("dDynOcc_vvm works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- c(0.4, 0.4, 0.1)
  probColonize <- c(0.4, 0.2, 0.1)
  p <- matrix(rep(0.8, 20), nrow = 4)

  probX <- dDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    for (t in 1:nyears) {
      if (end[t] - start[t] > 0) {
        numObs <- sum(x[t,start[t]:end[t]])
        if (numObs < 0) {
          print("Error in dDynamicOccupancy: numObs < 0 but nrep[t] > 0\n")
          stop("Error in dDynamicOccupancy: numObs < 0 but nrep[t] > 0\n")
        }
        ProbOccAndCount <- ProbOccNextTime *
            exp(sum(dbinom(x[t,start[t]:end[t]],
                           size = 1, prob = p[t,start[t]:end[t]], log = 1)))
        ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
        ProbCount <- ProbOccAndCount + ProbUnoccAndCount
        ProbOccGivenCount <- ProbOccAndCount / ProbCount
        ll <- ll + log(ProbCount)
        if (t < nyears)
            ProbOccNextTime <- ProbOccGivenCount * probPersist[t] +
                (1 - ProbOccGivenCount) * probColonize[t]
      } else {
        ## If there were no observations in year t,
        ## simply propagate probability of occupancy forward
        if (t < nyears)
            ProbOccNextTime <- ProbOccNextTime *  probPersist[t] +
                (1 - ProbOccNextTime) * probColonize[t]
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_vvm <- compileNimble(dDynOcc_vvm)
  CprobX <- CdDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_vvm(init, probPersist[1:3],
                              probColonize[1:3], p[1:4,1:5],
                              start[1:4], end[1:4])

  })

  m <- nimbleModel(code = nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_vvm(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_vvm <- compileNimble(rDynOcc_vvm)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_vvm(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)

# # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))

})


test_that("dDynOcc_vsm works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- c(0.4, 0.4, 0.1)
  probColonize <- 0.5
  p <- matrix(rep(0.8, 20), nrow = 4)

  probX <- dDynOcc_vsm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_vsm(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    if (nyears >= 1) {
      for (t in 1:nyears) {
        if (end[t] - start[t] + 1 > 0) {
          numObs <- sum(x[t,start[t]:end[t]])
          if (numObs < 0) {
            print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
            stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          }
          ProbOccAndCount <- ProbOccNextTime *
              exp(sum(dbinom(x[t,start[t]:end[t]],
                             size = 1, prob = p[t,start[t]:end[t]], log = 1)))
          ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
          ProbCount <- ProbOccAndCount + ProbUnoccAndCount
          ProbOccGivenCount <- ProbOccAndCount / ProbCount
          ll <- ll + log(ProbCount)
          if (t < nyears)
              ProbOccNextTime <- ProbOccGivenCount * probPersist[t] +
                  (1 - ProbOccGivenCount) * probColonize
        } else {
          ## If there were no observations in year t,
          ## simply propagate probability of occupancy forward
          if (t < nyears)
              ProbOccNextTime <- ProbOccNextTime *  probPersist[t] +
                  (1 - ProbOccNextTime) * probColonize
        }
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_vsm <- compileNimble(dDynOcc_vsm)
  CprobX <- CdDynOcc_vsm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_vsm(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_vsm(init, probPersist[1:3], probColonize, p[1:4,1:5],
                              start[1:4], end[1:4])

    init ~ dunif(0,1)

    for (i in 1:3) {
      probPersist[i] ~ dunif(0,1)
    }
    probColonize ~ dunif(0,1)

    for (i in 1:4) {
      for (j in 1:5) {
        p[i,j] ~ dunif(0,1)
      }
    }
  })

  m <- nimbleModel(nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_vsm(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_vsm <- compileNimble(rDynOcc_vsm)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_vsm(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)

# # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})

test_that("dDynOcc_svm works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- 0.5
  probColonize <- c(0.4, 0.4, 0.1)
  p <- matrix(rep(0.8, 20), nrow = 4)

  probX <- dDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
    if(nyears >= 1) {
      for(t in 1:nyears) {
        if (end[t] - start[t] + 1 > 0) {
          numObs <- sum(x[t,start[t]:end[t]])
          if (numObs < 0) {
            print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
            stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          }
          ProbOccAndCount <- ProbOccNextTime *
              exp(sum(dbinom(x[t,start[t]:end[t]],
                             size = 1, prob = p[t,start[t]:end[t]], log = 1)))
          ProbUnoccAndCount <- (1-ProbOccNextTime) * (numObs == 0)
          ProbCount <- ProbOccAndCount + ProbUnoccAndCount
          ProbOccGivenCount <- ProbOccAndCount / ProbCount
          ll <- ll + log(ProbCount)
          if(t < nyears)
              ProbOccNextTime <- ProbOccGivenCount * probPersist +
                  (1-ProbOccGivenCount) * probColonize[t]
        } else {
          ## If there were no observations in year t,
          ## simply propagate probability of occupancy forward
          if (t < nyears)
              ProbOccNextTime <- ProbOccNextTime * probPersist +
                  (1-ProbOccNextTime) * probColonize[t]
        }
      }
    }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_svm <- compileNimble(dDynOcc_svm)
  CprobX <- CdDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)



  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_svm(init, probPersist, probColonize[1:3], p[1:4,1:5],
                             start[1:4], end[1:4])

    init ~ dunif(0,1)

    probPersist ~ dunif(0,1)
    for (i in 1:3) {
      probColonize[i] ~ dunif(0,1)
    }


    for (i in 1:4) {
      for (j in 1:5) {
        p[i,j] ~ dunif(0,1)
      }
    }
  })

  m <- nimbleModel(nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_svm(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_svm <- compileNimble(rDynOcc_svm)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_svm(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)

#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("dDynOcc_ssm works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- 0.5
  probColonize <- 0.2
  p <- matrix(rep(0.8, 20), nrow = 4)

  probX <- dDynOcc_ssm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_ssm(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    for (t in 1:nyears) {
      if (end[t] - start[t] + 1 > 0) {
        numObs <- sum(x[t,start[t]:end[t]])
        if (numObs < 0) {
          print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
        }
        ProbOccAndCount <- ProbOccNextTime *
            exp(sum(dbinom(x[t,start[t]:end[t]],
                           size = 1, prob = p[t,start[t]:end[t]], log = 1)))
        ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
        ProbCount <- ProbOccAndCount + ProbUnoccAndCount
        ProbOccGivenCount <- ProbOccAndCount / ProbCount
        ll <- ll + log(ProbCount)
        if (t < nyears)
            ProbOccNextTime <- ProbOccGivenCount * probPersist +
                (1 - ProbOccGivenCount) * probColonize
      } else {
        ## If there were no observations in year t,
        ## simply propagate probability of occupancy forward
        if (t < nyears)
            ProbOccNextTime <- ProbOccNextTime * probPersist +
                (1 - ProbOccNextTime) * probColonize
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_ssm <- compileNimble(dDynOcc_ssm)
  CprobX <- CdDynOcc_ssm(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_ssm(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_ssm(init, probPersist, probColonize, p[1:4, 1:5],
                              start[1:4], end[1:4])

    init ~ dunif(0,1)

    probPersist ~ dunif(0,1)
    probColonize ~ dunif(0,1)

    for (i in 1:2) {
      for (j in 1:5) {
        p[i,j] ~ dunif(0,1)
      }
    }
  })

  m <- nimbleModel(code = nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_ssm(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_ssm <- compileNimble(rDynOcc_ssm)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_ssm(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)

#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("Case errors in compiled dDynOcc_**m work", {
  CdDynOcc_ssm <- compileNimble(dDynOcc_ssm)
  CdDynOcc_svm <- compileNimble(dDynOcc_svm)
  CdDynOcc_vsm <- compileNimble(dDynOcc_vsm)
  CdDynOcc_vvm <- compileNimble(dDynOcc_vvm)


  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  nrep <- c(5, 5)
  init <- 0.7
  probPersist <- 0.8
  probColonize <- 0.5
  p <- matrix(rep(0.8, 20), nrow = 4)

  expect_error(CdDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vsm(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  probPersist <- c(0.8, 0.5, 0.2)
  expect_error(CdDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  probColonize <- c(0.2, 0.5, 0.5)

  probPersist <- 0.3
  expect_error(CdDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  p <- p[1:2, 1:4]

  expect_error(CdDynOcc_svm(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  probPersist <- c(0.3, 0.5, 0.3)
  expect_error(CdDynOcc_vvm(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  probColonize <- 0.5
  expect_error(CdDynOcc_vsm(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  probPersist <- 0.3
  expect_error(CdDynOcc_ssm(x, init, probPersist, probColonize, p, start, end, log = FALSE))

})




test_that("dDynOcc_vvv works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- c(0.4, 0.4, 0.1)
  probColonize <- c(0.4, 0.2, 0.1)
  p <- rep(0.8, 4)

  probX <- dDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    for (t in 1:nyears) {
      if (end[t] - start[t] > 0) {
        numObs <- sum(x[t,start[t]:end[t]])
        if (numObs < 0) {
          print("Error in dDynamicOccupancy: numObs < 0 but nrep[t] > 0\n")
          stop("Error in dDynamicOccupancy: numObs < 0 but nrep[t] > 0\n")
        }
        ProbOccAndCount <- ProbOccNextTime *
            exp(sum(dbinom(x[t,start[t]:end[t]],
                           size = 1, prob = p[t], log = 1)))
        ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
        ProbCount <- ProbOccAndCount + ProbUnoccAndCount
        ProbOccGivenCount <- ProbOccAndCount / ProbCount
        ll <- ll + log(ProbCount)
        if (t < nyears)
            ProbOccNextTime <- ProbOccGivenCount * probPersist[t] +
                (1 - ProbOccGivenCount) * probColonize[t]
      } else {
        ## If there were no observations in year t,
        ## simply propagate probability of occupancy forward
        if (t < nyears)
            ProbOccNextTime <- ProbOccNextTime *  probPersist[t] +
                (1 - ProbOccNextTime) * probColonize[t]
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_vvv <- compileNimble(dDynOcc_vvv)
  CprobX <- CdDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_vvv(init = init, probPersist = probPersist[1:3],
                              probColonize = probColonize[1:3], p = p[1:4],
                              start = start[1:4], end = end[1:4])

    init ~ dunif(0,1)

    for (i in 1:3) {
      probPersist[i] ~ dunif(0,1)
      probColonize[i] ~ dunif(0,1)
    }

    for (i in 1:4) {
      p[i] ~ dunif(0,1)
    }
  })

  m <- nimbleModel(code = nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_vvv(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_vvv <- compileNimble(rDynOcc_vvv)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_vvv(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)
# # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("dDynOcc_vsv works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- c(0.4, 0.4, 0.1)
  probColonize <- 0.5
  p <- rep(0.8, 4)

  probX <- dDynOcc_vsv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_vsv(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    if (nyears >= 1) {
      for (t in 1:nyears) {
        if (end[t] - start[t] + 1 > 0) {
          numObs <- sum(x[t,start[t]:end[t]])
          if (numObs < 0) {
            print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
            stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          }
          ProbOccAndCount <- ProbOccNextTime *
              exp(sum(dbinom(x[t,start[t]:end[t]],
                             size = 1, prob = p[t], log = 1)))
          ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
          ProbCount <- ProbOccAndCount + ProbUnoccAndCount
          ProbOccGivenCount <- ProbOccAndCount / ProbCount
          ll <- ll + log(ProbCount)
          if (t < nyears)
              ProbOccNextTime <- ProbOccGivenCount * probPersist[t] +
                  (1 - ProbOccGivenCount) * probColonize
        } else {
          ## If there were no observations in year t,
          ## simply propagate probability of occupancy forward
          if (t < nyears)
              ProbOccNextTime <- ProbOccNextTime *  probPersist[t] +
                  (1 - ProbOccNextTime) * probColonize
        }
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_vsv <- compileNimble(dDynOcc_vsv)
  CprobX <- CdDynOcc_vsv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_vsv(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_vsv(init, probPersist[1:3], probColonize, p[1:4],
                              start[1:4], end[1:4])

    init ~ dunif(0,1)

    for (i in 1:3) {
      probPersist[i] ~ dunif(0,1)
    }
    probColonize ~ dunif(0,1)

    for (i in 1:4) {
      p[i] ~ dunif(0,1)
    }
  })

  m <- nimbleModel(nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_vsv(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_vsv <- compileNimble(rDynOcc_vsv)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_vsv(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)

#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})

test_that("dDynOcc_svv works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- 0.5
  probColonize <- c(0.4, 0.4, 0.1)
  p <- rep(0.8, 4)

  probX <- dDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
    if(nyears >= 1) {
      for(t in 1:nyears) {
        if (end[t] - start[t] + 1 > 0) {
          numObs <- sum(x[t,start[t]:end[t]])
          if (numObs < 0) {
            print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
            stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          }
          ProbOccAndCount <- ProbOccNextTime *
              exp(sum(dbinom(x[t,start[t]:end[t]],
                             size = 1, prob = p[t], log = 1)))
          ProbUnoccAndCount <- (1-ProbOccNextTime) * (numObs == 0)
          ProbCount <- ProbOccAndCount + ProbUnoccAndCount
          ProbOccGivenCount <- ProbOccAndCount / ProbCount
          ll <- ll + log(ProbCount)
          if(t < nyears)
              ProbOccNextTime <- ProbOccGivenCount * probPersist +
                  (1-ProbOccGivenCount) * probColonize[t]
        } else {
          ## If there were no observations in year t,
          ## simply propagate probability of occupancy forward
          if (t < nyears)
              ProbOccNextTime <- ProbOccNextTime * probPersist +
                  (1-ProbOccNextTime) * probColonize[t]
        }
      }
    }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_svv <- compileNimble(dDynOcc_svv)
  CprobX <- CdDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)



  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_svv(init, probPersist, probColonize[1:3], p[1:4],
                             start[1:4], end[1:4])

    init ~ dunif(0,1)

    probPersist ~ dunif(0,1)
    for (i in 1:3) {
      probColonize[i] ~ dunif(0,1)
    }


    for (i in 1:4) {
      p[i] ~ dunif(0,1)
    }
  })

  m <- nimbleModel(nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_svv(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_svv <- compileNimble(rDynOcc_svv)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_svv(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)

#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("dDynOcc_ssv works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- 0.5
  probColonize <- 0.2
  p <- rep(0.8, 4)

  probX <- dDynOcc_ssv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_ssv(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    for (t in 1:nyears) {
      if (end[t] - start[t] + 1 > 0) {
        numObs <- sum(x[t,start[t]:end[t]])
        if (numObs < 0) {
          print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
        }
        ProbOccAndCount <- ProbOccNextTime *
            exp(sum(dbinom(x[t,start[t]:end[t]],
                           size = 1, prob = p[t], log = 1)))
        ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
        ProbCount <- ProbOccAndCount + ProbUnoccAndCount
        ProbOccGivenCount <- ProbOccAndCount / ProbCount
        ll <- ll + log(ProbCount)
        if (t < nyears)
            ProbOccNextTime <- ProbOccGivenCount * probPersist +
                (1 - ProbOccGivenCount) * probColonize
      } else {
        ## If there were no observations in year t,
        ## simply propagate probability of occupancy forward
        if (t < nyears)
            ProbOccNextTime <- ProbOccNextTime * probPersist +
                (1 - ProbOccNextTime) * probColonize
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_ssv <- compileNimble(dDynOcc_ssv)
  CprobX <- CdDynOcc_ssv(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_ssv(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_ssv(init, probPersist, probColonize, p[1:4],
                              start[1:4], end[1:4])

    init ~ dunif(0,1)

    probPersist ~ dunif(0,1)
    probColonize ~ dunif(0,1)

    for (i in 1:4) {
      p[i] ~ dunif(0,1)
    }
  })

  m <- nimbleModel(code = nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_ssv(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_ssv <- compileNimble(rDynOcc_ssv)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_ssv(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)
#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("Case errors in compiled dDynOcc_**v work", {
  CdDynOcc_ssv <- compileNimble(dDynOcc_ssv)
  CdDynOcc_svv <- compileNimble(dDynOcc_svv)
  CdDynOcc_vsv <- compileNimble(dDynOcc_vsv)
  CdDynOcc_vvv <- compileNimble(dDynOcc_vvv)


  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  nrep <- c(5, 5)
  init <- 0.7
  probPersist <- 0.8
  probColonize <- 0.5
  p <- rep(0.8, 4)

  expect_error(CdDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vsv(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  probPersist <- c(0.8, 0.5, 0.2)
  expect_error(CdDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  probColonize <- c(0.2, 0.5, 0.5)

  probPersist <- 0.3
  expect_error(CdDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  p <- p[1:2]

  expect_error(CdDynOcc_svv(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  probPersist <- c(0.3, 0.5, 0.3)
  expect_error(CdDynOcc_vvv(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  probColonize <- 0.5
  expect_error(CdDynOcc_vsv(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  probPersist <- 0.3
  expect_error(CdDynOcc_ssv(x, init, probPersist, probColonize, p, start, end, log = FALSE))

})


context("Testing dDynOcc-related functions.")

test_that("dDynOcc_vvs works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- c(0.4, 0.4, 0.1)
  probColonize <- c(0.4, 0.2, 0.1)
  p <- 0.6

  probX <- dDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    for (t in 1:nyears) {
      if (end[t] - start[t] > 0) {
        numObs <- sum(x[t,start[t]:end[t]])
        if (numObs < 0) {
          print("Error in dDynamicOccupancy: numObs < 0 but nrep > 0\n")
          stop("Error in dDynamicOccupancy: numObs < 0 but nrep > 0\n")
        }
        ProbOccAndCount <- ProbOccNextTime *
            exp(sum(dbinom(x[t,start[t]:end[t]],
                           size = 1, prob = p, log = 1)))
        ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
        ProbCount <- ProbOccAndCount + ProbUnoccAndCount
        ProbOccGivenCount <- ProbOccAndCount / ProbCount
        ll <- ll + log(ProbCount)
        if (t < nyears)
            ProbOccNextTime <- ProbOccGivenCount * probPersist[t] +
                (1 - ProbOccGivenCount) * probColonize[t]
      } else {
        ## If there were no observations in year t,
        ## simply propagate probability of occupancy forward
        if (t < nyears)
            ProbOccNextTime <- ProbOccNextTime *  probPersist[t] +
                (1 - ProbOccNextTime) * probColonize[t]
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_vvs <- compileNimble(dDynOcc_vvs)
  CprobX <- CdDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_vvs(init = init, probPersist = probPersist[1:3],
                              probColonize = probColonize[1:3], p = p,
                              start = start[1:4], end = end[1:4])

    init ~ dunif(0,1)

    for (i in 1:3) {
      probPersist[i] ~ dunif(0,1)
      probColonize[i] ~ dunif(0,1)
    }

    p ~ dunif(0,1)
  })

  m <- nimbleModel(code = nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_vvs(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_vvs <- compileNimble(rDynOcc_vvs)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_vvs(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)
#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("dDynOcc_vss works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- c(0.4, 0.4, 0.1)
  probColonize <- 0.5
  p <- 0.6

  probX <- dDynOcc_vss(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_vss(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    if (nyears >= 1) {
      for (t in 1:nyears) {
        if (end[t] - start[t] + 1 > 0) {
          numObs <- sum(x[t,start[t]:end[t]])
          if (numObs < 0) {
            print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
            stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          }
          ProbOccAndCount <- ProbOccNextTime *
              exp(sum(dbinom(x[t,start[t]:end[t]],
                             size = 1, prob = p, log = 1)))
          ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
          ProbCount <- ProbOccAndCount + ProbUnoccAndCount
          ProbOccGivenCount <- ProbOccAndCount / ProbCount
          ll <- ll + log(ProbCount)
          if (t < nyears)
              ProbOccNextTime <- ProbOccGivenCount * probPersist[t] +
                  (1 - ProbOccGivenCount) * probColonize
        } else {
          ## If there were no observations in year t,
          ## simply propagate probability of occupancy forward
          if (t < nyears)
              ProbOccNextTime <- ProbOccNextTime *  probPersist[t] +
                  (1 - ProbOccNextTime) * probColonize
        }
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_vss <- compileNimble(dDynOcc_vss)
  CprobX <- CdDynOcc_vss(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_vss(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_vss(init, probPersist[1:3], probColonize, p,
                              start[1:4], end[1:4])

    init ~ dunif(0,1)

    for (i in 1:3) {
      probPersist[i] ~ dunif(0,1)
    }
    probColonize ~ dunif(0,1)

    p ~ dunif(0,1)
  })

  m <- nimbleModel(nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_vss(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_vss <- compileNimble(rDynOcc_vss)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_vss(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)
#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})

test_that("dDynOcc_svs works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- 0.5
  probColonize <- c(0.4, 0.4, 0.1)
  p <- 0.6

  probX <- dDynOcc_svs(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_svs(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
    if(nyears >= 1) {
      for(t in 1:nyears) {
        if (end[t] - start[t] + 1 > 0) {
          numObs <- sum(x[t,start[t]:end[t]])
          if (numObs < 0) {
            print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
            stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          }
          ProbOccAndCount <- ProbOccNextTime *
              exp(sum(dbinom(x[t,start[t]:end[t]],
                             size = 1, prob = p, log = 1)))
          ProbUnoccAndCount <- (1-ProbOccNextTime) * (numObs == 0)
          ProbCount <- ProbOccAndCount + ProbUnoccAndCount
          ProbOccGivenCount <- ProbOccAndCount / ProbCount
          ll <- ll + log(ProbCount)
          if(t < nyears)
              ProbOccNextTime <- ProbOccGivenCount * probPersist +
                  (1-ProbOccGivenCount) * probColonize[t]
        } else {
          ## If there were no observations in year t,
          ## simply propagate probability of occupancy forward
          if (t < nyears)
              ProbOccNextTime <- ProbOccNextTime * probPersist +
                  (1-ProbOccNextTime) * probColonize[t]
        }
      }
    }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_svs <- compileNimble(dDynOcc_svs)
  CprobX <- CdDynOcc_svs(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_svs(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)



  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_svs(init, probPersist, probColonize[1:3], p,
                             start[1:4], end[1:4])

    init ~ dunif(0,1)

    probPersist ~ dunif(0,1)
    for (i in 1:3) {
      probColonize[i] ~ dunif(0,1)
    }

    p ~ dunif(0,1)
  })

  m <- nimbleModel(nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))


  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_svs(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_svs <- compileNimble(rDynOcc_svs)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_svs(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)
#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("dDynOcc_sss works", {
  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)
  init <- 0.7
  probPersist <- 0.5
  probColonize <- 0.2
  p <- 0.8

  probX <- dDynOcc_sss(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  lProbX <- dDynOcc_sss(x, init, probPersist, probColonize, p, start, end, log = TRUE)

  ProbOccNextTime <- init
  ll <- 0
  nyears <- dim(x)[1]
  if (nyears >= 1) {
    for (t in 1:nyears) {
      if (end[t] - start[t] + 1 > 0) {
        numObs <- sum(x[t,start[t]:end[t]])
        if (numObs < 0) {
          print("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
          stop("Error in dDynamicOccupancy: numObs < 0 but number of obs in start/end > 0\n")
        }
        ProbOccAndCount <- ProbOccNextTime *
            exp(sum(dbinom(x[t,start[t]:end[t]],
                           size = 1, prob = p, log = 1)))
        ProbUnoccAndCount <- (1 - ProbOccNextTime) * (numObs == 0)
        ProbCount <- ProbOccAndCount + ProbUnoccAndCount
        ProbOccGivenCount <- ProbOccAndCount / ProbCount
        ll <- ll + log(ProbCount)
        if (t < nyears)
            ProbOccNextTime <- ProbOccGivenCount * probPersist +
                (1 - ProbOccGivenCount) * probColonize
      } else {
        ## If there were no observations in year t,
        ## simply propagate probability of occupancy forward
        if (t < nyears)
            ProbOccNextTime <- ProbOccNextTime * probPersist +
                (1 - ProbOccNextTime) * probColonize
      }
    }
  }

  lCorrectProbX <- ll

  expect_equal(exp(lCorrectProbX), probX)
  expect_equal(lCorrectProbX, lProbX)

  CdDynOcc_sss <- compileNimble(dDynOcc_sss)
  CprobX <- CdDynOcc_sss(x, init, probPersist, probColonize, p, start, end, log = FALSE)
  expect_equal(CprobX, probX)

  ClProbX <- CdDynOcc_sss(x, init, probPersist, probColonize, p, start, end, log = TRUE)
  expect_equal(ClProbX, lProbX)

  nc <- nimbleCode({

    x[1:4, 1:5] ~ dDynOcc_sss(init, probPersist, probColonize, p,
                              start[1:4], end[1:4])

    init ~ dunif(0,1)

    probPersist ~ dunif(0,1)
    probColonize ~ dunif(0,1)

    p ~ dunif(0,1)
  })

  m <- nimbleModel(code = nc, data = list(x = x),
                   constants = list(start = start, end = end),
                   inits = list(p = p, probPersist = probPersist,
                                init = init, probColonize = probColonize))

  m$calculate()
  MlProbX <- m$getLogProb("x")
  expect_equal(MlProbX, lProbX)

  cm <- compileNimble(m, showCompilerOutput = TRUE)
  cm$calculate()
  CMlProbX <- cm$getLogProb("x")
  expect_equal(CMlProbX, lProbX)

  set.seed(2468)
  cm$simulate('x')

  # Test simulation code
  set.seed(1)
  nSim <- 10
  xSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    xSim[i,,] <- rDynOcc_sss(1, init, probPersist, probColonize, p, start, end)
  set.seed(1)
  CrDynOcc_sss <- compileNimble(rDynOcc_sss)
  CxSim <- array(NA, dim = c(nSim, dim(x)))
  for(i in 1:nSim)
    CxSim[i,,] <- CrDynOcc_sss(1, init, probPersist, probColonize, p, start, end)
  expect_identical(xSim, CxSim)

  simNodes <- m$getDependencies(c('probPersist', 'probColonize', 'init', 'p'), self = FALSE)
  mxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    m$simulate(simNodes, includeData = TRUE)
    mxSim[i,,] <- m$x
  }
  expect_identical(mxSim, xSim)

  CmxSim <- array(NA, dim = c(nSim, dim(x)))
  set.seed(1)
  for(i in 1:nSim) {
    cm$simulate(simNodes, includeData = TRUE)
    CmxSim[i,,] <- cm$x
  }
  expect_identical(CmxSim, mxSim)
#   # Test imputing value for all NAs
#   xNA <- array(NA, dim = dim(x))
#   mNA <- nimbleModel(code = nc, data = list(x = xNA),
#                      constants = list(start = start, end = end),
#                      inits = list(p = p, probPersist = probPersist,
#                                   init = init, probColonize = probColonize))
#   mNAConf <- configureMCMC(mNA)
#   mNAConf$addMonitors('x')
#   mNA_MCMC <- buildMCMC(mNAConf)
#   cmNA <- compileNimble(mNA, mNA_MCMC)
#   set.seed(0)
#   cmNA$mNA_MCMC$run(10)
# # Did the imputed values come back?
#   expect_true(all(!is.na(as.matrix(cmNA$mNA_MCMC$mvSamples)[,"x[1]"])))
})


test_that("Case errors in compiled dDynOcc_** work", {
  CdDynOcc_sss <- compileNimble(dDynOcc_sss)
  CdDynOcc_svs <- compileNimble(dDynOcc_svs)
  CdDynOcc_vss <- compileNimble(dDynOcc_vss)
  CdDynOcc_vvs <- compileNimble(dDynOcc_vvs)


  x <- matrix(c(0,0,NA,0,
                1,1,1,0,
                0,0,0,0,
                0,0,1,0,
                0,0,0,NA), nrow = 4)
  nrep <- c(5, 5)
  init <- 0.7
  probPersist <- 0.8
  probColonize <- 0.5
  p <- 0.8
  start <- c(1,1,2,1)
  end <- c(5,5,5,4)

  expect_error(CdDynOcc_svs(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vss(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  probPersist <- c(0.8, 0.5, 0.2)
  expect_error(CdDynOcc_svs(x, init, probPersist, probColonize, p, start, end, log = FALSE))
  expect_error(CdDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = FALSE))

  probColonize <- c(0.2, 0.5, 0.5)

  probPersist <- 0.3
  expect_error(CdDynOcc_vvs(x, init, probPersist, probColonize, p, start, end, log = FALSE))

})