####################################################################################################################
###  Script to simulate spatial capture recapture data and to run CR, CR.dedge, and SCR models on each data set  ###
###  Jesse Whittington
###  Banff National Park
###  2015
####################################################################################################################

library(sp)
library(rgeos)
library(adehabitatHR)
library(coda)
library(jagsUI)
library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
options(stringsAsFactors=FALSE)
wd <- getwd()



### Similar to Gardner et al 2010 and Efford et al. 2013
# Base on sigma
# Track number of detections and number of animals for each simulated data set
# 10 x 10 grid of traps.
# Spacing (s) = sigma.                            # original sigma was 0.8 and spacing was 1 (ie each individual exposed to 2 different traps)
# D = 0.5/sigma2 and 1.0/sigma^2         Density at start
# g0 = 0.5 and 1.0
# State-Space:  buffer traps by 2.5 sigma.  detection probability = 0.04 when g0 = 1 for individuals with activity centres at the edge of the state space (at the nearest trap).
# Simulated 100 data sets per scenario
# Use 5 sampling occasions                                 # originally used.
# Simulate decrease with fixed apparent survival (phi = 0.8) and per capita recruitment (0.1) for T = 3 years
# jagsUI - go until Rhat reaches less than 1.1  (autojags, max iterations = 30000)


###########################################################
####  1.  Function to Simulate Data               #########
###########################################################

