## ----echo = FALSE-------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", tidy = FALSE)

## ----eval = FALSE-------------------------------------------------------------
#  install.packages('devtools')
#  
#  devtools::install_github("SMBC-NZP/MigConnectivity", build_vignettes = TRUE)
#  
#  library(MigConnectivity)

## ----warning=FALSE,message=FALSE,echo=FALSE-----------------------------------
library(MigConnectivity)

## -----------------------------------------------------------------------------
mlogitMat <- function(slope, dist) {
   preMat <- exp(-slope/mean(dist)*dist)
   diag(preMat) <- 0
   outMat <- t(apply(preMat, 1, function(x) x/(1 + sum(x))))
   diag(outMat) <- 1 - rowSums(outMat)
   return(outMat)
}


## -----------------------------------------------------------------------------
mlogitMC <- function(slope, 
                     MC.in, 
                     origin.dist,
                     target.dist,
                     origin.abund,
                     sample.size) {
  
    nBreeding <- nrow(origin.dist)
    nWintering <- nrow(target.dist)
    psi <- mlogitMat(slope, origin.dist)
    
  if (any(psi<0))
    
    return(5*slope^2)
    
  MC <- calcMC(origin.dist, 
               target.dist, 
               originRelAbund = origin.abund,
               psi, 
               sampleSize = sample.size)
  
  return((MC.in - MC)^2)
}

## -----------------------------------------------------------------------------
calcStrengthInd <- function(originDist, 
                            targetDist, 
                            locations, 
                            resamp=1000, 
                            verbose = 0) {
nInd <- dim(locations)[1]
originDist2 <- targetDist2 <- matrix(0, nInd, nInd)
for (i in 1:(nInd-1)) {
 for (j in (i+1):nInd) {
      originDist2[i,j] <- originDist[locations[i,1,1,1], locations[j,1,1,1]]
      targetDist2[i,j] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]]
      originDist2[j,i] <- originDist[locations[i,1,1,1], locations[j,1,1,1]]
      targetDist2[j,i] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]]
    }
  }
return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose))
} 

## -----------------------------------------------------------------------------
calcPsiMC <- function(originDist, 
                      targetDist, 
                      originRelAbund, 
                      locations, 
                      years = 1, 
                      months = 1, 
                      verbose=FALSE) {
  
  nOrigin <- nrow(originDist)
  nTarget <- nrow(targetDist)
  psiMat <- matrix(0, nOrigin, nTarget)
  nInd <- dim(locations)[1]
  nYears <- dim(locations)[3]
  nMonths <- dim(locations)[4]
for (i in 1:nInd) {
  if (i %% 1000 == 0 && verbose) #
      cat("Individual", i, "of", nInd, "\n")
    originMat <- locations[i, 1, years, months]
    targetMat <- locations[i, 2, years, months]
    bIndices <- which(!is.na(originMat))
    wIndices <- which(!is.na(targetMat))
  if (length(bIndices) && length(wIndices))
    for (bi in bIndices)
      for (wi in wIndices)
       psiMat[originMat[bi], targetMat[wi]] <- psiMat[originMat[bi], 
                                                      targetMat[wi]] + 1
  }
  psiMat <- apply(psiMat, 2, "/", rowSums(psiMat))
  MC <- calcMC(originDist, targetDist, psi = psiMat,
               originRelAbund = originRelAbund, sampleSize = nInd)
  return(list(psi=psiMat, MC=MC))
}

## -----------------------------------------------------------------------------
changeLocations <- function(animalLoc, 
                            breedingSiteTrans, 
                            winteringSiteTrans) {
  animalLoc[,1,,] <- breedingSiteTrans[animalLoc[,1,,]]
  animalLoc[,2,,] <- winteringSiteTrans[animalLoc[,2,,]]
  return(animalLoc)
}

## -----------------------------------------------------------------------------
simLocationError <- function(targetPoints, 
                             targetSites, 
              							 geoBias, 
              							 geoVCov, 
                             projection, 
              							 verbose = 0, 
              							 nSim = 1000) {
  
if(!inherits(targetPoints,"sf")){
  targetPoints <- sf::st_as_sf(targetPoints)}

nAnimals <- nrow(targetPoints)

geoBias2 <- matrix(rep(geoBias, nSim), nrow=nSim, byrow=TRUE)
target.sample <- rep(NA, nAnimals)
target.point.sample <- matrix(NA, nAnimals, 2)

for(i in 1:nAnimals){
  if (verbose > 0)
    cat('Animal', i, 'of', nAnimals)
  draws <- 0
  while (is.na(target.sample[i])) {
    draws <- draws + 1
    # Sample random point for each bird 
    # from parametric distribution of NB error
    ps <- MASS::mvrnorm(n=nSim, 
                  mu=cbind(sf::st_coordinates(targetPoints)[i,1],
                           sf::st_coordinates(targetPoints)[i,2]),
                  Sigma=geoVCov)+
                  geoBias2
    
    colnames(ps) <- c("Longitude","Latitude")
    
    point.sample <- sf::st_as_sf(
                  data.frame(ps),
                  coords = c("Longitude","Latitude"),
                  crs = sf::st_crs(targetPoints))
	  
    # filtered to stay in NB areas (land)
    
    target.sample0 <- unlist(unclass(sf::st_intersects(point.sample,
                                                       targetSites)))
    target.sample[i]<-target.sample0[!is.na(target.sample0)][1]
  }
	
  target.point.sample[i, ]<-sf::st_coordinates(point.sample[!is.na(target.sample0),])[1,]
  if (verbose > 0)
    cat(' ', draws, 'draws\n')
}
  return(list(targetSample = target.sample, 
              targetPointSample = target.point.sample))
}

## -----------------------------------------------------------------------------
toPoly <- function(siteCentroids,
                   projection.in = 4326,
                   projection.out = NA,
                   resolution = NA,
                   order.by.input = TRUE)
{
# This automatically sets the resolution so that all polygons touch and
# cover the entire surface. Alternatively the user can supply the 
# resolution of the cells in the units of the input projection (projection.in) 
  colnames(siteCentroids) <- c("Longitude","Latitude")
  
  if(is.na(resolution)){
     long <- unique(siteCentroids[,1])
     lat <- unique(siteCentroids[,2])
     if (length(long)>1)
       long.res <- abs(long[2]-long[1])
     else
       long.res <- 1
     if (length(lat)>1)
       lat.res <- abs(lat[2]-lat[1])
     else
       lat.res <- 1
     resolution <- c(long.res,lat.res)
  }
  centers <- sf::st_as_sf(data.frame(siteCentroids),
                          coords = c("Longitude","Latitude"),
                          crs = projection.in)
  polys <- sf::st_as_sf(sf::st_make_grid(centers,
                                 crs = projection.in, 
                                 cellsize = resolution, 
                                 offset = apply(siteCentroids, 2, min) - 
                                   resolution/2))
  

  if(order.by.input){
    # Reorders the polygons to match that of the input #
    orders <- suppressMessages(unlist(unclass(sf::st_intersects(centers,
                                                         polys,
                                                         sparse = TRUE))))
    polys <- polys[orders,] 
  }
  
  if(!is.na(projection.out)){
    polys <- sf::st_transform(polys, projection.out)
  }
  
  return(polys)
}

## -----------------------------------------------------------------------------
nBreeding <- 3 #number of breeding regions
nNonBreeding <- 3 #number of non-breeding regions

## -----------------------------------------------------------------------------
#distances between centroids of three breeding regions
breedDist <- matrix(c(0, 1, 2,
                      1, 0, 1,
                      2, 1, 0), nBreeding, nBreeding) 

#distances between centroids of three non-breeding regions
nonBreedDist <- matrix(c(0, 5, 10,
                         5, 0, 5,
                        10, 5, 0), nNonBreeding, nNonBreeding)

## ----echo=FALSE---------------------------------------------------------------
op <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(op))
par(bty="n")
plot(c(4,5,6),rep(1,3),
     ylim=c(-0.05,1.05),
     xlim=c(-1,11),
     pch=19,
     axes=FALSE,
     yaxt="n",
     xaxt="n",
     cex=2,
     ylab="",
     xlab="Non-breeding",
     main="Breeding")
points(c(0,5,10),rep(0,3),pch=15,cex=2)

## -----------------------------------------------------------------------------
#transition probabilities form each breeding to each non-breeding region
psi <- matrix(c(0.4, 0.35, 0.25,
                0.3, 0.4, 0.3,
                0.2, 0.3, 0.5), nBreeding, nNonBreeding, byrow = TRUE) 


## ----echo=FALSE---------------------------------------------------------------
plot(c(4,5,6),rep(1,3),
     ylim=c(-0.05,1.05),
     xlim=c(-1,11),
     pch=19,
     axes=FALSE,
     yaxt="n",
     xaxt="n",
     cex=2,
     ylab="",
     xlab="Non-breeding",
     main="Breeding")
points(c(0,5,10),rep(0,3),pch=15,cex=2)

segments(4,1,0,0,lwd=0.4*10)
segments(4,1,5,0,lwd=0.35*10)
segments(4,1,10,0,lwd=0.25*10)

