# set desired directory to save model outputs
setwd("C:/this_ist_just_an_example/results")

# define path to folder containing input data like capture histories, effort matrix
input_files = "C:/this_ist_just_an_example/input_files"

#### load required libraries ####
package.list=c("gdata",
               "readxl",
               "camtrapR",
               "optimbase",
               "plyr",
               "raster",
               "rgeos",
               "maptools",
               "stringr",
               "snowfall",
               "rjags",
               "rlecuyer",
               "abind")

for (package in package.list) {
  if (!require(package, character.only = T, quietly = T)) {
    install.packages(package)
    library(package, character.only = T)
  }
}

# load covariates table 
covtable <- read.csv(paste(input_files, "covariates.csv", sep = "/"))

cov.list <- colnames(covtable)[4:20]

# load detection histories for leeches.
# if there were no detections of a species then load empty matrix
detecthist1_l <- read.csv(paste(input_files, "spec_01_leech.csv", sep = "/"), row.names = 1)
detecthist2_l <- read.csv(paste(input_files, "spec_02_leech.csv", sep = "/"), row.names = 1)
detecthist3_l <- read.csv(paste(input_files, "spec_03_leech.csv", sep = "/"), row.names = 1)
detecthist4_l <- read.csv(paste(input_files, "spec_04_leech.csv", sep = "/"), row.names = 1)
detecthist5_l <- read.csv(paste(input_files, "spec_05_leech.csv", sep = "/"), row.names = 1)
detecthist6_l <- read.csv(paste(input_files, "spec_06_leech.csv", sep = "/"), row.names = 1)
detecthist7_l <- read.csv(paste(input_files, "spec_07_leech.csv", sep = "/"), row.names = 1)
detecthist8_l <- read.csv(paste(input_files, "spec_08_leech.csv", sep = "/"), row.names = 1)
detecthist9_l <- read.csv(paste(input_files, "EMPTY_DETECT.csv", sep = "/"), row.names = 1)
detecthist10_l <- read.csv(paste(input_files, "spec_10_leech.csv", sep = "/"), row.names = 1)
detecthist11_l <- read.csv(paste(input_files, "spec_11_leech.csv", sep = "/"), row.names = 1)
detecthist12_l <- read.csv(paste(input_files, "spec_12_leech.csv", sep = "/"), row.names = 1)
detecthist13_l <- read.csv(paste(input_files, "spec_13_leech.csv", sep = "/"), row.names = 1)
detecthist14_l <- read.csv(paste(input_files, "EMPTY_DETECT.csv", sep = "/"), row.names = 1)
detecthist15_l <- read.csv(paste(input_files, "EMPTY_DETECT.csv", sep = "/"), row.names = 1)
detecthist16_l <- read.csv(paste(input_files, "EMPTY_DETECT.csv", sep = "/"), row.names = 1)
detecthist17_l <- read.csv(paste(input_files, "spec_17_leech.csv", sep = "/"), row.names = 1)
detecthist18_l <- read.csv(paste(input_files, "spec_18_leech.csv", sep = "/"), row.names = 1)
detecthist19_l <- read.csv(paste(input_files, "spec_19_leech.csv", sep = "/"), row.names = 1)
detecthist20_l <- read.csv(paste(input_files, "EMPTY_DETECT.csv", sep = "/"), row.names = 1)
detecthist21_l <- read.csv(paste(input_files, "spec_21_leech.csv", sep = "/"), row.names = 1)
detecthist22_l <- read.csv(paste(input_files, "EMPTY_DETECT.csv", sep = "/"), row.names = 1)
detecthist23_l <- read.csv(paste(input_files, "spec_23_leech.csv", sep = "/"), row.names = 1)

