---
title: '''MigConnectivity'' package: isotope simulation'
author: "Michael T. Hallworth, Jeffrey A. Hostetler"
date: '`r Sys.Date()`'
output:
  rmarkdown::html_vignette: default
  rmarkdown::pdf_document: default
vignette: |
  %\usepackage[utf8]{inputenc}
  %\VignetteIndexEntry{'MigConnectivity' package: isotope simulation}
  %\VignetteEngine{knitr::rmarkdown}
---
```{r, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, 
                      comment = "#>", 
                      tidy = FALSE,
                      message = FALSE,
                      warning = FALSE)
```

Load required packages
```{r eval = FALSE}
library(terra)
library(MigConnectivity)
library(sf)
```

The following is a simulation that tests the how the spatial arrangement of target sites influences MC from stable-hydrogen isotopes. The following simulation is run using data generated within the code but we use the Ovenbird as an example species. 

### Ovenbird distribution
Read in the Ovenbird distribution and create a species distribution map from the abundance data. 
```{r eval = FALSE}
# read in raster layer
download.file(
"https://raw.githubusercontent.com/SMBC-NZP/MigConnectivity/master/data-raw/Spatial_Layers/bbsoven.txt",
              destfile = "bbsoven.txt")

OVENabund <- terra::rast("bbsoven.txt")

OVENdist <- OVENabund
OVENdist[OVENdist>0]<-1
OVENdist[OVENdist==0]<-NA

OVEN_single_poly <- terra::as.polygons(OVENdist)
terra::crs(OVEN_single_poly) <- MigConnectivity::projections$WGS84
terra::crs(OVENabund) <- MigConnectivity::projections$WGS84
```
To complete the simulation we need a template to ensure the raster resolution is the same as the assignment raster. To do this, we use the isotope data as our template. We grab the isotope base-map using the <code>getIsoMap</code> function. 
```{r eval = FALSE}
### get isotope raster for template
rasterTemplate <- getIsoMap()
```

### Crop template to distribution and mask with distribution
```{r eval = FALSE}
terra::crs(rasterTemplate) <- MigConnectivity::projections$WGS84
rasterTemplate <- terra::crop(terra::mask(rasterTemplate, 
                                          OVEN_single_poly),
                              OVEN_single_poly)

# rasterize the distribution for relative abundance so that raster
# dimensions and resolution match the isotope layer
relativeAbund <- terra::project(OVENabund,rasterTemplate)
relativeAbund  <- relativeAbund / 
  unlist(c(terra::global(relativeAbund, fun = "sum", na.rm = T)))
```

## Generate target sites
The simulation is focused on the effect of target sites on MC when using stable-hydrogen isotopes. The following code generates various target site layers used in the simulation.
```{r eval = FALSE}
# generate target sites
targetRanges <- vector('list',5)
# 3' latitude
targetRanges[[1]] <- terra::as.polygons(terra::rast(xmin = -180,
                                         xmax = -40,
                                         ymin = 25,
                                         ymax = 85,
                                         resolution = c(140,3)))

# 5' latitude
targetRanges[[2]] <- terra::as.polygons(terra::rast(xmin = -180,
                                             xmax = -40,
                                             ymin = 25,
                                             ymax = 85,
                                             resolution = c(140,5)))

# 10' latitude
targetRanges[[3]] <- terra::as.polygons(terra::rast(xmin = -180,
                                         xmax = -40,
                                         ymin = 25,
                                         ymax = 85,
                                         resolution = c(140,10)))

# 12 isotope units
featherIso <- (0.95*getIsoMap(period = "GrowingSeason")+(-17.57))
iso <- terra::crop(featherIso, terra::ext(c(-180,-40,25,85)))
isocut <- terra::classify(iso, 
                          rcl = seq(terra::global(iso, fun = "min", 
                                                  na.rm = TRUE)$min,
                                    terra::global(iso, fun = "max", 
                                                  na.rm = TRUE)$max,12))
targetRanges[[4]] <- terra::as.polygons(isocut)


# 12*2 isotope units
isocut <- terra::classify(iso, 
                          rcl = seq(terra::global(iso, fun = "min", 
                                                  na.rm = TRUE)$min, 
                                    terra::global(iso,fun = "max",
                                                  na.rm = TRUE)$max, 24))
targetRanges[[5]] <- terra::as.polygons(isocut)


TR <- lapply(targetRanges, sf::st_as_sf)
vTR <- lapply(TR, sf::st_make_valid)
vTR <- lapply(vTR, FUN = function(x){sf::st_set_crs(x, 4326)})

sf_use_s2(FALSE)
OVEN_single_poly <- sf::st_as_sf(OVEN_single_poly) %>% sf::st_make_valid()


# Keep only the targetSites that intersect with the OVEN polygon
targetRanges <- lapply(vTR, FUN = function(x){sf::st_intersection(x,OVEN_single_poly)})

targetRanges <- lapply(targetRanges, FUN = function(x){y <- sf::st_as_sf(x)
                                                       z <- sf::st_transform(y,4326)
                                                       return(z)})

for(i in 1:5){
  targetRanges[[i]]$Target <- 1:nrow(targetRanges[[i]])
}
targetRanges <- lapply(targetRanges, sf::st_make_valid)
sf_use_s2(TRUE)
```

#### Generate simulated data
```{r, eval = FALSE}
#Generate random breeding locations using the 10' target sites
Site1 <- st_as_sf(sf::st_sample(targetRanges[[3]][1:2,], size =100, type = "random"))
Site2 <- st_as_sf(sf::st_sample(targetRanges[[3]][2:3,], size =100, type = "random"))
Site3 <- st_as_sf(sf::st_sample(targetRanges[[3]][2:4,], size =100, type = "random"))

# Capture coordinates
capCoords <- array(NA,c(3,2))
capCoords[1,] <- cbind(-98.17,28.76)
capCoords[2,] <- cbind(-93.70,29.77)
capCoords[3,] <- cbind(-85.000,29.836)

featherIso <- (0.95*getIsoMap(period = "GrowingSeason")+(-17.57))

# Extract simulated data
iso_dat <- data.frame(Site = rep(1:3, each = 100),
                      xcoords = c(sf::st_coordinates(Site1)[,1], 
                                  sf::st_coordinates(Site2)[,1], 
                                  sf::st_coordinates(Site3)[,1]),
                      ycoords = c(sf::st_coordinates(Site1)[,2], 
                                  sf::st_coordinates(Site2)[,2],
                                  sf::st_coordinates(Site3)[,2]),
                      targetSite = unlist(unclass(sf::st_intersects(rbind(Site1,
                                                                          Site2,
                                                                          Site3),
                                        targetRanges[[3]]))),
                      featherIso = terra::extract(featherIso,
                                                  terra::vect(rbind(Site1,Site2,
                                                                    Site3))))

iso_dat <- iso_dat[complete.cases(iso_dat),]

# generate transition data from simulation
sim_psi <- table(iso_dat$Site, iso_dat$targetSite)

for(i in 1:nrow(sim_psi)){
  sim_psi[i,]<-sim_psi[i,]/sum(sim_psi[i,])
}

states <- rnaturalearth::ne_states(country = "United States of America", 
                                   returnclass = "sf")
originSites <- states[(states$woe_name %in% c("Texas","Louisiana","Florida")),]
originSites <- sf::st_transform(originSites, 4326)
                               
originDist <- distFromPos(st_coordinates(sf::st_centroid(originSites,
                                                            byid = TRUE)))
targetDistMC <- distFromPos(sf::st_coordinates(sf::st_centroid(targetRanges[[3]], 
                                                              byid = TRUE)))
```

## Start of simulations

*warning* this takes a long time to run. 

```{r, eval = FALSE}
originRelAbund <- rep(1/3, 3)
nTargetSetups <- 5
nSims <- 2            #SET LOW FOR EXAMPLE
# nSims <- 200        
nOriginSites = 3

targetPoints0 <- matrix(NA, 300, 2, dimnames = list(NULL, c("x", "y")))


a <- Sys.time()
output.sims <- vector("list", nTargetSetups)
for(target in 1:nTargetSetups){
  sim.output <- data.frame(targetSetup = rep(NA,nSims),
                           sim = rep(NA,nSims),
                           MC.generated = rep(NA,nSims),
                           MC.realized = rep(NA,nSims),
                           MC.est = rep(NA,nSims),
                           MC.low = rep(NA,nSims),
                           MC.high = rep(NA,nSims),
                           rM.realized = rep(NA,nSims),
                           rM.est = rep(NA,nSims),
                           rM.low = rep(NA,nSims),
                           rM.high = rep(NA,nSims))
  set.seed(9001)
  targetSites <- targetRanges[[target]]
  sf_use_s2(FALSE)
  targetSites <- sf::st_make_valid(targetSites)
  targetDist <- distFromPos(sf::st_coordinates(sf::st_centroid(targetSites, 
                                                               byid = TRUE)))

  # Extract simulated data
  iso_dat <- data.frame(Site = rep(1:3,each = 100),
                        xcoords = c(sf::st_coordinates(Site1)[,1],
                                    sf::st_coordinates(Site2)[,1],
                                    sf::st_coordinates(Site3)[,1]),
                        ycoords = c(sf::st_coordinates(Site1)[,2],
                                    sf::st_coordinates(Site2)[,2],
                                    sf::st_coordinates(Site3)[,2]),
                        targetSite = unlist(unclass(sf::st_intersects(rbind(Site1,
                                                                            Site2,
                                                                            Site3),
                                        targetSites))),
                        featherIso = terra::extract(featherIso,
                                                    terra::vect(rbind(Site1,
                                                                      Site2,
                                                                      Site3))))

  iso_dat <- iso_dat[complete.cases(iso_dat),]

  # generate transition data from simulation
  sim_psi <- prop.table(table(iso_dat$Site,factor(iso_dat$targetSite,
                                                  1:nrow(targetSites))), 1)
  
  sf_use_s2(TRUE)
  MC.generated <- calcMC(originDist = originDist,
                          targetDist = targetDist,
                          originRelAbund = originRelAbund,
                          psi = sim_psi)

  for (sim in 1:nSims) {
    cat("Simulation Run", sim, "of", nSims, "for target",target,"at", date(), "\n")
    sim_move <- simMove(rep(100, nOriginSites), originDist, targetDist, sim_psi, 1, 1)
    originAssignment <- sim_move$animalLoc[,1,1,1]
    targetAssignment <- sim_move$animalLoc[,2,1,1]
    sf_use_s2(FALSE)
    for (i in 1:300) {
      targetPoint1 <- sf::st_sample(targetSites[targetAssignment[i],], 
                               size = 1, type = "random", iter = 25)
      targetPoints0[i,] <- sf::st_coordinates(targetPoint1)
    }
    targetPoints <- sf::st_as_sf(as.data.frame(targetPoints0), crs = 4326,
                                 coords = c("x", "y"))
  
    # Extract simulated data
    iso_dat <- data.frame(Site = originAssignment,
                          xcoords = targetPoints0[,1],
                          ycoords = targetPoints0[,2],
                          targetSite = unlist(unclass(sf::st_intersects(targetPoints,
                                                                        targetSites))),
                          featherIso = terra::extract(featherIso,
                                                      terra::vect(targetPoints)))

  
    iso_dat <- iso_dat[complete.cases(iso_dat),]
  
    # generate transition data from simulation
    sim_psi0 <- table(iso_dat$Site,
                      factor(iso_dat$targetSite,
                             min(targetSites$Target):max(targetSites$Target)))
    sim_psi_realized <- prop.table(sim_psi0, 1)
  
    # get points ready for analysis
    nSite1 <- table(iso_dat$Site)[1]
    nSite2 <- table(iso_dat$Site)[2]
    nSite3 <- table(iso_dat$Site)[3]
  
    nTotal <- nSite1+nSite2+nSite3
  
    originCap <- array(NA, c(nTotal,2))
  
    wherecaught <- rep(paste0("Site", 1:3), c(nSite1, nSite2, nSite3))
  
    originCap[which(pmatch(wherecaught,"Site1",dup = TRUE)==1),1] <- capCoords[1,1]
    originCap[which(pmatch(wherecaught,"Site1",dup = TRUE)==1),2] <- capCoords[1,2]
  
    originCap[which(pmatch(wherecaught,"Site2",dup = TRUE)==1),1] <- capCoords[2,1]
    originCap[which(pmatch(wherecaught,"Site2",dup = TRUE)==1),2] <- capCoords[2,2]
  
    originCap[which(pmatch(wherecaught,"Site3",dup = TRUE)==1),1] <- capCoords[3,1]
    originCap[which(pmatch(wherecaught,"Site3",dup = TRUE)==1),2] <- capCoords[3,2]
  
    originPoints <- sf::st_as_sf(as.data.frame(originCap), crs = 4326, coords = 1:2)
    sf_use_s2(TRUE)

    MC.realized <- calcMC(originDist = originDist,
                          targetDist = targetDist,
                          originRelAbund = originRelAbund,
                          psi = sim_psi_realized,
                          sampleSize=nTotal)
  
  
  
    originPointDists <- distFromPos(originCap)
    targetPointDists <- distFromPos(cbind(iso_dat$xcoords, iso_dat$ycoords))
  
  
    simAssign <- isoAssign(isovalues = iso_dat$featherIso.GrowingSeasonD,
                           isoSTD = 12,
                           intercept = -17.57,
                           slope = 0.95,
                           odds = 0.67,
                           restrict2Likely = FALSE,
                           nSamples = 500,
                           sppShapefile = OVEN_single_poly,
                           relAbund = relativeAbund,
                           verbose = 2,
                           isoWeight = -0.7,
                           abundWeight = 0,
                           assignExtent = c(-179,-60,15,89),
                           element = "Hydrogen",
                           period = "GrowingSeason")
    psi <- estTransition(targetRaster = simAssign,
                        targetSites = targetSites,
                        originPoints = originPoints,
                        originSites = originSites, isRaster = TRUE,
                        nSamples = 100, nSim = 5, verbose = 0)
    rM <- estMantel(targetRaster = simAssign,
                    targetSites = targetSites,
                    originPoints = originPoints, 
                    isGL = FALSE, isTelemetry = FALSE,
                    originSites = originSites, isRaster = TRUE,
                    nBoot = 100, nSim = 5, verbose = 0, captured = "origin")
    simEst <- estStrength(originRelAbund = originRelAbund,
                          targetDist = targetDist,
                          psi = psi,
                          originDist = originDist,
                          nSamples = 100,
                          verbose = 0,
                          alpha = 0.05)
  
    sim.output$targetSetup[sim] <- target
    sim.output$sim[sim]<-sim
    sim.output$MC.generated[sim] <- MC.generated
    sim.output$MC.realized[sim] <- MC.realized
    sim.output$MC.est[sim] <- simEst$MC$mean
    sim.output$MC.low[sim] <- simEst$MC$bcCI[1]
    sim.output$MC.high[sim] <- simEst$MC$bcCI[2]
    sim.output$rM.realized[sim] <- calcMantel(originDist = originPointDists, 
                                              targetDist = targetPointDists)$pointCorr
    sim.output$rM.est[sim] <- rM$corr$mean
    sim.output$rM.low[sim] <- rM$corr$bcCI[1]
    sim.output$rM.high[sim] <- rM$corr$bcCI[2]
  
  }
  output.sims[[target]] <- sim.output
}
Sys.time()-a

output.sims <- do.call("rbind", output.sims)
#
```
Summarize the output 

```{r, eval = FALSE}
sim.output <- transform(output.sims,
                        MC.generated.cover = as.integer((MC.low <= MC.generated & 
                                                           MC.high >= MC.generated)),
                        MC.realized.cover = as.integer((MC.low <= MC.realized & 
                                                          MC.high >= MC.realized)),
                        MC.generated.error = MC.est - MC.generated, 
                        MC.realized.error = MC.est - MC.realized,
                        rM.cover = as.integer((rM.low <= rM.realized & 
                                                 rM.high >= rM.realized)),
                        rM.error = rM.est - rM.realized)

summary(sim.output)
# Examine results
aggregate(MC.generated.error ~ targetSetup, sim.output, mean)
aggregate(MC.generated.error ~ targetSetup, sim.output, function (x) mean(abs(x)))
aggregate(MC.generated.cover ~ targetSetup, sim.output, mean)
aggregate(MC.realized.error ~ targetSetup, sim.output, mean)
aggregate(MC.realized.error ~ targetSetup, sim.output, function (x) mean(abs(x)))
aggregate(MC.realized.cover ~ targetSetup, sim.output, mean)
aggregate(I(MC.realized - MC.generated) ~ targetSetup, sim.output, mean)
aggregate(rM.error ~ targetSetup, sim.output, mean)
aggregate(rM.error ~ targetSetup, sim.output, function (x) mean(abs(x)))
aggregate(rM.cover ~ targetSetup, sim.output, mean)
aggregate(MC.est ~ targetSetup, sim.output, mean)
```