segments(5,1,0,0,lwd=0.3*10, lty = 2)
segments(5,1,5,0,lwd=0.4*10, lty = 2)
segments(5,1,10,0,lwd=0.3*10, lty = 2)

segments(6,1,0,0,lwd=0.2*10, lty = 3)
segments(6,1,5,0,lwd=0.3*10, lty = 3)
segments(6,1,10,0,lwd=0.5*10, lty = 3)

## -----------------------------------------------------------------------------
#equal relative abundance among the three breeding regions, must sum to 1                      
relN <- rep(1/nBreeding, nBreeding) 

relN

# Define total sample size for psi data 
# for small sample corrected version of MC 

n <- 250


## -----------------------------------------------------------------------------
# Calculate the strength of migratory connectivity 
MC <- calcMC(originDist = breedDist,
             targetDist = nonBreedDist,
             originRelAbund = relN,
             psi = psi)

round(MC, 3)

MC <- calcMC(originDist = breedDist,
             targetDist = nonBreedDist,
             originRelAbund = relN, 
             psi = psi,
             sampleSize = n)

round(MC, 3)

## ----eval = FALSE-------------------------------------------------------------
#  # Load in projections
#  data("projections")
#  
#  # Define deployment locations (winter) #
#  captureLocations<-matrix(c(-77.93,18.04,   # Jamaica
#                             -80.94,25.13,   # Florida
#                             -66.86,17.97),  # Puerto Rico
#                              nrow=3, ncol=2, byrow = TRUE)
#  
#  colnames(captureLocations) <- c("Longitude","Latitude")
#  
#  # Convert capture locations into points #
#  
#  CapLocs <- sf::st_as_sf(data.frame(captureLocations),
#                        coords = c("Longitude","Latitude"),
#                        crs = 4326)
#  
#  # Project Capture locations #
#  
#  CapLocsM<-sf::st_transform(CapLocs, 'ESRI:54027')
#  
#  # Retrieve raw non-breeding locations from github
#  # First grab the identity of the bird so we can loop through the files
#  # For this example we are only interested in the error
#  # around non-breeding locations
#  # here we grab only the birds captured during the non-breeding season
#  # Using paste0 for vignette formatting purposes
#  
#  winterBirds <- dget(
#  "https://raw.githubusercontent.com/SMBC-NZP/MigConnectivity/master/data-raw/GL_NonBreedingFiles/winterBirds.txt"
#                      )
#  
#  # create empty list to store the location data #
#  Non_breeding_files <- vector('list',length(winterBirds))
#  
#  # Get raw location data from Github #
#  for(i in 1:length(winterBirds)){
#  Non_breeding_files[[i]] <- dget(paste0("https://raw.githubusercontent.com/",
#                                          "SMBC-NZP/MigConnectivity/master/data-raw/",
#                                          "GL_NonBreedingFiles/NonBreeding_",
#                                          winterBirds[i],".txt"))
#  }
#  
#  # Remove locations around spring Equinox and potential migration points
#  # same NB time frame as Hallworth et al. 2015
#  
#  # two steps because subset on shapefile doesn't like it in a single step
#  
#  Non_breeding_files <- lapply(Non_breeding_files,
#                        FUN = function(x){
#                        month <- as.numeric(format(x$Date,format = "%m"))
#                                 x[which(month != 3 & month != 4),]})
#  
#  
#  Jam <- c(1:9)   # locations w/in list of winterBirds captured in Jamaica
#  Fla <- c(10:12) # locations w/in list of winterBirds in Florida
#  PR <- c(13:16)  # locations w/in list of winterBirds in Puerto Rico
#  
#  # Turn the locations into shapefiles #
#  
#  NB_GL <- lapply(Non_breeding_files,
#                  FUN = function(x){
#                    sf::st_as_sf(data.frame(x),
#                                 coords = c("Longitude",
#                                            "Latitude"),
#                                 crs = 4326)})
#  
#  # Project into UTM projection #
#  
#  NB_GLmeters <- lapply(NB_GL,
#                        FUN = function(x){sf::st_transform(x,'ESRI:54027')})
#  
#  # Process to determine geolocator bias and variance-covariance in meters #
#  
#  # generate empty vector to store data #
#  LongError<-rep(NA,length(winterBirds))
#  LatError<-rep(NA,length(winterBirds))
#  
#  # Calculate the error in longitude derived
#  # from geolocators from the true capture location
#  
#  LongError[Jam] <- unlist(lapply(NB_GLmeters[Jam],
#                           FUN = function(x){mean(sf::st_coordinates(x)[,1]-
#                                                    sf::st_coordinates(CapLocsM)[1,1])}))
#  
#  LongError[Fla] <- unlist(lapply(NB_GLmeters[Fla],
#                           FUN = function(x){mean(sf::st_coordinates(x)[,1]-
#                                                    sf::st_coordinates(CapLocsM)[2,1])}))
#  
#  LongError[PR] <- unlist(lapply(NB_GLmeters[PR],
#                          FUN = function(x){mean(sf::st_coordinates(x)[,1]-
#                                                   sf::st_coordinates(CapLocsM)[3,1])}))
#  
#  # Calculate the error in latitude derived from
#  # geolocators from the true capture location
#  
#  LatError[Jam] <- unlist(lapply(NB_GLmeters[Jam],
#                          FUN = function(x){mean(sf::st_coordinates(x)[,2]-
#                                                    sf::st_coordinates(CapLocsM)[1,2])}))
#  
#  LatError[Fla] <- unlist(lapply(NB_GLmeters[Fla],
#                          FUN = function(x){mean(sf::st_coordinates(x)[,2]-
#                                                   sf::st_coordinates(CapLocsM)[2,2])}))
#  
#  LatError[PR] <- unlist(lapply(NB_GLmeters[PR],
#                          FUN = function(x){mean(sf::st_coordinates(x)[,2]-
#                                                    sf::st_coordinates(CapLocsM)[3,2])}))
#  
#  # Get co-variance matrix for error of
#  # known non-breeding deployment sites
#  
#  # lm does multivariate normal models if you give it a matrix dependent variable!
#  
#  geo.error.model <- lm(cbind(LongError,LatError) ~ 1)
#  
#  geo.bias <- coef(geo.error.model)
#  geo.vcov <- vcov(geo.error.model)

## -----------------------------------------------------------------------------
data(OVENdata) # Ovenbird 

names(OVENdata)

## ----eval = FALSE-------------------------------------------------------------
#  install.packages(c('shape', 'terra', 'maps'))

## ----message = FALSE, warning = FALSE, error=FALSE----------------------------
library(shape)
library(maps)
library(terra)

# Set the crs to WGS84
originSitesWGS84 <- sf::st_transform(OVENdata$originSites, 4326)
targetSitesWGS84 <- sf::st_transform(OVENdata$targetSites, 4326)

# Create a simple plot of the origin and and target sites 
par(mar=c(0,0,0,0))
plot(originSitesWGS84,
     xlim=c(terra::ext(targetSitesWGS84)[1],
            terra::ext(targetSitesWGS84)[2]),

     ylim=c(terra::ext(targetSitesWGS84)[3],
            terra::ext(originSitesWGS84)[4]))

plot(targetSitesWGS84,
     add = TRUE,
     col=c("gray70","gray35","gray10"))

map("world", add=TRUE)

## ----message=FALSE, warning = FALSE, error=FALSE------------------------------
M<-estMC(isGL=OVENdata$isGL, # Logical vector: light-level geolocator(T)/GPS(F)
         geoBias = OVENdata$geo.bias, # Light-level geolocator location bias
         geoVCov = OVENdata$geo.vcov, #Light-level geolocator covariance matrix
         targetDist = OVENdata$targetDist, # Target location distance matrix
         originDist = OVENdata$originDist, # Origin location distance matrix
         targetSites = OVENdata$targetSites, # Non-breeding / target sites
         originSites = OVENdata$originSites, # Breeding / origin sites 
         targetNames = OVENdata$targetNames, # Names of nonbreeding/target sites
         originNames = OVENdata$originNames, # Names of breeding/origin sites
         originPoints = OVENdata$originPoints, # Capture Locations 
         targetPoints = OVENdata$targetPoints, # Target locations from devices
         originRelAbund = OVENdata$originRelAbund, # Origin relative abundances
         resampleProjection = sf::st_crs(OVENdata$targetPoints), 
         verbose = 0,   # output options - see help ??estMC
         nSamples = 10) # This is set low for example 

M

## -----------------------------------------------------------------------------
str(M, max.level = 2)

