###########################################################################################
## Code for simulation in:
## "Dynamics of tree-cavity occupancy support the hypothesis that low use of 
##  avian-excavated cavities is driven by resource suitability"
## D.J. Zavala, K.L. Cockle, M.R. Gomez, C.A. Ferreyra, E.B. Bonaparte, 
## F.G. DiSallo, and G. Ferraz
## DOI: 10.5281/zenodo.14262784
###########################################################################################
## Simulates nest cavity use and estimates associated transition probabilities.
## Looks at two types of cavity (excavated vs non-excavated) and two types of user bird
## (excavating and secondary cavity-nesting)


# clean up
rm(list=ls())
# Load functions
source("SimCavsFunctions.R")

###########################################################################################
# 1. Define simulation parameters

nsim <- 3000             # number of simulations
ny <- 20                 # number of years
Nexc <- 200              # number of excavated cavities
Nnat <- round(1.6*Nexc)  # number of non-excavated cavities, in proportion to excavated
# vector of cavity types ("E" = excavated, "N" = non-excavated)
cvts <- c(rep("E",Nexc),rep("N",Nnat))

# define probabilities of use and cross-use (averaged from observations)
ne <- 0.01 # probability that natural cavity is taken by EXC (cross 1)
ns <- 0.46 # probability that natural cavity is taken by SCN
ee <- 0.48 # probability that excavated cavity is taken by EXC
es <- 0.13 # probability that excavated cavity is taken by SCN (cross 2)
useP <- c(ne,ns,ee,es) 

# number of pairs of each type of species
nEXC <- round(ee*Nexc+ne*Nnat)
nSCN <- round(ns*Nnat+es*Nexc)
birds <- c(nEXC,nSCN)

# setup probability of cavity loss
le0 <- 0.14 # excavated empty
le1 <- 0.09 # excavated occupied by EXC
le2 <- 0.08 # excavated occupied by SCN
ln0 <- 0.05 # non-excavated empty
ln1 <- 0.04 # non-excavated occupied by EXC
ln2 <- 0.02 # non-excavated occupied by SCN
elossvec <- c(le0,le1,le2) # later called eloss
nlossvec <- c(ln0,ln1,ln2) # later called nloss

###########################################################################################
# 2. Simulate nsim datasets and extract nsim transition matrices for both types of cavity
#    Load functions at end of script before running the commands bellow (see F1, F2, F3)

# setup empty arrays that will receive transition matrices
ETMarray <- NTMarray <- array(data=NA, dim=c(4,4,nsim))
freqvec_exc <- rep(NA,nsim)
freqvec_nat <- rep(NA,nsim)

set.seed(42)
for(k in 1:nsim) {
  #cat("\nsim",k,"\n")
  UseOut <- simCavUse(cvtp=cvts, birds=birds, ne=ne, ny=ny)
  CVU <- UseOut$STATES
  cvts <- UseOut$types
  CVUL<- simCavLoss(CVUse=CVU, cvtp=cvts, eloss=elossvec, nloss=nlossvec, ny=ny)
  #According to its origin (excavated 'E' , non-excavated 'N')
  ye <- CVU[which(cvts=="E"),]
  yn <- CVU[which(cvts=="N"),]
  ## get number of empty(1) per row
  yn0se <- apply(ye, 1, function(x) sum(x==0,na.rm=TRUE))
  yn0sn <- apply(yn, 1, function(x) sum(x==0,na.rm=TRUE))
  ## get number of empty(1) that occur before the last year of use
  yearly0e <- apply(ye, 1, function(x) sum(which(x==0) < max(which(x==1 | x==2)), na.rm=TRUE))
  yearly0n <- apply(yn, 1, function(x) sum(which(x==0) < max(which(x==1 | x==2)), na.rm=TRUE))
  ## get frequency of empty(1) that occur before last use
  yfearly0e <- sum(yearly0e)/sum(yn0se)
  yfearly0n <- sum(yearly0n)/sum(yn0sn)
  freqvec_exc[k]<-yfearly0e
  freqvec_nat[k]<-yfearly0n
  # save in Transition matrix
  TMs <- getTrans(CVUse=CVUL,cvtp=cvts)
  ETMarray[,,k] <- TMs[[1]]
  NTMarray[,,k] <- TMs[[2]]  
  
}

#Transition matrix
# to Excavated cavities
ExcT_mean <- apply(ETMarray,c(1,2),mean)
ExcT_sd <- apply(ETMarray,c(1,2),sd)
ExcT_025 <- apply(ETMarray,c(1,2), function(x) quantile(x,probs=0.025))
ExcT_975 <- apply(ETMarray,c(1,2), function(x) quantile(x,probs=0.975))
# to Non-excavated cavities
NatT_mean <- apply(NTMarray,c(1,2),mean)
NatT_sd <- apply(NTMarray,c(1,2),sd)
NatT_025 <- apply(NTMarray,c(1,2), function(x) quantile(x,probs=0.025))
NatT_975 <- apply(NTMarray,c(1,2), function(x) quantile(x,probs=0.975))

#Graphical representation "Simulated" results
plot(density(freqvec_exc), type='n', las = 1, yaxs="i", #xaxs="i",
     main = "", xlab = "",
     xlim = c(0,1),ylim = c(0,70))  
title(xlab = "Proportion of 'Empty' cavities before last occupation")
lines(density(freqvec_exc), col = "gray50", lwd = 2, lty = 5)
lines(density(freqvec_nat), col = "gray90", lwd = 2, lty = 5)