R/estConnectivity.R

Defines functions diffMantel diffMC getCMRexample estMantel estMC estMCisotope estMCGlGps estTransition estTransitionBoot estStrength

Documented in diffMantel diffMC estMantel estMC estStrength estTransition getCMRexample

###############################################################################
#' Estimate MC, migratory connectivity strength
#'
#' Resampling of uncertainty for MC (migratory connectivity strength)
#' from estimates of psi (transition probabilities) and/or relative abundance.
#' Psi estimates can come from an estMigConnectivity object, an RMark psi
#' matrix, MCMC samples, or other samples expressed in array form.
#' Abundance estimates for each origin site can be
#' either just point estimates (no uncertainty propagated) or MCMC samples.
#' Other inputs include distances between origin sites, distances between target
#' sites, and sample size used to estimate psi.
#'
#'
#' @param originDist Distances between the B origin sites. Symmetric B by B
#'  matrix
#' @param targetDist Distances between the W target sites. Symmetric W by W
#'  matrix
#' @param originRelAbund Relative abundance estimates at B origin sites. Either
#'  a numeric vector of length B that sums to 1, or an mcmc object (such as is
#'  produced by \code{\link{modelCountDataJAGS}}) or matrix with at least
#'  \code{nSamples} rows. If there are more than B columns, the relevant columns
#'  should be labeled "relN[1]" through "relN[B]"
#' @param psi Transition probabilities between B origin and W target sites.
#'  Either a matrix with B rows and W columns where rows sum to 1, an array with
#'  dimensions x, B, and W (with x samples of the transition probability matrix
#'  from another model), an 'estPsi' object (result of calling estTransition),
#'  or a MARK object with estimates of transition probabilities
#' @param sampleSize Total sample size of animals that psi will be estimated
#'  from. Should be the number of animals released in one of the origin sites
#'  and observed in one of the target sites (or vice-versa). Optional, but
#'  recommended, unless psi is an estPsi object (in which case this function can
#'  pull it from there)
#' @param originSites If \code{psi} is a MARK object, this must be a numeric
#'  vector indicating which sites are origin
#' @param targetSites If \code{psi} is a MARK object, this must be a numeric
#'  vector indicating which sites are target
#' @param originNames Optional. Vector of names for the origin sites. Mostly for
#'  internal use
#' @param targetNames Optional. Vector of names for the target sites. Mostly for
#'  internal use
#' @param nSamples Number of times to resample \code{psi} and/or
#'  \code{originRelAbund}. The purpose is to estimate sampling uncertainty;
#'  higher values here will do so with more precision
#' @param row0 If \code{originRelAbund} is an mcmc object or array, this can be
#'  set to 0 (default) or any greater integer to specify where to stop ignoring
#'  samples ("burn-in")
#' @param verbose 0 (default) to 2. 0 prints no output during run. 1 prints
#'  a progress update and summary every 100 samples. 2 prints a
#'  progress update and summary every sample
#' @param alpha Level for confidence/credible intervals provided. Default (0.05)
#'  gives 95 percent CI
#' @param approxSigTest Should function compute approximate one-sided
#'  significance tests (p-values) for MC from the resampling? Default is
#'  FALSE
#' @param sigConst Value to compare MC to in significance test. Default is 0
#' @param maintainLegacyOutput version 0.4.0 of \code{MigConnectivity}
#'  updated the structure of the estimates. If you have legacy code that refers
#'  to elements within an \code{estMigConnectivity} object (results of
#'  \code{estMC}), you can set this to TRUE to also keep the old structure.
#'  Defaults to FALSE
#' @param returnAllInput if TRUE (the default) the output includes all of the
#'  inputs. If FALSE, only the inputs currently used by another MigConnectivity
#'  function are included in the output.
#'
#' @return \code{estStrength} returns a list with the elements:
#' \describe{
#'   \item{\code{MC}}{List containing estimates of migratory connectivity
#'    strength:
#'    \itemize{
#'      \item{\code{sample}} \code{nSamples} sampled values for
#'       MC. Provided to allow the user to compute own summary statistics.
#'      \item{\code{mean}} Mean of \code{MC$sample}. Main estimate of MC,
#'       incorporating parametric uncertainty.
#'      \item{\code{se}} Standard error of MC, estimated from SD of
#'       \code{MC$sample}.
#'      \item{\code{simpleCI}} Default\code{1 - alpha} confidence interval for
#'       MC, estimated as \code{alpha/2} and \code{1 - alpha/2} quantiles of
#'       \code{MC$sample}.
#'      \item{\code{bcCI}} Bias-corrected \code{1 - alpha} confidence interval
#'       for MC. May be preferable to \code{MC$simpleCI} when \code{MC$mean} is
#'       the best estimate of MC. \code{MC$simpleCI} is preferred when
#'       \code{MC$median} is a better estimator. When \code{MC$mean==MC$median},
#'       these should be identical.  Estimated as the
#'       \code{pnorm(2 * z0 + qnorm(alpha / 2))} and
#'       \code{pnorm(2 * z0 + qnorm(1 - alpha / 2))} quantiles of \code{MC$sample},
#'       where z0 is the proportion of \code{MC$sample < MC$mean}.
#'      \item{\code{hpdCI}} \code{1 - alpha} credible interval for MC,
#'       estimated using the highest posterior density (HPD) method.
#'      \item{\code{median}} Median of MC, alternate point estimate also
#'       including parametric uncertainty.
#'      \item{\code{point}} Simple point estimate of MC, using the point
#'      estimates of \code{psi} and \code{originRelAbund} (usually the mean
#'      values), not accounting for sampling error.
#'      \item{\code{simpleP}} Approximate p-value for MC, estimated as the
#'      proportion of bootstrap iterations where MC < \code{sigConst} (or MC >
#'      \code{sigConst} if \code{pointMC < sigConst}).  Note that if the
#'      proportion is 0, a default value of 0.5 / \code{nSamples} is provided,
#'      but this is best interpreted as p < 1 / \code{nSamples}.  NULL when
#'      \code{approxSigTest==FALSE}.
#'      \item{\code{bcP}} Approximate bias-corrected p-value for MC, estimated as
#'      \code{pnorm(qnorm(simpleP) - 2 * z0)}, where z0 is the proportion of
#'      \code{sampleMC < meanMC}.  May be a better approximation of the p-value
#'      than \code{simpleP}, but many of the same limitations apply.  NULL when
#'      \code{approxSigTest==FALSE}.
#'    }
#'   }
#'   \item{\code{input}}{List containing the inputs to \code{estStrength}.}
#' }
#'
#' @export
#'
#' @seealso \code{\link{calcMC}}, \code{\link{estTransition}},
#'   \code{\link{estMC}}, \code{\link{estMantel}},
#'   \code{\link{plot.estMigConnectivity}}
#' @example inst/examples/estStrengthExamples.R
estStrength <- function(originDist, targetDist, originRelAbund, psi,
                        sampleSize = NULL,
                        originSites=NULL, targetSites=NULL,
                        originNames = NULL, targetNames = NULL,
                        nSamples = 1000, row0 = 0, verbose=0,
                        alpha = 0.05, approxSigTest = FALSE, sigConst = 0,
                        maintainLegacyOutput = FALSE,
                        returnAllInput = TRUE) {
  nOriginSites <- nrow(originDist)
  nTargetSites <- nrow(targetDist)
  absAbund <- !is.null(sampleSize)
  if (coda::is.mcmc(originRelAbund) || coda::is.mcmc.list(originRelAbund)) {
    originRelAbund <- as.matrix(originRelAbund)
  }
  if (is.matrix(originRelAbund) && all(dim(originRelAbund)>1)) {
    abundFixed <- FALSE
    if (dim(originRelAbund)[2]>nOriginSites)
      abundParams <- paste('relN[', 1:nOriginSites, ']', sep='')
    else if (dim(originRelAbund)[2]==nOriginSites)
      abundParams <- 1:nOriginSites
    else
      stop('Number of origin sites must be constant between distance matrix and abundance')
    if (dim(originRelAbund)[1] >= nSamples)
      abundRows <- round(seq(from = row0 + 1, to = dim(originRelAbund)[1],
                             length.out = nSamples))
    else
      stop("You need at least nSamples rows to originRelAbund")
    originRelAbund <- as.matrix(originRelAbund[abundRows, abundParams])
    abundBase <- colMeans(originRelAbund)
  }
  else {
    abundFixed <- TRUE
    if (length(originRelAbund)!=nOriginSites)
      stop('Number of origin sites must be constant between distance matrix and abundance')
    abundBase <- originRelAbund
  }
  if (is.matrix(psi)) {
    psiFixed <- TRUE
    psiVCV <- NULL
    if (nrow(psi)!=nOriginSites || ncol(psi)!=nTargetSites)
      stop('Size of psi matrix must be consistant with distance matrices')
    psiBase <- psi
    psiIn <- psi
  }
  else if (inherits(psi, "mark")) {
    psiFixed <- FALSE
    if (!is.numeric(originSites) || !is.numeric(targetSites))
      stop('Must specify which RMark Psi parameters represent origin and target sites')
    psiVCV <- psi$results$beta.vcv
    psiBase <- RMark::TransitionMatrix(RMark::get.real(psi, "Psi",
                                                       se=TRUE))[originSites,
                                                              targetSites]
    if (any(diag(psi$results$beta.vcv) < 0))
      stop("Can't sample model, negative beta variances")
    psiIn <- psi
  }
  else if (is.array(psi)) {
    if (length(dim(psi))!=3)
      stop('Psi should either be 2-(for fixed transition probabilities) or 3-dimensional array')
    psiFixed <- FALSE
    if (dim(psi)[2]!=nOriginSites || dim(psi)[3]!=nTargetSites)
      stop('Size of psi array must be consistent with distance matrices')
    psiBase <- apply(psi, 2:3, mean)
    psiVCV <- NULL
    if (dim(psi)[1]>=nSamples)
      psiSamples <- round(seq(from = 1, to = dim(psi)[1],
                              length.out = nSamples))
    else
      psiSamples <- sample.int(dim(psi)[1], nSamples, replace = TRUE)
    psiIn <- psi
  }
  else if (inherits(psi, "estPsi") || inherits(psi, "estMC")) {
    psiFixed <- FALSE
    psiIn <- psi
    psi <- psi$psi$sample
    if (dim(psi)[2]!=nOriginSites || dim(psi)[3]!=nTargetSites)
      stop('Size of psi array must be consistant with distance matrices')
    psiBase <- apply(psi, 2:3, mean)
    psiVCV <- NULL
    if (dim(psi)[1]>=nSamples)
      psiSamples <- 1:nSamples
    else
      psiSamples <- sample.int(dim(psi)[1], nSamples, replace = TRUE)
    if (is.null(sampleSize))
      sampleSize <- psiIn$input$sampleSize # or wherever we put it
  }
  pointMC <- ifelse(absAbund,
                    calcMC(originDist, targetDist, originRelAbund = abundBase,
                           psi = psiBase, sampleSize = sampleSize),
                    calcMC(originDist, targetDist, originRelAbund = abundBase,
                           psi = psiBase))
  sampleMC <- rep(NA, nSamples)
  psi.array <- array(NA, c(nSamples, nOriginSites, nTargetSites),
                     dimnames = list(NULL, originNames, targetNames))
  for (i in 1:nSamples) {
    if (verbose > 1 || verbose == 1 && i %% 100 == 0)
      cat("\tSample", i, "of", nSamples, "at", date(), "\n")
    # Generate random transition probability matrices
    if (psiFixed)
      psiNew <- psiBase
    else if (is.null(psiVCV))
      psiNew <- psi[psiSamples[i],,]
    else
      psiNew <- makePsiRand(psi, originSites, targetSites)
    psi.array[i, , ] <- psiNew
    # Point estimates of breeding densities
    if (abundFixed)
      abundNew <- abundBase
    else
      abundNew <- originRelAbund[i, abundParams]
    # Calculate MC for new psis and relative breeding densities
    sampleMC[i] <- ifelse(absAbund,
                          calcMC(originDist, targetDist, originRelAbund = abundNew,
                                 psi = psiNew, sampleSize = sampleSize),
                          calcMC(originDist, targetDist, originRelAbund = abundNew,
                                 psi = psiNew))
    if (verbose > 1 || verbose == 1 && i %% 100 == 0)
      cat(" MC mean:", mean(sampleMC, na.rm=TRUE),
          "SD:", sd(sampleMC, na.rm=TRUE),
          "low quantile:", quantile(sampleMC, alpha/2, na.rm=TRUE),
          "high quantile:", quantile(sampleMC, 1-alpha/2, na.rm=TRUE), "\n")
  }
  meanMC <- mean(sampleMC, na.rm=TRUE)
  medianMC <- median(sampleMC, na.rm=TRUE)
  seMC <- sd(sampleMC, na.rm=TRUE)
  # Calculate confidence intervals using quantiles of sampled MC
  simpleCI <- quantile(sampleMC, c(alpha/2, 1-alpha/2), na.rm=TRUE, type = 8,
                       names = FALSE)
  z0 <- qnorm(sum((sampleMC)<meanMC)/nSamples)
  bcCI <- quantile(sampleMC, pnorm(2*z0+qnorm(c(alpha/2, 1-alpha/2))),
                       na.rm=TRUE, names = FALSE)
  MC.mcmc <- coda::as.mcmc(sampleMC) # Ha!
  hpdCI <- as.vector(coda::HPDinterval(MC.mcmc, 1-alpha))
  if (!approxSigTest)
    simpleP <- bcP <- NULL
  else {
    if (pointMC > sigConst)
      simpleP <- sum(sampleMC < sigConst) / nSamples
    else
      simpleP <- sum(sampleMC > sigConst) / nSamples
    if (simpleP == 0)
      simpleP <- 0.5 / nSamples
    bcP <- pnorm(qnorm(simpleP) - 2 * z0)
    if (pointMC < sigConst)
      bcP <- 1 - bcP
  }
  if (returnAllInput) {
    input <- list(originDist = originDist, targetDist = targetDist,
                  originRelAbund = originRelAbund, psi = psiIn,
                  sampleSize = sampleSize, originSites = originSites,
                  targetSites = targetSites,
                  originNames = originNames,
                  targetNames = targetNames,
                  nSamples = nSamples, row0 = row0,
                  verbose = verbose, alpha = alpha,
                  approxSigTest = approxSigTest, sigConst = sigConst,
                  maintainLegacyOutput = maintainLegacyOutput,
                  returnAllInput = TRUE)
  }
  else {
    input <- list(sampleSize = sampleSize,
                  originNames = originNames, targetNames = targetNames,
                  alpha = alpha, returnAllInput = FALSE)
  }
  if (maintainLegacyOutput) {
    mc <- list(sampleMC=sampleMC, samplePsi = psi.array, pointPsi = psiBase,
                pointMC=pointMC, meanMC=meanMC,
                medianMC=medianMC, seMC=seMC, simpleCI=simpleCI,
                bcCI=bcCI, hpdCI=hpdCI, simpleP = simpleP, bcP = bcP,
                sampleCorr = NULL, pointCorr = NULL,
                meanCorr = NULL, medianCorr = NULL, seCorr=NULL,
                simpleCICorr=NULL, bcCICorr=NULL, inputSampleSize = sampleSize,
                alpha = alpha, sigConst = sigConst,
                MC = list(sample = sampleMC, mean = meanMC, se = seMC,
                          simpleCI = simpleCI, bcCI = bcCI, hpdCI = hpdCI,
                          median = medianMC, point = pointMC,
                          simpleP = simpleP, bcP = bcP),
                input = input)
  }
  else {
    mc <- list(MC = list(sample = sampleMC, mean = meanMC, se = seMC,
                          simpleCI = simpleCI, bcCI = bcCI, hpdCI = hpdCI,
                          median = medianMC, point = pointMC,
                          simpleP = simpleP, bcP = bcP),
               input = input)
  }
  class(mc) <- c("estMC", "estMigConnectivity")
  return(mc)
}