## ----echo = FALSE, warning = FALSE, error = FALSE, eval = FALSE---------------
#  plot(sf::st_geometry(wrld_simple),
#       xlim=sf::st_bbox(OVENdata$targetSites)[c(1,3)],
#       ylim=c(sf::st_bbox(OVENdata$targetSites)[2],
#              sf::st_bbox(OVENdata$originSites)[4]))
#  plot(sf::st_geometry(OVENdata$originSites[1,]),add=TRUE,lwd=1.75)
#  plot(sf::st_geometry(OVENdata$originSites[2,]),add=TRUE,lwd=1.75)
#  plot(sf::st_geometry(OVENdata$targetSites),add=TRUE,lwd=1.5,col=c("gray70","gray35","gray10"))
#  
#  legend(400000,-2105956,legend=paste("MC =",round(M$MC$mean,2), "\u00b1", round(M$MC$se,2)),bty="n",cex=1.8,bg="white",xjust=0)
#  
#  shape::Arrows(x0 = sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,1],
#                y0 = sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,2],
#                x1 = sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[2,]))[,1]+80000,
#                y1 = sf::st_bbox(OVENdata$targetSites[2,])[4]+150000,
#                arr.length = 0.3,
#                arr.adj = 0.5,
#                arr.lwd = 1,
#                arr.width = 0.4,
#                arr.type = "triangle",
#                lwd = M$psi$mean[1,2]*10,
#                lty = 1)
#  
#  shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,1],
#               sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,2],
#                sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[3,]))[,1],
#                sf::st_bbox(OVENdata$targetSites[3,])[4]+150000,
#                arr.length = 0.3,
#                arr.adj = 0.5,
#                arr.lwd = 1,
#                arr.width = 0.4,
#                arr.type = "triangle",
#                lwd = M$psi$mean[1,3]*10,
#                lty=1)
#  
#  shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,1],
#                sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,2],
#                sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[1,]))[,1],
#                sf::st_bbox(OVENdata$targetSites[1,])[4]+150000,
#                arr.length = 0.3,
#                arr.adj = 0.5,
#                arr.lwd = 1,
#                arr.width = 0.4,
#                arr.type = "triangle",
#                lwd = M$psi$mean[2,1]*10,
#                lty=1)
#  
#  shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,1],
#                sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,2],
#                sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[2,]))[,1]-80000,
#                sf::st_bbox(OVENdata$targetSites[2,])[4]+150000,
#                arr.length = 0.3,
#                arr.adj = 0.5,
#                arr.lwd = 1,
#                arr.width = 0.4,
#                arr.type = "triangle",
#                lwd = M$psi$mean[2,2]*10)
#  
#  box(which="plot")

## -----------------------------------------------------------------------------
set.seed(75)

# Parameters for simulations

nSeasons <- 2 # population with two discrete seasons - breeding / non-breeding
nYears <- 10 # Ten years of data
nMonths <- 4 # Four months within each season

nBreeding <- 100 # Number of populations
nWintering <- 100 # Number of populations 

## -----------------------------------------------------------------------------
breedingPos <- matrix(c(rep(seq(-99,-81,2), each=sqrt(nBreeding)),
                        rep(seq(49,31,-2), sqrt(nBreeding))), nBreeding, 2)
winteringPos <- matrix(c(rep(seq(-79,-61,2), each=sqrt(nWintering)),
                         rep(seq(9,-9,-2), sqrt(nWintering))), nWintering, 2)

# calculate distance between populations
breedDist <- distFromPos(breedingPos, 'ellipsoid') 

nonbreedDist <- distFromPos(winteringPos, 'ellipsoid') 

## -----------------------------------------------------------------------------
# Breeding Abundance
breedingN <- rep(500, nBreeding)
breedingRelN <- breedingN/sum(breedingN)

## ----eval = FALSE-------------------------------------------------------------
#  # Set up psi matrix
#  o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist,
#                target.dist = nonbreedDist, origin.abund = breedingRelN,
#                sample.size = sum(breedingN),
#                interval = c(1, 3), tol = .Machine$double.eps^0.5)#
#  
#  slope <- o$minimum
#  
#  psi <- mlogitMat(slope, breedDist)

## ----warning=FALSE,message=FALSE,echo=FALSE-----------------------------------
slope <- 1.92093332054914
psi <- mlogitMat(slope, breedDist)

## -----------------------------------------------------------------------------
# Baseline strength of migratory connectivity
MC <- calcMC(originDist = breedDist, 
             targetDist = nonbreedDist,
             psi = psi, 
             originRelAbund = breedingRelN,
             sampleSize = sum(breedingN))

# Show Results
MC

## -----------------------------------------------------------------------------
sims <- simMove(breedingAbund = breedingN, # Breeding relative abundance
               breedingDist = breedDist,  # Breeding distance
               winteringDist = nonbreedDist, # Non-breeding distance
               psi = psi, # Transition probabilities
               nYears = nYears, # Number of years
               nMonths = nMonths, # Number of months
               winMoveRate = 0, # winter movement rate
               sumMoveRate = 0, # breeding movement rate
               winDispRate = 0, # winter dispersal rate
               sumDispRate = 0, # summer disperal rate
               natalDispRate = 0, # natal dispersal rate
               breedDispRate = 0, # breeding dispersal rate
               verbose = 0) # verbose output


str(sims)

## -----------------------------------------------------------------------------
nScenarios1 <- length(samplePsis) # samplePsis - comes with MigConnectivity package

# Create vector of length nScenarios1
MC1 <- rep(NA, nScenarios1)

# Loop through the different senarios outlined above #
for (i in 1:nScenarios1) {
  MC1[i] <- calcMC(originDist = sampleOriginDist[[1]],
                 targetDist = sampleTargetDist[[1]],
                 psi = samplePsis[[i]],
                 originRelAbund = sampleOriginRelN[[1]])
}

# Give meaningful names to the MC1 vector
names(MC1) <- names(samplePsis)

# Print results 
round(MC1, 6)

## -----------------------------------------------------------------------------
nScenarios2 <- length(sampleOriginPos)

MC2 <- matrix(NA, nScenarios1, nScenarios2)

rownames(MC2) <- names(samplePsis)

colnames(MC2) <- names(sampleOriginPos)

for (i in 1:nScenarios1) {
  for (j in 1:nScenarios2) {
    MC2[i, j] <- calcMC(originDist = sampleOriginDist[[j]],
                      targetDist = sampleTargetDist[[j]],
                      psi = samplePsis[[i]],
                      originRelAbund = sampleOriginRelN[[1]])
  }
}

# Print results #
t(round(MC2, 4))

## -----------------------------------------------------------------------------

MC.diff2 <- apply(MC2, 2, "-", MC2[ , 1])

t(round(MC.diff2, 4))

## -----------------------------------------------------------------------------
nScenarios3 <- length(sampleOriginRelN)

MC3 <- matrix(NA, nScenarios1, nScenarios3)

rownames(MC3) <- names(samplePsis)
colnames(MC3) <- names(sampleOriginRelN)

for (i in 1:nScenarios1) {
for (j in 1) {
for (k in 1:nScenarios3) {
MC3[i, k] <- calcMC(originDist = sampleOriginDist[[j]],
                    targetDist = sampleTargetDist[[j]], 
                    psi = samplePsis[[i]],
                    originRelAbund = sampleOriginRelN[[k]])
}
}
}

t(round(MC3, 4)) # linear arrangement

## -----------------------------------------------------------------------------
MC4 <- matrix(NA, nScenarios1, nScenarios3)

rownames(MC4) <- names(samplePsis)
colnames(MC4) <- names(sampleOriginRelN)

for (i in 1:nScenarios1) {
for (j in nScenarios2) {
for (k in 1:nScenarios3) {
MC4[i, k] <- calcMC(originDist = sampleOriginDist[[j]],
                    targetDist = sampleTargetDist[[j]],
                    psi = samplePsis[[i]],
                    originRelAbund = sampleOriginRelN[[k]])
}
}
}

t(round(MC4, 4)) # grid arrangement

## -----------------------------------------------------------------------------
MC5 <- matrix(NA, nScenarios1, nScenarios3)

rownames(MC5) <- names(samplePsis)
colnames(MC5) <- names(sampleOriginRelN)

for (i in 1:nScenarios1) {
for (j in 6) {
for (k in 1:nScenarios3) {
MC5[i, k] <- calcMC(originDist = sampleOriginDist[[j]],
                    targetDist = sampleTargetDist[[j]], 
                    psi = samplePsis[[i]], 
                    originRelAbund = sampleOriginRelN[[k]])
}
}
}

t(round(MC5, 4)) # breeding grid, winter linear arrangement

