#
#  R Script to fit Occupancy model to Willow Tit Data
#  
#  Model Selection Approach: Cross-Validation and CPO
#


###
###  Load Covariates and Occupancy Data
###

wt.df=read.csv("wt.csv",header=TRUE)
head(wt.df)

n=200
Y=as.matrix(wt.df[1:n,1:3])
y=apply(wt.df[1:n,1:3],1,sum,na.rm=TRUE)
J=apply(!is.na(wt.df[1:n,1:3]),1,sum)

X=matrix(1,n,2)
X[,1]=scale(wt.df$elev[1:n])
X[,2]=scale(wt.df$forest[1:n])

cor(X)
pairs(X)

###
###  Fit Each Probit Occupancy Model using Parallelized C-V 
###
###  Issue this command in shell before starting R: export OMP_NUM_THREADS=1 
###

library(parallel) 
library(snowfall) 
library(rlecuyer) 
source("occ.aux.mcmc.R")

n.mcmc=160000
s2.beta=1.5^2

K=10
fold.idx.mat=matrix(1:n,,K)

cps=detectCores()
sfInit(parallel=TRUE, cpus=cps)
sfExportAll()
sfClusterSetupRNG()

cv.fcn <- function(k){
  score.vec=rep(0,4)
  fold.idx=fold.idx.mat[,k]
  y.oos=apply(Y[fold.idx,],1,sum,na.rm=TRUE)
  tmp.out.1=occ.aux.mcmc(Y[-fold.idx,],J[-fold.idx],X[-fold.idx,],y.oos,J[fold.idx],X[fold.idx,],n.mcmc,s2.beta,null.model=TRUE,no.print=TRUE)    
  score.vec[1]=tmp.out.1$score
  tmp.out.2=occ.aux.mcmc(Y[-fold.idx,],J[-fold.idx],X[-fold.idx,1],y.oos,J[fold.idx],X[fold.idx,1],n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE)    
  score.vec[2]=tmp.out.2$score
  tmp.out.3=occ.aux.mcmc(Y[-fold.idx,],J[-fold.idx],X[-fold.idx,2],y.oos,J[fold.idx],X[fold.idx,2],n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE)    
  score.vec[3]=tmp.out.3$score
  tmp.out.4=occ.aux.mcmc(Y[-fold.idx,],J[-fold.idx],X[-fold.idx,1:2],y.oos,J[fold.idx],X[fold.idx,1:2],n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE)    
  score.vec[4]=tmp.out.4$score
  score.vec
}

sfExport("Y","J","X","s2.beta","n.mcmc","cv.fcn","fold.idx.mat")
tmp.time=Sys.time()
score.list=sfClusterApplyLB(1:K,cv.fcn)
time.1=Sys.time()-tmp.time
sfStop()

time.1
score.cv.mat=matrix(unlist(score.list),,4,byrow=TRUE)
apply(score.cv.mat,2,sum)

###
###  Fit Each Probit Occupancy Model using Parallelized C-V 
###

tmp.time=proc.time()
cpo.mat=matrix(0,n,4)
cpo.mat[,1]=occ.aux.mcmc(Y,J,X,1,1,1,n.mcmc,s2.beta,null.model=TRUE,no.print=TRUE,no.pred=TRUE)$CPO.vec   
cpo.mat[,2]=occ.aux.mcmc(Y,J,X[,1],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,no.pred=TRUE)$CPO.vec   
cpo.mat[,3]=occ.aux.mcmc(Y,J,X[,2],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,no.pred=TRUE)$CPO.vec   
cpo.mat[,4]=occ.aux.mcmc(Y,J,X[,1:2],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,no.pred=TRUE)$CPO.vec   
time.2=proc.time()-tmp.time
time.2
-apply(log(cpo.mat),2,sum)

tmp.time=proc.time()
  occ.aux.mcmc(Y,J,X[,1:2],1,1,1,n.mcmc,s2.beta,null.model=FALSE,no.print=TRUE,no.pred=TRUE)$CPO.vec   
proc.time()-tmp.time