fnc.sim.data <- function (D1 = 0.2, K = 5, g0 = 0.1, sigma = 1, buffer.sigma = 5, phi = 0.8, R.n = 0.1, n.years=3, n.max=300) {
  # K = 5 occasions
  # n.max is the maximum number of individuals to include for the analysis.  This includes detected individuals and the augmented population
  # Set up Study Area
  buffer <- buffer.sigma * sigma
  trap.dx <- (1:10) * sigma
  trap.dx <- trap.dx - mean(trap.dx)   # Centre traps around 0
  traps <- cbind(rep(trap.dx, each=10), rep(trap.dx, 10))
  ntraps <- nrow(traps)
  x.lower <- min(traps[, 1] - buffer)
  x.upper <- max(traps[, 1] + buffer)
  y.lower <- min(traps[, 2] - buffer)
  y.upper <- max(traps[, 2] + buffer)
  study.area <- (y.upper - y.lower) * (x.upper - x.lower)
  #  Number of animals alive within state-space each year
  N <- D1 * study.area
  N.avail <- N * n.years * 2 # Create augmented population of available animals
  z <- matrix(0, nrow = N.avail, ncol=n.years)
  z[1:N ,1] <- 1
  N.ever.alive <- N
  N.prev <- sum(z[ ,1])
  for (t in 2:n.years){
    z[ ,t] <- rbinom(N.avail, 1, phi) * z[ , t-1 ]   # animals that live
    n.new <- rbinom(1, N.prev, R.n)
    i.new <- (N.ever.alive + 1): (N.ever.alive + n.new)
    z[i.new, t] <- 1
    N.ever.alive <- N.ever.alive + n.new
    N.prev <- sum(z[ ,t])
  }
  keep <- 1:N.ever.alive
  z <- z[keep, ]                                     # keep only the animals that were ever alive
  Nalive <- N.ever.alive
  sx <- runif(Nalive, x.lower, x.upper)
  sy <- runif(Nalive, y.lower, y.upper)
  S <- cbind(sx, sy)
  D <- sapply(1:ntraps, function(x) sqrt( (sx - traps[ x,1])^2 + (sy - traps[ x,2])^2  )  ) # Matrix: 1 row per animal, 1 column per trap
  alpha1 <- 1/(2 * sigma * sigma)
  probcap <- g0 * exp(-alpha1 * D * D)
  
  # Observation Process
  Y <- array(0, c(N.ever.alive, ntraps, n.years))
  Y.CR <- matrix(0, nrow = N.ever.alive, ncol = n.years)
  n.detect <- rep(0, n.years)
  n.both <- rep(0, n.years)
  n.trap.occ <- K * ntraps
  i.counter <- 0
  for (t in 1:n.years){
    i.counter <- i.counter + 1
    if (t == 1){
      while (n.detect[t] == 0){
        for (i in 1:N.ever.alive) {
          detect <- rbinom(n.trap.occ, 1, probcap[i, ]) * z[i, t]
          detect <- matrix(detect, nrow=K, ncol=ntraps, byrow=T)   # rows = occasions, columns = Traps
          Y[i, ,t] <- apply(detect, 2, sum)                        # number of detections (across occasions) for each trap
          Y.CR[i, t] <- sum(apply(detect, 1, max))                 # Number of occasions the animal was detected.
          #Y[i, ,t] <- rbinom(ntraps, K, probcap[i, ]) * z[i, t]
        }
        n.detect[t] <- sum(Y[ , ,t])
      }  # End of While for t = 1
    } else {  
      while (n.detect[t] == 0 | n.both[t] == 0){
        for (i in 1:N.ever.alive) {
          detect <- rbinom(n.trap.occ, 1, probcap[i, ]) * z[i, t]  #### Probability of trap specific capture is repeated K times.
          detect <- matrix(detect, nrow=K, ncol=ntraps,, byrow=T)
          Y[i, ,t] <- apply(detect, 2, sum)
          Y.CR[i, t] <- sum(apply(detect, 1, max))
          #Y[i, ,t] <- rbinom(ntraps, K, probcap[i, ]) * z[i, t]
        }
        n.detect[t] <- sum(Y[ , ,t])
        n.caps1 <- apply(Y[ , ,t-1], 1, sum)
        n.caps2 <- apply(Y[ , , t], 1, sum)
        both <- n.caps1 * n.caps2
        both <- ifelse(both > 0, 1, 0)
        n.both[t] <- sum(both)
        
      }  # End of While for t > 1
    }  # End of t > 1 
  }
  # 
  n <- n.max
  # SCR data
  obs0 <- array(0, c(n, ntraps, n.years))
  obs0[1:Nalive, , ] <- Y
  Y <- obs0
  # CR data
  obs0 <- matrix(0, nrow=n, ncol=n.years)
  obs0[1:Nalive,  ] <- Y.CR
  Y.CR <- obs0

  dimnames(Y) <- list(1:n, paste("trap", 1:ntraps, sep = ''), paste('year', 1:n.years, sep = ''))
  dimnames(Y.CR) <- list(1:n,paste('year', 1:n.years, sep = '') )
  
  ## Calculate Distance to Edge for CR models
  detect <- apply(Y, c(1,2), sum) # detections by animal and trap
  a1 <- data.frame(animal = rep(1:nrow(detect), ncol(detect)), trap=rep(1:ncol(detect), each=nrow(detect)), n=as.vector(detect)) # down column1, down column2, etc
  a1$detect <- ifelse(a1$n > 0, 1, 0)
  a2 <- filter(a1, n > 0)
  traps.tmp <- data.frame(trap = 1:nrow(traps), x=traps[,1], y=traps[ ,2])
  a3 <- merge(a2, traps.tmp, by='trap', all.x=T, all.y=F)
  detect <- a3
  hr.centre <- data.frame(animal = 1:nrow(Y.CR), x=NA, y=NA) 
  for (i in 1:nrow(Y.CR)){
    b <- filter(detect, animal == i)
    if (nrow(b) > 0) {
      # Composite multi-year Home Range Centre
      i.hull <- chull(cbind(b$x, b$y))   # convex hull
      xy2 <- b[i.hull,]
      if (nrow(xy2) <= 2) {
        x.mean <- mean(xy2$x, na.rm=T)
        y.mean <- mean(xy2$y, na.rm=T)
      } else
      {
        b1 <- Polygon(cbind(c(xy2$x, xy2$x[1]), c(xy2$y, xy2$y[1])))
        b2 <- Polygons(list(b1), ID='1')
        b3 <- SpatialPolygons(list(b2),proj4string=CRS("+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m"))
        b4 <- gCentroid(b3)
        b5 <- coordinates(b4)
        x.mean <- b5[1,1]
        y.mean <- b5[1,2] 
      }
      hr.centre[i, 2:3] <- c(x.mean, y.mean)
    } # End of n.detect > 0
  }  # End of animal
  traps.sp <- SpatialPoints(cbind(traps.tmp$x, traps.tmp$y),proj4string=CRS("+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m"))
  traps.sp <- SpatialPointsDataFrame(traps.sp, data.frame(id = rep(1, nrow(traps.tmp))))  #site.xy)
  a <- mcp(traps.sp, perc=100)
  edge <- as(a, "SpatialLines")
  hr.good <- ifelse(is.finite(hr.centre[ ,2]), 1, 0)
  hr.centre$x <- ifelse(hr.good == 1, hr.centre$x, 0)
  hr.centre$y <- ifelse(hr.good == 1, hr.centre$y, 0)
  hr.centre.sp <- SpatialPoints(cbind(hr.centre$x, hr.centre$y), proj4string=CRS("+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m")) 
  hr.centre.sp <- SpatialPointsDataFrame(hr.centre.sp, hr.centre)
  #plot(traps.sp, col='blue', pch=19)
  #points(hr.centre$x, hr.centre$y, col='red', pch=17)
  dedge <- gDistance(hr.centre.sp, edge, byid=T) 
  dedge[hr.good == 0] <- NA
  dedge <- as.vector(dedge)
  list(y.scr = Y, y.cr=Y.CR, traps = traps, dedge = dedge, x.range = c(x.lower, x.upper), y.range = c(y.lower, y.upper), D1 = D1, N=N, Nalive = Nalive, g0 = g0, alpha1 = alpha1, sigma = sigma, phi=phi, R.n=R.n,K = K)
}