## ----eval=FALSE---------------------------------------------------------------
#  #Run functions and parameters above first
#  set.seed(75)
#  
#  scenarios14 <- c("Base",
#                   "Breeding10",
#                   "Wintering10",
#                   "Breeding10Wintering10",
#                   "Breeding4",
#                   "Wintering4",
#                   "Breeding4Wintering10",
#                   "Breeding10Wintering4",
#                   "Breeding4Wintering4",
#                   "CentroidSampleBreeding10",
#                   "CentroidSampleBreeding10Wintering10",
#                   "CentroidSampleBreeding10Wintering4",
#                   "CentroidSampleBreeding4",
#                   "CentroidSampleBreeding4Wintering10",
#                   "CentroidSampleBreeding4Wintering4")
#  
#  # Each element is for a scenario (see above 1-8), transferring from natural
#  # breeding populations to defined ones
#  breedingSiteTrans14 <- list(1:nBreeding,
#                             rep(1:10, each=10),
#                             1:nBreeding,
#                             rep(1:10, each=10),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             1:nBreeding,
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             rep(1:10, each=10),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             rep(1:10, each=10),
#                             rep(1:10, each=10),
#                             rep(1:10, each=10),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
#  
#  # Same for non-breeding populations
#  winteringSiteTrans14 <- list(1:nWintering,
#                               1:nWintering,
#                               rep(1:10, each=10),
#                               rep(1:10, each=10),
#                               1:nWintering,
#                               c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                               rep(1:10, each=10),
#                               c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                               c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                               1:nWintering,
#                               rep(1:10, each=10),
#                               c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                               1:nWintering,
#                               rep(1:10, each=10),
#                               c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
#  
#  # Examine the transfers in matrix form
#  #lapply(breedingSiteTrans14, matrix, nrow=10, ncol=10)
#  #lapply(winteringSiteTrans14, matrix, nrow=10, ncol=10)
#  
#  #positions of the human defined populations
#  breedingPos14 <- list(breedingPos,
#                        rowsum(breedingPos, rep(1:10, each=10))/10,
#                        breedingPos,
#                        rowsum(breedingPos, rep(1:10, each=10))/10,
#                        rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                              rep(3:4, 5, each=5)))/25,
#                        breedingPos,
#                        rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                              rep(3:4, 5, each=5)))/25,
#                        rowsum(breedingPos, rep(1:10, each=10))/10,
#                        rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                              rep(3:4, 5, each=5)))/25,
#                        rowsum(breedingPos, rep(1:10, each=10))/10,
#                        rowsum(breedingPos, rep(1:10, each=10))/10,
#                        rowsum(breedingPos, rep(1:10, each=10))/10,
#                        rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                              rep(3:4, 5, each=5)))/25,
#                        rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                              rep(3:4, 5, each=5)))/25,
#                        rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                              rep(3:4, 5, each=5)))/25)
#  
#  winteringPos14 <- list(winteringPos,
#                         winteringPos,
#                         rowsum(winteringPos, rep(1:10, each=10))/10,
#                         rowsum(winteringPos, rep(1:10, each=10))/10,
#                         winteringPos,
#                         rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                                rep(3:4, 5, each=5)))/25,
#                         rowsum(winteringPos,rep(1:10, each=10))/10,
#                         rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                                rep(3:4, 5, each=5)))/25,
#                         rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                                rep(3:4, 5, each=5)))/25,
#                         winteringPos,
#                         rowsum(winteringPos, rep(1:10, each=10))/10,
#                         rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                                rep(3:4, 5, each=5)))/25,
#                         winteringPos,
#                         rowsum(winteringPos, rep(1:10, each=10))/10,
#                         rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                                rep(3:4, 5, each=5)))/25)
#  
#  
#  # Calculate distances between defined breeding populations
#  breedDist14 <- lapply(breedingPos14, distFromPos, surface = 'ellipsoid')
#  
#  # Calculate distances between defined non-breeding populations
#  nonbreedDist14 <- lapply(winteringPos14, distFromPos, surface = 'ellipsoid')
#  
#  # Numbers of defined populations
#  nBreeding14 <- c(100, 10, 100, 10, 4, 100, 4, 10, 4, 10, 10, 10, 4, 4, 4)
#  
#  nWintering14 <- c(100, 100, 10, 10, 100, 4, 10, 4, 4, 100, 10, 4, 100, 10, 4)
#  
#  # Relative abundance by scenario and breeding population
#  breedingRelN14 <- lapply(nBreeding14, function(x) rep(1/x, x))
#  
#  MC14 <- 0.25 # Same for all scenarios in this set
#  
#  nSample14 <- 1000 # Total number sampled per simulation
#  
#  # How many sampled from each natural population (sampling scenarios separate
#  # from definition scenarios)
#  
#  sampleBreeding14 <- list(round(breedingRelN14[[1]]*nSample14),
#                           c(rep(0, 22), round(breedingRelN14[[5]][1]*nSample14),
#                             rep(0, 4), round(breedingRelN14[[5]][2]*nSample14),
#                             rep(0, 44), round(breedingRelN14[[5]][3]*nSample14),
#                             rep(0, 4), round(breedingRelN14[[5]][4]*nSample14),
#                             rep(0, 22)),
#                           rep(c(rep(0, 4),
#                               rep(round(breedingRelN14[[2]][1]*nSample14/2), 2),
#                               rep(0, 4)), 10))
#  
#  # for the baseline use the simulation from above, sims
#  
#  animalLoc14base <- sims$animalLoc
#  
#  
#  # Number of scenarios and number of simulations to run
#  nScenarios14 <- length(breedingSiteTrans14)
#  nSims14 <- 100
#  
#  # Connections between scenarios and sampling regimes
#  scenarioToSampleMap14 <- c(rep(1, 9), rep(3, 3), rep(2, 3))
#  
#  # Set up data structures for storing results
#  animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill
#  
#  compare14 <- data.frame(Scenario = c("True", scenarios14),
#                          MC = c(MC14, rep(NA, nScenarios14)),
#                          Mantel = c(MC14, rep(NA, nScenarios14)))
#  
#  compare14.array <- array(NA, c(nSims14, nScenarios14, 2), dimnames =
#                             list(1:nSims14,
#                             scenarios14,
#                             c("MC", "Mantel")))
#  
#  results14 <- vector("list", nScenarios14)
#  
#  # Run simulations
#  for (sim in 1:nSims14) {
#    cat("Simulation", sim, "of", nSims14, '\n')
#    sim14 <- lapply(sampleBreeding14,
#                    simMove,
#                    breedingDist = breedDist14[[1]],
#                    winteringDist=nonbreedDist14[[1]],
#                    psi=psi,
#                    nYears=nYears,
#                    nMonths=nMonths)
#  
#    for (i in 1:nScenarios14) {
#     cat("\tScenario", i, "\n")
#     animalLoc14[[i]]<-changeLocations(
#       sim14[[scenarioToSampleMap14[i]]]$animalLoc,
#                                       breedingSiteTrans14[[i]],
#                                       winteringSiteTrans14[[i]])
#  
#     results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]],
#                                 animalLoc14[[i]],
#                                 originRelAbund = breedingRelN14[[i]],
#                                 verbose = FALSE)
#  
#     compare14.array[sim, i, 'MC'] <- results14[[i]]$MC
#  
#     compare14.array[sim,i,'Mantel']<-calcStrengthInd(breedDist14[[1]],
#                                                      nonbreedDist14[[1]],
#                                     sim14[[scenarioToSampleMap14[i]]]$animalLoc,
#                                                      resamp=0)$correlation
#    }
#  }
#  
#  # Compute means for each scenario
#  compare14$MC[1:nScenarios14 + 1] <- apply(compare14.array[,,'MC'], 2, mean)
#  compare14$Mantel[1:nScenarios14 + 1] <- apply(compare14.array[,,'Mantel'], 2,
#                                                mean)
#  
#  compare14 <- transform(compare14, MC.diff=MC - MC[1],
#                         Mantel.diff=Mantel - Mantel[1])
#  compare14
#  

