# R code for Table 1 of Efford & Fletcher 2025

# repeat simulations of Howe et al. 2022, Scenario 2
# with both static and stochastic densities
# originally with 6 grids; repeated also with 18

#-------------------------------------------------------------------------------
# start by defining functions runone(), runall() and some summary functions

runone <- function (D, i, distribution = 'poisson') {
	ngrids <- length(D)
	pops <- sim.popn(D         = D, 
					 core      = trps, 
					 buffer    = popbuffer, 
					 nsessions = ngrids)
	ch <- sim.capthist(trps, 
					   popn       = pops, 
					   detectfn   = 'HN', 
					   detectpar  = list(g0 = g0, sigma = sigma), 
					   noccasions = noccasions, 
					   nsessions  = ngrids, 
					   renumber   = FALSE)
	ni <- sapply(ch, nrow)
	chT <- rbind(ch, verify = FALSE)       # pool grids  ('mash')
	attr(chT, 'n.mash') <- ni
	
	fit <- secr.fit(
		chT, 
		mask     = msk, 
		detectfn = 'HN', 
		trace    = FALSE,
		details  = list(distribution = distribution))
	
	pred <- predict(fit)   # automatically unmashes to per-grid density 
	chat <- function(ni, np) {
		n <- sum(ni)
		J <- length(ni)
		J/(n * (J-np)) * sum((ni - n/J)^2)
	}
	adjpred <- adjustVarD(pred, chat = chat(ni,1))
	cat('Completed fit ', i, ' at ', format(Sys.time(), "%H:%M:%S %d %b %Y"), '\n')
	list(
		ni        = ni, 
		varration = var(ni)/mean(ni), 
		pred      = pred, 
		adjpred   = adjpred
	)
}

runall <- function (ngrids = c(6,18), type = c('static','stochastic'), nrepl = 1000) {
	simone <- function (ngrids, type) {
		# run one scenario
		out <- list()
		for (i in 1:nrepl) {
			if (type == 'static')
				D <- rep(c(6,18,12), length.out = ngrids) * 1e-4
			else
				D <- sample (c(6,18,12), size = ngrids, replace = TRUE) * 1e-4
			out[[i]] <- runone (D, i)
		}
		out
	}
	# run each scenario
    scen <- expand.grid(ngrids, type)
    mapply(simone, ngrids = scen[,1], type = scen[,2], SIMPLIFY = FALSE)
}
#-------------------------------------------------------------------------------

# summary functions
RB <- function(x, type = 'pred', true = mean(c(6,18,12)*1e-4)) {
	(x[[type]]['D','estimate'] - true) / true
}
rse <- function(x, type = 'pred') {
	x[[type]]['D','SE.estimate'] / x[[type]]['D','estimate']
}
COV <- function(x, type = 'pred', true = mean(c(6,18,12)*1e-4)) {
	x[[type]]['D','lcl'] < true && x[[type]]['D','ucl'] > true
}
sumout <- function(out) {
	rb <- sapply(out, RB, 'pred')
	rse <- sapply(out, rse, 'pred')
	rseadj <- sapply(out, rse, 'adjpred')
	vr <- sapply(out, '[[', 'varration')
	se <- function(x) sd(x, na.rm = TRUE) / sqrt(sum(!is.na(x)))
	chat <- sapply(out, function(x) x$adjpred['D', 'c-hat'])
	data.frame (
		nvalid = sum(!is.na(rb)), 
		RB = mean(rb),
		seRB = se(rb),
		RSE = mean(rse),
		seRSE = se(rse),
		COV = mean(sapply(out, COV, 'pred')),
		varration = mean(vr),
		sevarration = se(vr),
		chat = mean(chat),
		sechat = se(chat),
		RSEadj = mean(rseadj),
		seRSEadj = se(rseadj),
		COVadj = mean(sapply(out, COV, 'adjpred'))
	)
}

bracket <- function (df, col = c(2,4,7,9,11), dec = 4) {
	# formats output dataframe for table, enclosing SE in parentheses
	df <- round(df,dec)
	for (j in col) {
		df[,j] <- paste0(df[,j], ' (', df[,j+1], ')')
	}
	df <- df[,-(col+1)]
	data.frame(IHPP = rep(c('static','stochastic'), each = 2), 
			   ngrids = rep(c(6,18), 2),
			   df)
}

#-------------------------------------------------------------------------------

library(secr)

g0         <- 0.3
sigma      <- 1500
noccasions <- 6
popbuffer  <- 15000
buffer     <- 8000
trps       <- make.grid(5,8, detector = 'proximity', spacing = 2000)
msk        <- make.mask(trps, buffer = buffer)

set.seed(123)
setNumThreads(24)   # adjust as appropriate

# sims is a list with four components, one for each scenario. 
# Each component is a list with the output of 1000 replicates
sims  <- runall(nrepl = 1000)

# restore from file
# sims <- readRDS('matrix/spatialrep2.RDS')

# summarise each scenario
tmp <- do.call(rbind, lapply(sims, sumout))

# place SE in parentheses
tab <- bracket(tmp)

# Latex code to generate Table 1 (some final editing needed)
knitr::kable(tab[,-3], format = 'latex', row.names = FALSE)

#-------------------------------------------------------------------------------