# Test data.
b <- fnc.sim.data(D1 = 0.2, K = 5, g0 = 0.1, sigma = 1, buffer.sigma = 5, phi = 0.8, R.n = 0.1, n.years=3, n.max=200)



#####################################################################
###  2.  Write CR, CR.dedge, and SCR models to be run in JAGS     ###
#####################################################################

sink("model_CR_open.txt")
cat("
  model {
    alpha0 ~ dunif(-10, 10)
    logit(p.mean) <- alpha0
    for(jj in 1:T){             # T = 3 years
      gamma[jj] ~ dunif(0, 1)   #recruitment - conditional on number of available animals
      N[jj] <- sum(z[1:M, jj])
      Recr[jj] <- sum(R[1:M, jj])
    }
    for(jjj in 1:(T-1)){             # T = 3 years
      phi[jjj] ~ dunif(0, 1)         #survival  - use if WANT survival to vary by year
    }
    for (i in 1:M){                      # Loop over individuals
      z[i,1] ~ dbern(gamma[1])                # probability of inclusion in year 1.  Gamma could be calculated later to represent prob of entry given ever alive
      ncaps[i, 1] <- z[i, 1]
      available[i, 1] <- (1 - z[i,1])
      R[i,1] <- z[i,1]                          # New individuals in Year 1
      for (t in 2:T){
        mu[i,t] <- (phi[t-1] * z[i, t-1]) + (gamma[t] * available[i, t-1])
        z[i,t] ~ dbern(mu[i,t])
        ncaps[i,t] <- sum(z[i,1:t])
        available[i,t] <- 1 - step(ncaps[i,t] - 1)
        R[i, t] <- z[i,t] * available[i,t-1]    # New individual in Year t
      }
      alive[i] <- step(ncaps[i, T] - 1)
      for(k in 1:T) {
        p[i,k] <-   p.mean * z[i, k] 
        n.detect[i,k] ~ dbinom(p[i,k], n.checks) 
      }
    }
    Nalive <- sum(alive[1:M])
  }
  ",fill=TRUE)
sink()


sink("model_CR_dedge.txt")
cat("
  model {
    alpha0 ~ dunif(-10, 10)
    a.dedge ~ dunif(-10,10) 
    #  shape.dedge ~ dunif(0.001, 20)         # Set to equal 1
    rate.dedge ~ dunif(0.001, 20)             # Set shape to equal 1 (higher density near zero)
    for(t in 1:T){                            # T = 3 years
      gamma[t] ~ dunif(0, 1)                  # Inclusion probability: - conditional on number of available animals for recruitment
      N[t] <- sum(z[1:M, t])
      Recr[t] <- sum(R[1:M, t])
    }
    for(t in 1:(T-1)){                        # T = 3 years
      phi[t] ~ dunif(0, 1)                    # survival  - use if WANT survival to vary by year
    }
    for (i in 1:M){                           # Loop over individuals
      dedge[i] ~ dgamma(1, rate.dedge)        # dgamma(shape.dedge, rate.dedge)
      z[i,1] ~ dbern(gamma[1])                # probability of inclusion in year 1.
      q[i, 1] <- (1 - z[i,1])
      R[i,1] <- z[i,1]                        # New individuals in Year 1
      for (t in 2:T){
        mu[i,t] <- (phi[t-1] * z[i, t-1]) + (gamma[t] * prod(q[i, 1:(t-1)]) )  # available[i, t-1])
        z[i,t] ~ dbern(mu[i,t])
        q[i, t] <- (1 - z[i,t])
        R[i, t] <- (1-z[i,t-1]) * z[i,t]       # available[i,t-1]    # New individual in Year t
      }
      N.yrs.alive[i] <- sum(z[i, 1:T])
      alive[i] <- 1 - equals(N.yrs.alive[i], 0)
      logit(p.mean[i]) <- alpha0 + a.dedge * dedge[i]
      for(t in 1:T) {
        p[i, t] <- p.mean[i] * z[i, t] 
        n.detect[i,t] ~ dbinom(p[i,t], n.checks) 
      }
    }
    Nalive <- sum(alive[1:M]) 
  }
  ",fill=TRUE)
sink()

sink("model_SCR_open.txt")
cat("
  model {
    sigma ~ dunif(0, 15)                       # sigma 
    inv.sigma2 <- 1 / (2 * sigma * sigma)
    alpha0 ~ dunif(-10,10) 
    logit(g0) <- alpha0
    for(jj in 1:T){                            # T = 3 years
      gamma[jj] ~ dunif(0, 1)                  # recruitment - conditional on number of available animals
      N[jj] <- sum(z[1:M, jj])
      Recr[jj] <- sum(R[1:M, jj])
    }
    for(jjj in 1:(T-1)){                       # T = 3 years
      phi[jjj] ~ dunif(0, 1)                   # survival  - use if WANT survival to vary by year
    }
    
    for (i in 1:M){                            # Loop over individuals
      SX[i]~ dunif(x.lower, x.upper)           # sets the prior on the X and Y- coordinates as uniform
      SY[i]~ dunif(y.lower, y.upper)           # over the given range
      z[i,1] ~ dbern(gamma[1])                 # probability of inclusion in year 1.  Gamma could be calculated later to represent prob of entry given ever alive
      ncaps[i, 1] <- z[i, 1]
      available[i, 1] <- (1 - z[i,1])
      R[i,1] <- z[i,1]                          # New individuals in Year 1
      for (t in 2:T){
        mu[i,t] <- (phi[t-1] * z[i, t-1]) + (gamma[t] * available[i, t-1])
        z[i,t] ~ dbern(mu[i,t])
        ncaps[i,t] <- sum(z[i,1:t])
        available[i,t] <- 1 - step(ncaps[i,t] - 1)
        R[i, t] <- z[i,t] * available[i,t-1]    # New individual in Year t
      }
      alive[i] <- step(ncaps[i, T] - 1)
      for (j in 1:ntraps){
        Dist2[i,j] <- pow(SX[i] - x[j], 2) + pow(SY[i] - y[j], 2)
        pmean[i,j] <- g0 * exp(-Dist2[i,j] * inv.sigma2) 
        for (k in 1:T){
          p[i,j,k] <- pmean[i, j] * z[i, k]
          n.detect[i,j,k] ~ dbin(p[i,j,k], n.checks)
        }
      }  
    }
    Nalive <- sum(alive[1:M]) 
  }
  ",fill=TRUE)
sink()


########################################################
###  3.  Function to Save Results from Models       ####
########################################################

fnc.CI <- function(model, method='CR', scenario=1, D1=0.5, g0=0.1){   # model is an mcmc.list - result from jagsUI
  a.samples <- model$samples
  a1 <- lapply(a.samples, function(x) as.data.frame(x))
  a2 <- rbind_all(a1)
  a2$Lambda.N1 <- a2$'N[2]'/a2$'N[1]'
  a2$Lambda.N2 <- a2$'N[3]'/a2$'N[2]'
  a2$Lambda <- exp( (log(a2$Lambda.N1) + log(a2$Lambda.N2)) * 0.5)
  a2$R1 <- a2$'Recr[2]' / a2$'N[1]'
  a2$R2 <- a2$'Recr[3]' / a2$'N[2]'
  a2$Lambda.R1 <- a2$'phi[1]' + a2$R1
  a2$Lambda.R2 <- a2$'phi[2]' + a2$R2
  a2$Lambda.R <- exp( (log(a2$Lambda.R1) + log(a2$Lambda.R2)) * 0.5)
  m.CI <- data.frame(HPDinterval(as.mcmc(a2)))
  m.CI$Parameter <- row.names(m.CI)
  m.med <- apply(as.matrix(a2), 2, median)
  m.mean <- apply(as.matrix(a2), 2, mean)
  m.max <- apply(as.matrix(a2), 2, max)
  m.med <- data.frame(Parameter= names(m.med), Median = m.med, Mean = m.mean, Max = m.max, scenario = scenario, D1 = D1, g0 = g0)
  m.med <- merge(m.med, m.CI, by='Parameter')
  a1 <- model$summary
  a1 <- data.frame(Parameter = rownames(a1), Rhat=a1[ ,colnames(a1) == 'Rhat'])
  df <- merge(m.med, a1, by='Parameter', all.x=T, all.y=F)
  df$n.samples <- model$mcmc.info$n.samples
  df$Method <- method
  return(df)  
}


######################################################
### 4.  Run Simulations                           ####
######################################################

# NOTE:  The script is currently set to simulate 5 dataset for each scenario and run 3 mcmc chains with a burn in of 10 and then another 10 sampling iterations.
#        In the manuscript we simulated 100 datasets for each scenario and ran autojags with an initial burnin of 4000, niterations = 5000, and max iterations = 10000
#        The simulations in the manuscript took many, many days to run spread across many, many computers.
wd <- getwd()
wd.save <- paste0(wd, '/.RData')
n.reps <- 5                                     # NUMBER OF SIMULATED DATA SETS
nb <- 10
ni <- 20

# Set up data.frame of scenarios
a <- data.frame(expand.grid(data.set = 1:n.reps, D1=c(0.5, 1.0), g0= c(0.1, 0.5)))
a <- arrange(a, data.set, D1, g0)
a1 <- data.frame( scenario=1:nrow(a), a)

tbl.nmax <- data.frame(D1 = c(0.5, 1.0), 
                       n.max= c(175, 300),        # Maximum number of individuals (including augmented population) for SCR models - set this after trial and error looking at Max Nalive
                       n.max.cr = c(140, 240))    # Maximum number of individuals for CR models (usually fewer than required for SCR)
a1 <- merge(a1, tbl.nmax, by=c('D1'))
a1 <- arrange(a1, data.set, D1, g0)
a1$n.detect <- 0
a1$n.sims <- 0
scenarios <- a1

# Create lists to hold data and model results.
sim.data.list <- vector(mode='list', nrow(scenarios))
result.list.SCR <- vector(mode='list', nrow(scenarios))
result.list.CR <- vector(mode = 'list', nrow(scenarios))
result.list.CR.dedge <- vector(mode = 'list', nrow(scenarios))
times <- data.frame(scenarios, time.start = c(now(), rep(NA, nrow(scenarios) -1)), time.end = c(now(), rep(NA, nrow(scenarios) -1)) )

scen <- scenarios                                 # Temporary data.frame with a shorter name

# Generate Data
for (i in 1:nrow(scenarios)){  
  sim.data.list[[i]] <- fnc.sim.data(D1 = scen$D1[i], K = 5, g0 = scen$g0[i], sigma = 1.0, buffer = 2.5, phi=0.8, R.n = 0.1, n.years = 3, n.max=scen$n.max[i])
}


for (i in 1:nrow(scenarios)){
  print(i)
  print(now())
  times$time.start[i] <- now()
  df <- sim.data.list[[i]]
  y.scr <- df$y.scr
  zst <- apply(y.scr, c(1, 3), max)                               # Initial values for alive/dead
  zst[zst > 0] <- 1
  nyr <- ncol(zst)
  if (ncol(zst) > 2){                                             # Fill in gap years:  eg. animal alive in years 1 and 3 but not detected in year 2
    for (t in 2:(nyr-1)){
      z.prev <- rowSums(data.frame(zst[ ,1:(t-1)]))
      z.prev[z.prev > 0] <- 1
      z.post <- rowSums(data.frame(zst[ ,(t+1):nyr]))
      z.post[z.post > 0] <- 1
      zst[ ,t] <- ifelse(z.prev + z.post == 2, 1, zst[ ,t])
    }  
  }
  n.detect <- apply(zst, 1, max)
  scenarios$n.detect[i] <- sum(n.detect)
  detect <- apply(df$y.scr, c(1,2), sum)                          # Number of detections per trap
  detect[detect > 0] <- 1
  s.xy <- matrix(NA, nrow=nrow(y.scr), 2)                         # Calculate home range centre
  for (ii in 1:nrow(s.xy)){
    if (sum(detect[ii, ]) > 0){
      s.xy[ii, 1] <- sum(detect[ii, ] * df$traps[ ,1]) / sum(detect[ii, ])
      s.xy[ii, 2] <- sum(detect[ii, ] * df$traps[ ,2]) / sum(detect[ii, ])
    } else {
      s.xy[ii, 1] <- runif(1, df$x.range[1], df$x.range[2])
      s.xy[ii, 2] <- runif(1, df$y.range[1], df$y.range[2])
    }
  }
  
  # Spatial Capture Recapture
  data <- list(M = nrow(df$y.scr), n.detect = df$y.scr, ntraps=nrow(df$traps), x = df$traps[,1], y = df$traps[ ,2], n.checks = 5, T=3,
               x.lower = df$x.range[1], x.upper = df$x.range[2], y.lower = df$y.range[1], y.upper = df$y.range[2])
  
  inits <- function(i){list( z=zst, gamma=runif(3, 0.05, 0.2), phi=runif(2, 0.8, 0.95), SX=s.xy[ ,1], SY=s.xy[ ,2] ,
                             sigma=1, alpha0 = runif(1, -2.3, -2.1) )}
  params <- c("sigma", "alpha0",  'N', 'Recr', "phi", "Nalive" )
  model <- "model_SCR_open.txt"
  #m <- autojags(data = data, inits = inits, parameters.to.save = params, model.file = model, n.chains = 3, iter.increment = 5000, n.burnin = 4000, parallel = TRUE, Rhat.limit = 1.1, max.iter = 10000)
  m <- jags(data = data, inits = inits, parameters.to.save = params, model.file = model, n.chains = 3, n.burnin = nb, n.iter=ni, parallel = TRUE)  
  result.list.SCR[[i]] <- fnc.CI(m, method='SCR', scenario=i, D1=scen$D1[i], g0 = scen$g0[i])
  rm(m, data, params, model, inits)
  
  ## Capture-Recapture - Distance to edge
  y.cr <- df$y.cr
  y.cr <- y.cr[1:scen$n.max.cr[i], ]    # require a smaller population of augmented individuals for CR
  zst <- zst[1:scen$n.max.cr[i], ]
  
  data <- list(M = nrow(y.cr), n.detect = y.cr, n.checks = 5, T=3, dedge=df$dedge)
  inits <- function(i){list( z=zst, gamma=runif(3, 0.05, 0.2), phi=runif(2, 0.88, 0.92), alpha0 = runif(1, -2.3, -2.1 ), 
                             a.dedge=runif(1, 0, 5),   rate.dedge=runif(1, 1, 3)) }   
  params <- c('alpha0', 'a.dedge', 'shape.dedge', 'rate.dedge',  'N', 'Recr', "phi" , "Nalive")
  model <- "model_CR_dedge.txt"
  #m <- autojags(data = data, inits = inits, parameters.to.save = params, model.file = model, n.chains = 3, iter.increment = 5000, n.burnin = 4000, parallel = TRUE, Rhat.limit = 1.1, max.iter = 10000)
  m <- jags(data = data, inits = inits, parameters.to.save = params, model.file = model, n.chains = 3, n.burnin = nb, n.iter=ni, parallel = TRUE)  
  result.list.CR.dedge[[i]] <- fnc.CI(m, method='CR.dedge', scenario=i, D1=scen$D1[i], g0 = scen$g0[i])
  rm(m, data, params, model, inits)
  
  ## Capture-Recapture
  data <- list(M = nrow(y.cr), n.detect = y.cr, n.checks = 5, T=3)
  inits <- function(i){list( z=zst, gamma=runif(3, 0.05, 0.2), phi=runif(2, 0.88, 0.92), alpha0 = runif(1, -2.3, -2.1 )) } # shape.dedge=runif(1, 1, 3),
  params <- c('alpha0', 'N', 'Recr', "phi" , "Nalive")
  model <- "model_CR_open.txt"
  #m <- autojags(data = data, inits = inits, parameters.to.save = params, model.file = model, n.chains = 3, iter.increment = 5000, n.burnin = 4000, parallel = TRUE, Rhat.limit = 1.1, max.iter = 10000)
  m <- jags(data = data, inits = inits, parameters.to.save = params, model.file = model, n.chains = 3, n.burnin = nb, n.iter=ni, parallel = TRUE)  
  result.list.CR[[i]] <- fnc.CI(m, method='CR', scenario=i, D1=scen$D1[i], g0 = scen$g0[i])
  rm(m, data, params, model, inits)
  rm(df, y.scr, y.cr, detect, zst, s.xy)
  
  times$time.end[i] <- now()
  times$dhours <- (times$time.end - times$time.start) / dhours(1)
  write.csv(times, paste0('Times ', as.Date(now()), '.csv'), row.names=F, na='')   # Handy to look at to determine how long simulations will take to run
  save.image(wd.save)
}


##########################################################
###  4.  Append model results and a couple quick plots ###
##########################################################

i.combine <- list(result.list.SCR, result.list.CR.dedge, result.list.CR)
a <- lapply(i.combine,function(x) rbind_all(x))
a <- rbind_all(a)
a <- mutate(a, scenario2 = ifelse(Method == 'SCR', scenario - 0.2, ifelse(Method == 'CR', scenario + 0.2, scenario) ))
a <- mutate(a, D.f = ifelse(D1 == 0.5, 'Density = 0.5', 'Density = 1.0'),
            g0.f = paste('g0 =', g0), D.g0 = paste('D =', D1, '\ng0 = ', g0))
a <- mutate(a, Parameter2 = gsub('1', '', Parameter), Parameter2 = gsub('2', '', Parameter2),
            Parameter2 = gsub('[', '', Parameter2, fixed=T), Parameter2 = gsub(']', '', Parameter2, fixed=T))
results.all <- a


### Boxplots  and Plot Median and 95% HPD for each simulated dataset
a <- results.all
# Apparent survival
a1 <- filter(a, Parameter %in% c('phi[1]', 'phi[2]'))
ggplot(a1, aes(Method, Median)) + geom_hline(yintercept=0.8, linetype=2) +  geom_boxplot(fill='lightblue', alpha=0.5) +facet_grid( D.f ~ g0.f) + ggtitle('Apparent Survival')
ggplot(a1, aes(scenario2, Median, ymin=lower, ymax=upper, col=Method)) + geom_hline(yintercept = 0.8, linetype=2) + geom_linerange() +  geom_point() + facet_grid(D.g0~Parameter) + ggtitle('Apparent Survival')

# Recruitment
a1 <- filter(a, Parameter %in% c('R1', 'R2'))
ggplot(a1, aes(Method, Median)) + geom_hline(yintercept=0.1, linetype=2) +  geom_boxplot(fill='lightblue', alpha=0.5) + facet_grid( D.f ~ g0.f) + ggtitle('Per Capita Recruitment')
ggplot(a1, aes(scenario2, Median, ymin=lower, ymax=upper, col=Method)) + geom_hline(yintercept = 0.1, linetype=2) + geom_linerange() +  geom_point() + facet_grid(D.g0~Parameter) + ggtitle('Per Capita Recruitment')

# Growth Rate
a1 <- filter(a, Parameter %in% c('Lambda.R1', 'Lambda.R2'))
ggplot(a1, aes(Method, Median)) + geom_hline(yintercept = 0.9, linetype=2) + geom_boxplot(fill = 'lightblue', alpha=0.5) + facet_grid( D.f ~ g0.f) + ggtitle('Growth Rate')
ggplot(a1, aes(scenario2, Median, ymin=lower, ymax=upper, col=Method)) + geom_hline(yintercept = 0.9, linetype=2) + geom_linerange() +  geom_point() + facet_grid(D.g0~Parameter) + ggtitle('Growth Rate')