## ----eval = FALSE-------------------------------------------------------------
#  set.seed(75)
#  
#  # Transfer between true populations and researcher defined ones (only for
#  # breeding, as not messing with winter populations here)
#  
#  breedingSiteTrans15 <- list(1:100,
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              1:100,
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
#  
#  nScenarios15 <- length(breedingSiteTrans15)
#  
#  nSims15 <- 100
#  
#  # Basing positions of researcher defined breeding populations on above
#  breedingPos15 <- list(breedingPos,
#                        breedingPos14[[5]],
#                        breedingPos14[[5]],
#                        breedingPos,
#                        breedingPos14[[5]],
#                        breedingPos14[[5]])
#  
#  winteringPos15 <- rep(list(winteringPos), nScenarios15)
#  
#  breedDist15 <- lapply(breedingPos15, distFromPos)
#  
#  nonbreedDist15 <- lapply(winteringPos15, distFromPos)
#  
#  nBreeding15 <- rep(c(100, 4, 4), 2)
#  
#  nWintering15 <- rep(100, nScenarios15)
#  
#  # Highest abundance in lower right corner, lowest in top left
#  # Making symmetrical
#  
#  breedingN15base <- rep(NA, 100)
#  for (i in 1:10) #row
#    for (j in 1:10)  #column
#      breedingN15base[i+10*(j-1)] <- 500 + 850*i*j
#  
#  sum(breedingN15base)
#  
#  # For researcher defined populations
#  breedingN15 <- lapply(breedingSiteTrans15, rowsum, x=breedingN15base)
#  
#  breedingRelN15 <- lapply(breedingN15, "/", sum(breedingN15base))
#  
#  nSample15 <- 1000 # Total number sampled per simulation
#  
#  # Number sampled per natural population
#  sampleBreeding15 <- list(round(breedingRelN15[[1]]*nSample15),
#                           c(rep(0, 22), round(breedingRelN15[[3]][1]*nSample15),
#                             rep(0, 4), round(breedingRelN15[[3]][2]*nSample15),
#                             rep(0, 44), round(breedingRelN15[[3]][3]*nSample15),
#                             rep(0, 4), round(breedingRelN15[[3]][4]*nSample15),
#                             rep(0, 22)),
#                           round(breedingRelN15[[1]]*nSample15)[100:1],
#                           c(rep(0, 22), round(breedingRelN15[[3]][1]*nSample15),
#                             rep(0, 4), round(breedingRelN15[[3]][2]*nSample15),
#                             rep(0, 44), round(breedingRelN15[[3]][3]*nSample15),
#                             rep(0, 4), round(breedingRelN15[[3]][4]*nSample15),
#                             rep(0, 22))[100:1])
#  
#  # Set up psi matrix
#  o15 <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist15[[1]],
#                  target.dist = nonbreedDist15[[1]],
#                  origin.abund = breedingN15[[1]]/sum(breedingN15[[1]]),
#                  sample.size = sum(breedingN15[[1]]),
#                  interval = c(0,10),
#                  tol = .Machine$double.eps^0.5)
#  
#  slope15 <- o15$minimum
#  
#  psi15 <- mlogitMat(slope15, breedDist15[[1]])
#  
#  # Baseline strength of migratory connectivity
#  MC15 <- calcMC(originDist = breedDist15[[1]],
#                 targetDist = nonbreedDist15[[1]],
#                 psi = psi15,
#                 originRelAbund = breedingN15[[1]]/sum(breedingN15[[1]]),
#                 sampleSize = sum(breedingN15[[1]]))
#  
#  # Run sampling regimes
#  scenarioToSampleMap15 <- c(1, 1, 2, 3, 3, 4)
#  
#  animalLoc15 <- vector("list", nScenarios15)
#  
#  results15 <- vector("list", nScenarios15)
#  
#  compare15 <- data.frame(Scenario = c("True",
#                                       "Base",
#                                       "Breeding4",
#                                       "CentroidSampleBreeding4",
#                                       "BiasedSample",
#                                       "BiasedSampleBreeding4",
#                                       "BiasedCentroidSampleBreeding4"),
#                          MC = c(MC15, rep(NA, nScenarios15)),
#                          Mantel = c(MC15, rep(NA, nScenarios15)))
#  
#  compare15.array <- array(NA, c(nSims15, nScenarios15, 2),
#                           dimnames = list(1:nSims15,
#                                      c("Base", "Breeding4",
#                                        "CentroidSampleBreeding4",
#                                        "BiasedSample", "BiasedSampleBreeding4",
#                                        "BiasedCentroidSampleBreeding4"),
#                                      c("MC", "Mantel")))
#  
#  
#  for (sim in 1:nSims15) {
#    cat("Simulation", sim, "of", nSims15, '\n')
#  
#    sim15 <- lapply(sampleBreeding15,
#                    simMove,
#                    breedingDist = breedDist15[[1]],
#                    winteringDist=nonbreedDist15[[1]],
#                    psi=psi15,
#                    nYears=nYears,
#                    nMonths=nMonths)
#  
#  for (i in 1:nScenarios15) {
#    cat("\tScenario", i, "\n")
#      animalLoc15[[i]] <- changeLocations(
#                          animalLoc=sim15[[scenarioToSampleMap15[i]]]$animalLoc,
#                          breedingSiteTrans = breedingSiteTrans15[[i]],
#                          winteringSiteTrans = 1:nWintering15[i])
#  
#      results15[[i]] <- calcPsiMC(originDist = breedDist15[[i]],
#                                  targetDist = nonbreedDist15[[i]],
#                                  originRelAbund = breedingRelN15[[i]],
#                                  locations = animalLoc15[[i]],
#                                  verbose = FALSE)
#  
#  compare15.array[sim, i, 'MC'] <- results15[[i]]$MC
#  
#  compare15.array[sim, i, 'Mantel'] <- calcStrengthInd(breedDist15[[1]],
#                                                       nonbreedDist15[[1]],
#                                     sim15[[scenarioToSampleMap15[i]]]$animalLoc,
#                                                      resamp=0)$correlation
#    }
#  }
#  
#  compare15$MC[1:nScenarios15+1] <- apply(compare15.array[,,'MC'], 2, mean,
#                                          na.rm = TRUE)
#  
#  compare15$Mantel[1:nScenarios15+1] <- apply(compare15.array[,,'Mantel'], 2,
#                                              mean, na.rm = TRUE)
#  compare15 <- transform(compare15, MC.diff=MC - MC[1],
#                         Mantel.diff=Mantel - Mantel[1],
#                         MC.prop=MC/MC[1],
#                         Mantel.prop=Mantel/Mantel[1])
#  
#  
#  compare15a <- as.matrix(compare15[1 + 1:nScenarios15,
#                                    c("MC", "MC.diff", "Mantel", "Mantel.diff")])
#  
#  rownames(compare15a) <- compare15$Scenario[1 + 1:nScenarios15]

## ----eval = FALSE-------------------------------------------------------------
#  set.seed(75)
#  
#  # Transfer between true populations and researcher defined ones
#  # (only for breeding, as not messing with winter populations here)
#  breedingSiteTrans16 <- list(1:100, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:100,
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
#  
#  #lapply(breedingSiteTrans16, matrix, nrow=10, ncol=10)
#  
#  nScenarios16 <- length(breedingSiteTrans16)
#  
#  nSims16 <- 100
#  
#  # Basing positions of researcher defined breeding populations on above
#  breedingPos16 <- breedingPos15
#  winteringPos16 <- winteringPos15
#  
#  breedDist16 <- breedDist15
#  nonbreedDist16 <- nonbreedDist15
#  nBreeding16 <- nBreeding15
#  nWintering16 <- rep(100, nScenarios16)
#  
#  # Highest abundance in lower right corner, lowest in top left
#  # In fact basing on distance from top left population
#  
#  breedingN16base <- breedingN15base
#  breedingN16 <- breedingN15
#  breedingRelN16 <- lapply(breedingN16, "/", sum(breedingN16base))
#  
#  # Set up psi matrix
#  # Each quadrant of breeding range has different MC
#  
#  MC.levels16 <- seq(0.15, 0.6, 0.15)
#  
#  nLevels16 <- 4
#  
#  psi16 <- matrix(NA, nBreeding16[1], nWintering16[1])
#  
#  
#  for (i in 1:nLevels16) {
#    cat("MC", MC.levels16[i])
#    # Find a psi matrix that produces the given MC (for whole species)
#    o16a <- optimize(mlogitMC, MC.in = MC.levels16[i],
#                     origin.dist = breedDist16[[1]],
#                     target.dist = nonbreedDist16[[1]],
#                     origin.abund = breedingN16[[1]]/sum(breedingN16[[1]]),
#                     sample.size = sum(breedingN16[[1]]),
#                     interval=c(0,10), tol=.Machine$double.eps^0.5)
#  
#    slope16a <- o16a$minimum
#  
#    cat(" slope", slope16a, "\n")
#  
#    psi16a <- mlogitMat(slope16a, breedDist16[[1]])
#  
#    # Then use the rows of that psi matrix only for the one breeding quadrant
#    rows <- 50*(i %/% 3) + rep(1:5, 5) + rep(seq(0, 40, 10), each=5) +
#      ((i-1) %% 2) * 5
#    psi16[rows, ] <- psi16a[rows, ]
#  }
#  
#  
#  # Baseline strength of migratory connectivity
#  MC16 <- calcMC(originDist = breedDist16[[1]],
#                 targetDist = nonbreedDist16[[1]],
#                 psi = psi16,
#                 originRelAbund = breedingN16[[1]]/sum(breedingN16[[1]]),
#                 sampleSize = sum(breedingN16[[1]]))
#  
#  
#  # Set up sampling regimes (different number than number of scenarios)
#  nSample16 <- 1000
#  sampleBreeding16 <- list(round(breedingRelN16[[1]]*nSample16),
#                           c(rep(0, 22), round(breedingRelN16[[3]][1]*nSample16),
#                             rep(0, 4), round(breedingRelN16[[3]][2]*nSample16),
#                             rep(0, 44), round(breedingRelN16[[3]][3]*nSample16),
#                             rep(0, 4), round(breedingRelN16[[3]][4]*nSample16),
#                             rep(0, 22)),
#                           round(breedingRelN16[[1]]*nSample16)[100:1],
#                           c(rep(0, 22), round(breedingRelN16[[3]][1]*nSample16),
#                             rep(0, 4), round(breedingRelN16[[3]][2]*nSample16),
#                             rep(0, 44), round(breedingRelN16[[3]][3]*nSample16),
#                             rep(0, 4), round(breedingRelN16[[3]][4]*nSample16),
#                             rep(0, 22))[100:1])
#  
#  sapply(sampleBreeding16, sum)
#  
#  # Run sampling regimes
#  scenarioToSampleMap16 <- c(1, 1, 2, 3, 3, 4)
#  animalLoc16 <- vector("list", nScenarios16)
#  results16 <- vector("list", nScenarios16)
#  compare16 <- data.frame(Scenario = c("True",
#                                       "Base",
#                                       "Breeding4",
#                                       "CentroidSampleBreeding4",
#                                       "BiasedSample",
#                                       "BiasedSampleBreeding4",
#                                       "BiasedCentroidSampleBreeding4"),
#                          MC = c(MC16, rep(NA, nScenarios16)),
#                          Mantel = c(MC16, rep(NA, nScenarios16)))
#  
#  compare16.array <- array(NA, c(nSims16, nScenarios16, 2),
#                           dimnames = list(1:nSims16,
#                                           c("Base", "Breeding4",
#                                             "CentroidSampleBreeding4",
#                                             "BiasedSample",
#                                             "BiasedSampleBreeding4",
#                                             "BiasedCentroidSampleBreeding4"),
#                                           c("MC", "Mantel")))
#  
#  set.seed(80)
#  for (sim in 1:nSims16) {
#    cat("Simulation", sim, "of", nSims16, '\n')
#  
#    sim16 <- lapply(sampleBreeding16, simMove, breedingDist = breedDist16[[1]],
#                    winteringDist=nonbreedDist16[[1]], psi=psi16, nYears=nYears,
#                    nMonths=nMonths)
#  
#  for (i in 1:nScenarios16) {
#    cat("\tScenario", i, "\n")
#    animalLoc16[[i]]<-changeLocations(sim16[[scenarioToSampleMap16[i]]]$animalLoc,
#                                        breedingSiteTrans16[[i]],
#                                        1:nWintering16[i])
#  
#    results16[[i]] <- calcPsiMC(originDist = breedDist16[[i]],
#                                targetDist = nonbreedDist16[[i]],
#                                locations = animalLoc16[[i]],
#                                originRelAbund = breedingRelN16[[i]],
#                                verbose = FALSE)
#  
#  compare16.array[sim, i, 'MC'] <- results16[[i]]$MC
#  
#  compare16.array[sim, i, 'Mantel'] <- calcStrengthInd(breedDist16[[1]],
#                                                           nonbreedDist16[[1]],
#                                     sim16[[scenarioToSampleMap16[i]]]$animalLoc,
#                                                           resamp=0)$correlation
#    }
#  }
#  
#  compare16$MC[1:nScenarios16+1] <- apply(compare16.array[,,'MC'], 2, mean, na.rm = TRUE)
#  
#  compare16$Mantel[1:nScenarios16+1] <- apply(compare16.array[,,'Mantel'], 2,
#                                              mean, na.rm = TRUE)
#  
#  compare16 <- transform(compare16, MC.diff=MC - MC[1],
#                         Mantel.diff=Mantel - Mantel[1])
#  
#  compare16a <- as.matrix(compare16[1 + 1:nScenarios16, c('MC','MC.diff',
#                                                          "Mantel",
#                                                          "Mantel.diff")])
#  
#  rownames(compare16a) <- compare16$Scenario[1 + 1:nScenarios16]