estTransitionJAGS <- function (banded, reencountered,
                               originAssignment = NULL, targetAssignment = NULL,
                               alpha = 0.05, nSamples = 1000, verbose=0,
                               originNames = NULL, targetNames = NULL,
                               nThin = 1, nBurnin = 5000, nChains = 3,
                               fixedZero = NULL, psiPrior = NULL,
                               returnAllInput = TRUE,
                               originPoints = NULL, targetPoints = NULL,
                               originSites = NULL, targetSites = NULL) {
  if (is.null(originAssignment) && !is.null(originPoints) &&
      !is.null(originSites)) {
    if (inherits(originSites,"SpatialPolygonsDataFrame")){
      originSites <- sf::st_as_sf(originSites)
    }
    if (!inherits(originPoints, "sf")){
      originPoints <- sf::st_as_sf(originPoints)
    }
    originAssignment <- suppressMessages(unclass(sf::st_intersects(x = originPoints,
                                                                   y = originSites,
                                                                   sparse = TRUE)))
  }
  if (is.null(targetAssignment) && !is.null(targetPoints) &&
      !is.null(targetSites)) {
    if (inherits(targetSites,"SpatialPolygonsDataFrame")){
      targetSites <- sf::st_as_sf(targetSites)
    }
    if (!inherits(targetPoints, "sf")){
      targetPoints <- sf::st_as_sf(targetPoints)
    }
    targetAssignment <- suppressMessages(unclass(sf::st_intersects(x = targetPoints,
                                                                   y = targetSites,
                                                                   sparse = TRUE)))
  }
  if (is.array(originAssignment) && length(dim(originAssignment))>1) {
    if (length(dim(originAssignment))>2 || !all(originAssignment %in% 0:1))
      stop("originAssigment should either be a vector of assignments or a matrix of 0s and 1s")
    originAssignment <- apply(originAssignment, 1, which.max)
  }
  if (is.array(targetAssignment) && length(dim(targetAssignment))>1) {
    if (length(dim(targetAssignment))>2 || !all(targetAssignment %in% 0:1))
      stop("targetAssigment should either be a vector of assignments or a matrix of 0s and 1s")
    targetAssignment <- apply(targetAssignment, 1, which.max)
  }
  jags.inits <- vector("list", nChains)
  if (is.null(banded)) {
    if (is.null(originAssignment) || is.null(targetAssignment)) {
      stop("If running estTransition with MCMC, need to provide banding and/or telemetry data")
    }
    nDim <- 0
    nTargetSites <- max(length(unique(targetAssignment)),
                        length(targetNames))
    nOriginSites <- max(length(unique(originAssignment)),
                        length(originNames))
    if (is.null(originNames))
      originNames <- LETTERS[1:nOriginSites]
    if (is.null(targetNames))
      targetNames <- as.character(1:nTargetSites)

    jags.data <- list(npop = nOriginSites, ndest = nTargetSites)
    sampleSize <- 0
    # Parameters to monitor
    params <- c("psi")
    for (i in 1:nChains)
      jags.inits[[i]] <- list(m0 = matrix(runif(nOriginSites * nTargetSites),
                                          nOriginSites, nTargetSites))
  }
  else {
    nDim <- length(dim(reencountered))
    if (is.null(originNames))
      originNames <- dimnames(reencountered)[[1]]
    if (is.null(targetNames))
      targetNames <- dimnames(reencountered)[[nDim]]
    nTargetSites <- dim(reencountered)[nDim]
    nOriginSites <- dim(reencountered)[1]
    if (nDim == 2) {
      nAges <- 1
      if (length(banded)!=nOriginSites)
        stop('Number of origin sites must be constant between reencountered and banded')
      nfound <- apply(reencountered, 1, sum)
      sampleSize <- sum(nfound)
      tmat <- cbind(reencountered, banded - nfound)
      # Data
      jags.data <- list(recmat = tmat, npop = nOriginSites,
                        ndest = nTargetSites, nreleased = banded)
      for (i in 1:nChains)
        jags.inits[[i]] <- list(m0 =  matrix(runif(nOriginSites * nTargetSites),
                                             nOriginSites, nTargetSites),
                                r = rbeta(nTargetSites, 1, 1))
    }
    else {
      nAges <- dim(banded)[2]
      if (dim(banded)[1]!=nOriginSites)
        stop('Number of origin sites must be consistant between reencountered and banded')
      if (dim(reencountered)[2]!=nAges)
        stop('Number of ages must be consistant between banded and reencountered')
      nfound <- apply(reencountered, 1:2, sum)
      sampleSize <- sum(nfound)
      tmat <- array(NA, c(nOriginSites, nAges, nTargetSites + 1))
      tmat[ , , 1:nTargetSites] <- reencountered
      tmat[ , , 1 + nTargetSites] <- banded - nfound
      # Data
      jags.data <- list(recmat = tmat, npop = nOriginSites, nages = nAges,
                        ndest = nTargetSites, nreleased = banded)
      for (i in 1:nChains)
        jags.inits[[i]] <- list(m0 = matrix(runif(nOriginSites * nTargetSites),
                                            nOriginSites, nTargetSites),
                                r = matrix(rbeta(nTargetSites * nAges, 1, 1),
                                           nAges, nTargetSites))
    }
    # Parameters to monitor
    params <- c("psi", "r")
  }
  if (is.null(psiPrior)) {
    psiPrior <- matrix(1, nOriginSites, nTargetSites)
  }
  # else {
  #   z <- which(psiPrior==0)
  #   if (length(z)>0) {
  #     psiFixed <- matrix(NA, nOriginSites, nTargetSites)
  #     psiFixed[z] <- 0
  #     jags.data$m0 <- psiFixed
  #     for (i in 1:nChains){
  #       jags.inits[[i]]$m0[z] <- NA
  #     }
  #   }
  # }
  jags.data$psiPrior <- psiPrior
  if (!is.null(originAssignment)) {
    telmat <- table(factor(originAssignment, levels = 1:nOriginSites),
                    factor(targetAssignment, levels = 1:nTargetSites))
    jags.data$telmat <- telmat
    jags.data$ntel <- as.vector(table(factor(originAssignment,
                                             levels = 1:nOriginSites)))
    sampleSize <- sampleSize + sum(jags.data$ntel)
  }

  if (!is.null(fixedZero)) {
    psiFixed <- matrix(NA, nOriginSites, nTargetSites)
    for (i in 1:nrow(fixedZero)) {
      psiFixed[fixedZero[i, 1], fixedZero[i, 2]] <- 0
      for (j in 1:nChains)
        jags.inits[[j]]$m0[fixedZero[i, 1], fixedZero[i, 2]] <- NA
    }
    jags.data$m0 <- psiFixed
  }
  file <- tempfile(fileext = ".txt")
  sink(file = file)
  cat("
model{
  # psi prior
  for(i in 1:npop){
    for(k in 1:ndest){
      m0[i, k] ~ dgamma(psiPrior[i, k], 1)
      psi[i, k] <- m0[i, k] / sum(m0[i, 1:ndest])
    } #k
  }#i")
  if (!is.null(banded)){
    if (nAges==1) {
      cat("
  # model for recoveries with known number of ringed
  for (i in 1:npop){
    for(k in 1:ndest){
      p[i, k] <- psi[i, k] * r[k]
    }
    p[i, (ndest+1)] <- 1 - sum(p[i, 1:ndest])
  }

  for(k in 1:ndest) {
    r[k] ~ dunif(0, 1)
  }
  for(i in 1:npop){
    recmat[i, 1:(ndest+1)] ~ dmulti(p[i, 1:(ndest+1)], nreleased[i])
  }")
    }
    else {
      cat("
  for (j in 1:nages) {
        for(k in 1:ndest) {
          r[j, k] ~ dunif(0, 1)
        }
  }
  for (i in 1:npop){
    for (j in 1:nages) {
      for(k in 1:ndest){
        p[i, j, k] <- psi[i, k] * r[j, k]
      }
      p[i, j, (ndest+1)] <- 1 - sum(p[i, j, 1:ndest])
    }
  }
  for(i in 1:npop){
    for (j in 1:nages) {
      recmat[i, j, 1:(ndest+1)] ~ dmulti(p[i, j, 1:(ndest+1)], nreleased[i, j])
    }
  }")
    }
  }
  if (!is.null(originAssignment)) {
    cat("
  for(i in 1:npop){
    telmat[i, 1:ndest] ~ dmulti(psi[i, 1:ndest], ntel[i])
  }
")
  }
  cat("}")
  sink()
  # print(nAges)
  # print(file)
  # print(jags.data)
  # print(rowSums(jags.data$recmat))
  # print(jags.inits)
  out <- R2jags::jags(data = jags.data, inits = jags.inits, params,
                      file,
                      n.chains = nChains, n.thin = nThin,
                      n.iter = nBurnin + ceiling(nSamples * nThin / nChains),
                      n.burnin = nBurnin, DIC = FALSE,
                      progress.bar = ifelse(verbose==0, 'none', 'text'))

  if (nChains > 1) {
    maxRhat <- max(out$BUGSoutput$summary[, "Rhat"], na.rm = TRUE)
    if (verbose>0) {
      if (maxRhat < 1.1)
        cat("Successful convergence based on Rhat values (all < 1.1).\n")
      else
        cat("**WARNING** Rhat values indicate convergence failure.\n")
    }
  }
  file.remove(file)
  psi <- out$BUGSoutput$sims.list$psi
  dimnames(psi) <- list(NULL, originNames, targetNames)
  meanPsi <- out$BUGSoutput$mean$psi
  simpleCIPsi <- apply(psi, 2:3, quantile, probs = c(alpha/2, 1-alpha/2),
                       na.rm=TRUE, names = FALSE)
  psi.matrix <- array(c(psi), c(dim(psi)[[1]], nOriginSites * nTargetSites),
                      list(NULL, paste(rep(originNames, nTargetSites),
                                       rep(targetNames, each = nOriginSites),
                                       sep = "#")))
  psi.mcmc <- coda::as.mcmc(psi.matrix)
  hpdCI <- coda::HPDinterval(psi.mcmc, 1-alpha)
  hpdCI <- array(hpdCI, c(nOriginSites, nTargetSites, 2),
                 list(originNames, targetNames, c("lower", "upper")))
  hpdCI <- aperm(hpdCI, c(3, 1, 2))
  bcCIPsi <- array(NA, dim = c(2, nOriginSites, nTargetSites),
                   dimnames = list(NULL, originNames, targetNames))
  for (i in 1:nOriginSites) {
    for (j in 1:nTargetSites) {
      psi.z0 <- qnorm(sum(psi[, i, j] < meanPsi[i, j], na.rm = TRUE) /
                        length(which(!is.na(psi[, i, j]))))
      bcCIPsi[ , i, j] <- quantile(psi[, i, j],
                                   pnorm(2 * psi.z0 + qnorm(c(alpha/2, 1-alpha/2))),
                                   na.rm=TRUE, names = FALSE)
    }
  }
  results <- list(psi = list(sample = psi, mean = meanPsi,
                             se = out$BUGSoutput$sd$psi,
                             simpleCI = simpleCIPsi, bcCI = bcCIPsi,
                             hpdCI = hpdCI,
                             median = out$BUGSoutput$median$psi))
  if (!is.null(out$BUGSoutput$sims.list$r)) {
    if (nAges == 1) {
      bcCIr <- array(NA, dim = c(2, nTargetSites),
                     dimnames = list(NULL, targetNames))
      for (j in 1:nTargetSites) {
        psi.z0 <- qnorm(sum(out$BUGSoutput$sims.list$r[ ,j] <
                              out$BUGSoutput$mean$r[j], na.rm = TRUE) /
                          length(which(!is.na(out$BUGSoutput$sims.list$r[ ,j]))))
        bcCIr[ , j] <- quantile(out$BUGSoutput$sims.list$r[ ,j],
                                pnorm(2 * psi.z0 + qnorm(c(alpha/2, 1-alpha/2))),
                                na.rm=TRUE, names = FALSE)
      }
      colnames(out$BUGSoutput$sims.list$r) <- names(out$BUGSoutput$mean$r) <-
        names(out$BUGSoutput$sd$r) <- names(out$BUGSoutput$median$r) <-
        targetNames
      results$r <- list(sample = out$BUGSoutput$sims.list$r,
                       mean = out$BUGSoutput$mean$r,
                       se = out$BUGSoutput$sd$r,
                       simpleCI = apply(out$BUGSoutput$sims.list$r, 2, quantile,
                                        probs = c(alpha/2, 1-alpha/2), type = 8),
                       bcCI = bcCIr, median = out$BUGSoutput$median$r)
    }
    else {
      bcCIr <- array(NA, dim = c(2, nAges, nTargetSites),
                     dimnames = list(NULL, NULL, targetNames))
      for (i in 1:nAges) {
        for (j in 1:nTargetSites) {
          psi.z0 <- qnorm(sum(out$BUGSoutput$sims.list$r[ , i, j] <
                              out$BUGSoutput$mean$r[i, j], na.rm = TRUE) /
                          length(which(!is.na(out$BUGSoutput$sims.list$r[,i,j]))))
          bcCIr[,i,j] <- quantile(out$BUGSoutput$sims.list$r[ , i, j],
                                  pnorm(2 * psi.z0 + qnorm(c(alpha/2,1-alpha/2))),
                                  na.rm=TRUE, names = FALSE)
        }
      }
      dimnames(out$BUGSoutput$sims.list$r)[[3]]<-colnames(out$BUGSoutput$mean$r) <-
        colnames(out$BUGSoutput$sd$r) <- colnames(out$BUGSoutput$median$r) <-
        targetNames
      results$r <- list(sample = out$BUGSoutput$sims.list$r,
                       mean = out$BUGSoutput$mean$r,
                       se = out$BUGSoutput$sd$r,
                       simpleCI = apply(out$BUGSoutput$sims.list$r, 2:3, quantile,
                                        probs = c(alpha/2, 1-alpha/2), type = 8),
                       bcCI = bcCIr, median = out$BUGSoutput$median$r)
    }
  }
  if (returnAllInput) {
    results$input <- list(banded = banded, reencountered = reencountered,
                          originAssignment = originAssignment,
                          targetAssignment = targetAssignment,
                          sampleSize = sampleSize, alpha = alpha,
                          nSamples = nSamples, verbose=verbose,
                          originNames = originNames, targetNames = targetNames,
                          nThin = nThin, nBurnin = nBurnin, nChains = nChains,
                          fixedZero = fixedZero, psiPrior = psiPrior,
                          method = "MCMC", returnAllInput = TRUE)
  }
  else {
    results$input <- list(sampleSize = sampleSize, alpha = alpha,
                          originNames = originNames, targetNames = targetNames,
                          method = "MCMC", returnAllInput = FALSE)
  }
  results$BUGSoutput <- out$BUGSoutput
  return(results)
}

estTransitionBoot <- function(originSites = NULL,
                              targetSites = NULL,
                              originPoints = NULL,
                              targetPoints = NULL,
                              originAssignment = NULL,
                              targetAssignment = NULL,
                              originNames = NULL,
                              targetNames = NULL,
                              nBoot = 1000,
                              isGL = FALSE,
                              isTelemetry = FALSE,
                              isRaster = FALSE,
                              isProb = FALSE,
                              captured = "origin",
                              geoBias = NULL,
                              geoVCov = NULL,
                              geoBiasOrigin = geoBias,
                              geoVCovOrigin = geoVCov,
                              targetRaster = NULL,
                              originRaster = NULL,
                              verbose = 0,
                              alpha = 0.05,
                              resampleProjection = 'ESRI:102010',#MigConnectivity::projections$EquidistConic,
                              nSim = ifelse(any(isRaster), 10, 1000),
                              maxTries = 300,
                              dataOverlapSetting = "dummy",
                              fixedZero = NULL,
                              targetRelAbund = NULL,
                              banded = NULL,
                              reencountered = NULL,
                              method = "bootstrap",
                              m = NULL,
                              returnAllInput = TRUE) {
  # Input checking and assignment
  if (any(captured != "origin" & captured != "target" & captured != "neither")){
    stop("captured should be 'origin', 'target', 'neither', or a vector of those options")}
  if (!(verbose %in% 0:3)){
    stop("verbose should be integer 0-3 for level of output during bootstrap: 0 = none, 1 = every 10, 2 = every run, 3 = number of draws")}
  if (length(geoBias)!=2 && any(isGL & (captured == "origin" | captured == "neither"))){
    stop("geoBias should be vector of length 2 (expected bias in longitude and latitude of targetPoints, in resampleProjection units, default meters)")}
  if (!isTRUE(all.equal(dim(geoVCov), c(2, 2), check.attributes = FALSE)) &&
      any(isGL & (captured == "origin" | captured == "neither"))){
    stop("geoVCov should be 2x2 matrix (expected variance/covariance in longitude and latitude of targetPoints, in resampleProjection units, default meters)")}
  if ((is.null(originPoints) && is.null(originRaster) && is.null(originSites)) &&
      is.null(originAssignment) && is.null(banded)){
    stop("Need to define either originAssignment, originSites, originRaster, originPoints, or banded")}
  if ((is.null(targetPoints) && is.null(targetRaster) &&
       is.null(targetSites)) && is.null(targetAssignment) && is.null(reencountered)){
    stop("Need to define either targetAssignment, targetSites, targetRaster, targetPoints, or reencountered")}
  if ((is.null(banded) && !is.null(reencountered) ||
       !is.null(banded)) && is.null(reencountered)){
    stop("Need to define both banded and reencountered")}
  if(inherits(originSites,"SpatialPolygonsDataFrame")){
    originSites <- sf::st_as_sf(originSites)}
  if(inherits(targetSites,"SpatialPolygonsDataFrame")){
    targetSites <- sf::st_as_sf(targetSites)}

  targetStats <- assignRasterStats(targetRaster)
  targetPointsAssigned <- targetStats$PointsAssigned
  targetSingleCell <- targetStats$SingleCell
  targetRasterXYZ <- targetStats$RasterXYZ
  targetRasterXYZcrs <- targetStats$RasterXYZcrs
  targetRaster <- targetStats$Raster

  originStats <- assignRasterStats(originRaster)
  originPointsAssigned <- originStats$PointsAssigned
  originSingleCell <- originStats$SingleCell
  originRasterXYZ <- originStats$RasterXYZ
  originRasterXYZcrs <- originStats$RasterXYZcrs
  originRaster <- originStats$Raster

  if (dataOverlapSetting != "dummy") {
    if (verbose > 0)
      cat("Configuring data overlap settings\n")
    temp <- reassignInds(dataOverlapSetting = dataOverlapSetting,
                         originPoints = originPoints,
                         targetPoints = targetPoints,
                         originAssignment = originAssignment,
                         targetAssignment = targetAssignment,
                         isGL = isGL, isTelemetry = isTelemetry,
                         isRaster = isRaster, isProb = isProb,
                         captured = captured,
                         originRasterXYZ = originRasterXYZ,
                         originSingleCell = originSingleCell,
                         targetRasterXYZ = targetRasterXYZ,
                         targetSingleCell = targetSingleCell,
                         targetSites = targetSites, originSites = originSites)
    originPoints <- temp$originPoints; targetPoints <- temp$targetPoints
    originAssignment <- temp$originAssignment
    targetAssignment <- temp$targetAssignment
    isGL <- temp$isGL; isTelemetry <- temp$isTelemetry
    isRaster <- temp$isRaster; isProb <- temp$isProb
    originRasterXYZ <- temp$originRasterXYZ
    if (!is.null(originRasterXYZ)){
      colnames(originRasterXYZ) <- c("x", "y", paste0("lyr.", 1:length(isRaster)))
      originRaster <- terra::rast(originRasterXYZ, crs = originRasterXYZcrs,
                                  extent = terra::ext(originRaster), type = "xyz")
    }
    originSingleCell <- temp$originSingleCell
    targetRasterXYZ <- temp$targetRasterXYZ
    if (!is.null(targetRasterXYZ)){
      colnames(targetRasterXYZ) <- c("x", "y", paste0("lyr.", 1:length(isRaster)))
      targetRaster <- terra::rast(targetRasterXYZ, crs = targetRasterXYZcrs,
                                  extent = terra::ext(targetRaster), type = "xyz")
    }
    targetSingleCell <- temp$targetSingleCell
  }
  if (any(isProb & (captured != "target")) && (is.null(targetAssignment) || length(dim(targetAssignment))!=2)){
    stop("With probability assignment (isProb==TRUE) animals captured at origin, targetAssignment must be a [number of animals] by [number of target sites] matrix")}
  if (any(isProb & captured != "origin") && (is.null(originAssignment) || length(dim(originAssignment))!=2)){
    stop("With probability assignment (isProb==TRUE) animals captured at target, originAssignment must be a [number of animals] by [number of origin sites] matrix")}

  if (is.null(targetPoints) && is.null(originPoints) &&
      is.null(targetAssignment) && is.null(originAssignment) &&
      is.null(targetRaster) && is.null(originRaster))
    nAnimals <- 0
  else
    nAnimals <- max(nrow(targetPoints), nrow(originPoints), length(isGL),
                  length(isTelemetry), length(isRaster), length(isProb),
                  min(length(targetAssignment), dim(targetAssignment)[1]),
                  min(length(originAssignment), dim(originAssignment)[1]),
                  ifelse(is.null(targetRaster), 0,
                         ifelse(targetPointsAssigned, dim(targetSingleCell)[3],
                         dim(targetRasterXYZ)[2] - 2)),
                  ifelse(is.null(originRaster), 0,
                         ifelse(originPointsAssigned, dim(originSingleCell)[3],
                                dim(originRasterXYZ)[2] - 2)),
                  length(captured))
  nAnimalsTotal <- nAnimals + sum(banded) #+ sum(reencountered)#
  isCMR <- c(rep(FALSE, nAnimals), rep(TRUE, nAnimalsTotal - nAnimals))
  if (length(isGL)==1){
    isGL <- c(rep(isGL, nAnimals), rep(FALSE, nAnimalsTotal - nAnimals))
  }
  else
    isGL <- c(isGL, rep(FALSE, nAnimalsTotal - nAnimals))
  if (length(isTelemetry)==1){
    isTelemetry <- c(rep(isTelemetry, nAnimals),
                     rep(FALSE, nAnimalsTotal - nAnimals))
  }
  else
    isTelemetry <- c(isTelemetry, rep(FALSE, nAnimalsTotal - nAnimals))
  if (length(isRaster)==1){
    isRaster <- c(rep(isRaster, nAnimals), rep(FALSE, nAnimalsTotal - nAnimals))
  }
  else
    isRaster <- c(isRaster, rep(FALSE, nAnimalsTotal - nAnimals))
  if (length(isProb)==1){
    isProb <- c(rep(isProb, nAnimals), rep(FALSE, nAnimalsTotal - nAnimals))
  }
  else
    isProb <- c(isProb, rep(FALSE, nAnimalsTotal - nAnimals))
  if (length(captured)==1){captured <- rep(captured, nAnimals)}

  isCMR <- c(rep(FALSE, nAnimals), rep(TRUE, nAnimalsTotal - nAnimals))
  if (!is.null(banded)) {
    captured <- c(captured, rep("origin", nAnimalsTotal - nAnimals)) #sum(banded)
  }
  if (nAnimals > 0)
    if (any(!isGL[1:nAnimals] & !isTelemetry[1:nAnimals] & !isRaster[1:nAnimals] &
            !isProb[1:nAnimals]))
      stop("For each individual animal (not in banded) one of the following must be set to TRUE:
           isGL, isTelemetry, isRaster, or isProb")
  if (method=="m-out-of-n-bootstrap" && is.null(m))
    m <- ceiling(nAnimalsTotal / 4) # don't know if this is a good default or not!
  else if (method == "bootstrap")
    m <- nAnimalsTotal

  if (any(isRaster & !isProb)) {
    conversion <- rasterToProb(originSites = originSites,
                               targetSites = targetSites,
                               originAssignment = originAssignment,
                               targetAssignment = targetAssignment,
                               isRaster = isRaster, isProb = isProb,
                               captured = captured, targetRaster = targetRaster,
                               originRaster = originRaster)
    if (!conversion$failTarget)
      targetAssignment <- conversion$targetAssignment
    if (!conversion$failOrigin)
      originAssignment <- conversion$originAssignment
    isRaster <- conversion$isRaster
    isProb <- conversion$isProb
  }

  if (!is.null(originPoints) && !is.null(originSites))
    if(!identical(sf::st_crs(originPoints), sf::st_crs(originSites)))
      # project if needed
      originPoints <- sf::st_transform(originPoints, sf::st_crs(originSites))
  if (!is.null(targetPoints) && !is.null(targetSites))
    if(!identical(sf::st_crs(targetPoints), sf::st_crs(targetSites)))
      # project if needed
      targetPoints <- sf::st_transform(targetPoints, sf::st_crs(targetSites))

  # IF originAssignment is NULL - we need to generate originAssignments from
  # the data provided
  if (is.null(originAssignment)){
    if (verbose > 0)
      cat("Creating originAssignment\n")
    # if geolocator, telemetry and captured in origin then simply get the origin site
    if (all(isGL | isTelemetry | captured != "target") && !is.null(originPoints)){
      originAssignment <- suppressMessages(unclass(sf::st_intersects(x = originPoints,
                                                          y = originSites,
                                                          sparse = TRUE)))
    # if raster and not captured in origin sites then determine the origin site
    }
    else if (all(isRaster & captured != "origin")) {
      # if isRaster == TRUE and captured != origin
      # WEIGHTED XY COORDIANTES FROM THE RASTER
      # get geographically weighted median value
      xyOriginRast <- apply(originRasterXYZ[,3:ncol(originRasterXYZ)],
                            MARGIN = 2,
                            FUN = function(x){
                        # xy <-cbind(weighted.mean(originRasterXYZ[,1], w = x, na.rm = TRUE),
                        #           weighted.mean(originRasterXYZ[,2], w = x, na.rm = TRUE))
                  #select the cell with the highest posterior probability #
                              xy <- cbind(originRasterXYZ[which.max(x)[1],1],
                                          originRasterXYZ[which.max(x)[1],2])
                        return(xy)})
      # returns a point estimate for each bird - turn it into a sf object
      xyOriginRast <- t(xyOriginRast)
      colnames(xyOriginRast) <- c("x","y")
      # right now the assignment CRS is WGS84 - should be the same as the origin raster

      #cat("--originRasterXYZcrs -- \n")
      originAssignRast <- sf::st_as_sf(data.frame(xyOriginRast),
                                       coords = c("x","y"),
                                       crs = originRasterXYZcrs)
                                       # crs = sf::st_crs(originRasterXYZcrs))
                                      # crs = 4326)
      # transform to match originSites
      originAssignRast <- sf::st_transform(originAssignRast, sf::st_crs(originSites))
      originAssignment <- suppressMessages(unclass(sf::st_intersects(x = originAssignRast,
                                            y = originSites,
                                            sparse = TRUE)))
    }   # originAssignment <- what # need point assignment for raster (mean location?)
    else if (!is.null(originPoints))
      # originAssignment <- what # points over where we have them, raster assignment otherwise
      originAssignment <- suppressMessages(unclass(sf::st_intersects(x = originPoints,
                                                                 y = originSites,
                                                                sparse = TRUE)))
    else
      originAssignment <- NULL
    if (!is.null(originAssignment)) {
      originAssignment[lengths(originAssignment)==0] <- NA
      if (any(lengths(originAssignment)>1)){
        warning("Overlapping originSites may cause issues\n")
        originAssignment <- lapply(originAssignment, function (x) x[1])
      }
      originAssignment <- array(unlist(originAssignment))
      duds <- is.na(originAssignment) & captured[1:nAnimals] == "origin"
      if (any(duds)){
        if (verbose > 0)
          cat("Not all origin capture locations are within originSites. Assigning to closest site\n")
        warning("Not all origin capture locations are within originSites. Assigning to closest site.\n",
                "Affects animals: ", paste(which(duds), collapse = ","))
        originAssignment[duds] <-
          sf::st_nearest_feature(x = originPoints[duds,],
                                 y = originSites)

      }
    }
    if (!is.null(reencountered)) {
      originAssignment <- array(c(originAssignment,
                                  rep(1:length(banded), banded)))
    }
  }
  else if (!is.null(reencountered)) {
    if (is.array(originAssignment) && length(dim(originAssignment))>1){
      originAssignment <- rbind(originAssignment,
                                array(0, c(nAnimalsTotal - nAnimals,
                                           dim(originAssignment)[2])))
      place <- nAnimals
      for (i in 1:length(banded)) {
        originAssignment[place + 1:banded[i], i] <- 1
        place <- place + banded[i]
      }
    }
    else {
      nOriginSites <- length(banded)
      originAssignment <- array(c(originAssignment,
                                rep(1:nOriginSites, banded)))
    }
  }

  if (is.null(targetAssignment)){
    if (verbose > 0)
      cat("Creating targetAssignment\n")
    if (all(isGL | isTelemetry | captured != "origin")) {
      targetAssignment <- suppressMessages(unclass(sf::st_intersects(x = targetPoints,
                                                    y = targetSites,
                                                    sparse = TRUE)))
    }
    else if (all(isRaster & captured != "target")){
      #targetAssignment <- what # need point assignment for raster (mean location?)
      xyTargetRast <- apply(targetRasterXYZ[,3:ncol(targetRasterXYZ)],
                          MARGIN = 2,
                          FUN = function(x){
                            #xy <-cbind(Hmisc::wtd.quantile(targetRasterXYZ[,1],probs = 0.5, weight = x, na.rm = TRUE),
                            #           Hmisc::wtd.quantile(targetRasterXYZ[,2],probs = 0.5, weight = x, na.rm = TRUE))
                            # xy <- cbind(weighted.mean(targetRasterXYZ[,1], w = x, na.rm = TRUE),
                            #            weighted.mean(targetRasterXYZ[,2], w = x, na.rm = TRUE))
                  # Select cell with the maximum posterior probability #
                            xy <- cbind(targetRasterXYZ[which.max(x)[1],1],
                                        targetRasterXYZ[which.max(x)[1],2])
                            return(xy)})
      # returns a point estimate for each bird - turn it into a sf object
      xyTargetRast <- t(xyTargetRast)
      colnames(xyTargetRast) <- c("x","y")
      # right now the assignment CRS is WGS84 - should be the same as the origin raster
      targetAssignRast <- sf::st_as_sf(data.frame(xyTargetRast),
                                       coords = c("x","y"),
                                       crs = targetRasterXYZcrs)
                                      # crs = sf::st_crs(targetRasterXYZcrs))
                                       #crs = 4326)
      # transform to match originSites
      targetSites_wgs <- sf::st_transform(targetSites, 4326)
      #targetAssignRast <- sf::st_transform(targetAssignRast, sf::st_crs(targetSites))
      #targetAssignment <- suppressMessages(array(unlist(unclass(sf::st_intersects(x = targetAssignRast,
      #                                                           y = targetSites_wgs,
      #                                                           sparse = TRUE)))))
      # Check which points are in target sites
      targetAssignment <- suppressMessages(unclass(sf::st_intersects(x = targetAssignRast,
                                                      y = targetSites_wgs,
                                                      sparse = TRUE)))

    # NEED TO ADD A CATCH HERE TO ASSIGN THE MAX_prob to CLOSEST TARGET REGION
    # IF INTERSECTS IS (empty)
    }
    else if (!is.null(targetPoints))
   #   targetAssignment <- what # points over where we have them, raster assignment otherwise
      targetAssignment <- suppressMessages(unclass(sf::st_intersects(x = targetPoints,
                                                        y = targetSites,
                                                        sparse = TRUE)))
    else
      targetAssignment <- NULL
   if (!is.null(targetAssignment)) {
     targetAssignment[lengths(targetAssignment)==0] <- NA
     if (any(lengths(targetAssignment)>1)){
       warning("Overlapping targetSites may cause issues\n")
       targetAssignment <- lapply(targetAssignment, function(x) x[1])
     }
     targetAssignment <- array(unlist(targetAssignment))
     duds <- is.na(targetAssignment) & captured[1:nAnimals] == "target"
     if (any(duds)){
       if (verbose > 0)
         cat("Not all target capture locations are within targetSites. Assigning to closest site\n")
       warning("Not all target capture locations are within targetSites. Assigning to closest site.\n",
               "Affects animals: ", paste(which(duds), collapse = ","))
       targetAssignment[duds] <-
         sf::st_nearest_feature(x = targetPoints[duds,],
                                y = targetSites)

     }
   }
   if (!is.null(reencountered)) {
     nTargetSites <- dim(reencountered)[2]
     nOriginSites <- dim(reencountered)[1]
     for (j in 1:nOriginSites)
       targetAssignment <- array(c(targetAssignment,
                                   rep(c(1:nTargetSites, NA),
                                       c(reencountered[j, ],
                                         banded[j] - sum(reencountered[j, ])))))
   }
  }
  else if (!is.null(reencountered)) {
    if (is.array(targetAssignment) && length(dim(targetAssignment))>1){
      nTargetSites <- dim(reencountered)[2]
      nOriginSites <- dim(reencountered)[1]
      targetAssignment <- rbind(targetAssignment,
                                array(0,
                                      c(nAnimalsTotal - nAnimals, nTargetSites)))
      place <- nAnimals
      for (j in 1:nOriginSites) {
        for (i in 1:nTargetSites) {
          targetAssignment[place + 1:reencountered[j, i], i] <- 1
          place <- place + reencountered[j, i]
        }
        targetAssignment[place + 1:(banded[j] - sum(reencountered[j, ])), ] <- NA
        place <- place + banded[j] - sum(reencountered[j, ])
      }
    }
    else {
      nTargetSites <- dim(reencountered)[2]
      nOriginSites <- dim(reencountered)[1]
      for (j in 1:nOriginSites)
        targetAssignment <- array(c(targetAssignment,
                                    rep(c(1:nTargetSites, NA),
                                        c(reencountered[j, ],
                                          banded[j] - sum(reencountered[j, ])))))
    }
  }

  if(is.null(originSites)){
    if (is.array(originAssignment) && length(dim(originAssignment))>1){
      nOriginSites <- ncol(originAssignment)
    }
    else {
      nOriginSites <- length(unique(originAssignment))
    }
  }
  else{
    nOriginSites <- nrow(originSites)
  }

  if(is.null(targetSites)){
    if (is.array(targetAssignment) && length(dim(targetAssignment))>1){
      nTargetSites <- ncol(targetAssignment)
    }
    else {
      nTargetSites <- max(length(unique(targetAssignment[!is.na(targetAssignment)])),
                          length(targetNames),
                          dim(reencountered)[2])
    }
  }
  else {
    nTargetSites <- nrow(targetSites)
  }
  # if (length(targetPoints)!=nAnimals &&
  #     dim(targetAssignment)[length(dim(targetAssignment))]!=nAnimals ||
  #     nrow(originAssignment)!=nAnimals)
  #   stop("isGL should be the same length as originAssignment/originPoints and targetPoints/targetAssignment (number of animals)")
  # if (any(is.na(originAssignment)))
  #   stop("NAs in origin sites (make sure all points fall within polygons)")
  if(!is.null(originPoints) && is.na(sf::st_crs(originPoints)))
    stop('Coordinate system definition needed for originPoints')
  if(!is.null(originSites) && is.na(sf::st_crs(originSites)))
    stop('Coordinate system definition needed for originSites')
  if(!is.null(targetPoints) && is.na(sf::st_crs(targetPoints))){
    stop('Coordinate system definition needed for targetPoints')
  }
  if(!is.null(targetSites) && is.na(sf::st_crs(targetSites))){
    stop('Coordinate system definition needed for targetSites')
  }
  if(!is.null(targetPoints)){
    targetPoints <- sf::st_transform(targetPoints, crs = resampleProjection)}
  if(!is.null(targetSites)){
    targetSites <- sf::st_transform(targetSites, crs = resampleProjection)}
  if(!is.null(originPoints)){
    originPoints <- sf::st_transform(originPoints, crs = resampleProjection)}
  if(!is.null(originSites)){
    originSites <- sf::st_transform(originSites, crs = resampleProjection)}
  # Names - this isn't going to work as currently written for sf objects #
  if (is.null(targetNames)){
    # if (is.null(targetSites[[1]]) || anyDuplicated(targetSites[[1]])){
    #   targetNames <- as.character(1:nTargetSites)}else{
    targetNames <- as.character(1:nTargetSites)
  }

  if (is.null(originNames)){
#    if (is.null(originSites[[1]]) || anyDuplicated(originSites[[1]])){
    if (nOriginSites > 26){
      originNames <- as.character(1:nOriginSites)
    }else{
      originNames <- LETTERS[1:nOriginSites]}
    # }else{
    #   originNames <- 1:nOriginSites
    # }
  }

  targetPointsInSites <- FALSE

  if (targetPointsAssigned && !is.null(targetSites) && any(isRaster)) {
    if (verbose > 0){
      cat('Checking if single cell target points in targetSites, may take a moment\n')}
    targetPointSample2 <- apply(targetSingleCell,
                                FUN = function(x){sf::st_as_sf(data.frame(x),
                                                               coords = c("Longitude", "Latitude"),
                                                               crs = 4326)},
                                MARGIN = 3)
    if(!sf::st_crs(targetSites)==sf::st_crs(targetPointSample2[[1]])){
      targetPointSample2 <- sapply(targetPointSample2, sf::st_transform, crs = resampleProjection)
    }

    targetCon <- sapply(targetPointSample2, FUN = function(z){
      suppressMessages(as.numeric(unclass(sf::st_intersects(x = z, y = targetSites,
                                           sparse = TRUE))))})

    if (!any(is.na(targetCon)))
      targetPointsInSites <- TRUE
    else if (verbose > 0)
      cat('Single cell target points supplied, but some points (proportion',
          format(sum(is.na(targetCon))/length(targetCon), digits = 2), ') not in targetSites\n')
  }else {targetCon <- NULL}

  originPointsInSites <- FALSE
  if (originPointsAssigned && !is.null(originSites) && any(isRaster)) {
    if (verbose > 0)
      cat('Checking if single cell origin points in originSites, may take a moment\n')
    nSamples <- dim(originSingleCell)[1]
    originPointSample2 <- apply(originSingleCell,
                                FUN = function(x){sf::st_as_sf(data.frame(x),
                                                               coords = c("Longitude", "Latitude"),
                                                               crs = 4326)},
                                MARGIN = 3)
    if(!sf::st_crs(originSites)==sf::st_crs(originPointSample2[[1]])){
      originPointSample2 <- sapply(originPointSample2, sf::st_transform,
                                   crs = resampleProjection)
    }

    originCon <- sapply(originPointSample2, FUN = function(z){
      suppressMessages(as.numeric(unclass(sf::st_intersects(x = z, y = originSites,
                                           sparse = TRUE))))})
    if (!any(is.na(originCon)))
      originPointsInSites <- TRUE
    else if (verbose > 0){
      cat('Single cell origin points supplied, but some points (proportion',
          sum(is.na(originCon))/length(originCon), ') not in originSites\n')
    }
  }else{
    originCon <- NULL
  }

  # if (!is.null(targetPoints) && nAnimals < nAnimalsTotal) {
  #   if (verbose > 0)
  #     cat("Filling in dummy target points for CMR data")
  #   dummyPoint <- targetPoints[1, ]
  #   for (i in 1:(nAnimalsTotal - nAnimals)) {
  #     targetPoints <- rbind(targetPoints,
  #                           dummyPoint)
  #   }
  # }
  # if (!is.null(originPoints) && nAnimals < nAnimalsTotal) {
  #   if (verbose > 0)
  #     cat("Filling in dummy origin points for CMR data")
  #   dummyPoint <- originPoints[1, ]
  #   for (i in 1:(nAnimalsTotal - nAnimals)) {
  #     originPoints <- rbind(originPoints,
  #                           dummyPoint)
  #   }
  # }

  if (!is.null(targetRelAbund) && any(captured=="target")) {
    if (coda::is.mcmc(targetRelAbund) || coda::is.mcmc.list(targetRelAbund)) {
      targetRelAbund <- as.matrix(targetRelAbund)
    }
    if (is.matrix(targetRelAbund) && dim(targetRelAbund)[1]>1) {
      abundFixed <- FALSE
      if (dim(targetRelAbund)[2]>nTargetSites)
        abundParams <- paste('relN[', 1:nTargetSites, ']', sep='')
      else if (dim(targetRelAbund)[2]==nTargetSites)
        abundParams <- 1:nTargetSites
      else
        stop('Number of target sites must be constant between sites and abundance')
      if (dim(targetRelAbund)[1] >= nBoot)
        abundRows <- round(seq(from = 1, to = dim(targetRelAbund)[1],
                                      length.out = nBoot))
      else
        stop("You need at least nSamples rows to targetRelAbund")
      targetRelAbund <- as.matrix(targetRelAbund[abundRows, abundParams])
      abundBase <- colMeans(targetRelAbund)
    }
    else {
      abundFixed <- TRUE
      if (length(targetRelAbund)!=nTargetSites)
        stop('Number of target sites must be constant between sites and abundance')
      abundBase <- targetRelAbund
      targetRelAbund <- matrix(targetRelAbund, nBoot, nTargetSites, TRUE)
    }
    weights <- array(0, c(nBoot, nAnimalsTotal))
    if (length(dim(targetAssignment))==2) {
      ta <- apply(targetAssignment, 1, which.max)
      if (is.list(ta)) {
        ta[lengths(ta)==0] <- NA
        ta <- unlist(ta)
      }
    }
    else
      ta <- targetAssignment
    nOriginAnimals <- sum(captured != "target")
    nTargetAnimals <- rep(NA, nTargetSites)
    for (i in 1:nTargetSites) {
      nTargetAnimals[i] <- sum(captured=="target" & ta==i)
      if (nTargetAnimals[i] > 0)
        weights[ , captured=="target" & ta==i] <- targetRelAbund[ , i] /
          nTargetAnimals[i] * (nAnimalsTotal - nOriginAnimals) / nAnimalsTotal
    }
    if (nOriginAnimals>0) {
      if (all(nTargetAnimals>0))
        weights[ , captured!="target"] <- 1/nAnimalsTotal
      else {
        t0 <- which(nTargetAnimals>0)
        weights[ , captured!="target"] <- 1/nAnimalsTotal +
          rowSums(targetRelAbund[ , t0]) *
          sum(nTargetAnimals) / nAnimalsTotal / nOriginAnimals
        # weights[captured!="target" & ta %in% t0] <- (1 - sum(weights)) /
        #   sum(captured!="target" & ta %in% t0)
        if (sum(captured!="target" & ta %in% t0) == 0)
          warning("Not all target sites have likely data. Estimates will be biased.")
      }
    }
    else if (any(nTargetAnimals==0)) {
      warning("Not all target sites have data. Estimates will be biased.")
    }
  }
  else {
    weights <- NULL
    if (any(captured!="origin"))
      warning("Unless target site data were collected in proportion to abundance ",
              "at those sites, estimates will be biased. We recommend including ",
              "an estimate of target site relative abundance with the ",
              "targetRelAbund argument, to allow resampling in proportion to ",
              "abundance.")
  }

  sites.array <- psi.array <- array(0, c(nBoot, nOriginSites, nTargetSites),
                                    dimnames = list(1:nBoot, originNames,
                                                    targetNames))
  if (!is.null(banded))
    r.array <- array(0, c(nBoot, nTargetSites),
                     dimnames = list(1:nBoot, targetNames))
  else
    r.array <- NULL
  if (is.null(dim(originAssignment))){
    originAssignment <- array(originAssignment)}
  if (is.null(dim(targetAssignment))){
    targetAssignment <- array(targetAssignment)}

  if (is.null(fixedZero)) {
    nFixedZero <- 0
  }
  else {
    nFixedZero <- nrow(fixedZero)
  }
  countFailed <- 0
  failed <- FALSE
  if (length(dim(originAssignment))==2){
    pointOriginAssignment <- apply(originAssignment, 1, which.max)
  }
  else{
    pointOriginAssignment <- as.vector(originAssignment)
  }
  if (length(pointOriginAssignment) > nAnimals){
    if (nAnimals==0)
      pointOriginAssignment <- NULL
    else
      pointOriginAssignment <- pointOriginAssignment[1:nAnimals]
  }
  if (length(dim(targetAssignment))==2){
    pointTargetAssignment <- apply(targetAssignment, 1, which.max)
  }
  else{
    pointTargetAssignment <- as.vector(targetAssignment)
  }
  if (length(pointTargetAssignment) > nAnimals) {
    if (nAnimals==0)
      pointTargetAssignment <- NULL
    else
      pointTargetAssignment <- pointTargetAssignment[1:nAnimals]
  }
  if (length(pointTargetAssignment) == length(pointOriginAssignment)) {
   #pointSites <- table(pointOriginAssignment, pointTargetAssignment)
   psi_r <- calcTransition(banded, reencountered,
                           originAssignment = pointOriginAssignment,
                           targetAssignment = pointTargetAssignment,
                           originNames = originNames,
                           targetNames = targetNames,
                           method = "BFGS")
   pointPsi <- psi_r$psi
   point_r <- psi_r$r
  }
  else {
    pointPsi <- NULL
    point_r <- NULL
  }
  boot <- 1
  if (verbose > 0)
    cat("Starting bootstrap\n")
  while (boot <= nBoot) {
    if (verbose > 1 || verbose == 1 && boot %% 100 == 0)
      cat("Bootstrap Run", boot, "of", nBoot, "at", date(), "\n")
    # Make sure have animals from every origin site
    origin.sample <- c() # Start with zero origin sites
    while (length(unique(origin.sample)) < nOriginSites) { #2
      # Sample individual animals with replacement
      animal.sample <- sample.int(m, replace=TRUE, prob = weights[boot,])
      if (any(captured[animal.sample]!='origin')) {
        if (length(dim(originAssignment))==2)
          assignment <- originAssignment[animal.sample, , drop = FALSE]
        else
          assignment <- originAssignment[animal.sample, drop = FALSE]
        oSamp <- locSample(isGL = (isGL[animal.sample] & captured[animal.sample]!='origin'),
                           isRaster = (isRaster[animal.sample] & captured[animal.sample]!='origin'),
                           isProb = (isProb[animal.sample] & captured[animal.sample]!='origin'),
                           isTelemetry = (isTelemetry[animal.sample] |
                                            isCMR[animal.sample] |
                                            captured[animal.sample]=='origin'),
                           geoBias = geoBiasOrigin,
                           geoVCov = geoVCovOrigin,
                           points = originPoints[animal.sample, ],
                           matvals = originRasterXYZ[, c(1:2, animal.sample + 2)],
                           matvals_crs = originRasterXYZcrs,
                           singleCell = originSingleCell[,,animal.sample],
                           overlap1 = originCon[,animal.sample],
                           pointsInSites = originPointsInSites,
                           assignment = assignment,
                           sites = originSites,
                           resampleProjection = resampleProjection,
                           nSim = nSim,
                           maxTries = maxTries)
        if (!is.null(oSamp$notfind)) {
          oSamp$notfind$Animal <- animal.sample[oSamp$notfind$Animal]
          notfind <- unique(oSamp$notfind)
          stop('maxTries (',maxTries,') reached during origin location sampling, exiting. ',
               'Animal(s) where location sampling failed to fall in sites:\n',
               paste(utils::capture.output(print(notfind, row.names = FALSE)), collapse = "\n"),
               '\nExamine originSites',
               ifelse(any(notfind$isGL),
                      ', geoBiasOrigin, geoVcovOrigin, originPoints', ''),
               ifelse(any(notfind$isRaster), ', originRaster', ''),
               ifelse(any(notfind$isTelemetry), ', originPoints/captured', ''),
               ', and resampleProjection to determine why sampled points fell outside sites.')
        }
        origin.sample <- oSamp$site.sample
        if (verbose > 2)
          cat(' ', oSamp$draws, 'origin draw(s) (of length', nSim, 'and of', maxTries, 'possible).\n')
      }
      else {
        # Get origin population for each animal sampled
        if (length(dim(originAssignment))==2){
          origin.sample <- apply(originAssignment[animal.sample, ], 1, which.max)
          if (is.list(origin.sample)) {
            origin.sample[lengths(origin.sample)==0] <- NA
            origin.sample <- unlist(origin.sample)
          }
        }
        else
          origin.sample <- originAssignment[animal.sample]
      }
    }
    if (any(captured[animal.sample]!="target")) {
      if (length(dim(targetAssignment))==2)
        assignment <- targetAssignment[animal.sample, , drop = FALSE]
      else
        assignment <- targetAssignment[animal.sample, drop = FALSE]
      tSamp <- locSample(isGL = (isGL[animal.sample] & captured[animal.sample] != "target"),
                         isRaster = (isRaster[animal.sample] & captured[animal.sample] != "target"),
                         isProb = (isProb[animal.sample] & captured[animal.sample] != "target"),
                         isTelemetry = (isTelemetry[animal.sample] |
                                          isCMR[animal.sample] |
                                          captured[animal.sample] == "target"),
                         geoBias = geoBias, geoVCov = geoVCov,
                         points = targetPoints[animal.sample, ],
                         matvals = targetRasterXYZ[, c(1:2, animal.sample + 2)],
                         matvals_crs = targetRasterXYZcrs,
                         singleCell = targetSingleCell[,,animal.sample],
                         pointsInSites = targetPointsInSites,
                         overlap1 = targetCon[, animal.sample],
                         sites = targetSites,
                         assignment = assignment,
                         resampleProjection = resampleProjection, nSim = nSim,
                         maxTries = maxTries)
      if (!is.null(tSamp$notfind)) {
        tSamp$notfind$Animal <- animal.sample[tSamp$notfind$Animal]
        notfind <- unique(tSamp$notfind)
        stop('maxTries (',maxTries,') reached during target location sampling, exiting. ',
             'Animal(s) where location sampling failed to fall in sites:\n',
             paste(utils::capture.output(print(notfind, row.names = FALSE)), collapse = "\n"),
             '\nExamine targetSites',
             ifelse(any(notfind$isGL),
                    ', geoBiasOrigin, geoVcovOrigin, targetPoints', ''),
             ifelse(any(notfind$isRaster), ', targetRaster', ''),
             ifelse(any(notfind$isTelemetry), ', targetPoints/captured', ''),
             ', and resampleProjection to determine why sampled points fell outside sites.')
      }
      target.sample <- tSamp$site.sample
      target.sample[target.sample==0] <- NA
      #target.point.sample <- tSamp$target.point.sample
      if (verbose > 2)
        cat(' ', tSamp$draws, 'target draw(s) (of length', nSim, 'and of', maxTries, 'possible).\n')
    }
    else {
      # Get target population for each animal sampled
      if (length(dim(targetAssignment))==2){
        target.sample <- apply(targetAssignment[animal.sample, ], 1, which.max)
        if (is.list(target.sample)) {
          target.sample[lengths(target.sample)==0] <- NA
          target.sample <- unlist(target.sample)
        }
      }
      else
        target.sample <- targetAssignment[animal.sample]
    }
    # Now that we have breeding and non-breeding site for point, add to transition count matrix
    sites <- table(origin.sample,
                   target.sample,
                   useNA = "no")
    sites.array[boot, as.integer(rownames(sites)),
                as.integer(colnames(sites))] <- sites
    if (nFixedZero > 0) {
      for (i in 1:nFixedZero) {
        if (sites.array[boot, fixedZero[i, 1], fixedZero[i, 2]] > 0) {
          failed <- TRUE
          countFailed <- countFailed + 1
          if (countFailed > nBoot * 100)
            stop("estTransition stopped because getting very high number of bootstrap runs:\n",
                 countFailed, " where animals use transitions fixed to zero.\n",
                 "You should examine fixedZero and the data to make sure those ",
                 "transition probabilities are really zero")
          sites.array[boot,,] <- 0
          break
        }
      }
      if (failed) {
        failed <- FALSE
        next
      }
    }
    if (is.null(banded)) {
      banded.sample <- NULL
      reencountered.sample <- NULL
    }
    else {
      reencountered.sample <- table(factor(origin.sample[isCMR[animal.sample]],
                                           levels = 1:nOriginSites),
                                    factor(target.sample[isCMR[animal.sample]],
                                           levels = 1:nTargetSites),
                                    useNA = "no")
      banded.sample <- table(factor(origin.sample[isCMR[animal.sample]],
                                    levels = 1:nOriginSites))
    }
    # Use new function that allows for CMR data
    psi_r <- calcTransition(banded.sample, reencountered.sample,
                            originAssignment = origin.sample[!isCMR[animal.sample]],
                            targetAssignment = target.sample[!isCMR[animal.sample]],
                            originNames = originNames,
                            targetNames = targetNames,
                            method = "BFGS")
    if (any(is.na(psi_r$psi)) || any(psi_r$psi < 0) || any(psi_r$psi > 1)) {
      if (verbose > 2)
        cat(' Bootstrap estimate producing nonsense psi; drawing again\n')
      next
    }
    psi.array[boot, , ] <- psi_r$psi
    if (!is.null(banded))
      r.array[boot, ] <- psi_r$r
    boot <- boot + 1
  }
  if (countFailed > 0)
    warning(countFailed, " bootstrap ",ifelse(countFailed>1, "runs", "run"),
            " failed due to animals using transitions fixed at zero with ",
            nBoot, " successful. If this ratio is high, you should examine ",
            "fixedZero and the data to make sure those transition ",
            "probabilities really are zero\n")
  if (method=="bootstrap") {
    meanPsi <- apply(psi.array, 2:3, mean)
    medianPsi <- apply(psi.array, 2:3, median)
    sePsi <- apply(psi.array, 2:3, sd)
    simpleCIPsi <- apply(psi.array, 2:3, quantile, probs = c(alpha/2, 1-alpha/2),
                         na.rm=TRUE, names = FALSE)
    psi.matrix <- array(c(psi.array), c(nBoot, nOriginSites * nTargetSites),
                        list(NULL, paste(rep(originNames, nTargetSites),
                                         rep(targetNames, each = nOriginSites),
                                         sep = "#")))
    psi.mcmc <- coda::as.mcmc(psi.matrix)
    hpdCI <- coda::HPDinterval(psi.mcmc, 1-alpha)
    hpdCI <- array(hpdCI, c(nOriginSites, nTargetSites, 2),
                   list(originNames, targetNames, c("lower", "upper")))
    hpdCI <- aperm(hpdCI, c(3, 1, 2))
    bcCIPsi <- array(NA, dim = c(2, nOriginSites, nTargetSites),
                     dimnames = list(NULL, originNames, targetNames))
    for (i in 1:nOriginSites) {
      for (j in 1:nTargetSites) {
        psi.z0 <- qnorm(sum(psi.array[, i, j] < meanPsi[i, j], na.rm = TRUE) /
                          length(which(!is.na(psi.array[, i, j]))))
        bcCIPsi[ , i, j] <- quantile(psi.array[, i, j],
                                     pnorm(2 * psi.z0 + qnorm(c(alpha/2, 1-alpha/2))),
                                     na.rm=TRUE, names = FALSE)
      }
    }
    if (!is.null(r.array)){
      mean.r <- apply(r.array, 2, mean)
      median.r <- apply(r.array, 2, median)
      se.r <- apply(r.array, 2, sd)
      simpleCIr <- apply(r.array, 2, quantile, probs = c(alpha/2, 1-alpha/2),
                         na.rm=TRUE, names = FALSE)
      r.mcmc <- coda::as.mcmc(r.array)
      hpdCIr <- coda::HPDinterval(r.mcmc, 1-alpha)
      hpdCIr <- aperm(hpdCIr, c(2, 1))
      bcCIr <- array(NA, dim = c(2, nTargetSites),
                     dimnames = list(c("lower", "upper"), targetNames))
      for (j in 1:nTargetSites) {
        r.z0 <- qnorm(sum(r.array[ ,j] < mean.r[j], na.rm = TRUE) /
                          length(which(!is.na(r.array[ ,j]))))
        bcCIr[ , j] <- quantile(r.array[ ,j],
                                pnorm(2 * r.z0 + qnorm(c(alpha/2, 1-alpha/2))),
                                na.rm=TRUE, names = FALSE)
      }
      r <- list(sample = r.array, mean = mean.r,
                se = se.r, simpleCI = simpleCIr,
                bcCI = bcCIr, hpdCI = hpdCIr, median = median.r)
    }
    else
      r <- NULL
  }
  else {
    meanPsi <- apply(psi.array, 2:3, mean) * m / nAnimals
    medianPsi <- apply(psi.array, 2:3, median)
    sePsi <- apply(psi.array, 2:3, sd) * sqrt(m / nAnimals)
    simpleCIPsi <- apply(psi.array, 2:3, quantile, probs = c(alpha/2, 1-alpha/2),
                         na.rm=TRUE, names = FALSE)
  }
  if (returnAllInput) {
    input <-list(sampleSize = nAnimals, originSites = originSites,
                 targetSites = targetSites,
                 originPoints = originPoints,
                 targetPoints = targetPoints,
                 originAssignment = originAssignment,
                 targetAssignment = targetAssignment,
                 originNames = originNames,
                 targetNames = targetNames,
                 nSamples = nBoot,
                 isGL=isGL, isTelemetry = isTelemetry,
                 isRaster = isRaster, isProb = isProb,
                 captured = captured,
                 geoBias=geoBias, geoVCov=geoVCov,
                 geoBiasOrigin = geoBiasOrigin,
                 geoVCovOrigin = geoVCovOrigin,
                 targetRaster = targetRaster,
                 originRaster = originRaster,
                 verbose = verbose,
                 alpha = alpha,
                 resampleProjection = resampleProjection,
                 nSim = nSim, maxTries = maxTries,
                 dataOverlapSetting = dataOverlapSetting,
                 fixedZero = fixedZero,
                 targetRelAbund = targetRelAbund,
                 banded = banded,
                 reencountered = reencountered,
                 method = method,
                 m = m,
                 returnAllInput = TRUE)
  }
  else {
    input <-list(sampleSize = nAnimals,
                 originNames = originNames,
                 targetNames = targetNames,
                 alpha = alpha,
                 method = method,
                 returnAllInput = FALSE)
  }
  return(list(psi = list(sample = psi.array, mean = meanPsi, se = sePsi,
                         simpleCI = simpleCIPsi, bcCI = bcCIPsi, hpdCI = hpdCI,
                         median = medianPsi, point = pointPsi),
              r = r,
              input = input,
              BUGSoutput = NULL))
}

#' Estimate psi (transition probabilities between locations in two phases of
#' the annual cycle)
#'
#' Estimation and resampling of uncertainty for psi (transition probabilities
#' between origin sites in one phase of the annual cycle and target sites in
#' another for migratory animals). Data can be from any combination of
#' geolocators (GL), telemetry/GPS, intrinsic markers such as isotopes and
#' genetics, and band/ring reencounter data.
#'
#' @param originSites A polygon spatial layer (sf - MULTIPOLYGON)
#'  defining the geographic representation of sites in the
#'  origin season.
#' @param targetSites A polygon spatial layer (sf - MULTIPOLYGON)
#'  defining the geographic representation of sites in the
#'  target season.
#' @param originPoints A \code{sf} or \code{SpatialPoints} object, with number
#'  of rows or length being the number of animals tracked. Each point indicates
#'  the origin location of an animal (or point estimate of same, for GL animals
#'  released on target sites). Note that to simplify input of multiple
#'  data types both between and for the same animal, if origin points are
#'  provided for any animal, they must be provided for all except banding data
#'  (can be dummy values), unless \code{dataOverlapSetting} is set to "none".
#' @param targetPoints For GL or telemetry data, a \code{sf} or
#'  \code{SpatialPoints} object, with length or number of rows number of animals
#'  tracked. Each point indicates the point estimate location of an animal in
#'  the target season. Note that to simplify input of multiple
#'  data types both between and for the same animal, if target points are
#'  provided for any animal, they must be provided for all except banding data
#'  (can be dummy values), unless \code{dataOverlapSetting} is set to "none".
#' @param originAssignment Assignment of animals to origin season sites. Either
#'  an integer vector with length number of animals tracked or a matrix of
#'  probabilities with number of animals tracked rows and number of origin sites
#'  columns (and rows summing to 1). The latter only applies to animals released
#'  in the target sites where there is uncertainty about their origin site, for
#'  example from genetic population estimates from the rubias package.
#'  Optional, but some combination of these inputs should be defined. Note that
#'  if \code{originAssignment} is a probability table, animals with known origin
#'  sites can have 1 in that column and 0s in all others. Also note that if
#'  \code{method} is "MCMC", anything in \code{originAssignment} and
#'  \code{targetAssignment} will be assumed to represent animals tracked via
#'  telemetry, with known origin and target sites.
#' @param targetAssignment Assignment of animals to target season sites. Either
#'  an integer vector with length number of animals tracked or a matrix of
#'  probabilities with number of animals tracked rows and number of target sites
#'  columns (and rows summing to 1). The latter only applies to animals released
#'  in the origin sites where there is uncertainty about their target site, for
#'  example from genetic population estimates from the rubias package.
#'  Optional, but some combination of these inputs needs to be defined. Note
#'  that if \code{targetAssignment} is a probability table, animals with known
#'  target sites can have 1 in that column and 0s in all others.
#' @param originNames Optional, but recommended to keep track. Vector of names
#'  for the origin sites. If not provided, the function will either try to get
#'  these from another input or provide default names (capital letters).
#' @param targetNames Optional, but recommended to keep track. Vector of names
#'  for the target sites. If not provided, the function will either try to get
#'  these from another input or provide default names (numbers).
#' @param nSamples Number of post-burn-in MCMC samples to store (\code{method}
#'  == "MCMC") OR number of bootstrap runs for \code{method}
#'  == "bootstrap". In the latter case, animals are sampled with replacement
#'  for each. For all, the purpose is to estimate sampling uncertainty.
#' @param isGL Indicates whether or which animals were tracked with geolocators.
#'  Should be either single TRUE or FALSE value, or vector with length of
#'  number of animals tracked, with TRUE or FALSE for each animal in data
#'  (except those in \code{banded}, which are handled separately). For
#'  TRUE animals, the model applies \code{geoBias} and \code{geoVCov} to
#'  \code{targetPoints} where \code{captured} == "origin" or "neither" and
#'  \code{geoBiasOrigin} and \code{geoVCovOrigin} to
#'  \code{originPoints} where \code{captured} == "target" or "neither".
#'  Geolocator data should be entered as \code{originPoints} and
#'  \code{targetPoints}.
#' @param isTelemetry Indicates whether or which animals were tracked with
#'  telemetry/GPS (no location uncertainty on either end).
#'  Should be either single TRUE or FALSE value, or vector with length of
#'  number of animals tracked, with TRUE or FALSE for each animal in data
#'  (except those in \code{banded}, which are handled separately).
#'  Telemetry data can be entered as points or using the \code{targetAssignment}
#'  and \code{originAssignment} arguments.
#' @param isRaster Indicates whether or which animals were tracked with
#'  intrinsic markers (e.g., genetics or isotopes), with location uncertainty
#'  expressed as a raster of probabilities by grid cells, either in
#'  \code{targetRaster} or \code{originRaster}. Should be either single TRUE or
#'  FALSE value, or vector with length of number of animals tracked, with TRUE
#'  or FALSE for each animal in data (except those in \code{banded}, which are
#'  handled separately).
#' @param isProb Indicates whether or which animals were tracked with
#'  intrinsic markers (e.g., genetics or isotopes), with location uncertainty
#'  expressed as a probability table, either in \code{targetAssignment} or
#'  \code{originAssignment}. Should be either single TRUE or FALSE value, or
#'  vector with length of number of animals tracked, with TRUE or FALSE for each
#'  animal in data (except those in \code{banded}, which are handled separately).
#' @param captured Indicates whether or which animals were captured in the
#'  origin sites, the target sites, or neither (another phase of the annual
#'  cycle). Location uncertainty will only be applied where the animal was not
#'  captured. So this doesn't matter for telemetry data, and is assumed to be
#'  "origin" for band return data. Should be either single "origin" (default),
#'  "target", or "neither" value, or a character vector with length of number of
#'  animals tracked, with "origin", "target", or "neither" for each animal.
#' @param geoBias For GL data, vector of length 2 indicating expected bias
#'  in longitude and latitude of \code{targetPoints}, in
#'  \code{resampleProjection} units (default meters).
#' @param geoVCov For GL data, 2x2 matrix with expected variance/covariance
#'    in longitude and latitude of \code{targetPoints}, in
#'    \code{resampleProjection} units (default meters).
#' @param geoBiasOrigin For GL data where \code{captured}!="origin", vector of
#'  length 2 indicating expected bias in longitude and latitude of
#'  \code{originPoints}, in \code{resampleProjection} units (default meters).
#' @param geoVCovOrigin For GL data where \code{captured}!="origin", 2x2 matrix
#'  with expected variance/covariance in longitude and latitude of
#'  \code{targetPoints}, in \code{resampleProjection} units (default meters).
#' @param targetRaster For intrinsic tracking data, the results of
#'  \code{isoAssign} or a similar function of class \code{intrinsicAssign} or
#'  class \code{RasterBrick}/\code{RasterStack}, for example from the package
#'  \code{assignR}. In any case, it expresses location uncertainty on target
#'  range, through a raster of probabilities by grid cells.
#' @param originRaster For intrinsic tracking data, the results of
#'  \code{isoAssign} or a similar function of class \code{intrinsicAssign} or
#'  class \code{RasterBrick}/\code{RasterStack}, for example from the package
#'  \code{assignR}. In any case, it expresses location uncertainty on origin
#'  range, through a raster of probabilities by grid cells.
#' @param banded For band return data, a vector or matrix of the number of
#'  released animals from each origin site (including those never reencountered
#'  in a target site). If a matrix, the second dimension is taken as the number
#'  of age classes of released animals; the model estimates reencounter
#'  probability by age class but assumes transition probabilities are the same.
#'  Note that this age model is currently implemented only for \code{method}
#'  set to "MCMC", and only when banding data is analyzed alone (no telemetry
#'  data).
#' @param reencountered For band return data, either a matrix with B rows and W
#'  columns or a B x [number of ages] x W array. Number of animals reencountered
#'  on each target site (by age class banded as) by origin site they came from.
#' @param verbose 0 (default) to 3. 0 prints no output during run (except on
#'  convergence for \code{method} set to "MCMC"). 1 prints an update every 100
#'  samples or bootstraps (or a status bar for "MCMC").  2 prints an update
#'  every sample or bootstrap. 3 also prints the number of draws (for
#'  tuning \code{nSim}).
#' @param alpha Level for confidence/credible intervals provided. Default (0.05)
#'  gives 95 percent CI.
#' @param resampleProjection Projection when sampling from location uncertainty.
#'  Default is Equidistant Conic. The default setting preserves distances
#'  around latitude = 0 and longitude = 0. Other projections may work well,
#'  depending on the location of sites. Ignored unless data are entered using
#'  sites and points and/or rasters.
#' @param nSim Tuning parameter for GL or intrinsic data. Affects only the
#'  speed; 1000 seems to work well with our GL data and 10 for our intrinsic
#'  data, but your results may vary. For data combinations, we put the default
#'  higher (5000) to allow for more data conflicts. Should be integer > 0.
#'  Ignored when \code{method} is "MCMC".
#' @param maxTries Maximum number of times to run a single GL/intrinsic
#'  bootstrap before exiting with an error. Default is 300; you may want to make
#'  a little higher if your \code{nSim} is low and \code{nSamples} is high. Set
#'  to NULL to never exit. This parameter was added to prevent setups where some
#'  sample points never land on target sites from running indefinitely.
#' @param nBurnin For \code{method} set to "MCMC", \code{estTransition} runs a
#'  \code{JAGS} multinomial non-Markovian transitions model, for which it needs
#'  the number of burn-in samples before beginning to store results. Default
#'  5000.
#' @param nChains For \code{method} set to "MCMC", \code{estTransition} runs a
#'  \code{JAGS} multinomial non-Markovian transitions model, for which it needs
#'  the number of MCMC chains (to test for convergence). Default 3.
#' @param nThin For \code{method} set to "MCMC", \code{estTransition} runs a
#'  \code{JAGS} multinomial non-Markovian transitions model, for which it needs
#'  the thinning rate. Default 1.
#' @param dataOverlapSetting When there is more than one type of data, this
#'  setting allows the user some flexibility for clarifying which type(s) of
#'  data apply to which animals. Setting "dummy" (the default) indicates that
#'  there are dummy values within each dataset for the animals that isGL,
#'  isTelemetry, etc. don't have that data type (FALSE values). If no animals
#'  have a data type, no dummy values are required. If no animals have more than
#'  one type of data, the user can simplify processing their data by choosing
#'  setting "none" here. In this case, there should be no dummy values, and only
#'  the animals with a type of data should be included in that dataset. The
#'  third setting ("named") is not yet implemented, but will eventually allow
#'  another way to allow animals with more than one type of data with named
#'  animals linking records. When there is only one type of data, it is fastest
#'  to leave this on the default. Note that banding data entered through
#'  \code{banded} and \code{reencountered} are assumed to have no
#'  overlap with other data types, so none of this applies to those.
#' @param fixedZero When the user has a priori reasons to believe one or more
#'  transition probabilities are zero, they can indicate those here, and the
#'  model will keep them fixed at zero. This argument should be a matrix with
#'  two columns (for row and column of the transition probability matrix) and
#'  number of transitions being fixed to zero rows. For MCMC modeling,
#'  substantial evidence that a transition fixed to zero isn't zero may
#'  cause an error. For bootstrap modeling, a warning
#'  will come up if any bootstrap runs generate the transition fixed to zero,
#'  and the function will quit with an error if a very large number of runs do
#'  (> 10 * nSamples). Fixing transitions to zero may also slow down the
#'  bootstrap model somewhat.
#' @param targetRelAbund When some/all data have location error at origin sites
#'  (i.e., GL, raster, or probability table data with captured = "target" or
#'  "none"), unless the data were collected in proportion to abundance at target
#'  sites, simulation work indicates substantial bias in transition probability
#'  estimates can result. However, if these data are resampled in proportion to
#'  target site abundance, this bias is removed. This argument allows the user
#'  to provide an estimate of relative abundance at the target sites. Either
#'  a numeric vector of length [number target sites] that sums to 1, or an mcmc
#'  object (such as is produced by \code{\link{modelCountDataJAGS}}) or matrix
#'  with at least \code{nSamples} rows. If there are more than [number target
#'  sites] columns, the relevant columns should be labeled "relN[1]" through
#'  "relN[number target sites]".
#' @param method This important setting lets the user choose the estimation
#'  method used: bootstrap or MCMC (Markov chain Monte Carlo). Bootstrap (the
#'  default) now works with any and all types of data, whereas MCMC currently
#'  only works with banding and telemetry data (enter telemetry data for MCMC
#'  using \code{originAssignment} and \code{targetAssignment}, not
#'  \code{originPoints} and \code{targetPoints}). However, MCMC is
#'  usually faster (and may be a bit more accurate). The third option,
#'  "m-out-of-n-bootstrap", is still under development and should be left alone.
#' @param m We read that the m-out-of-n-bootstrap method may improve the
#'  coverage of confidence intervals for parameters on or near a boundary (0 or
#'  1 in this case). So we're testing that out. This still under development and
#'  not for the end user. In the m-out-of-n-bootstrap, m is the number of
#'  samples taken each time (less than the true sample size, n). If the
#'  "m-out-of-n-bootstrap" is chosen under \code{method} but this is left blank,
#'  currently the default is n/4, rounded up (no idea if that is reasonable).
#' @param psiPrior matrix with same dimensions as psi. Only relevant when
#'  \code{method} is "MCMC". Each row provides a Dirichlet
#'  (https://en.wikipedia.org/wiki/Dirichlet_distribution) prior on the
#'  transition probabilities from that origin site. The default (NULL) supplies
#'  Dirichlet parameters of all 1s, which is a standard uninformative Dirichlet
#'  prior. Setting these to other positive numbers is useful when you think a
#'  priori that certain transitions are unlikely, but don't want to rule them
#'  out altogether using \code{fixedZero}.
#' @param returnAllInput if TRUE (the default) the output includes all of the
#'  inputs. If FALSE, only the inputs currently used by another MigConnectivity
#'  function are included in the output. Switch this if you're worried about
#'  computer memory (and the output will be much slimmer).
#'
#' @return \code{estTransition} returns a list with the elements:
#' \describe{
#'   \item{\code{psi}}{List containing estimates of transition probabilities:
#'   \itemize{
#'    \item{\code{sample}} Array of sampled values for psi. \code{nSamples} x
#'      [number of origin sites] x [number of target sites]. Provided to allow
#'      the user to compute own summary statistics.
#'    \item{\code{mean}} Main estimate of psi matrix. [number of origin sites]
#'      x [number of target sites].
#'    \item{\code{se}} Standard error of psi, estimated from SD of
#'      \code{psi$sample}.
#'    \item{\code{simpleCI}} \code{1 - alpha} confidence interval for psi,
#'      estimated as \code{alpha/2} and \code{1 - alpha/2} quantiles of
#'      \code{psi$sample}.
#'    \item{\code{bcCI}} Bias-corrected \code{1 - alpha} confidence interval
#'      for psi. May be preferable to \code{simpleCI} when \code{mean} is the
#'      best estimate of psi. \code{simpleCI} is preferred when
#'      \code{median} is a better estimator. When the mean and median are equal,
#'      these should be identical.  Estimated as the
#'      \code{pnorm(2 * z0 + qnorm(alpha / 2))} and
#'      \code{pnorm(2 * z0 + qnorm(1 - alpha / 2))} quantiles of \code{sample},
#'      where z0 is the proportion of \code{sample < mean}.
#'    \item{\code{hpdCI}} \code{1 - alpha} credible interval for psi,
#'      estimated using the highest posterior density (HPD) method.
#'    \item{\code{median}} Median estimate of psi matrix.
#'    \item{\code{point}} Simple point estimate of psi matrix, not accounting
#'      for sampling error.
#'   }
#'   }
#'   \item{\code{r}}{List containing estimates of reencounter probabilities at
#'    each target site. NULL except when using direct band/ring reencounter
#'    data.}
#'   \item{\code{input}}{List containing the inputs to \code{estTransition}.}
#'   \item{\code{BUGSoutput}}{List containing \code{R2jags} output. Only present
#'    when using \code{method} of "MCMC".}
#' }
#'
#' @export
#'
#' @seealso \code{\link{estStrength}}, \code{\link{plot.estMigConnectivity}},
#'   \code{\link{estMC}}, \code{\link{estMantel}}
#'
#' @example inst/examples/estTransitionExamples.R
estTransition <- function(originSites = NULL, targetSites = NULL,
                          originPoints = NULL, targetPoints = NULL,
                          originAssignment = NULL, targetAssignment = NULL,
                          originNames = NULL, targetNames = NULL,
                          nSamples = 1000, isGL = FALSE, isTelemetry = FALSE,
                          isRaster = FALSE, isProb = FALSE,
                          captured = "origin", geoBias = NULL, geoVCov = NULL,
                          geoBiasOrigin = geoBias, geoVCovOrigin = geoVCov,
                          targetRaster = NULL, originRaster = NULL,
                          banded = NULL, reencountered = NULL,
                          verbose = 0, alpha = 0.05,
                          resampleProjection = 'ESRI:102010',
                          nSim = ifelse(any(isRaster & isGL) ||
                                          any(isRaster & isProb) ||
                                          any(isGL & isProb), 5000,
                                        ifelse(any(isGL), 1000,
                                               ifelse(any(isRaster), 10, 1))),
                          maxTries = 300,
                          nBurnin = 5000, nChains = 3, nThin = 1,
                          dataOverlapSetting = c("dummy", "none", "named"),
                          fixedZero = NULL,
                          targetRelAbund = NULL,
                          method = c("bootstrap", "MCMC",
                                     "m-out-of-n-bootstrap"),
                          m = NULL,
                          psiPrior = NULL,
                          returnAllInput = TRUE) {
  dataOverlapSetting <- match.arg(dataOverlapSetting)
  method <- match.arg(method)
  if (method != "MCMC") {
    psi <- estTransitionBoot(isGL=isGL, isTelemetry = isTelemetry,
                             isRaster = isRaster, isProb = isProb,
                             geoBias=geoBias, geoVCov=geoVCov,
                             geoBiasOrigin = geoBiasOrigin,
                             geoVCovOrigin=geoVCovOrigin,
                             targetPoints=targetPoints, targetSites=targetSites,
                             targetAssignment=targetAssignment,
                             originPoints=originPoints, originSites=originSites,
                             originAssignment=originAssignment,
                             originNames=originNames, targetNames=targetNames,
                             targetRaster = targetRaster,
                             originRaster = originRaster,
                             captured = captured,
                             nBoot = nSamples, verbose=verbose,
                             nSim = nSim, alpha = alpha,
                             resampleProjection = resampleProjection,
                             maxTries = maxTries,
                             dataOverlapSetting = dataOverlapSetting,
                             fixedZero = fixedZero,
                             targetRelAbund = targetRelAbund,
                             method = method, m = m,
                             banded = banded, reencountered = reencountered,
                             returnAllInput = returnAllInput)
  }
  else {
    psi <- estTransitionJAGS(banded = banded, reencountered = reencountered,
                             originAssignment = originAssignment,
                             targetAssignment = targetAssignment,
                             alpha = alpha,
                             nSamples = nSamples, verbose = verbose,
                             originNames = originNames,
                             targetNames = targetNames,
                             nBurnin = nBurnin, nThin = nThin,
                             nChains = nChains, fixedZero = fixedZero,
                             psiPrior = psiPrior,
                             returnAllInput = returnAllInput,
                             originPoints = originPoints,
                             targetPoints = targetPoints,
                             originSites = originSites,
                             targetSites = targetSites)
  }
  class(psi) <- c("estPsi", "estMigConnectivity")
  return(psi)
}

###############################################################################
# Resampling of uncertainty from geolocators and/or GPS data
###############################################################################
estMCGlGps <- function(originDist, targetDist, originRelAbund, isGL,
                       geoBias, geoVCov,
                       targetPoints, targetSites,
                       sampleSize = NULL,
                       targetAssignment=NULL,
                       originPoints=NULL, originSites=NULL,
                       originAssignment=NULL, originNames=NULL,
                       targetNames=NULL, nBoot = 1000, verbose=0,
                       nSim = 1000, calcCorr=TRUE, alpha = 0.05,
                       approxSigTest = FALSE, sigConst = 0,
            resampleProjection = 'ESRI:102010',
                       maxTries = 300,
                       row0=0,
                       maintainLegacyOutput = FALSE) {

  # Input checking and assignment
  if (!(verbose %in% 0:3))
    stop("verbose should be integer 0-3 for level of output during bootstrap: 0 = none, 1 = every 10, 2 = every run, 3 = number of draws")
  if (length(geoBias)!=2 && any(isGL))
    stop("geoBias should be vector of length 2 (expected bias in longitude and latitude of targetPoints, in resampleProjection units, default meters)")
  if (!isTRUE(all.equal(dim(geoVCov), c(2, 2), check.attributes = FALSE)) && any(isGL))
    stop("geoVCov should be 2x2 matrix (expected variance/covariance in longitude and latitude of targetPoints, in resampleProjection units, default meters)")
  if ((is.null(originPoints) || is.null(originSites)) &&
      is.null(originAssignment))
    stop("Need to define either originAssignment or originSites and originPoints")
  if (calcCorr && is.null(originPoints))
    stop('If calcCorr is TRUE, need to define originPoints')
  if ((is.null(targetPoints) || is.null(targetSites)) &&
      is.null(targetAssignment))
    stop("Need to define either targetAssignment or targetSites and targetPoints")
  nAnimals <- max(nrow(targetPoints), length(targetAssignment))
  if (length(isGL)==1)
    isGL <- rep(isGL, nAnimals)
  if(inherits(originSites, "SpatialPolygonsDataFrame") |
     inherits(originSites, "SpatialPolygons")){
    originSites <- sf::st_as_sf(originSites)}
    if(inherits(targetSites, "SpatialPolygonsDataFrame") |
       inherits(targetSites, "SpatialPolygons")){
         targetSites <- sf::st_as_sf(targetSites)}
  if (is.null(originAssignment)) {
    originAssignment <- suppressMessages(unclass(sf::st_intersects(x = originPoints,
                                                                   y = originSites,
                                                                   sparse = TRUE)))
    originAssignment[lengths(originAssignment)==0] <- NA
    if (any(lengths(originAssignment)>1)){
      stop("Overlapping originSites not allowed\n")
    }
    originAssignment <- unlist(originAssignment)
  }
  if (is.null(targetAssignment)) {
    targetAssignment <- suppressMessages(unclass(sf::st_intersects(x = targetPoints,
                                                                   y = targetSites,
                                                                   sparse = TRUE)))
    targetAssignment[lengths(targetAssignment)==0] <- NA
    if (any(lengths(targetAssignment)>1))
      stop("Overlapping targetSites not allowed\n")
    targetAssignment <- unlist(targetAssignment)
  }
  nOriginSites <- ifelse(is.null(originSites), nrow(originDist),
                         nrow(originSites))
  nTargetSites <- ifelse(is.null(targetSites), nrow(targetDist),
                         nrow(targetSites))
  if (nrow(targetPoints)!=nAnimals && length(targetAssignment)!=nAnimals ||
      length(originAssignment)!=nAnimals)
    stop("isGL should be the same length as originAssignment/originPoints and targetPoints/targetAssignment (number of animals)")
  # if (any(is.na(originAssignment)))
  #   stop("NAs in origin sites (make sure all points fall within polygons)")
  if (length(originRelAbund)!=nOriginSites || sum(originRelAbund)!=1)
    stop('originRelAbund should be vector with length number of origin sites that sums to 1')
  if(!is.null(originPoints))
    if(is.na(sf::st_crs(originPoints)) || !is.null(originSites) & is.na(sf::st_crs(originSites))) {
      stop('Coordinate system definition needed for originSites & originPoints')
    }
  if(!is.null(targetSites) & is.na(sf::st_crs(targetSites)) || !is.null(targetPoints) & is.na(sf::st_crs(targetPoints))){
    stop('Coordinate system definition needed for targetSites & targetPoints')
  }
  if(!is.null(targetPoints)){targetPoints <- sf::st_transform(targetPoints, resampleProjection)}
  if(!is.null(targetSites)){targetSites <- sf::st_transform(targetSites, resampleProjection)}
  if (is.null(targetNames))
    targetNames <- as.character(1:nTargetSites)
  if (is.null(originNames)) {
    originNames <- LETTERS[1:nOriginSites]
  }
  if (is.null(sampleSize))
    sampleSize <- nAnimals
  if (!identical(dim(originDist),rep(nOriginSites,2)) ||
      !identical(dim(targetDist),rep(nTargetSites,2)))
    stop('Distance matrices should be square with same number of sites of each type as assignments/points (with distances in meters)')
  sites.array <- psi.array <- array(0, c(nBoot, nOriginSites, nTargetSites),
                                    dimnames = list(1:nBoot, originNames,
                                                    targetNames))


  MC <- corr <- rep(NA, nBoot)

  # determine the number of animals from the input data
  nAnimals <- length(originAssignment)

  # Point estimate of MC
  pointSites <- array(0, c(nOriginSites, nTargetSites),
                       dimnames = list(originNames, targetNames))

  for(i in 1:nAnimals)
      pointSites[originAssignment[i], targetAssignment[i]] <-
      pointSites[originAssignment[i], targetAssignment[i]] + 1

  pointPsi <- prop.table(pointSites, 1)

  if (any(is.na(pointPsi)))
    pointMC <- NA
  else
    pointMC <- calcMC(originDist, targetDist, originRelAbund, pointPsi,
                    sampleSize = sampleSize)

  if (calcCorr) {
    pointMantel <- calcMantel(targetPoints = targetPoints,
                              originPoints = originPoints)
    targetDist1 <- pointMantel$targetDist
    originDistStart <- pointMantel$originDist
    pointCorr <- pointMantel$pointCorr
  }

  boot <- 1
  while (boot <= nBoot) {
    if (verbose > 1 || verbose == 1 && boot %% 100 == 0)
      cat("Bootstrap Run", boot, "of", nBoot, "at", date(), "\n")
    # Make sure have birds from every origin site
    origin.sample <- 'Filler' # Start with one origin site
    while (length(unique(origin.sample)) < nOriginSites) { #2
      # Sample individual animals with replacement
      animal.sample <- sample.int(nAnimals, replace=TRUE)
      # Get origin population for each animal sampled
      origin.sample <- originAssignment[animal.sample]
    }
    tSamp <- targetSample(isGL = isGL, geoBias = geoBias, geoVCov = geoVCov,
                          targetPoints = targetPoints, animal.sample = animal.sample,
                          targetSites = targetSites, targetAssignment = targetAssignment,
                          resampleProjection = resampleProjection, nSim = nSim,
                          maxTries = maxTries)
    target.sample <- tSamp$target.sample
    target.point.sample <- tSamp$target.point.sample
    if (verbose > 2)
      cat(' ', tSamp$draws, 'draw(s) (of length', nSim, 'and of', maxTries, 'possible).\n')
    # Now that we have breeding and non-breeding site for point, add to transition count matrix
    sites <- table(origin.sample, target.sample)
    sites.array[boot, as.integer(rownames(sites)), as.integer(colnames(sites))] <- sites
    # Create psi matrix as proportion of those from each breeding site that went to each NB site
    psi.array[boot, , ] <- prop.table(sites.array[boot, , ], 1)
    # Calculate MC from that psi matrix
    MC[boot] <- calcMC(originDist, targetDist, originRelAbund,
                       psi.array[boot, , ], sampleSize)
    if (verbose > 1 || verbose == 1 && boot %% 10 == 0)
      cat(" MC mean:", mean(MC, na.rm=TRUE), "SD:", sd(MC, na.rm=TRUE),
          "low quantile:", quantile(MC, alpha/2, na.rm=TRUE),
          "high quantile:", quantile(MC, 1-alpha/2, na.rm=TRUE), "\n")
    if (calcCorr) {
      colnames(target.point.sample) <- c("x", "y")
      target.point.sample <- sf::st_as_sf(data.frame(target.point.sample),
                                          coords = c("x","y"),
                                          crs = resampleProjection)
      originDist1 <- originDistStart[animal.sample, animal.sample]
      corr[boot] <- calcMantel(originDist = originDist1,
                               targetPoints = target.point.sample)$pointCorr
      if (verbose > 1 || verbose == 1 && boot %% 10 == 0)
        cat(" Correlation mean:", mean(corr, na.rm=TRUE), "SD:", sd(corr, na.rm=TRUE),
            "low quantile:", quantile(corr, alpha/2, na.rm=TRUE),
            "high quantile:", quantile(corr, 1-alpha/2, na.rm=TRUE), "\n")
    }
    if (!is.na(MC[boot]))
      boot <- boot + 1
  }
  MC.z0 <- qnorm(sum(MC<mean(MC, na.rm = TRUE), na.rm = TRUE)/length(which(!is.na(MC))))
  bcCI <- quantile(MC, pnorm(2*MC.z0+qnorm(c(alpha/2, 1-alpha/2))),
                       na.rm=TRUE, type = 8, names = FALSE)
  MC.mcmc <- coda::as.mcmc(MC) # Ha!
  hpdCI <- as.vector(coda::HPDinterval(MC.mcmc, 1-alpha))
  if (!approxSigTest)
    simpleP <- bcP <- NULL
  else {
    if (pointMC > sigConst)
      simpleP <- sum(MC < sigConst) / nBoot
    else
      simpleP <- sum(MC > sigConst) / nBoot
    if (simpleP == 0)
      simpleP <- 0.5 / nBoot
    bcP <- pnorm(qnorm(simpleP) - 2 * MC.z0)
    if (pointMC < sigConst)
      bcP <- 1 - bcP
  }
  if (calcCorr) {
    meanCorr <- mean(corr, na.rm=TRUE)
    medianCorr <- median(corr, na.rm=TRUE)
    seCorr <- sd(corr, na.rm=TRUE)
    simpleCICorr <- quantile(corr, c(alpha/2, 1-alpha/2), na.rm=TRUE,
                             names = FALSE)
    corr.z0 <- qnorm(sum((corr)<meanCorr)/nBoot)
    bcCICorr <- quantile(corr, pnorm(2*corr.z0+qnorm(c(alpha/2, 1-alpha/2))),
                           na.rm=TRUE, names = FALSE)
  } else
    pointCorr <- meanCorr <- medianCorr <- seCorr <- simpleCICorr <- bcCICorr <- NULL
  meanMC <- mean(MC, na.rm=TRUE)
  medianMC <- median(MC, na.rm=TRUE)
  seMC <- sd(MC, na.rm=TRUE)
  simpleCI <- quantile(MC, c(alpha/2, 1-alpha/2), na.rm=TRUE, names = FALSE)
  meanPsi <- apply(psi.array, 2:3, mean)
  medianPsi <- apply(psi.array, 2:3, median)
  sePsi <- apply(psi.array, 2:3, sd)
  simpleCIPsi <- apply(psi.array, 2:3, quantile, probs = c(alpha/2, 1-alpha/2),
                       na.rm=TRUE, names = FALSE)
  bcCIPsi <- array(NA, dim = c(2, nOriginSites, nTargetSites),
                   dimnames = list(NULL, originNames, targetNames))
  for (i in 1:nOriginSites) {
    for (j in 1:nTargetSites) {
      psi.z0 <- qnorm(sum(psi.array[, i, j] < meanPsi[i, j], na.rm = TRUE) /
                        length(which(!is.na(psi.array[, i, j]))))
      bcCIPsi[ , i, j] <- quantile(psi.array[, i, j],
                                   pnorm(2 * psi.z0 + qnorm(c(alpha/2, 1-alpha/2))),
                                   na.rm=TRUE, type = 8, names = FALSE)
    }
  }
  if (maintainLegacyOutput) {
    return(list(sampleMC=MC, samplePsi = psi.array, pointPsi = pointPsi,
                pointMC=pointMC, meanMC=meanMC,
                medianMC=medianMC, seMC=seMC, simpleCI=simpleCI,
                bcCI=bcCI, hpdCI=hpdCI, simpleP = simpleP, bcP = bcP,
                sampleCorr = corr, pointCorr = pointCorr,
                meanCorr = meanCorr, medianCorr = medianCorr, seCorr=seCorr,
                simpleCICorr=simpleCICorr, bcCICorr=bcCICorr,
                inputSampleSize = sampleSize,
                alpha = alpha, sigConst = sigConst,
                psi = list(sample = psi.array, mean = meanPsi, se = sePsi,
                           simpleCI = simpleCIPsi, bcCI = bcCIPsi,
                           median = medianPsi, point = pointPsi),
                MC = list(sample = MC, mean = meanMC, se = seMC,
                          simpleCI = simpleCI, bcCI = bcCI, hpdCI = hpdCI,
                          median = medianMC, point = pointMC,
                          simpleP = simpleP, bcP = bcP),
                corr = list(sample = corr, mean = meanCorr, se = seCorr,
                            simpleCI = simpleCICorr, bcCI = bcCICorr,
                            median = medianCorr, point = pointCorr),
                input = list(originDist = originDist, targetDist = targetDist,
                             originRelAbund = originRelAbund,
                             sampleSize = sampleSize, originSites = originSites,
                             targetSites = targetSites,
                             originPoints = originPoints,
                             targetPoints = targetPoints,
                             originAssignment = originAssignment,
                             targetAssignment = targetAssignment,
                             originNames = originNames,
                             targetNames = targetNames,
                             nSamples = nBoot, nSim = nSim,
                             isGL=isGL, geoBias=geoBias, geoVCov=geoVCov,
                             row0 = row0,
                             verbose = verbose, calcCorr = calcCorr,
                             alpha = alpha,
                             approxSigTest = approxSigTest, sigConst = sigConst,
                             resampleProjection = resampleProjection,
                             maxTries = maxTries, maintainLegacyOutput = TRUE)))
  }
  else {
    return(list(psi = list(sample = psi.array, mean = meanPsi, se = sePsi,
                           simpleCI = simpleCIPsi, bcCI = bcCIPsi,
                           median = medianPsi, point = pointPsi),
                MC = list(sample = MC, mean = meanMC, se = seMC,
                          simpleCI = simpleCI, bcCI = bcCI, hpdCI = hpdCI,
                          median = medianMC, point = pointMC,
                          simpleP = simpleP, bcP = bcP),
                corr = list(sample = corr, mean = meanCorr, se = seCorr,
                            simpleCI = simpleCICorr, bcCI = bcCICorr,
                            median = medianCorr, point = pointCorr),
                input = list(originDist = originDist, targetDist = targetDist,
                             originRelAbund = originRelAbund,
                             sampleSize = sampleSize, originSites = originSites,
                             targetSites = targetSites,
                             originPoints = originPoints,
                             targetPoints = targetPoints,
                             originAssignment = originAssignment,
                             targetAssignment = targetAssignment,
                             originNames = originNames,
                             targetNames = targetNames,
                             nSamples = nBoot, nSim = nSim,
                             isGL=isGL, geoBias=geoBias, geoVCov=geoVCov,
                             row0 = row0,
                             verbose = verbose, calcCorr = calcCorr,
                             alpha = alpha,
                             approxSigTest = approxSigTest, sigConst = sigConst,
                             resampleProjection = resampleProjection,
                             maxTries = maxTries,
                             maintainLegacyOutput = FALSE)))
  }
}

###############################################################################
#
#
# Resampling of uncertainty from intrinsic markers (Stable-Isotopes)
#
#
###############################################################################
estMCisotope <- function(targetDist=NULL,
                         originRelAbund,
                         targetIntrinsic,
                         targetSites = NULL,
                         sampleSize = NULL,
                         originPoints=NULL,
                         originSites=NULL,
                         originDist = NULL,
                         originAssignment=NULL,
                         originNames=NULL,
                         targetNames=NULL,
                         nBoot = 1000,
                         verbose=0,
                         nSim = NULL,
                         calcCorr=TRUE,
                         alpha = 0.05,
                         approxSigTest = FALSE,
                         sigConst = 0,
                         resampleProjection = sf::st_crs(4326),# MigConnectivity::projections$WGS84,
                         maxTries = 300,
                         maintainLegacyOutput = FALSE) {

  # Input checking and assignment
  if (!(verbose %in% 0:3)){
    stop("verbose should be integer 0-3 for level of output during bootstrap: 0 = none, 1 = every 10, 2 = every run, 3 = number of draws")}

  if ((is.null(originPoints) || is.null(originSites)) && is.null(originAssignment)){
    stop("Need to define either originAssignment or originSites and originPoints")}

  if (calcCorr && is.null(originPoints)){
    stop('If calcCorr is TRUE, need to define originPoints')}

  if (is.null(targetIntrinsic)){
    stop("Need to define targetIntrinsic")
  }

  if (is.null(targetSites))
    targetSites <- targetIntrinsic$targetSites
  if(!inherits(targetSites, "sf")){
    targetSites <- sf::st_as_sf(targetSites)
    targetSites <- sf::st_make_valid(targetSites)
    targetSites <- sf::st_union(targetSites)
  }

  if (is.null(targetDist))
    targetDist <- distFromPos(sf::st_coordinates(sf::st_centroid(targetSites)))


  if (!inherits(targetIntrinsic, 'isoAssign'))
    stop("targetIntrinsic should be output of isoAssign when isIntrinsic == TRUE")

  pointsAssigned <- !(is.null(targetIntrinsic$SingleCell) ||
                        all(is.na(targetIntrinsic$SingleCell)))

  if(inherits(originSites, "SpatialPolygonsDataFrame")){
    originSites <- sf::st_as_sf(originSites)
    originSites <- sf::st_union(originSites)
  }

  if (is.null(originAssignment))
    originAssignment <-
      suppressMessages(as.numeric(unclass(sf::st_intersects(x = originPoints,
                                                            y = originSites,
                                                            sparse = TRUE))))

  nAnimals <- ifelse(pointsAssigned, dim(targetIntrinsic$SingleCell)[3],
                     dim(targetIntrinsic$probassign)[3])
  targetSites <- sf::st_transform(targetSites, resampleProjection)

  if (length(originAssignment)!=nAnimals)
    stop("originAssignment/originPoints should be the same length as targetIntrinsic (number of animals)")

  pointsInSites <- FALSE
  if (pointsAssigned && !is.null(targetSites)) {
    nSamples <- dim(targetIntrinsic$SingleCell)[1]
    targCon <- array(NA, c(nSamples, nAnimals))

    for(i in 1:nSamples) {
      temptargCon <- sf::st_as_sf(data.frame(t(targetIntrinsic$SingleCell[i,,])),
                                  coords = c("Longitude","Latitude"),
                                  crs = sf::st_crs(targetSites))

      targCon[i, ] <- suppressMessages(as.numeric(unclass(sf::st_intersects(x = temptargCon, y = targetSites, sparse = TRUE))))
    }

    if (!any(is.na(targCon)))
      pointsInSites <- TRUE
    else if (verbose > 0)
      cat('Single cell assignment points supplied, but some points (proportion',
        sum(is.na(targCon))/length(targCon), ') not in targetSites\n')
  }
  else
    targCon <- NULL
  nOriginSites <- length(unique(originAssignment))
  nTargetSites <- ifelse(is.null(targetSites), nrow(targetDist), nrow(targetSites))
  # cat("nSites figured\n")

  if (any(is.na(originAssignment)))
    stop("NAs in origin sites (make sure all points fall within polygons)")

  if (length(originRelAbund)!=nOriginSites || sum(originRelAbund)!=1)
    stop('originRelAbund should be vector with length number of origin sites that sums to 1')

  if(!is.null(originPoints))
    if(is.na(sf::st_crs(originPoints)) || is.na(sf::st_crs(originSites))) {
      stop('Coordinate system definition needed for originSites & originPoints')
    }



  if (is.null(targetNames))
    # Need to be able to conserve the names from input
    targetNames <- 1:nrow(targetSites)

  if (is.null(originNames))
    # need to conserve names from inputs
    originNames <- 1:nrow(originSites)
  # cat("names figured", targetNames, originNames, "\n")

  if (is.null(sampleSize))
    sampleSize <- nAnimals

  # cat(nOriginSites, "\n")
  # cat(nTargetSites, "\n")
  if (any(dim(originDist)!=rep(nOriginSites,2)) ||
      any(dim(targetDist)!=rep(nTargetSites,2)))
    stop('Distance matrices should be square with same number of sites of each type as assignments/points (with distances in meters)')


  if (calcCorr) {
    originPoints2 <- sf::st_transform(originPoints, 4326)
    originDistStart <- distFromPos(sf::st_coordinates(originPoints2))
  }



  sites.array <- psi.array <- array(0, c(nBoot, nOriginSites, nTargetSites),
                                    dimnames = list(1:nBoot, originNames,
                                                    targetNames))


  MC <- corr <- rep(NA, nBoot)


  boot <- 1
  while (boot <= nBoot) {
    if (verbose > 1 || verbose == 1 && boot %% 100 == 0)
      cat("Bootstrap Run", boot, "of", nBoot, "at", date(), "\n")
    # Make sure have birds from every origin site
    origin.sample <- 'Filler' # Start with one origin site
    while (length(unique(origin.sample)) < nOriginSites) { #2
      # Sample individual animals with replacement
      animal.sample <- sample.int(nAnimals, replace=TRUE)
      # Get origin points for those animals
      if (calcCorr)
        origin.point.sample <- originPoints[animal.sample,]
      # Get origin population for each animal sampled
      origin.sample <- originAssignment[animal.sample]
    }
    if (verbose > 2)
      cat("Origin sample complete, making target sample with",
          ifelse(pointsAssigned, "points assigned,", "no points assigned,"),
          ifelse(pointsInSites, "all points in sites\n", "not all points in sites\n"))
    # Resample from points for each animal
    tSamp <- targetSampleIsotope(targetIntrinsic = targetIntrinsic,
                                 animal.sample = animal.sample,
                                 targetSites = targetSites,
                                 resampleProjection = resampleProjection, nSim = nSim,
                                 maxTries = maxTries, pointsAssigned = pointsAssigned,
                                 targCon = targCon, pointsInSites = pointsInSites)

    target.sample <- tSamp$target.sample
    target.point.sample <- tSamp$target.point.sample
    if (verbose > 2 & !pointsInSites)
      cat(' ', tSamp$draws, 'draw(s) (of length', nSim, 'and of', maxTries, 'possible).\n')
    # Now that we have breeding and non-breeding site for point, add to transition count matrix
    sites <- table(origin.sample, target.sample)
    sites.array[boot, as.integer(rownames(sites)), as.integer(colnames(sites))] <- sites
    # Create psi matrix as proportion of those from each breeding site that went to each NB site
    psi.array[boot, , ] <- prop.table(sites.array[boot, , ], 1)
    # Calculate MC from that psi matrix
    MC[boot] <- calcMC(originDist, targetDist, originRelAbund,
                       psi.array[boot, , ], sampleSize)
    if (verbose > 1 || verbose == 1 && boot %% 10 == 0)
      cat(" MC mean:", mean(MC, na.rm=TRUE), "SD:", sd(MC, na.rm=TRUE),
          "low quantile:", quantile(MC, alpha/2, na.rm=TRUE),
          "high quantile:", quantile(MC, 1-alpha/2, na.rm=TRUE), "\n")
    if (calcCorr) {
      originDist1 <- originDistStart[animal.sample, animal.sample]
      colnames(target.point.sample) <- c("x","y")
      target.point.sample <- sf::st_as_sf(data.frame(target.point.sample),
                                          coords = c("x","y"),
                                          crs = resampleProjection)
      #target.point.sample <- sf::st_sfc(t.p.s.geom, crs = resampleProjection)
      #target.point.sample <- sp::SpatialPoints(target.point.sample,sp::CRS(resampleProjection))
      corr[boot] <- calcMantel(originDist = originDist1, targetPoints = target.point.sample)$pointCorr
      if (verbose > 1 || verbose == 1 && boot %% 10 == 0)
        cat(" Correlation mean:", mean(corr, na.rm=TRUE), "SD:", sd(corr, na.rm=TRUE),
            "low quantile:", quantile(corr, alpha/2, na.rm=TRUE),
            "high quantile:", quantile(corr, 1-alpha/2, na.rm=TRUE), "\n")
    }
    if (!is.na(MC[boot]))
      boot <- boot + 1
  }
  MC.z0 <- qnorm(sum(MC<mean(MC, na.rm = TRUE), na.rm = TRUE) /
                   length(which(!is.na(MC))))
  bcCI <- quantile(MC, pnorm(2*MC.z0+qnorm(c(alpha/2, 1-alpha/2))),
                   na.rm=TRUE, type = 8, names = FALSE)
  MC.mcmc <- coda::as.mcmc(MC) # Ha!
  hpdCI <- as.vector(coda::HPDinterval(MC.mcmc, 1-alpha))
  if (!approxSigTest)
    simpleP <- bcP <- NULL
  else {
    if (mean(MC, na.rm=TRUE) > sigConst)
      simpleP <- sum(MC < sigConst) / nBoot
    else
      simpleP <- sum(MC > sigConst) / nBoot
    if (simpleP == 0)
      simpleP <- 0.5 / nBoot
    bcP <- pnorm(qnorm(simpleP) - 2 * MC.z0)
    if (mean(MC, na.rm=TRUE) < sigConst)
      bcP <- 1 - bcP
  }
  if (calcCorr) {
    meanCorr <- mean(corr, na.rm=TRUE)
    medianCorr <- median(corr, na.rm=TRUE)
    seCorr <- sd(corr, na.rm=TRUE)
    simpleCICorr <- quantile(corr, c(alpha/2, 1-alpha/2), na.rm=TRUE, type = 8,
                             names = FALSE)
    corr.z0 <- qnorm(sum((corr)<meanCorr)/nBoot)
    bcCICorr <- quantile(corr, pnorm(2*corr.z0+qnorm(c(alpha/2, 1-alpha/2))),
                         na.rm=TRUE, type = 8, names = FALSE)
    pointCorr <- NULL
  } else
    pointCorr <- meanCorr <- medianCorr <- seCorr <- simpleCICorr <- bcCICorr <- NULL
  meanMC <- mean(MC, na.rm=TRUE)
  medianMC <- median(MC, na.rm=TRUE)
  seMC <- sd(MC, na.rm=TRUE)
  simpleCI <- quantile(MC, c(alpha/2, 1-alpha/2), na.rm=TRUE, type = 8,
                       names = FALSE)
  meanPsi <- apply(psi.array, 2:3, mean)
  medianPsi <- apply(psi.array, 2:3, median)
  sePsi <- apply(psi.array, 2:3, sd)
  simpleCIPsi <- apply(psi.array, 2:3, quantile, probs = c(alpha/2, 1-alpha/2),
                       na.rm=TRUE, type = 8, names = FALSE)
  bcCIPsi <- array(NA, dim = c(2, nOriginSites, nTargetSites),
                   dimnames = list(NULL, originNames, targetNames))
  for (i in 1:nOriginSites) {
    for (j in 1:nTargetSites) {
      psi.z0 <- qnorm(sum(psi.array[, i, j] < meanPsi[i, j], na.rm = TRUE) /
                        length(which(!is.na(psi.array[, i, j]))))
      bcCIPsi[ , i, j] <- quantile(psi.array[, i, j],
                                   pnorm(2 * psi.z0 + qnorm(c(alpha/2, 1-alpha/2))),
                                   na.rm=TRUE, type = 8, names = FALSE)
    }
  }
  if (maintainLegacyOutput) {
    return(list(sampleMC = MC, samplePsi = psi.array, pointPsi = NULL,
                pointMC = NULL, meanMC = meanMC,
                medianMC = medianMC, seMC = seMC, simpleCI = simpleCI,
                bcCI = bcCI, hpdCI = hpdCI, simpleP = simpleP, bcP = bcP,
                sampleCorr = corr, pointCorr = pointCorr,
                meanCorr = meanCorr, medianCorr = medianCorr, seCorr=seCorr,
                simpleCICorr=simpleCICorr, bcCICorr=bcCICorr,
                inputSampleSize = sampleSize,
                alpha = alpha, sigConst = sigConst,
                psi = list(sample = psi.array, mean = meanPsi, se = sePsi,
                           simpleCI = simpleCIPsi, bcCI = bcCIPsi,
                           median = medianPsi, point = NULL),
                MC = list(sample = MC, mean = meanMC, se = seMC,
                          simpleCI = simpleCI, bcCI = bcCI, hpdCI = hpdCI,
                          median = medianMC, point = NULL,
                          simpleP = simpleP, bcP = bcP),
                corr = list(sample = corr, mean = meanCorr, se = seCorr,
                            simpleCI = simpleCICorr, bcCI = bcCICorr,
                            median = medianCorr, point = pointCorr),
                input = list(originDist = originDist, targetDist = targetDist,
                             originRelAbund = originRelAbund,
                             sampleSize = sampleSize, originSites = originSites,
                             targetSites = targetSites,
                             originPoints = originPoints,
                             originAssignment = originAssignment,
                             originNames = originNames,
                             targetNames = targetNames,
                             nSamples = nBoot, nSim = nSim,
                             verbose = verbose, calcCorr = calcCorr,
                             alpha = alpha,
                             approxSigTest = approxSigTest, sigConst = sigConst,
                             resampleProjection = resampleProjection,
                             maxTries = maxTries,
                             targetIntrinsic = targetIntrinsic,
                             isIntrinsic = TRUE,
                             maintainLegacyOutput = TRUE)))
  }
  else {
    return(list(psi = list(sample = psi.array, mean = meanPsi, se = sePsi,
                           simpleCI = simpleCIPsi, bcCI = bcCIPsi,
                           median = medianPsi, point = NULL),
                MC = list(sample = MC, mean = meanMC, se = seMC,
                          simpleCI = simpleCI, bcCI = bcCI, hpdCI = hpdCI,
                          median = medianMC, point = NULL,
                          simpleP = simpleP, bcP = bcP),
                corr = list(sample = corr, mean = meanCorr, se = seCorr,
                            simpleCI = simpleCICorr, bcCI = bcCICorr,
                            median = medianCorr, point = pointCorr),
                input = list(originDist = originDist, targetDist = targetDist,
                             originRelAbund = originRelAbund,
                             sampleSize = sampleSize, originSites = originSites,
                             targetSites = targetSites,
                             originPoints = originPoints,
                             originAssignment = originAssignment,
                             originNames = originNames,
                             targetNames = targetNames,
                             nSamples = nBoot, nSim = nSim,
                             verbose = verbose, calcCorr = calcCorr,
                             alpha = alpha,
                             approxSigTest = approxSigTest, sigConst = sigConst,
                             resampleProjection = resampleProjection,
                             maxTries = maxTries,
                             targetIntrinsic = targetIntrinsic,
                             isIntrinsic = TRUE,
                             maintainLegacyOutput = FALSE)))
  }
}

###############################################################################
#' Estimate migratory connectivity
#'
#' Resampling of uncertainty for migratory connectivity strength (MC)
#' and transition probabilities (psi) from RMark psi matrix estimates or
#' samples of psi and/or JAGS
#' relative abundance MCMC samples OR SpatialPoints geolocators and/or GPS
#' data OR intrinsic markers such as isotopes. NOTE: active development of this
#' function is ending. We suggest users estimate psi with
#' \code{\link{estTransition}}, MC with \code{\link{estStrength}}, and Mantel
#' correlations (rM) with \code{\link{estMantel}}.
#'
#' @param originDist Distances between the B origin sites.  Symmetric B by B
#'  matrix
#' @param targetDist Distances between the W target sites.  Symmetric W by W
#'  matrix.  Optional for intrinsic data
#' @param originRelAbund Relative abundance estimates at B origin sites. Either
#'  a numeric vector of length B that sums to 1 or an mcmc object with
#'  \code{nSamples} rows  and columns including 'relN[1]' through 'relN[B]'.
#'  Currently, an mcmc object doesn't work with geolocator, GPS, or intrinsic
#'  data
#' @param psi Transition probabilities between B origin and W target sites.
#'  Either a matrix with B rows and W columns where rows sum to 1, an array with
#'  dimensions x, B, and W (with x samples of the transition probability matrix
#'  from another model), or a MARK object with estimates of transition
#'  probabilities.  If you are estimating MC from GPS, geolocator, or intrinsic
#'  data, leave this as NULL
#' @param sampleSize Total sample size of animals that psi will be estimated
#'  from. Should be the number of animals released in one of the origin sites
#'  and observed in one of the target sites.  Optional, but recommended, unless
#'  you are estimating MC from GPS, geolocator, intrinsic, or direct band return
#'  data (in which case the function can calculate it for you)
#' @param originSites If \code{psi} is a MARK object, this must be a numeric
#'  vector indicating which sites are origin.  If using GPS, geolocator, or
#'  intrinsic data, this can be the geographic definition of sites in the
#'  release season
#' @param targetSites If \code{psi} is a MARK object, this must be a numeric
#'  vector indicating which sites are target.  If using GPS, geolocator, or
#'  intrinsic data, this must be the geographic definition of sites in the
#'  non-release season.  Optional for intrinsic data; if left out, the function
#'  will use the \code{targetSites} defined in \code{targetIntrinsic}
#' @param originPoints A \code{POINT} sf object, with length number of
#'    animals tracked.  Each point indicates the release location of an animal
#' @param targetPoints For GL or GPS data, a \code{POINT} sf object, with
#'    length number of animals tracked.  Each point indicates the point estimate
#'    location in the non-release season
#' @param originAssignment Assignment of \code{originPoints} to release season
#'    sites. Integer vector with length number of animals tracked. Optional,
#'    but if using GL or GPS data, either \code{originAssignment} or
#'    \code{originSites} and \code{originPoints} should be defined
#' @param targetAssignment Optional. Point estimate assignment of
#'    \code{targetPoints} to non-release season sites. Integer vector with
#'    length number of animals tracked
#' @param originNames Optional but recommended. Vector of names for the release season sites
#' @param targetNames Optional but recommended. Vector of names for the non-release season
#'    sites
#' @param nSamples Number of times to resample \code{psi} and/or
#'  \code{originRelAbund} OR number of post-burn-in MCMC samples to store (band
#'  data) OR number of times to resample \code{targetPoints}
#'  for intrinsic data OR number of bootstrap runs for GL or GPS data. In
#'  the last two cases, animals are sampled with replacement for each. For all,
#'  the purpose is to estimate sampling uncertainty
#' @param nSim Tuning parameter for GL or intrinsic data. Affects only the
#'    speed; 1000 seems to work well with our GL data and 10 for our intrinsic
#'    data, but your results may vary.  Should be integer > 0
#' @param isGL Indicates whether or which animals were tracked with geolocators.
#'    Should be either single TRUE or FALSE value, or vector with length of
#'    number of animals tracked, with TRUE for animals in
#'    \code{targetPoints} with geolocators and FALSE for animals with GPS
#' @param geoBias For GL data, vector of length 2 indicating expected bias
#'    in longitude and latitude of \code{targetPoints}, in
#'    \code{resampleProjection} units (default meters)
#' @param geoVCov For GL data, 2x2 matrix with expected variance/covariance
#'    in longitude and latitude of \code{targetPoints}, in
#'    \code{resampleProjection} units (default meters)
#' @param row0 If \code{originRelAbund} is an mcmc object, this can be set
#'  to 0 (default) or any greater integer to specify where to stop ignoring
#'  samples ("burn-in")
#' @param verbose 0 (default) to 3. 0 prints no output during run. 1 prints
#'  a line every 100 samples or bootstraps and a summary every 10.  2 prints a
#'  line and summary every sample or bootstrap. 3 also prints the number of
#'  draws (for tuning nSim for GL/intrinsic data only)
#' @param calcCorr In addition to MC, should function also estimate Mantel
#'    correlation between release and non-release locations (GPS or GL data
#'    only)?  Default is FALSE
#' @param alpha Level for confidence/credible intervals provided
#' @param approxSigTest Should function compute approximate one-sided
#'    significance tests (p-values) for MC from the bootstrap?  Default is
#'    FALSE
#' @param sigConst Value to compare MC to in significance test.
#'    Default is 0
#' @param resampleProjection Projection when sampling from geolocator
#'    bias/error. This projection needs units = m. Default is Equidistant
#'    Conic. The default setting preserves distances around latitude = 0 and
#'    longitude = 0. Other projections may work well, depending on the location
#'    of \code{targetSites}.  Ignored unless data are geolocator or GPS
#' @param maxTries Maximum number of times to run a single GL/intrinsic
#'    bootstrap before exiting with an error.  Default is 300.  Set to NULL to
#'    never stop.  This parameter was added to prevent GL setups where some
#'    sample points never land on target sites from running indefinitely
#' @param targetIntrinsic For intrinsic tracking data, the results of
#'    \code{isoAssign} or a similar function, of class \code{intrinsicAssign}
#' @param isIntrinsic Logical indicating whether the animals are tracked via
#'  intrinsic marker (e.g. isotopes) or not.  Currently estMC will only estimate
#'  connectivity for all intrinsically marked animals or all extrinsic (e.g.,
#'  bands, GL, or GPS), so isIntrinsic should be a single TRUE or FALSE
#' @param maintainLegacyOutput version 0.4.0 of \code{MigConnectivity}
#'  updated the structure of the estimates. If you have legacy code that refers
#'  to elements within a \code{estMigConnectivity} object, you can set this
#'  to TRUE to also keep the old structure. Defaults to FALSE
#'
#' @return NOTE: Starting with version 0.4.0 of \code{MigConnectivity}, we've
#' updated the structure of \code{MigConnectivityEstimate} objects. Below we
#' describe the updated structure. If parameter \code{maintainLegacyOutput} is
#' set to TRUE, the list will start with the old structure: \code{sampleMC},
#' \code{samplePsi}, \code{pointPsi}, \code{pointMC}, \code{meanMC},
#' \code{medianMC}, \code{seMC}, \code{simpleCI}, \code{bcCI}, \code{hpdCI},
#' \code{simpleP}, \code{bcP}, \code{sampleCorr}, \code{pointCorr},
#' \code{meanCorr, medianCorr, seCorr, simpleCICorr, bcCICorr},
#' \code{inputSampleSize}, \code{alpha}, and \code{sigConst}.
#'
#' \code{estMC} returns a list with the elements:
#' \describe{
#'   \item{\code{psi}}{List containing estimates of transition probabilities:
#'   \itemize{
#'    \item{\code{sample}} Array of sampled values for psi. \code{nSamples} x
#'      [number of origin sites] x [number of target sites]. Provided to allow
#'      the user to compute own summary statistics.
#'    \item{\code{mean}} Main estimate of psi matrix. [number of origin sites]
#'      x [number of target sites].
#'    \item{\code{se}} Standard error of psi, estimated from SD of
#'      \code{psi$sample}.
#'    \item{\code{simpleCI}} \code{1 - alpha} confidence interval for psi,
#'      estimated as \code{alpha/2} and \code{1 - alpha/2} quantiles of
#'      \code{psi$sample}.
#'    \item{\code{bcCI}} Bias-corrected \code{1 - alpha} confidence interval
#'      for psi.  Preferable to \code{simpleCI} when \code{mean} is the
#'      best estimate of psi. \code{simpleCI} is preferred when
#'      \code{median} is a better estimator. When \code{meanMC==medianMC},
#'      these should be identical.  Estimated as the
#'      \code{pnorm(2 * z0 + qnorm(alpha / 2))} and
#'      \code{pnorm(2 * z0 + qnorm(1 - alpha / 2))} quantiles of \code{sample},
#'      where z0 is the proportion of \code{sample < mean}.
#'    \item{\code{median}} Median estimate of psi matrix.
#'    \item{\code{point}} Simple point estimate of psi matrix, not accounting
#'      for sampling error. NULL when \code{isIntrinsic == TRUE}.
#'   }
#'   }
#'   \item{\code{MC}}{List containing estimates of migratory connectivity
#'    strength:
#'    \itemize{
#'      \item{\code{sample}} \code{nSamples} sampled values for
#'       MC. Provided to allow the user to compute own summary statistics.
#'      \item{\code{mean}} Mean of \code{MC$sample}. Main estimate of MC,
#'       incorporating parametric uncertainty.
#'      \item{\code{se}} Standard error of MC, estimated from SD of
#'       \code{MC$sample}.
#'      \item{\code{simpleCI}} Default\code{1 - alpha} confidence interval for
#'       MC, estimated as \code{alpha/2} and \code{1 - alpha/2} quantiles of
#'       \code{MC$sample}.
#'      \item{\code{bcCI}} Bias-corrected \code{1 - alpha} confidence interval
#'       for MC.  Preferable to \code{MC$simpleCI} when \code{MC$mean} is the
#'       best estimate of MC. \code{MC$simpleCI} is preferred when
#'       \code{MC$median} is a better estimator. When \code{MC$mean==MC$median},
#'       these should be identical.  Estimated as the
#'       \code{pnorm(2 * z0 + qnorm(alpha / 2))} and
#'       \code{pnorm(2 * z0 + qnorm(1 - alpha / 2))} quantiles of \code{MC$sample},
#'       where z0 is the proportion of \code{MC$sample < MC$mean}.
#'      \item{\code{hpdCI}} \code{1 - alpha} credible interval for MC,
#'       estimated using the highest posterior density (HPD) method.
#'      \item{\code{median}} Median of MC, alternate estimator also including
#'       parametric uncertainty.
#'      \item{\code{point}} Simple point estimate of MC, using the point
#'      estimates of \code{psi} and \code{originRelAbund}, not accounting
#'      for sampling error. NULL when \code{isIntrinsic == TRUE}.
#'      \item{\code{simpleP}} Approximate p-value for MC, estimated as the
#'      proportion of bootstrap iterations where MC < \code{sigConst} (or MC >
#'      \code{sigConst} if \code{pointMC < sigConst}).  Note that if the
#'      proportion is 0, a default value of 0.5 / \code{nSamples} is provided,
#'      but this is best interpreted as p < 1 / \code{nSamples}.  NULL when
#'      \code{approxSigTest==FALSE}.
#'      \item{\code{bcP}} Approximate bias-corrected p-value for MC, estimated as
#'      \code{pnorm(qnorm(simpleP) - 2 * z0)}, where z0 is the proportion of
#'      \code{sampleMC < meanMC}.  May be a better approximation of the p-value
#'      than \code{simpleP}, but many of the same limitations apply.  NULL when
#'      \code{approxSigTest==FALSE}.
#'    }
#'   }
#'   \item{\code{corr}}{List containing estimates of rM, an alternate measure of
#'    migratory connectivity strength. NULL when \code{calcCorr==FALSE} or
#'    \code{!is.null(psi)}:
#'    \itemize{
#'     \item{\code{sample}} \code{nBoot} sampled values for continuous
#'      correlation. Provided to allow the user to compute own summary
#'      statistics.
#'     \item{\code{mean, se, simpleCI, bcCI, median, point}} Summary
#'      statistics for continuous correlation bootstraps.
#'    }
#'   }
#'   \item{\code{input}}{List containing the inputs to \code{estMC}, or at least
#'    the relevant ones, such as sampleSize.}
#' }
#' @example inst/examples/estMCExamples.R
#' @seealso \code{\link{estStrength}}, \code{\link{estTransition}},
#'   \code{\link{estMantel}}, \code{\link{calcMC}}, \code{\link{projections}},
#'   \code{\link{isoAssign}}, \code{\link{plot.estMigConnectivity}}
#' @export
#' @references
#' Cohen, E. B., J. A. Hostetler, M. T. Hallworth, C. S. Rushing, T. S. Sillett,
#' and P. P. Marra. 2018. Quantifying the strength of migratory connectivity.
#' Methods in Ecology and Evolution 9: 513 - 524.
#' \doi{10.1111/2041-210X.12916}
#'
#' Cohen, E. B., C. S. Rushing, F. R. Moore, M. T. Hallworth, J. A. Hostetler,
#' M. Gutierrez Ramirez, and P. P. Marra. 2019. The strength of
#' migratory connectivity for birds en route to breeding through the Gulf of
#' Mexico. Ecography 42: 658–669.
#' \doi{10.1111/ecog.03974}

estMC <- function(originDist, targetDist = NULL, originRelAbund, psi = NULL,
                  sampleSize = NULL,
                  originSites = NULL, targetSites = NULL,
                  originPoints = NULL, targetPoints = NULL,
                  originAssignment = NULL, targetAssignment = NULL,
                  originNames = NULL, targetNames = NULL,
                  nSamples = 1000, nSim = ifelse(isTRUE(isIntrinsic), 10, 1000),
                  isGL = FALSE, geoBias = NULL, geoVCov = NULL, row0 = 0,
                  verbose = 0, calcCorr = FALSE, alpha = 0.05,
                  approxSigTest = FALSE, sigConst = 0,
                  resampleProjection = 'ESRI:102010',
                  maxTries = 300, targetIntrinsic = NULL,
                  isIntrinsic = FALSE,
                  maintainLegacyOutput = FALSE) {
  if (is.null(psi)) {
    if(isIntrinsic) {
      mc <- estMCisotope(targetDist = targetDist,
                         originRelAbund = originRelAbund,
                         targetIntrinsic = targetIntrinsic,
                         targetSites = targetSites, sampleSize = sampleSize,
                         # targetAssignment = targetAssignment,
                         originPoints=originPoints, originSites=originSites,
                         originDist = originDist,
                         originAssignment = originAssignment,
                         originNames = originNames, targetNames = targetNames,
                         nBoot = nSamples, verbose = verbose, nSim = nSim,
                         calcCorr=calcCorr, alpha = alpha,
                         approxSigTest = approxSigTest, sigConst = sigConst,
                         maxTries = maxTries,
                         maintainLegacyOutput = maintainLegacyOutput)
    }
    else #if (is.null(nReleased)) {
      mc <- estMCGlGps(isGL=isGL, geoBias=geoBias, geoVCov=geoVCov,
                       originRelAbund=originRelAbund, sampleSize = sampleSize,
                       originDist=originDist,
                       targetDist=targetDist,
                       targetPoints=targetPoints, targetSites=targetSites,
                       targetAssignment=targetAssignment,
                       originPoints=originPoints, originSites=originSites,
                       originAssignment=originAssignment,
                       originNames=originNames, targetNames=targetNames,
                       nBoot = nSamples, verbose=verbose,
                       nSim = nSim, calcCorr=calcCorr, alpha = alpha,
                       approxSigTest = approxSigTest, sigConst = sigConst,
                       resampleProjection = resampleProjection,
                       maxTries = maxTries,
                       maintainLegacyOutput = maintainLegacyOutput)
  }
  else {
    mc <- estStrength(originRelAbund = originRelAbund,
                        sampleSize = sampleSize, psi = psi,
                        originDist = originDist, targetDist = targetDist,
                        originSites=originSites, targetSites=targetSites,
                        nSamples = nSamples, row0 = row0, verbose=verbose,
                        alpha = alpha, approxSigTest = approxSigTest,
                        sigConst = sigConst, originNames = originNames,
                        targetNames = targetNames,
                        maintainLegacyOutput = maintainLegacyOutput)
  }
  if (calcCorr && is.null(psi))
    class(mc) <- c("estMC", "estPsi", "estMantel", "estMigConnectivity")
  else if (is.null(psi))
    class(mc) <- c("estMC", "estPsi", "estMigConnectivity")
  else
    class(mc) <- c("estMC", "estMigConnectivity")
  return(mc)
}

#' Estimate Mantel correlation (rM) from geolocator, GPS, and/or raster data.
#'
#' Resampling of uncertainty for migratory connectivity strength, as quantified
#' by Mantel correlation (rM), from geolocators, GPS, and/or raster (e.g.,
#' genoscape or isotope) data.
#'
#' @param targetPoints A \code{POINTS} from sf
#'  object, with length number of animals tracked.  Each point indicates the
#'  point estimate location in the non-release season.
#' @param originPoints A \code{POINTS} from sf
#'  object, with length number of animals tracked.  Each point indicates the
#'  release location of an animal.
#' @param isGL Indicates whether or which animals were tracked with geolocators
#'  Should be either single TRUE or FALSE value, or vector with length of
#'  number of animals tracked, with TRUE for animals in  \code{targetPoints}
#'  with geolocators and FALSE for animals without.
#' @param geoBias For GL data, vector of length 2 indicating expected bias
#'  in longitude and latitude of \code{targetPoints}, in
#'  \code{resampleProjection} units (default meters).
#' @param geoVCov For GL data, 2x2 matrix with expected variance/covariance
#'  in longitude and latitude of \code{targetPoints}, in
#'  \code{resampleProjection} units (default meters).
#' @param targetSites A \code{SpatialPolygons}, \code{SpatialPolygonsDataFrame},
#'  or \code{POLYGONS} sf object indicating valid target location(s). Not
#'  needed unless you want to mask out certain areas (e.g. water) and
#'  \code{captured} is "origin" or you want to use a weighted bootstrap based on
#'  \code{targetRelAbund} for animals captured on the target side.
#' @param nBoot Number of bootstrap runs. Animals are sampled with replacement
#'  for each, to estimate sampling uncertainty.
#' @param nSim Tuning parameter for GL or raster data. Affects only the speed;
#'  1000 seems to work well with our GL data.  Should be integer > 0.
#' @param verbose 0 (default) to 3. 0 prints no output during run. 1 prints
#'  a line every 100 bootstraps.  2 prints a line every bootstrap.
#'  3 also prints the number of draws (for tuning nSim only).
#' @param alpha Level for confidence/credible intervals provided.
#' @param resampleProjection Projection when sampling from geolocator
#'  bias/error. This projection needs units = m. Default is Equidistant
#'  Conic. The default setting preserves distances around latitude = 0 and
#'  longitude = 0. Other projections may work well, depending on the location
#'  of \code{targetPoints}.
#' @param maxTries Maximum number of times to run a single GL bootstrap before
#'  exiting with an error.  Default is 300.  Set to NULL to never stop.  This
#'  parameter was added to prevent GL setups where some sample points never
#'  land on target sites from running indefinitely.
#' @param maintainLegacyOutput version 0.4.0 of \code{MigConnectivity}
#'  updated the structure of the estimates. If you have legacy code that refers
#'  to elements within a \code{estMigConnectivity} object, you can set this
#'  to TRUE to also keep the old structure. Defaults to FALSE.
#' @param originSites A \code{SpatialPolygons}, \code{SpatialPolygonsDataFrame},
#'  or \code{POLYGONS} sf object indicating valid origin location(s). Not
#'  needed unless you want to mask out certain areas (e.g. water) and
#'  \code{captured} is "target" or you want to use a weighted bootstrap based on
#'  \code{originRelAbund} for animals captured on the origin side.
#' @param isTelemetry Indicates whether or which animals were tracked with
#'  telemetry/GPS (no location uncertainty on either end).
#'  Should be either single TRUE or FALSE value, or vector with length of
#'  number of animals tracked, with TRUE or FALSE for each animal in data.
#' @param isRaster Indicates whether or which animals were tracked with
#'  intrinsic markers (e.g., genetics or isotopes), with location uncertainty
#'  expressed as a raster of probabilities by grid cells, either in
#'  \code{targetRaster} or \code{originRaster}. Should be either single TRUE or
#'  FALSE value, or vector with length of number of animals tracked, with TRUE
#'  or FALSE for each animal in data.
#' @param captured Indicates whether or which animals were captured in the
#'  origin sites, the target sites, or neither (another phase of the annual
#'  cycle). Location uncertainty will only be applied where the animal was not
#'  captured. So this doesn't matter for telemetry data. Should be either single
#'  "origin" (default), "target", or "neither" value, or a character vector with
#'  length of number of animals tracked, with "origin", "target", or "neither"
#'  for each animal.
#' @param geoBiasOrigin For GL data where \code{captured}!="origin", vector of
#'  length 2 indicating expected bias in longitude and latitude of
#'  \code{originPoints}, in \code{resampleProjection} units (default meters).
#' @param geoVCovOrigin For GL data where \code{captured}!="origin", 2x2 matrix
#'  with expected variance/covariance in longitude and latitude of
#'  \code{targetPoints}, in \code{resampleProjection} units (default meters).
#' @param targetRaster For intrinsic tracking data, the results of
#'  \code{isoAssign} or a similar function of class \code{intrinsicAssign} or
#'  class \code{RasterBrick}/\code{RasterStack}, for example from the package
#'  \code{assignR}. In any case, it expresses location uncertainty on target
#'  range, through a raster of probabilities by grid cells.
#' @param originRaster For intrinsic tracking data, the results of
#'  \code{isoAssign} or a similar function of class \code{intrinsicAssign} or
#'  class \code{RasterBrick}/\code{RasterStack}, for example from the package
#'  \code{assignR}. In any case, it expresses location uncertainty on origin
#'  range, through a raster of probabilities by grid cells.
#' @param dataOverlapSetting When there is more than one type of data, this
#'  setting allows the user some flexibility for clarifying which type(s) of
#'  data apply to which animals. Setting "dummy" (the default) indicates that
#'  there are dummy values within each dataset for the animals that isGL,
#'  isTelemetry, etc. don't have that data type (FALSE values). If no animals
#'  have a data type, no dummy values are required. If no animals have more than
#'  one type of data, the user can simplify processing their data by choosing
#'  setting "none" here. In this case, there should be no dummy values, and only
#'  the animals with a type of data should be included in that dataset. The
#'  third setting ("named") is not yet implemented, but will eventually allow
#'  another way to allow animals with more than one type of data with named
#'  animals linking records. When there is only one type of data, it is fastest
#'  to leave this on the default.
#' @param originRelAbund the proportion of the total abundance in each of B
#'  \code{originSites}. Used to set up the bootstrap to be weighted by relative
#'  abundance (for animals captured on the origin side). Either a numeric vector
#'  of length B that sums to 1, or an mcmc object (such as is produced by
#'  \code{\link{modelCountDataJAGS}}) or matrix with at least B columns.
#'  If there are more than B columns, the relevant columns should be
#'  labeled "relN[1]" through "relN[B]". Optional, but if you don't set it and
#'  at least some animals are captured on the origin side, there's potential for
#'  rM to be biased (if sampling isn't proportional to abundance).
#' @param targetRelAbund the proportion of the total abundance in each of W
#'  \code{targetSites}. Used to set up the bootstrap to be weighted by relative
#'  abundance (for animals captured on the target side). Either a numeric vector
#'  of length W that sums to 1, or an mcmc object (such as is produced by
#'  \code{\link{modelCountDataJAGS}}) or matrix with at least W columns.
#'  If there are more than W columns, the relevant columns should be
#'  labeled "relN[1]" through "relN[W]". Optional, but if you don't set it and
#'  at least some animals are captured on the target side, there's potential for
#'  rM to be biased (if sampling isn't proportional to abundance).
#'
#' @return \code{estMantel} returns a list with elements:
#' \describe{
#'   \item{\code{corr}}{List containing estimates of rM:
#'    \itemize{
#'     \item{\code{sample}} \code{nBoot} sampled values for Mantel
#'      correlation. Provided to allow the user to compute own summary
#'      statistics.
#'     \item{\code{mean, se, simpleCI, bcCI, median, point}} Summary
#'      statistics for Mantel correlation bootstraps.
#'    }
#'   }
#'   \item{\code{input}}{List containing the inputs to \code{estMantel}}
#' }
#' @export
#'
#' @examples
#' \donttest{
#' data('OVENdata')
#' rM1 <- estMantel(isGL=OVENdata$isGL,#Logical vector: light-level GL(T)/GPS(F)
#'                  geoBias = OVENdata$geo.bias, # Geolocator location bias
#'                  geoVCov = OVENdata$geo.vcov, # Location covariance matrix
#'                  targetSites = OVENdata$targetSites,#Nonbreeding/target sites
#'                  originPoints = OVENdata$originPoints, # Capture Locations
#'                  targetPoints = OVENdata$targetPoints, # Target locations
#'                  verbose = 1,   # output options
#'                  nBoot = 10, # This is set low for example
#'                  resampleProjection = sf::st_crs(OVENdata$targetSites))
#' rM1
#' str(rM1, max.level = 2)
#' }
#' @seealso \code{\link{estMC}}
#'
#' @references
#' Cohen, E. B., J. A. Hostetler, M. T. Hallworth, C. S. Rushing, T. S. Sillett,
#' and P. P. Marra. 2018. Quantifying the strength of migratory connectivity.
#' Methods in Ecology and Evolution 9: 513 - 524.
#' \doi{10.1111/2041-210X.12916}

estMantel <- function(targetPoints = NULL, originPoints = NULL, isGL,
                      geoBias = NULL, geoVCov = NULL, targetSites = NULL,
                      nBoot = 1000,
                      nSim = ifelse(any(isRaster & isGL), 5000,
                                    ifelse(any(isGL), 1000,
                                           ifelse(any(isRaster), 10, 1))),
                      verbose=0, alpha = 0.05,
                      resampleProjection = 'ESRI:102010',
                      maxTries = 300, maintainLegacyOutput = FALSE,
                      originSites = NULL, isTelemetry = !isGL, isRaster = FALSE,
                      captured = "origin", geoBiasOrigin = geoBias,
                      geoVCovOrigin = geoVCov,
                      targetRaster = NULL, originRaster = NULL,
                      dataOverlapSetting = c("dummy", "none", "named"),
                      originRelAbund = NULL, targetRelAbund = NULL) {
  dataOverlapSetting <- match.arg(dataOverlapSetting)
  # double check that spatial data coming in is in sf format #
  if (inherits(targetPoints, "SpatialPoints"))
    targetPoints <- sf::st_as_sf(targetPoints)
  if (inherits(originPoints, "SpatialPoints"))
    originPoints <- sf::st_as_sf(originPoints)

  # Input checking and assignment
  if (!(verbose %in% 0:3))
    stop("verbose should be integer 0-3 for level of output during bootstrap: 0 = none, 1 = every 10, 2 = every run, 3 = every animal")
  if (length(geoBias)!=2 && any(isGL) && any(captured != "target"))
    stop("geoBias should be vector of length 2 (expected bias in longitude and latitude of targetPoints, in meters)")
  if (!isTRUE(all.equal(dim(geoVCov), c(2, 2), check.attributes = FALSE)) &&
      any(isGL) && any(captured != "target"))
    stop("geoVCov should be 2x2 matrix (expected variance/covariance in longitude and latitude of targetPoints, in meters)")

  targetStats <- assignRasterStats(targetRaster)
  targetPointsAssigned <- targetStats$PointsAssigned
  targetSingleCell <- targetStats$SingleCell
  targetRasterXYZ <- targetStats$RasterXYZ
  targetRasterXYZcrs <- targetStats$RasterXYZcrs

  originStats <- assignRasterStats(originRaster)
  originPointsAssigned <- originStats$PointsAssigned
  originSingleCell <- originStats$SingleCell
  originRasterXYZ <- originStats$RasterXYZ
  originRasterXYZcrs <- originStats$RasterXYZcrs

  if (is.null(targetPoints)) {
    if (any(isGL) || any(isTelemetry) || !all(captured=="origin") ||
        is.null(targetRaster))
      stop("Need to define targetPoints unless all data is raster captured on origin")
  }
  else {
    if(is.na(sf::st_crs(targetPoints))){
      stop('Coordinate system definition needed for targetPoints')
    }
    targetPoints <- sf::st_transform(targetPoints, crs = resampleProjection)
  }

  if (is.null(originPoints)) {
    if (any(isGL) || any(isTelemetry) || any(captured=="origin") ||
        is.null(originRaster))
      stop("Need to define originPoints unless all data is raster captured on target")
  }
  else if (is.na(sf::st_crs(originPoints))) {
    stop('Coordinate system definition needed for originPoints')
  }
  if (dataOverlapSetting != "dummy") {
    temp <- reassignInds(dataOverlapSetting = dataOverlapSetting,
                         originPoints = originPoints,
                         targetPoints = targetPoints,
                         isGL = isGL, isTelemetry = isTelemetry,
                         isRaster = isRaster, isProb = FALSE,
                         captured = captured,
                         originRasterXYZ = originRasterXYZ,
                         #originRasterXYZcrs = originRasterXYZcrs,
                         originSingleCell = originSingleCell,
                         targetRasterXYZ = targetRasterXYZ,
                         #targetRasterXYZcrs = targetRasterXYZcrs,
                         targetSingleCell = targetSingleCell)
    originPoints <- temp$originPoints; targetPoints <- temp$targetPoints
    # originAssignment <- temp$originAssignment
    # targetAssignment <- temp$targetAssignment
    isGL <- temp$isGL; isTelemetry <- temp$isTelemetry
    isRaster <- temp$isRaster; isProb <- temp$isProb
    originRasterXYZ <- temp$originRasterXYZ
    originSingleCell <- temp$originSingleCell
    targetRasterXYZ <- temp$targetRasterXYZ
    targetSingleCell <- temp$targetSingleCell
  }

  nAnimals <- max(nrow(targetPoints), nrow(originPoints), length(isGL),
                  length(isTelemetry), length(isRaster),
                  ifelse(is.null(targetRaster), 0,
                         ifelse(targetPointsAssigned, dim(targetSingleCell)[3],
                                dim(targetRasterXYZ)[2] - 2)),
                  ifelse(is.null(originRaster), 0,
                         ifelse(originPointsAssigned, dim(originSingleCell)[3],
                                dim(originRasterXYZ)[2] - 2)),
                  length(captured))

  if (length(isGL)==1)
    isGL <- rep(isGL, nAnimals)
  if (length(isTelemetry)==1)
    isTelemetry <- rep(isTelemetry, nAnimals)
  if (length(isRaster)==1)
    isRaster <- rep(isRaster, nAnimals)
  if (length(captured)==1)
    captured <- rep(captured, nAnimals)

  if(!is.null(targetSites)){
    #if sp object covert to sf #
    if("SpatialPolygonsDataFrame" %in% class(targetSites)){
      targetSites <- sf::st_as_sf(targetSites)
    }

    if(is.na(sf::st_crs(targetSites))){
      stop('Coordinate system definition needed for targetSites')
    }

    targetSites <- sf::st_transform(targetSites, resampleProjection)
  }
  targetPointsInSites <- FALSE

  if (targetPointsAssigned && !is.null(targetSites) && any(isRaster)) {
    if (verbose > 0){
      cat('Checking if single cell target points in targetSites, may take a moment\n')}
    targetPointSample2 <- apply(targetSingleCell,
                                FUN = function(x)
                                  sf::st_as_sf(data.frame(x),
                                               coords = c("Longitude",
                                                          "Latitude"),
                                               crs = 4326),
                                MARGIN = 3)
    if(!sf::st_crs(targetSites)==sf::st_crs(targetPointSample2[[1]])){
      targetPointSample2 <- sapply(targetPointSample2, sf::st_transform,
                                   crs = resampleProjection)
    }

    targetCon <- sapply(targetPointSample2, FUN = function(z){
      suppressMessages(as.numeric(unclass(sf::st_intersects(x = z,
                                                            y = targetSites,
                                                            sparse = TRUE))))})

    if (!any(is.na(targetCon)))
      targetPointsInSites <- TRUE
    else if (verbose > 0)
      cat('Single cell target points supplied, but some points (proportion',
          format(sum(is.na(targetCon))/length(targetCon), digits = 2),
          ') not in targetSites\n')
  }
  else {
    targetCon <- NULL
  }

  originPointsInSites <- FALSE
  if (originPointsAssigned && !is.null(originSites) && any(isRaster)) {
    if (verbose > 0)
      cat('Checking if single cell origin points in originSites, may take a moment\n')
    nSamples <- dim(originSingleCell)[1]
    originPointSample2 <- apply(originSingleCell,
                                FUN = function(x){sf::st_as_sf(data.frame(x),
                                                               coords = c("Longitude", "Latitude"),
                                                               crs = 4326)},
                                MARGIN = 3)
    if(!sf::st_crs(originSites)==sf::st_crs(originPointSample2[[1]])){
      originPointSample2 <- sapply(originPointSample2, sf::st_transform,
                                   crs = resampleProjection)
    }

    originCon <- sapply(originPointSample2, FUN = function(z){
      suppressMessages(as.numeric(unclass(sf::st_intersects(x = z, y = originSites,
                                                            sparse = TRUE))))})

    if (!any(is.na(originCon)))
      originPointsInSites <- TRUE
    else if (verbose > 0){
      cat('Single cell origin points supplied, but some points (proportion',
          sum(is.na(originCon))/length(originCon), ') not in originSites\n')
    }
  }else{
    originCon <- NULL
  }

  weights <- array(0, c(nBoot, nAnimals))
  # if (is.null(originRelAbund) && any(captured!="target") ||
  #     is.null(targetRelAbund) && any(captured!="origin")) {
  #   warning("rM can be biased if you don't weight by abundance")
  # }
  if (is.null(originSites) && !is.null(originRelAbund)) {
    stop("If you want to set origin relative abundances, you need to define the origin sites")
  }
  if (is.null(targetSites) && !is.null(targetRelAbund)) {
    stop("If you want to set target relative abundances, you need to define the target sites")
  }
  if (is.null(originRelAbund) && is.null(targetRelAbund))
    weights[,] <- 1/nAnimals
  else if (!is.null(originRelAbund) && any(captured!="target")) {
    nOriginSites <- nrow(originSites)
    if (coda::is.mcmc(originRelAbund) || coda::is.mcmc.list(originRelAbund)) {
      originRelAbund <- as.matrix(originRelAbund)
    }
    if (is.matrix(originRelAbund) && dim(originRelAbund)[1]>1) {
      abundFixed <- FALSE
      if (dim(originRelAbund)[2]>nOriginSites)
        abundParams <- paste('relN[', 1:nOriginSites, ']', sep='')
      else if (dim(originRelAbund)[2]==nOriginSites)
        abundParams <- 1:nOriginSites
      else
        stop('Number of origin sites must be constant between sites and abundance')
      if (dim(originRelAbund)[1] >= nBoot)
        abundRows <- round(seq(from = 1, to = dim(originRelAbund)[1],
                               length.out = nBoot))
      else
        abundRows <- sample.int(n = dim(originRelAbund)[1], replace = TRUE,
                                size = nBoot)
      originRelAbund <- as.matrix(originRelAbund[abundRows, abundParams])
      abundBase <- colMeans(originRelAbund)
    }
    else {
      abundFixed <- TRUE
      if (length(originRelAbund)!=nOriginSites)
        stop('Number of origin sites must be constant between sites and abundance')
      abundBase <- originRelAbund
      originRelAbund <- matrix(originRelAbund, nBoot, nOriginSites, TRUE)
    }
    if (verbose > 0)
      cat("Creating origin assignment\n")
    # if geolocator, telemetry and captured in origin then simply get the origin site
    if (!is.null(originPoints)){
      if(!identical(sf::st_crs(originPoints),sf::st_crs(originSites))){
        # project if needed
        originSites <- sf::st_transform(originSites, sf::st_crs(originPoints))
      }
      originAssignment <- suppressMessages(unclass(sf::st_intersects(x = originPoints,
                                                                     y = originSites,
                                                                     sparse = TRUE)))
    }
    else
      originAssignment <- NULL
    if (!is.null(originAssignment)) {
      originAssignment[lengths(originAssignment)==0] <- NA
      if (any(lengths(originAssignment)>1)){
        warning("Overlapping originSites may cause issues\n")
        originAssignment <- lapply(originAssignment, function (x) x[1])
      }
      originAssignment <- array(unlist(originAssignment))
      duds <- is.na(originAssignment) & captured[1:nAnimals] != "target"
      if (any(duds)){
        if (verbose > 0)
          cat("Not all origin capture locations are within originSites. Assigning to closest site\n")
        warning("Not all origin capture locations are within originSites. Assigning to closest site.\n",
                "Affects animals: ", paste(which(duds), collapse = ","))
        originAssignment[duds] <-
          sf::st_nearest_feature(x = originPoints[duds,],
                                 y = originSites)

      }
      nOriginAnimals <- rep(NA, nOriginSites)
      for (i in 1:nOriginSites) {
        nOriginAnimals[i] <- sum(originAssignment==i & captured!="target", na.rm = TRUE)
        if (nOriginAnimals[i] > 0)
          weights[ , originAssignment==i] <- originRelAbund[, i]/nOriginAnimals[i]
      }
    }
    else{
      warning("Can't assign animals to origin sites; using unweighted sampling",
              immediate. = TRUE)
      weights[ , captured!="target"] <- 1/sum(captured!="target")
    }
  }
  if (!is.null(targetRelAbund) && any(captured=="target")) {
    nTargetSites <- nrow(targetSites)
    if (coda::is.mcmc(targetRelAbund) || coda::is.mcmc.list(targetRelAbund)) {
      targetRelAbund <- as.matrix(targetRelAbund)
    }
    if (is.matrix(targetRelAbund) && dim(targetRelAbund)[1]>1) {
      abundFixedT <- FALSE
      if (dim(targetRelAbund)[2]>nTargetSites)
        abundParams <- paste('relN[', 1:nTargetSites, ']', sep='')
      else if (dim(targetRelAbund)[2]==nTargetSites)
        abundParams <- 1:nTargetSites
      else
        stop('Number of target sites must be constant between sites and abundance')
      if (dim(targetRelAbund)[1] >= nBoot)
        abundRows <- round(seq(from = 1, to = dim(targetRelAbund)[1],
                               length.out = nBoot))
      else
        abundRows <- sample.int(n = dim(targetRelAbund)[1], replace = TRUE,
                                size = nBoot)
      targetRelAbund <- as.matrix(targetRelAbund[abundRows, abundParams])
    }
    else {
      abundFixedT <- TRUE
      if (length(targetRelAbund)!=nTargetSites)
        stop('Number of target sites must be constant between sites and abundance')
      targetRelAbund <- matrix(targetRelAbund, nBoot, nTargetSites, TRUE)
    }
    if (verbose > 0)
      cat("Creating target assignment\n")
    # if geolocator, telemetry and captured in target then simply get the target site
    if (!is.null(targetPoints)){
      if(!identical(sf::st_crs(targetPoints),sf::st_crs(targetSites))){
        # project if needed
        targetSites <- sf::st_transform(targetSites, sf::st_crs(targetPoints))
      }
      targetAssignment <- suppressMessages(unclass(sf::st_intersects(x = targetPoints,
                                                                     y = targetSites,
                                                                     sparse = TRUE)))
    }
    else
      targetAssignment <- NULL
    if (!is.null(targetAssignment)) {
      targetAssignment[lengths(targetAssignment)==0] <- NA
      if (any(lengths(targetAssignment)>1)){
        warning("Overlapping targetSites may cause issues\n")
        targetAssignment <- lapply(targetAssignment, function (x) x[1])
      }
      targetAssignment <- array(unlist(targetAssignment))
      duds <- is.na(targetAssignment) & captured[1:nAnimals] == "target"
      if (any(duds)){
        if (verbose > 0)
          cat("Not all target capture locations are within targetSites. Assigning to closest site\n")
        warning("Not all target capture locations are within targetSites. Assigning to closest site.\n",
                "Affects animals: ", paste(which(duds), collapse = ","))
        targetAssignment[duds] <-
          sf::st_nearest_feature(x = targetPoints[duds,],
                                 y = targetSites)

      }
      nTargetAnimals <- rep(NA, nTargetSites)
      for (i in 1:nTargetSites) {
        nTargetAnimals[i] <- sum(targetAssignment==i & captured=="target", na.rm = TRUE)
        if (nTargetAnimals[i] > 0)
          weights[ , targetAssignment==i] <- targetRelAbund[, i]/nTargetAnimals[i]
      }
    }
    else{
      warning("Can't assign animals to targetsites; using unweighted sampling",
              immediate. = TRUE)
      weights[ , captured=="target"] <- 1/sum(captured=="target")
    }
  }

  corr <- rep(NA, nBoot)
  if (!is.null(targetPoints) && !is.null(originPoints)){
    pointMantel <- calcMantel(targetPoints = targetPoints,
                              originPoints = originPoints)
    targetDistStart <- pointMantel$targetDist
    originDistStart <- pointMantel$originDist
    pointCorr <- pointMantel$pointCorr
  }
  else {
    pointCorr <- NA
    if (!is.null(targetPoints)){
      targetPoints2 <- sf::st_transform(targetPoints,4326)
      targetDistStart <- distFromPos(sf::st_coordinates(targetPoints2))
    }
    else { # Can't both be null
      originPoints2 <- sf::st_transform(originPoints,4326)

      originDistStart <- distFromPos(sf::st_coordinates(originPoints2))
    }
  }

  for (boot in 1:nBoot) {
    if (verbose > 1 || verbose == 1 && boot %% 100 == 0)
      cat("Bootstrap Run", boot, "of", nBoot, "at", date(), "\n")
    # Sample individual animals with replacement
    animal.sample <- sample.int(nAnimals, replace=TRUE, prob = weights[boot, ])
    if (any(captured[animal.sample]!='origin')) {
      oSamp <- locSample(isGL = (isGL[animal.sample] & captured[animal.sample]!='origin'),
                         isRaster = (isRaster[animal.sample] & captured[animal.sample]!='origin'),
                         isProb = rep(FALSE, nAnimals),
                         isTelemetry = (isTelemetry[animal.sample] | captured[animal.sample]=='origin'),
                         geoBias = geoBiasOrigin,
                         geoVCov = geoVCovOrigin,
                         points = originPoints[animal.sample, ],
                         matvals = originRasterXYZ[, c(1:2, animal.sample + 2)],
                         matvals_crs = originRasterXYZcrs,
                         singleCell = originSingleCell[,,animal.sample],
                         overlap1 = originCon[,animal.sample],
                         pointsInSites = originPointsInSites,
                         assignment = NULL,
                         sites = originSites,
                         resampleProjection = resampleProjection,
                         nSim = nSim,
                         maxTries = maxTries)
      if (!is.null(oSamp$notfind)) {
        oSamp$notfind$Animal <- animal.sample[oSamp$notfind$Animal]
        notfind <- unique(oSamp$notfind)
        stop('maxTries (',maxTries,') reached during origin location sampling, exiting. ',
             'Animal(s) where location sampling failed to fall in sites:\n',
             paste(utils::capture.output(print(notfind, row.names = FALSE)), collapse = "\n"),
             '\nExamine originSites',
             ifelse(any(notfind$isGL),
                    ', geoBiasOrigin, geoVcovOrigin, originPoints', ''),
             ifelse(any(notfind$isRaster), ', originRaster', ''),
             ifelse(any(notfind$isTelemetry), ', originPoints/captured', ''),
             ', and resampleProjection to determine why sampled points fell outside sites.')
      }
      origin.point.sample <- oSamp$point.sample
      origin.point.sample <- sf::st_as_sf(data.frame(origin.point.sample),
                                          coords = c("x","y"),
                                          crs = resampleProjection)
      if (verbose > 2)
        cat(' ', oSamp$draws, 'origin draw(s) (of length', nSim, 'and of', maxTries, 'possible).\n')
    }
    else {
      # Get origin point for each animal sampled
      origin.point.sample <- originPoints[animal.sample, ]
    }

    if (any(captured[animal.sample]!="target")) {
      tSamp <- locSample(isGL = (isGL[animal.sample] &
                                   captured[animal.sample] != "target"),
                         isRaster = (isRaster[animal.sample] &
                                       captured[animal.sample] != "target"),
                         isProb = rep(FALSE, nAnimals),
                         isTelemetry = (isTelemetry[animal.sample] |
                                          captured[animal.sample] == "target"),
                         geoBias = geoBias, geoVCov = geoVCov,
                         points = targetPoints[animal.sample, ],
                         matvals = targetRasterXYZ[, c(1:2, animal.sample + 2)],
                         matvals_crs = targetRasterXYZcrs,
                         singleCell = targetSingleCell[,,animal.sample],
                         pointsInSites = targetPointsInSites,
                         overlap1 = targetCon[, animal.sample],
                         sites = targetSites,
                         resampleProjection = resampleProjection, nSim = nSim,
                         maxTries = maxTries)
      if (!is.null(tSamp$notfind)) {
        tSamp$notfind$Animal <- animal.sample[tSamp$notfind$Animal]
        notfind <- unique(tSamp$notfind)
        stop('maxTries (',maxTries,') reached during target location sampling, exiting. ',
             'Animal(s) where location sampling failed to fall in sites:\n',
             paste(utils::capture.output(print(notfind, row.names = FALSE)), collapse = "\n"),
             '\nExamine targetSites',
             ifelse(any(notfind$isGL),
                    ', geoBiasOrigin, geoVcovOrigin, targetPoints', ''),
             ifelse(any(notfind$isRaster), ', targetRaster', ''),
             ifelse(any(notfind$isTelemetry), ', targetPoints/captured', ''),
             ', and resampleProjection to determine why sampled points fell outside sites.')
      }
      target.point.sample <- tSamp$point.sample
      target.point.sample <- sf::st_as_sf(data.frame(target.point.sample),
                                          coords = c("x","y"),
                                          crs = resampleProjection)
      if (verbose > 2){
        cat(' ', tSamp$draws, 'target draw(s) (of length', nSim, 'and of', maxTries, 'possible).\n')
      }
    }
    else {
      # Get target point for each animal sampled
      target.point.sample <- targetPoints[animal.sample, ]
    }

    if (all(captured[animal.sample]=="target")) {
      targetDist1 <- targetDistStart[animal.sample, animal.sample]
      corr[boot] <- calcMantel(targetDist = targetDist1,
                               originPoints = origin.point.sample)$pointCorr
    }
    else if (all(captured[animal.sample]=="origin")) {
      originDist1 <- originDistStart[animal.sample, animal.sample]
      corr[boot] <- calcMantel(originDist = originDist1,
                               targetPoints = target.point.sample)$pointCorr
    }
    else {
      corr[boot] <- calcMantel(originPoints = origin.point.sample,
                               targetPoints = target.point.sample)$pointCorr
    }

    if (verbose > 1 || verbose == 1 && boot %% 10 == 0)
      cat(" Correlation mean:", mean(corr, na.rm=TRUE), "SD:", sd(corr, na.rm=TRUE),
          "low quantile:", quantile(corr, alpha/2, na.rm=TRUE),
          "high quantile:", quantile(corr, 1-alpha/2, na.rm=TRUE), "\n")
  }
  meanCorr <- mean(corr, na.rm=TRUE)
  medianCorr <- median(corr, na.rm=TRUE)
  seCorr <- sd(corr, na.rm=TRUE)
  simpleCICorr <- quantile(corr, c(alpha/2, 1-alpha/2), na.rm=TRUE, type = 8,
                           names = FALSE)
  corr.z0 <- qnorm(sum((corr)<meanCorr)/nBoot)
  bcCICorr <- quantile(corr, pnorm(2*corr.z0+qnorm(c(alpha/2, 1-alpha/2))),
                       na.rm=TRUE, names = FALSE)
  if (maintainLegacyOutput)
    return(structure(list(sampleCorr = corr, pointCorr = pointCorr,
                          meanCorr = meanCorr, medianCorr = medianCorr,
                          seCorr=seCorr, simpleCICorr=simpleCICorr,
                          bcCICorr=bcCICorr, alpha = alpha,
                          corr = list(sample = corr, mean = meanCorr,
                                      se = seCorr, simpleCI = simpleCICorr,
                                      bcCI = bcCICorr, median = medianCorr,
                                      point = pointCorr),
                          input = list(targetPoints = targetPoints,
                                       originPoints = originPoints, isGL = isGL,
                                       geoBias = geoBias, geoVCov = geoVCov,
                                       targetSites = targetSites,
                                       targetSites = targetSites, nBoot = nBoot,
                                       nSim = nSim, verbose = verbose,
                                       alpha = alpha,
                                       resampleProjection = resampleProjection,
                                       maxTries = maxTries,
                                       maintainLegacyOutput = TRUE)),
                     class = c("estMantel", "estMigConnectivity")))
  else
    return(structure(list(corr = list(sample = corr, mean = meanCorr,
                                      se = seCorr, simpleCI = simpleCICorr,
                                      bcCI = bcCICorr, median = medianCorr,
                                      point = pointCorr),
                          input = list(targetPoints = targetPoints,
                                       originPoints = originPoints, isGL = isGL,
                                       geoBias = geoBias, geoVCov = geoVCov,
                                       targetSites = targetSites, nBoot = nBoot,
                                       nSim = nSim, verbose = verbose,
                                       alpha = alpha,
                                       resampleProjection = resampleProjection,
                                       maxTries = maxTries,
                                       maintainLegacyOutput = FALSE,
                                       originSites = originSites,
                                       isTelemetry = isTelemetry,
                                       isRaster = isRaster, captured = captured,
                                       geoBiasOrigin = geoBiasOrigin,
                                       geoVCovOrigin = geoVCovOrigin,
                                       targetRaster = targetRaster,
                                       originRaster = originRaster,
                                       dataOverlapSetting = dataOverlapSetting,
                                       originRelAbund = originRelAbund,
                                       targetRelAbund = targetRelAbund)),
                     class = c("estMantel", "estMigConnectivity")))
}

###############################################################################
#' Grab (from https://github.com/SMBC-NZP/MigConnectivity) example RMark
#' transition probability estimates obtained from simulated data
#'
#' Get a dataset containing RMark transition probability estimates from
#' simulated mark-recapture-recovery data from Cohen et al. (2014).  These all
#' represent the intermediate scenario for all settings (moderate connectivity,
#' low re-encounter, 100,000 banded in each breeding area).  Each estimate can
#' be used in \code{estMC} function to estimate MC with uncertainty.  Requires
#' internet connection.
#'
#' @param number Integer 1 - 100, which simulation and RMark estimate you want
#'
#' @return RMark object
#' @export
#' @seealso \code{\link{estMC}}
getCMRexample <- function(number = 1) {
  obj.name <- paste0('psiB.enc2.band100.', number)
  file.name <- paste0('out_', obj.name, '.rds')
  url1 <- paste0('https://github.com/SMBC-NZP/MigConnectivity/blob/master/data-raw/', file.name, '?raw=true')
  temp <- paste(tempdir(), file.name, sep = '/')
  utils::download.file(url1, temp, mode = 'wb')
  fm <- readRDS(temp)
  unlink(temp)
  return(fm)
}


#' Pairwise differences between two or more independent MC estimates
#'
#' Estimates mean (and median) differences in MC, and includes measures of
#' uncertainty (SE and CI).  For those measures of uncertainty to be accurate,
#' only apply this function to MC estimates where all data sources are
#' independent (e.g., different species).
#'
#' @param estimates List of at least two MC estimates, provided by the estMC
#'    function. If this is a named list (recommended), the function will use
#'    these names in labeling the differences.
#' @param nSamples A positive integer, number of samples (with replacement)
#'    to draw from each pair of MC estimates (default 100000).  If set to NULL,
#'    compares all MC samples from each pair.
#' @param alpha Level for confidence/credible intervals provided.
#' @param returnSamples Should the function return all the sampled differences?
#'    Defaults to FALSE to reduce storage requirements. Change to TRUE to
#'    compute your own summary statistics.
#'
#' @return  \code{diffMC} returns a list with elements:
#' \describe{
#'   \item{\code{meanDiff, medianDiff}}{Vectors with mean and medians of sampled
#'      differences for each pairwise comparison. Estimates of difference
#'      between MC values incorporating parametric uncertainty.}
#'   \item{\code{seDiff}}{Vector with standard errors of MC differences for each
#'      pairwise comparison, estimated from SD of sampled differences.}
#'   \item{\code{simpleCI}}{Matrix of \code{1 - alpha} confidence intervals for
#'      MC differences, estimated as \code{alpha/2} and \code{1 - alpha/2}
#'      quantiles of \code{sampleMC}.}
#'   \item{\code{bcCI}}{Matrix of bias-corrected \code{1 - alpha} confidence
#'      intervals for MC differences for each pairwise comparison. Preferable
#'      to \code{simpleCI} when \code{meanDiff} is the best estimate of the MC
#'      difference. \code{simpleCI} is preferred when
#'      \code{medianDiff} is a better estimator. When \code{meanDiff==medianDiff},
#'      these should be identical.  Estimated as the
#'      \code{pnorm(2 * z0 + qnorm(alpha / 2))} and
#'      \code{pnorm(2 * z0 + qnorm(1 - alpha / 2))} quantiles of sampled
#'      differences, where z0 is the proportion of \code{sampleDiff < meanDiff}.}
#   \item{\code{hpdCI}}{Matrix of \code{1 - alpha} credible intervals for MC
#      differences for each pairwise comparison, etimated using the highest
#      posterior density (HPD) method.}
#'   \item{\code{sampleDiff}}{Only provided if \code{returnSamples} is TRUE.
#'      List of sampled values for each pairwise MC difference.}
#' }
#' @export
#'
#' @references
#' Cohen, E. B., C. S. Rushing, F. R. Moore, M. T. Hallworth, J. A. Hostetler,
#' M. Gutierrez Ramirez, and P. P. Marra. 2019. The strength of
#' migratory connectivity for birds en route to breeding through the Gulf of Mexico.
#'
#' @examples
#' \donttest{
#' data('OVENdata')
#' ovenPsi <- estTransition(isGL = OVENdata$isGL, #Logical vector:light-level GL(T)
#'                  isTelemetry = !OVENdata$isGL,
#'                  geoBias = OVENdata$geo.bias, # Light-level GL location bias
#'                  geoVCov = OVENdata$geo.vcov, # Location covariance matrix
#'                  targetSites = OVENdata$targetSites, # Non-breeding target sites
#'                  originSites = OVENdata$originSites, # Breeding origin sites
#'                  originPoints = OVENdata$originPoints, # Capture Locations
#'                  targetPoints = OVENdata$targetPoints, # Device target locations
#'                  verbose = 0,   # output options
#'                  nSamples = 100, # This is set low for example
#'                  resampleProjection = sf::st_crs(OVENdata$targetSites))
#' ovenEst <- estStrength(targetDist = OVENdata$targetDist, # targetSites distance matrix
#'                  originDist = OVENdata$originDist, # originSites distance matrix
#'                  originRelAbund = OVENdata$originRelAbund,#Origin relative abund
#'                  psi = ovenPsi,
#'                  verbose = 1,   # output options
#'                  nSamples = 1000)
#' fm <- getCMRexample()
#' originPos13 <- matrix(c(rep(seq(-99, -81, 2), each = 10),
#'                         rep(seq(49, 31, -2), 10)), 100, 2)
#' targetPos13 <- matrix(c(rep(seq(-79, -61, 2), each = 10),
#'                         rep(seq(9, -9, -2), 10)), 100, 2)
#' originPosCMR <- rowsum(originPos13, c(rep(1:2, 5, each = 5),
#'                                       rep(3:4, 5, each = 5))) / 25
#' targetPosCMR <- rowsum(targetPos13, c(rep(1:2, 5, each = 5),
#'                                       rep(3:4, 5, each = 5))) / 25
#' originDist <- distFromPos(originPosCMR, 'ellipsoid')
#' targetDist <- distFromPos(targetPosCMR, 'ellipsoid')
#' originRelAbundTrue <- rep(0.25, 4)

#' theorEst <- estStrength(originRelAbund = originRelAbundTrue, psi = fm,
#'                   originDist = originDist, targetDist = targetDist,
#'                   originSites = 5:8, targetSites = c(3,2,1,4),
#'                   nSamples = 1000, verbose = 0,
#'                   sampleSize = length(grep("[2-5]", fm$data$data$ch)))
#' ovenEst
#' theorEst
#' diff1 <- diffMC(estimates = list(Ovenbird = ovenEst, Theorybird = theorEst),
#'                 nSamples = 10000, returnSamples = TRUE)
#'
#' }
diffMC <- function(estimates, nSamples = 100000, alpha = 0.05,
                   returnSamples = FALSE) {
  nEst <- length(estimates)
  nComparisons <- choose(nEst, 2)
  for (i in 1:nEst)
    if (is.null(estimates[[i]]$MC))
      estimates[[i]]$MC <- list(sample = estimates[[i]]$sampleMC)
  nSamplesEst <- sapply(estimates, function(x) length(x$MC$sample))
  diffSamples <- vector('list', nComparisons)
  if (is.null(names(estimates)))
    names(estimates) <- 1:nEst
  comparisons <- matrix(c(sequence(1:(nEst - 1)), rep(2:nEst, 1:(nEst - 1))),
                        nComparisons, 2)
  for (i in 1:nComparisons) {
    if (is.null(nSamples))
      diffSamples[[i]] <- rep(estimates[[comparisons[i, 1]]]$MC$sample,
                              times = nSamplesEst[comparisons[i, 2]]) -
        rep(estimates[[comparisons[i, 2]]]$MC$sample,
            each = nSamplesEst[comparisons[i, 1]])
    else
      diffSamples[[i]] <- sample(estimates[[comparisons[i, 1]]]$MC$sample,
                                 nSamples, replace = TRUE) -
        sample(estimates[[comparisons[i, 2]]]$MC$sample, nSamples, replace = TRUE)
  }
  meanDiff <- sapply(diffSamples, mean, na.rm=TRUE)
  medianDiff <- sapply(diffSamples, median, na.rm=TRUE)
  seDiff <- sapply(diffSamples, sd, na.rm=TRUE)
  simpleCI <- sapply(diffSamples, quantile, c(alpha/2, 1-alpha/2), na.rm=TRUE,
                     names = FALSE)
  diff.z0 <- sapply(diffSamples,
                    function(MC) qnorm(sum(MC<mean(MC, na.rm = TRUE),
                                           na.rm = TRUE)/length(which(!is.na(MC)))))
  bcCI <- mapply(function(MC, z0)
    quantile(MC, pnorm(2*z0+ qnorm(c(alpha/2, 1-alpha/2))),
                       na.rm=TRUE, names = FALSE), diffSamples, diff.z0)
  names(diffSamples) <- names(meanDiff) <- paste(names(estimates[comparisons[,1]]),
                                                 '-', names(estimates[comparisons[,2]]))
  names(medianDiff) <- names(seDiff) <- names(diffSamples)
  colnames(simpleCI) <- colnames(bcCI) <- names(diffSamples) #colnames(hpdCI) <-
  sampleDiff <- ifelse(returnSamples, diffSamples, NA)
  return(structure(list(meanDiff = meanDiff, medianDiff = medianDiff,
                        seDiff = seDiff, simpleCI = simpleCI, bcCI = bcCI,
                        sampleDiff = sampleDiff, alpha = alpha),
                   class = c('diffMC', 'diffMigConnectivity')))
}

#' Pairwise differences between two or more independent Mantel correlation
#' estimates
#'
#' Estimates mean (and median) differences in Mantel correlations (rM), and
#' includes measures of uncertainty (SE and CI).  For those measures of
#' uncertainty to be accurate, only apply this function to rM estimates where
#' all data sources are independent (e.g., different species).
#'
#' @param estimates List of at least two Mantel correlation estimates, provided
#'    by either the estMC or the estMantel functions. If this is a named list
#'    (recommended), the function will use these names in labeling the
#'    differences.
#' @param nSamples A positive integer, number of samples (with replacement)
#'    to draw from each pair of MC estimates (default 100000).  If set to NULL,
#'    compares all Mantel correlation samples from each pair.
#' @param alpha Level for confidence/credible intervals provided.
#' @param returnSamples Should the function return all the sampled differences?
#'    Defaults to FALSE to reduce storage requirements. Change to TRUE to
#'    compute your own summary statistics.
#'
#' @return  \code{diffMantel} returns a list with elements:
#' \describe{
#'   \item{\code{meanDiff, medianDiff}}{Vectors with mean and medians of sampled
#'      differences for each pairwise comparison. Estimates of difference
#'      between rM values incorporating parametric uncertainty.}
#'   \item{\code{seDiff}}{Vector with standard errors of rM differences for each
#'      pairwise comparison, estimated from SD of sampled differences.}
#'   \item{\code{simpleCI}}{Matrix of \code{1 - alpha} confidence intervals for
#'      rM differences, estimated as \code{alpha/2} and \code{1 - alpha/2}
#'      quantiles of \code{sampleCorr}.}
#'   \item{\code{bcCI}}{Matrix of bias-corrected \code{1 - alpha} confidence
#'      intervals for rM differences for each pairwise comparison. Preferable
#'      to \code{simpleCI} when \code{meanDiff} is the best estimate of the rM
#'      difference. \code{simpleCI} is preferred when
#'      \code{medianDiff} is a better estimator. When \code{meanDiff==medianDiff},
#'      these should be identical.  Estimated as the
#'      \code{pnorm(2 * z0 + qnorm(alpha / 2))} and
#'      \code{pnorm(2 * z0 + qnorm(1 - alpha / 2))} quantiles of sampled
#'      differences, where z0 is the proportion of \code{sampleDiff < meanDiff}.}
#   \item{\code{hpdCI}}{Matrix of \code{1 - alpha} credible intervals for rM
#      differences for each pairwise comparison, etimated using the highest
#      posterior density (HPD) method.}
#'   \item{\code{sampleDiff}}{Only provided if \code{returnSamples} is TRUE.
#'      List of sampled values for each pairwise rM difference.}
#' }
#' @export
#'
#' @references
#' Cohen, E. B., C. S. Rushing, F. R. Moore, M. T. Hallworth, J. A. Hostetler,
#' M. Gutierrez Ramirez, and P. P. Marra. 2019. The strength of
#' migratory connectivity for birds en route to breeding through the Gulf of Mexico.
#'
# @examples
diffMantel <- function(estimates, nSamples = 100000, alpha = 0.05,
                       returnSamples = FALSE) {
  nEst <- length(estimates)
  nComparisons <- choose(nEst, 2)
  for (i in 1:nEst)
    if (is.null(estimates[[i]]$corr))
      estimates[[i]]$corr <- list(sample = estimates[[i]]$sampleCorr)
  nSamplesEst <- sapply(estimates, function(x) length(x$corr$sample))
  diffSamples <- vector('list', nComparisons)
  if (is.null(names(estimates)))
    names(estimates) <- 1:nEst
  comparisons <- matrix(c(sequence(1:(nEst - 1)), rep(2:nEst, 1:(nEst - 1))),
                        nComparisons, 2)
  for (i in 1:nComparisons) {
    if (is.null(nSamples))
      diffSamples[[i]] <- rep(estimates[[comparisons[i, 1]]]$corr$sample,
                              times = nSamplesEst[comparisons[i, 2]]) -
        rep(estimates[[comparisons[i, 2]]]$corr$sample,
            each = nSamplesEst[comparisons[i, 1]])
    else
      diffSamples[[i]] <- sample(estimates[[comparisons[i, 1]]]$corr$sample,
                                 nSamples, replace = TRUE) -
        sample(estimates[[comparisons[i, 2]]]$corr$sample, nSamples,
               replace = TRUE)
  }
  meanDiff <- sapply(diffSamples, mean, na.rm=TRUE)
  medianDiff <- sapply(diffSamples, median, na.rm=TRUE)
  seDiff <- sapply(diffSamples, sd, na.rm=TRUE)
  simpleCI <- sapply(diffSamples, quantile, c(alpha/2, 1-alpha/2), na.rm=TRUE,
                     names = FALSE)
  diff.z0 <- sapply(diffSamples,
                    function(MC) qnorm(sum(MC<mean(MC, na.rm = TRUE), na.rm = TRUE)/
                                         length(which(!is.na(MC)))))
  bcCI <- mapply(function(MC, z0) quantile(MC,
                                           pnorm(2*z0+qnorm(c(alpha/2, 1-alpha/2))),
                       na.rm=TRUE, names = FALSE), diffSamples, diff.z0)
  names(diffSamples) <- names(meanDiff) <- paste(names(estimates[comparisons[,1]]),
                                                 '-', names(estimates[comparisons[,2]]))
  names(medianDiff) <- names(seDiff) <- names(diffSamples)
  colnames(simpleCI) <- colnames(bcCI) <- names(diffSamples)
  sampleDiff <- ifelse(returnSamples, diffSamples, NA)
  return(structure(list(meanDiff = meanDiff, medianDiff = medianDiff,
                        seDiff = seDiff, simpleCI = simpleCI, bcCI = bcCI,
                        sampleDiff = sampleDiff, alpha = alpha),
                   class = c('diffMantel', 'diffMigConnectivity')))
}

#' @rdname estTransition
#' @export
estPsi <- estTransition

#' @rdname estMantel
#' @export
estCorr <- estMantel

#' @rdname diffMC
#' @export
diffStrength <- diffMC

#' @rdname diffMantel
#' @export
diffCorr <- diffMantel
SMBC-NZP/MigConnectivity documentation built on March 26, 2024, 4:22 p.m.