# load detection histories for camera traps
# if there were ne detections of a species then load empty matrix
detecthist1_c <- read.csv(paste(input_files, "spec_01__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist2_c <- read.csv(paste(input_files, "spec_02__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist3_c <- read.csv(paste(input_files, "spec_03__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist4_c <- read.csv(paste(input_files, "spec_04__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist5_c <- read.csv(paste(input_files, "spec_05__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist6_c <- read.csv(paste(input_files, "spec_06__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist7_c <- read.csv(paste(input_files, "spec_07__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist8_c <- read.csv(paste(input_files, "spec_08__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist9_c <- read.csv(paste(input_files, "spec_09__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist10_c <- read.csv(paste(input_files, "spec_10__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist11_c <- read.csv(paste(input_files, "spec_11__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist12_c <- read.csv(paste(input_files, "spec_12__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist13_c <- read.csv(paste(input_files, "spec_13__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist14_c <- read.csv(paste(input_files, "spec_14__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist15_c <- read.csv(paste(input_files, "spec_15__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist16_c <- read.csv(paste(input_files, "spec_16__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist17_c <- read.csv(paste(input_files, "spec_17__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist18_c <- read.csv(paste(input_files, "spec_18__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist19_c <- read.csv(paste(input_files, "spec_19__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist20_c <- read.csv(paste(input_files, "spec_20__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist21_c <- read.csv(paste(input_files, "spec_21__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist22_c <- read.csv(paste(input_files, "spec_22__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)
detecthist23_c <- read.csv(paste(input_files, "spec_23__detection_history__with_effort__10_days__2019-01-15.csv", sep = "/"), row.names=1)

# combine leech and camera trap detection histories
detecthist1 <- cbind(detecthist1_c,detecthist1_l)
detecthist2 <- cbind(detecthist2_c,detecthist2_l)
detecthist3 <- cbind(detecthist3_c,detecthist3_l)
detecthist4 <- cbind(detecthist4_c,detecthist4_l)
detecthist5 <- cbind(detecthist5_c,detecthist5_l)
detecthist6 <- cbind(detecthist6_c,detecthist6_l)
detecthist7 <- cbind(detecthist7_c,detecthist7_l)
detecthist8 <- cbind(detecthist8_c,detecthist8_l)
detecthist9 <- cbind(detecthist9_c,detecthist9_l)
detecthist10 <- cbind(detecthist10_c,detecthist10_l)
detecthist11 <- cbind(detecthist11_c,detecthist11_l)
detecthist12 <- cbind(detecthist12_c,detecthist12_l)
detecthist13 <- cbind(detecthist13_c,detecthist13_l)
detecthist14 <- cbind(detecthist14_c,detecthist14_l)
detecthist15 <- cbind(detecthist15_c,detecthist15_l)
detecthist16 <- cbind(detecthist16_c,detecthist16_l)
detecthist17 <- cbind(detecthist17_c,detecthist17_l)
detecthist18 <- cbind(detecthist18_c,detecthist18_l)
detecthist19 <- cbind(detecthist19_c,detecthist19_l)
detecthist20 <- cbind(detecthist20_c,detecthist20_l)
detecthist21 <- cbind(detecthist21_c,detecthist21_l)
detecthist22 <- cbind(detecthist22_c,detecthist22_l)
detecthist23 <- cbind(detecthist23_c,detecthist23_l)

# combine all detection histories into 3 dimensional array and
# rearrange it to correct dimensions
site.det <- abind(detecthist1, detecthist2, detecthist3, detecthist4, 
                  detecthist5, detecthist6, detecthist7, detecthist8,
                  detecthist9, detecthist10, detecthist11, 
                  detecthist12, detecthist13, detecthist14, detecthist15, 
                  detecthist16, detecthist17, detecthist18, detecthist19,
                  detecthist20, detecthist21, detecthist22,detecthist23, along = 3 )

det <- aperm(site.det, c(3, 1, 2))

# create the "method" covariate
# 1 = CT
# 2 = brown
# 3 = tiger
method <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 2, 3)

# load the effort table for camera traps and leeches
# create the expanded effort tables
# scale the effort 
CT_effort <- read.csv(paste(input_files, "effort_matrix_ct.csv", sep = "/"), row.names = 1)
CT_effort <- as.matrix(CT_effort)
CT_effort[is.na(CT_effort)] <- 0

leech_effort <- read.csv(paste(input_files, "effort_matrix_leech.csv", sep = "/"), row.names=1)
leech_effort <- as.matrix(leech_effort)
leech_effort[is.na(leech_effort)] <- 0
# number of stations, number of occasions take dim from det
leech_effort_expand <- matrix(0, nrow = dim(det)[2], ncol = dim(det)[3])
ct_effort_expand <- matrix(0, nrow = dim(det)[2], ncol = dim(det)[3])

leech_effort_expand[, (ncol(CT_effort) + 1) : (ncol(CT_effort) + ncol(leech_effort))] <- leech_effort
ct_effort_expand[, 1 : ncol(CT_effort)] <- CT_effort

eff <- leech_effort_expand
eff[,  1 : ncol(CT_effort)] <- ct_effort_expand[, 1 : ncol(CT_effort)]

leech_efftot <- scale(as.vector(leech_effort_expand))
leech_efftot <- matrix(leech_efftot, dim(leech_effort_expand))
CT_efftot <- scale(as.vector(ct_effort_expand))
CT_efftot <- matrix(CT_efftot, dim(ct_effort_expand))

# assign covariates from covariate table
elevation = covtable$gps_elevation
village_density = covtable$village_density
city_distance_linear = covtable$city_distance_linear
roads_average_LCP = covtable$roads_average_LCP_3_years

# list of stations with covariate values
stationlist = covtable$station

# change NA to zero
eff[is.na(eff)] <- 0
det[is.na(det)] <- 0

# binary info of trapping effort
eff.abs = eff
eff.abs[eff.abs > 1] <- 1

# number of camera stations
num.stations <- nrow(covtable)

# input data matrix for the occupancy model and delete occasion names for efftot and effabs ####
occ.in <- list(y= det,
               method = method,
               leech_efftot = leech_efftot,
               CT_efftot = CT_efftot,
               efftot=eff,
               effabs=eff.abs,
               J=num.stations,
               M=23,  # number of species           
               ELE = as.vector(scale(elevation)),
               elevation_unscaled = as.vector(elevation),
               VD = as.vector(scale(village_density)),
               village_density_unscaled = as.vector(village_density),
               CDL = as.vector(scale(log(city_distance_linear))),
               city_distance_log_unscaled = as.vector(log(city_distance_linear)),
               LCP = as.vector(scale(roads_average_LCP)),
               roads_average_LCP_unscaled = as.vector(roads_average_LCP),
               maxocc=length(method))  

# remove attributes that aren't needed
attr(occ.in$y, "dimnames") <- NULL

# save the occupancy input file so that we don't have to run
# everything above here again unless we change the input
dput(occ.in, 'OccData_mammal_community_RERUN.R')

# load the occupancy input file
occ.in <- dget('OccData_mammal_community_RERUN.R')

# set which parameters from the model we want to save
params<-c("R3", "new.R3",
          "mu.a", "sigma.a",
          "mu.b1", "sigma.b1",
          "mu.b2", "sigma.b2",
          "mu.b3", "sigma.b3",
          "mu.b4", "sigma.b4",
          "mu.bp", "sigma.bp",
          "alpha", "b.eff",
          "beta1", "beta2", 
          "beta3", "beta4", 
          "mu.p", "b.leech",
          "b.CT", "res",
          "R1")
        
# load JAGS model file 
modelFile2 <- paste("com_occ_model_FULL.txt", sep = "")

# create initial values
zin <- apply(occ.in$y, 1 : 2, sum)
zin[zin > 1] <- 1

inits2 <- function(){
  list(mu.a = runif(1, -4.5, -1.5), 
       sig.a = runif(1, 0.0, 2.5),
       mu.bp = runif(3, 0, 2),
       tau.bp = runif(3, 0.5, 1), 
       mu.b1 = runif(1),
       tau.b1 = runif(1, 0.5, 1),
       mu.b2 = runif(1),
       tau.b2 = runif(1, 0.5, 1),
       mu.b3 = runif(1),
       tau.b3 = runif(1, 0.5, 1),
       mu.b4 = runif(1),
       tau.b4 = runif(1, 0.5, 1),
       b.eff = runif(1), 
       z = zin,
       b.leech = runif(1, 0, 0.5),
       b.CT = runif(1)
       )
  }

# set model parameters
n.iter <- 250000
thin <- 1
n.burnin <- 50000

# create wrapper function for model to run chains in parallel
wrapper2 <- function(a){
  mod <- jags.model(modelFile2, occ.in, inits2, n.chain = 1, n.adapt = 500)
  out <- coda.samples(mod, params, n.iter, thin)
  return(out)
  }

# run model in parallel
# setting for parallel computing with the snowflake package; 
# cpus = 3 defines the number of chains started
sfInit( parallel = TRUE, cpus = 3)
sfLibrary(rjags)
sfExportAll()
sfClusterSetupRNG()
time.2 <- system.time(
  (out.sp = sfLapply(1 : 3, wrapper2))
)
sfStop()

# save model output
saveRDS(out.sp, "model_rerun_NOTHIN_RESIDUALS.rds")


######Residual Analysis#######
out.f <- mcmc.list(mcmc(out.sp[[1]][[1]]), mcmc(out.sp[[2]][[1]]), mcmc(out.sp[[3]][[1]]))

saveRDS(out.f, "model_rerun_NOTHIN_outf.rds")

# generate summary statistics of model output
sm <- summary(window(out.f, start = 50000))

saveRDS(sm, "model_rerun_NOTHIN_sm.rds")

# test for model convergence
gd <- gelman.diag(window(out.f, start = 50000), multivariate = FALSE)
gd[["psrf"]]

# get BP value
outlist <- mcmc(cbind(out.f[[1]][, "R3"], out.f[[2]][, "R3"], out.f[[3]][, "R3"]))
outlist2 <- mcmc(c(out.f[[1]][, "new.R3"], out.f[[2]][, "new.R3"], out.f[[3]][, "new.R3"]))

BP_value <- sum(outlist > outlist2) / length(outlist) ##0.60
BP_sd <- sd(outlist > outlist2) / length(outlist) ##1.36e-5

# create and save table with summary statistics, gd, and BP values
out <- cbind(sm$statistics, sm$quantiles, gd[["psrf"]], BP_value, BP_sd)
write.csv(out, paste("model_rerun_NOTHIN_NOBLOCK.csv", sep = ""))