## ----eval = FALSE-------------------------------------------------------------
#  # Load in projections
#  data("projections")
#  
#  set.seed(75)
#  capLocs14 <- lapply(breedingPos14,
#                      FUN = function(x){
#                        colnames(x)<-c("Longitude","Latitude")
#                        sf::st_as_sf(data.frame(x),
#                                     coords = c("Longitude","Latitude"),
#                                     crs = 4326)})
#  
#  targLocs14 <- lapply(winteringPos14,
#                       FUN = function(x){
#                        colnames(x)<-c("Longitude","Latitude")
#                        sf::st_as_sf(data.frame(x),
#                                     coords = c("Longitude","Latitude"),
#                                     crs = 4326)})
#  
#  
#  # Relative abundance by scenario and breeding population
#  breedingN14base <- rep(25, nBreeding14[1])
#  breedingN14 <- lapply(breedingSiteTrans14, rowsum, x=breedingN14base)
#  breedingRelN14 <- lapply(breedingN14, "/", sum(breedingN14base))
#  
#  
#  MC14 <- 0.25
#  
#  nSample14.1 <- 100 # Total number sampled per simulation
#  
#  # How many sampled from each natural population (sampling scenarios separate
#  # from definition scenarios)
#  sampleBreeding14.1 <- list(round(breedingRelN14[[1]]*nSample14.1),
#                           c(rep(0,22),round(breedingRelN14[[5]][1]*nSample14.1),
#                             rep(0,4), round(breedingRelN14[[5]][2]*nSample14.1),
#                             rep(0,44),round(breedingRelN14[[5]][3]*nSample14.1),
#                             rep(0,4), round(breedingRelN14[[5]][4]*nSample14.1),
#                             rep(0, 22)),
#                             rep(c(rep(0, 4),
#                             rep(round(breedingRelN14[[2]][1]*nSample14.1/2), 2),
#                             rep(0, 4)), 10))
#  
#  lapply(sampleBreeding14.1, matrix, nrow=10, ncol=10)
#  
#  # Number of simulations to run
#  nSims14 <- 100
#  nSimsLarge14 <- 250 #0riginally 2500
#  nYears <- 1
#  nMonths <- 1
#  
#  # Set up data structures for storing results
#  animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill
#  sim14 <- vector("list", nSimsLarge14) #making an empty list to fill
#  
#  compare14.1.array <- array(NA, c(nSimsLarge14, 2),
#                             dimnames = list(1:nSimsLarge14,
#                                             c("MC", "Mantel")))
#  
#  results14 <- vector("list", nScenarios14)
#  
#  # Run simulations
#  set.seed(7)
#  
#  system.time(for (sim in 1:nSimsLarge14) {
#  
#    cat("Simulation", sim, "of", nSimsLarge14, '\n')
#  
#    sim14[[sim]] <- lapply(sampleBreeding14.1,
#                           simMove,
#                           breedingDist = breedDist14[[1]],
#                           winteringDist=nonbreedDist14[[1]],
#                           psi=psi,
#                           nYears=nYears,
#                           nMonths=nMonths)
#  
#  for (i in 13) {
#   animalLoc14[[i]] <- changeLocations(sim14[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                       breedingSiteTrans14[[i]],
#                                       winteringSiteTrans14[[i]])
#  
#   results14[[i]] <- calcPsiMC(breedDist14[[i]],
#                               nonbreedDist14[[i]],
#                               animalLoc14[[i]],
#                               originRelAbund = breedingRelN14[[i]],
#                               verbose = FALSE)
#  
#    compare14.1.array[sim,'MC'] <- results14[[i]]$MC
#  
#    compare14.1.array[sim,'Mantel'] <- calcStrengthInd(breedDist14[[1]],
#                                                    nonbreedDist14[[1]],
#                              sim14[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                                    resamp=0)$correlation
#    }
#  })
#  
#  means14 <- apply(compare14.1.array, 2, mean)
#  
#  vars14 <- apply(compare14.1.array, 2, var)
#  
#  rmse14 <- apply(compare14.1.array, 2, function(x) sqrt(mean((x - MC14)^2)))
#  
#  # Set up data structures for storing estimation results
#  est14.array <- array(NA, c(nSims14, 2),
#                       dimnames = list(1:nSims14,c("MC", "Mantel")))
#  
#  var14.array <- array(NA, c(nSims14, 2),
#                       dimnames =list(1:nSims14,c("MC", "Mantel")))
#  
#  ci14.array <- array(NA, c(nSims14, 2, 2),
#                      dimnames = list(1:nSims14,
#                                      c("MC", "Mantel"),
#                                      c('lower', 'upper')))
#  
#  #making an empty list to fill
#  
#  animalLoc14 <- vector("list", nSims14)
#  
#  results14 <- vector("list", nSims14)
#  
#  # Run estimations
#  set.seed(567)
#  
#  sim14.sub <- sim14[sample.int(nSimsLarge14, nSims14, TRUE)]
#  
#  for (sim in 1:nSims14) {
#    cat("Estimation", sim, "of", nSims14, '\n')
#  
#   for (i in 13) {#:nScenarios14) {
#    animalLoc14[[sim]] <-
#      changeLocations(sim14.sub[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                          breedingSiteTrans14[[i]],
#                                          winteringSiteTrans14[[i]])
#  
#    results14[[sim]] <- estMC(originDist = breedDist14[[i]],
#                              targetDist = nonbreedDist14[[i]],
#                              originRelAbund = breedingRelN14[[i]],
#                originPoints = capLocs14[[i]][animalLoc14[[sim]][, 1, 1, 1], ],
#                targetPoints = targLocs14[[i]][animalLoc14[[sim]][, 2, 1, 1], ],
#                              originAssignment = animalLoc14[[sim]][, 1, 1, 1],
#                              targetAssignment = animalLoc14[[sim]][, 2, 1, 1],
#                              nSamples = 1000,
#                              verbose = 1,
#                              calcCorr = TRUE,
#                              geoBias = c(0, 0),
#                              geoVCov = matrix(0, 2, 2))
#  
#      est14.array[sim, 'MC'] <- results14[[sim]]$MC$mean
#      est14.array[sim, 'Mantel'] <- results14[[sim]]$corr$mean
#      var14.array[sim, 'MC'] <- results14[[sim]]$MC$se ^ 2
#      var14.array[sim, 'Mantel'] <- results14[[sim]]$corr$se ^ 2
#      ci14.array[sim, 'MC', ] <- results14[[sim]]$MC$bcCI
#      ci14.array[sim, 'Mantel', ] <- results14[[sim]]$corr$bcCI
#    }
#  }
#  
#  # Crude point estimates (bootstrap means are better)
#  pointEsts14.1 <- t(sapply(results14, function(x) c(x$MC$point, x$corr$point)))
#  
#  # Summarize results
#  summary(var14.array)
#  vars14
#  summary(est14.array)
#  summary(pointEsts14.1)
#  means14
#  colMeans(pointEsts14.1) - MC14
#  sqrt(colMeans((pointEsts14.1 - MC14)^2))
#  summary(ci14.array[, "MC", 'lower'] <= MC14 &
#            ci14.array[, "MC", 'upper'] >= MC14)
#  summary(ci14.array[, "Mantel", 'lower'] <= MC14 &
#            ci14.array[, "Mantel", 'upper'] >= MC14)
#  
#  # Plot histograms of bootstrap estimation results
#  est.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2, each = nSims14),
#                       Quantity = rep(c('Mean', 'Variance'), each = 2 * nSims14),
#                       sim = rep(1:nSims14, 4),
#                       Estimate = c(est14.array, var14.array))
#  
#  trues.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2),
#                         Quantity = rep(c('Mean', 'Variance'), each = 2),
#                         Value = c(MC14, MC14, vars14))
#  
#  library(ggplot2)
#  g.est <- ggplot(est.df, aes(Estimate)) + geom_histogram(bins = 15) +
#    facet_grid(Parameter ~ Quantity, scales = 'free_x') +
#    geom_vline(aes(xintercept = Value), data = trues.df) +
#    theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3))
#  
#  g.est
#  
#  
#  # Summary table of bootstrap estimation results
#  qualities14 <- data.frame(Parameter = rep(c("MC", 'rM'), 5),
#                            Quantity = rep(rep(c('Mean', 'Variance'), c(2, 2)),
#                                           3, 10),
#                            Measure = rep(c('Bias', 'RMSE', "Coverage"),
#                                          c(4, 4, 2)),
#                            Value = c(colMeans(est14.array - MC14),
#                                      colMeans(var14.array) - vars14,
#                                      sqrt(colMeans((est14.array - MC14)^2)),
#                                      sqrt(mean((var14.array[,1]-vars14[1])^2)),
#                                      sqrt(mean((var14.array[,2]-vars14[2])^2)),
#                                      mean(ci14.array[, "MC", 'lower'] <= MC14 &
#                                           ci14.array[, "MC", 'upper'] >= MC14),
#                                      mean(ci14.array[,"Mantel",'lower']<=MC14 &
#                                           ci14.array[,"Mantel",'upper']>=MC14)))
#  
#  format(qualities14, digits = 2, scientific = FALSE)
#  

## ----eval = FALSE-------------------------------------------------------------
#  # Simulation using error associated with light-level geolocators
#  # Assign geolocator bias / variance co-variance matrix
#  
#  geoBias <- c(-10000, 50000)
#  geoVCov <- matrix(c(1.2e+8, 0, 0, 1.2e+9), 2, 2)
#  sqrt(diag(geoVCov))
#  
#  targLocs14.2 <- lapply(targLocs14,
#                         FUN = function(x){sf::st_transform(x, 6933)})
#  
#  # convert wintering locations to polygons using the helper function
#  WinteringPolys <- lapply(winteringPos14,
#                           toPoly,
#                           projection.in = 4326,
#                           projection.out = 6933)
#  
#  winteringPolys <- lapply(winteringPos14,
#                           toPoly,
#                           projection.in = 4326,
#                           projection.out = 4326)
#  
#  winteringPolys2 <- lapply(winteringPos14,
#                            toPoly,
#                            projection.in = 4326,
#                            projection.out = 'ESRI:54027')
#  
#  
#  # Simulate capture and non-breeding locations of 100 individuals #
#  
#  nSample14.2 <- 100 # Total number sampled per simulation
#  
#  # How many sampled from each natural population (sampling scenarios separate
#  # from definition scenarios)
#  sampleBreeding14.2 <- list(round(breedingRelN14[[1]]*nSample14.2),
#                           c(rep(0,22),round(breedingRelN14[[5]][1]*nSample14.2),
#                             rep(0, 4),round(breedingRelN14[[5]][2]*nSample14.2),
#                             rep(0,44),round(breedingRelN14[[5]][3]*nSample14.2),
#                             rep(0, 4),round(breedingRelN14[[5]][4]*nSample14.2),
#                             rep(0, 22)),
#                           rep(c(rep(0, 4),
#                                 rep(round(breedingRelN14[[2]][1]*nSample14.2/2),
#                                     2),
#                                 rep(0, 4)), 10))
#  
#  lapply(sampleBreeding14.2, matrix, nrow=10, ncol=10)
#  
#  # Number of simulations to run
#  nSims14.2 <- 10
#  nSimsLarge14.2 <- 25
#  nYears <- 1
#  nMonths <- 1
#  
#  # Set up data structures for storing results
#  
#  #making an empty list to fill
#  animalLoc14 <- vector("list", nScenarios14)
#  sim14.2 <- vector("list", nSimsLarge14.2)
#  
#  compare14.2.array <- array(NA, c(nSimsLarge14.2, 2),
#                             dimnames = list(1:nSimsLarge14.2,
#                                             c("MC", "Mantel")))
#  
#  results14 <- vector("list", nScenarios14)
#  results14.2 <- vector("list", nSimsLarge14.2)
#  
#  # Run simulations
#  set.seed(7)
#  
#  system.time(for (sim in 1:nSimsLarge14.2) {
#    cat("Simulation", sim, "of", nSimsLarge14.2, '\n')
#    sim14.2[[sim]] <- lapply(sampleBreeding14.2,
#                             simMove,
#                             breedingDist = breedDist14[[1]],
#                             winteringDist=nonbreedDist14[[1]],
#                             psi=psi,
#                             nYears=nYears,
#                             nMonths=nMonths)
#  
#  for (i in 13) { # only run one scenario
#    animalLoc14.0 <- sim14.2[[sim]][[scenarioToSampleMap14[i]]]$animalLoc
#    animalLoc14[[i]] <- changeLocations(animalLoc14.0,
#                                        breedingSiteTrans14[[i]],
#                                        winteringSiteTrans14[[i]])
#  
#      breedingTruePoints14 <- capLocs14[[i]][animalLoc14[[i]][,1,1,1], ]
#  
#      winteringTruePoints14 <- targLocs14.2[[i]][animalLoc14.0[,2,1,1], ]
#  
#      results14.2[[sim]] <- simLocationError(winteringTruePoints14,
#                                             WinteringPolys[[i]],
#                                             geoBias,
#                                             geoVCov,
#                                             projection = sf::st_crs(winteringTruePoints14))
#  
#      animalLoc14.2 <- animalLoc14[[i]]
#  
#      animalLoc14.2[ , 2, 1, 1] <- results14.2[[sim]]$targetSample
#  
#      results14[[i]] <- calcPsiMC(breedDist14[[i]],
#                                  nonbreedDist14[[i]],
#                                  animalLoc14.2,
#                                  originRelAbund = breedingRelN14[[i]],
#                                  verbose = FALSE)
#  
#      compare14.2.array[sim, 'MC'] <- results14[[i]]$MC
#  
#      originDist1 <- distFromPos(sf::st_coordinates(breedingTruePoints14))
#  
#      target.point.sample <- sf::st_as_sf(
#        data.frame(results14.2[[sim]]$targetPointSample),
#        coords = c(1,2),
#        crs = 6933)
#  
#      target.point.sample2 <- sf::st_transform(target.point.sample, 4326)
#  
#      targetDist1 <-distFromPos(sf::st_coordinates(target.point.sample2))
#  
#      compare14.2.array[sim, 'Mantel'] <- calcMantel(originDist = originDist1,
#                                                     targetDist = targetDist1)$pointCorr
#    }
#  })
#  
#  means14.2 <- apply(compare14.2.array, 2, mean)
#  vars14.2 <- apply(compare14.2.array, 2, var)
#  rmse14.2 <- apply(compare14.2.array, 2, function(x) sqrt(mean((x - MC14)^2)))
#  
#  # Set up data structures for storing estimation results
#  est14.2.array <- array(NA, c(nSims14.2, 2),
#                         dimnames = list(1:nSims14.2,
#                                         c("MC", "Mantel")))
#  
#  var14.2.array <- array(NA, c(nSims14.2, 2),
#                         dimnames = list(1:nSims14.2,
#                                         c("MC", "Mantel")))
#  
#  ci14.2.array <- array(NA, c(nSims14.2, 2, 2),
#                        dimnames = list(1:nSims14.2,
#                                        c("MC", "Mantel"),
#                                        c('lower', 'upper')))
#  
#  animalLoc14 <- vector("list", nSims14.2) #making an empty list to fill
#  results14 <- vector("list", nSims14.2)
#  
#  # Run estimations
#  
#  set.seed(567)
#  sim14.2.sub <- sim14.2[sample.int(nSimsLarge14.2, nSims14.2, TRUE)]
#  
#  for (sim in 1:nSims14.2) {
#    cat("Estimation", sim, "of", nSims14.2, '\n')
#   for (i in 13) {
#    points0 <- sf::st_as_sf(data.frame(results14.2[[sim]]$targetPointSample),
#                            coords=c("X1","X2"),
#                            crs = 6933)
#  
#    points1 <- sf::st_transform(points0, crs = "ESRI:54027")
#  
#    animalLoc14[[sim]] <- changeLocations(
#      sim14.2.sub[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                          breedingSiteTrans14[[i]],
#                                          winteringSiteTrans14[[i]])
#  
#    results14[[sim]] <- estMC(originDist = breedDist14[[i]],
#                              targetDist = nonbreedDist14[[i]],
#                              originRelAbund = breedingRelN14[[i]],
#                  originPoints = capLocs14[[i]][animalLoc14[[sim]][, 1, 1, 1], ],
#                              targetPoints = points1,
#                              originAssignment = animalLoc14[[sim]][, 1, 1, 1],
#                              targetAssignment = results14.2[[sim]]$targetSample,
#                              nSamples = 1000,
#                              verbose = 0,
#                              calcCorr = TRUE,
#                              geoBias = geoBias,
#                              geoVCov = geoVCov,
#                              isGL = TRUE,
#                              targetSites = winteringPolys2[[i]],
#                              maintainLegacyOutput = TRUE)
#  
#      est14.2.array[sim, 'MC'] <- results14[[sim]]$meanMC
#  
#      est14.2.array[sim, 'Mantel'] <- results14[[sim]]$meanCorr
#  
#      var14.2.array[sim, 'MC'] <- results14[[sim]]$seMC ^ 2
#  
#      var14.2.array[sim, 'Mantel'] <- results14[[sim]]$seCorr ^ 2
#  
#      ci14.2.array[sim, 'MC', ] <- results14[[sim]]$bcCI
#  
#      ci14.2.array[sim, 'Mantel', ] <- results14[[sim]]$bcCICorr
#  
#    }
#  }
#  
#  pointEsts14.2 <- t(sapply(results14, function(x) c(x$pointMC, x$pointCorr)))
#  
#  # Summarize estimation results
#  summary(var14.2.array)
#  vars14.2
#  summary(est14.2.array)
#  summary(pointEsts14.2)
#  colMeans(pointEsts14.2) - MC14
#  sqrt(colMeans((pointEsts14.2 - MC14)^2))
#  means14.2
#  mean(ci14.2.array[, "MC", 'lower'] <= MC14 &
#         ci14.2.array[, "MC", 'upper'] >= MC14)
#  mean(ci14.2.array[, "Mantel", 'lower'] <= MC14 &
#         ci14.2.array[, "Mantel", 'upper'] >= MC14)
#  sqrt(vars14.2) / MC14
#  summary(sqrt(var14.2.array)/est14.2.array)
#  
#  # Plot histograms of bootstrap estimation results
#  est14.2.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2, each = nSims14.2),
#                       Quantity = rep(c('Mean', 'Variance'),
#                                      each = 2 * nSims14.2),
#                       sim = rep(1:nSims14.2, 4),
#                       Estimate = c(est14.2.array, var14.2.array))
#  
#  trues14.2.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2),
#                         Quantity = rep(c('Mean', 'Variance'), each = 2),
#                         Value = c(MC14, MC14, vars14.2))
#  
#  g.est <- ggplot(est14.2.df, aes(Estimate)) + geom_histogram(bins = 15) +
#    facet_grid(Parameter ~ Quantity, scales = 'free_x') +
#    geom_vline(aes(xintercept = Value), data = trues14.2.df) +
#    theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3))
#  
#  g.est
#  
#  
#  # Summary table of bootstrap estimation results
#  qualities14.2 <- data.frame(Parameter = rep(c("MC", 'rM'), 5),
#                              Quantity = rep(rep(c('Mean', 'Variance'),
#                                                 c(2, 2)), 3, 10),
#                              Measure = rep(c('Bias', 'RMSE', "Coverage"),
#                                            c(4, 4, 2)),
#                              Value = c(colMeans(est14.2.array - MC14),
#                                        colMeans(var14.2.array) - vars14.2,
#                                        sqrt(colMeans((est14.2.array - MC14)^2)),
#                                        sqrt(mean((var14.2.array[,1] -
#                                                     vars14[1])^2)),
#                                        sqrt(mean((var14.2.array[,2] -
#                                                     vars14[2])^2)),
#                                        mean(ci14.2.array[,"MC",'lower']<=MC14 &
#                                             ci14.2.array[,"MC",'upper']>=MC14),
#                                mean(ci14.2.array[, "Mantel", 'lower'] <= MC14 &
#                                     ci14.2.array[, "Mantel", 'upper'] >= MC14)))
#  
#  format(qualities14.2, digits = 2, scientific = FALSE)
#  
#  
#  
#  # Make summary table and figures of results with (Example 9) and without
#  # (Example 8) location error
#  
#  qualities.all <- rbind(data.frame(DataType = 'GPS', qualities14),
#                         data.frame(DataType = 'GL', qualities14.2))
#  
#  
#  est14.all.df <- rbind(data.frame(DataType = 'GPS', est.df),
#                         data.frame(DataType = 'GL', est14.2.df))
#  
#  trues14.all.df <- data.frame(DataType = rep(c("GPS", "GL"), 2, each = 2),
#                               Parameter = rep(c("MC", 'rM'), 4),
#                         Quantity = rep(c('Mean', 'Variance'), each = 4),
#                         Value = c(rep(MC14, 4), vars14, vars14.2))
#  
#  g.est <- ggplot(est14.all.df, aes(Estimate)) + geom_histogram(bins = 12) +
#    facet_grid(DataType + Parameter ~ Quantity, scales = 'free_x') +
#    geom_vline(aes(xintercept = Value), data = trues14.all.df) +
#    theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3))
#  
#  g.est
#  
#  

## ----eval = FALSE-------------------------------------------------------------
#  ## Simulation ----
#  nBreeding <- 100
#  nWintering <- 100
#  breedingPos <- matrix(c(rep(seq(-99, -81, 2), each=sqrt(nBreeding)),
#                          rep(seq(49, 31, -2), sqrt(nBreeding))), nBreeding, 2)
#  winteringPos <- matrix(c(rep(seq(-79, -61, 2), each=sqrt(nWintering)),
#                           rep(seq(9, -9, -2), sqrt(nWintering))), nWintering, 2)
#  head(breedingPos)
#  tail(breedingPos)
#  head(winteringPos)
#  tail(winteringPos)
#  
#  breedDist <- distFromPos(breedingPos, 'ellipsoid')
#  nonbreedDist <- distFromPos(winteringPos, 'ellipsoid')
#  
#  # Breeding Abundance
#  breedingN <- rep(5000, nBreeding)
#  breedingRelN <- breedingN/sum(breedingN)
#  
#  # Set up psi matrix
#  o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist,
#                target.dist = nonbreedDist, origin.abund = breedingRelN,
#                sample.size = sum(breedingN),
#                interval = c(0, 10), tol = .Machine$double.eps^0.5)
#  
#  slope <- o$minimum
#  psi <- mlogitMat(slope, breedDist)
#  
#  # Baseline strength of migratory connectivity
#  MC <- calcMC(breedDist, nonbreedDist, breedingRelN, psi, sum(breedingN))
#  MC
#  
#  # Other basic simulation parameters
#  
#  ## Dispersal simulations---
#  set.seed(1516)
#  
#  nYears <- 15
#  
#  nMonths <- 4 # Each season
#  
#  Drates <- c(0.02, 0.04, 0.08, 0.16, 0.32, 0.64)    #rates of dispersal
#  
#  birdLocDisp <- vector('list', length(Drates))
#  
#  Disp.df  <- data.frame(Year=rep(1:nYears, length(Drates)),
#                         Rate=rep(Drates, each = nYears), MC = NA)
#  
#  for(i in 1:length(Drates)){
#    cat('Dispersal Rate', Drates[i], '\n')
#    birdLocDisp[[i]] <- simMove(breedingN,
#                                breedDist,
#                                nonbreedDist,
#                                psi,
#                                nYears,
#                                nMonths,
#                                sumDispRate = Drates[i])
#    for(j in 1:nYears){
#      cat('\tYear', j, '\n')
#      temp.results <- calcPsiMC(breedDist,
#                                nonbreedDist,
#                                breedingRelN,
#                                birdLocDisp[[i]]$animalLoc,
#                                years=j,
#                                verbose = FALSE)
#  
#      Disp.df$MC[j + (i - 1) * nYears] <- temp.results$MC
#    }
#  } # end i loop
#  
#  Disp.df$Year <- Disp.df$Year - 1 #just run once!
#  
#  data.frame(Disp.df, roundMC = round(Disp.df$MC, 2),
#             nearZero = Disp.df$MC < 0.01)
#  
#  
#  # Convert dispersal rates to probabilities of dispersing at least certain distance
#  threshold <- 1000
#  probFarDisp <- matrix(NA, nBreeding, length(Drates), dimnames = list(NULL, Drates))
#  
#  for (i in 1:length(Drates)) {
#    for (k in 1:nBreeding) {
#      probFarDisp[k, i] <-
#        sum(birdLocDisp[[i]]$natalDispMat[k, which(breedDist[k, ]>= threshold)])
#    }
#  }
#  
#  summary(probFarDisp)
#  

