##########################################################################################################################################
##########################################################################################################################################
##########################################################################################################################################
# ------ S2 Script to run the four estimators
##########################################################################################################################################
##########################################################################################################################################
##########################################################################################################################################

# First run the simulation using your parameters
sim_res <- pedisim(Npopstart, Npoplong, 
					sexratiostart, sexratiolong, 
					agestart, 
					sexmatumale, sexmatufemale, 
					fec, fecvar, 
					repromaleperc, reprofemaleperc, 
					multimaleperc, multifemaleperc, femalepartner, 
					lifeexp, 
					mortA, mortY, mortN, 
					ngen, thresh,
					silent)


# Then the four estimators are using two objects:
# "true_pop" which is the full population for a given generation (here the last one)
# "samp_pop" which is a subset of this population

# Keep the last generation of the simulation
true_popV <- sim_res$pedigrees[[length(sim_res$pedigrees)]]

# Sample a given percentage of this population
percsamp <- 30
popsampleV <- true_popV[sample(1:nrow(true_popV), percsamp*nrow(true_popV), replace=F),]



############################################################
# Close-Kin Mark-Recapture Bravington et al. 2016 
# CKMR estimator
############################################################

CKMR_O <- function(samp_pop, true_pop) {

#--------------------------------------------------------------------------------#
#----# Number of adults sampled 
na <- nrow(subset(samp_pop, samp_pop$ageclass %in% c("A", "Y")))

#----# Number of juveniles sampled 
nj <- nrow(subset(samp_pop, samp_pop$ageclass=="N"))

#----# Number of parent-offspring pairs 
H <- sum(subset(samp_pop, samp_pop$ageclass=="N")$father %in% subset(samp_pop, samp_pop$sex=="M" & samp_pop$ageclass %in% c("A","Y"))$ID) +
sum(subset(samp_pop, samp_pop$ageclass %in% c("A","Y"))$father %in% subset(samp_pop, samp_pop$sex=="M" & samp_pop$ageclass %in% c("A","Y"))$ID) +
sum(subset(samp_pop, samp_pop$ageclass=="N")$mother %in% subset(samp_pop, samp_pop$sex=="F" & samp_pop$ageclass %in% c("A","Y"))$ID) +
sum(subset(samp_pop, samp_pop$ageclass %in% c("A","Y"))$mother %in% subset(samp_pop, samp_pop$sex=="F" & samp_pop$ageclass %in% c("A","Y"))$ID)

if (H==0) {
return(c(NA, NA)) 
} else {

#--------------------------------------------------------------------------------#
#----# Compute the estimator
N_A <- 2*nj*na/H

#----# Put the estimation and the true pop size 
CKMR <- data.frame(N_estim=N_A, N_true_pop = nrow(true_pop[true_pop$ageclass %in%c("A", "Y"),]))
return(CKMR)

}
}



############################################################
# Using pedigree reconstruction to estimate population size: genotypes are more than individually unique marks - Creel and Rosenblat 2013
# CRE estimator
############################################################

CRE_O <- function(samp_pop, true_pop) {

#--------------------------------------------------------------------------------#
#----# Ns Number of individuals sampled
Ns <- nrow(samp_pop)

#--------------------------------------------------------------------------------#
#----# Number of breeders sampled 
BsID <- c()

for (a in unique(samp_pop$age)) {

# I subset each generation
temp <- subset(samp_pop, samp_pop$age==a)

# I search for father or mother in the older generations, if a parent is parent of several offspring of the same age, I count it only once
BsID <- c(BsID, c(unique(temp$father[temp$father %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID]),
unique(temp$mother[temp$mother %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID])))

}

# I count the number of breeders sampled
Bs <- length(BsID)

if (Bs==0) {
return(c(NA, NA)) 
} else {

#--------------------------------------------------------------------------------#
#----# Nin Number of individuals not directly sampled, but whose presence inferred by pedigree reconstruction (i.e., individuals not sampled but for which we have identified an offspring and a parent)

Nintemp <- c()

for (a in unique(samp_pop$age)) {

# I subset each generation
temp <- subset(samp_pop, samp_pop$age==a)

# I first identify the fathers
fatherid <- temp$father[temp$father %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID]
# Then I subset the generation to keep only the offspring for which I've identified a father
fathersub <- temp[temp$father %in% fatherid,]
# Then I check the number of offspring for which I've also identified the mother
fathersubmother <- fathersub$mother[fathersub$mother %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID]
# Then the Nin is the number of individuals for which I've identified a father minus the number of individuals for which I've also identified a mother - I count as many inferred mother as offspring
Ninfather <- c(nrow(fathersub)-length(fathersubmother))

# I do the same thing for the mothers this time
motherid <- temp$mother[temp$mother %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID]
mothersub <- temp[temp$mother %in% motherid,]
mothersubfather <- mothersub$father[mothersub$father %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID]
Ninmother <- c(nrow(mothersub)-length(mothersubfather))

# I store the Nin
Nintemp <- c(Nintemp, Ninfather, Ninmother)

}

# Sum all the Nin for all the generation
Nin <- sum(Nintemp)

#--------------------------------------------------------------------------------#
#----# Compute the CRE
N_total <- Ns + 2*Nin - c(c(Nin*Bs)/c(Ns+Nin))

#----# Put the estimation and the true pop size 
CRE <- data.frame(N_estim=N_total, N_true_pop = nrow(true_pop[true_pop$ageclass %in%c("A", "Y"),]))
return(CRE)

}
}



############################################################
# Population size estimates based on the frequency of genetically assigned parent–offspring pairs within a subsample - Muller et al. 2020
# gCMR estimator
############################################################

gCMR_O <- function(samp_pop, true_pop) {

#----# Adult males number
nAm <- nrow(subset(samp_pop, samp_pop$sex=="M" & samp_pop$ageclass %in% c("A","Y")))

#----# Adult females number
nAf <- nrow(subset(samp_pop, samp_pop$sex=="F" & samp_pop$ageclass %in% c("A","Y")))

#----# Juveniles number
nJ <- nrow(subset(samp_pop, samp_pop$ageclass=="N"))

#----# Father-Offspring pairs number (2 blocs, the F-O pairs for the newborns, and the F-O for the adults and subadults)
nFOP <- sum(subset(samp_pop, samp_pop$ageclass=="N")$father %in% subset(samp_pop, samp_pop$sex=="M" & samp_pop$ageclass %in% c("A","Y"))$ID) +
sum(subset(samp_pop, samp_pop$ageclass %in% c("A","Y"))$father %in% subset(samp_pop, samp_pop$sex=="M" & samp_pop$ageclass %in% c("A","Y"))$ID)

#----# Mother-Offspring pairs number (2 blocs, the M-O pairs for the newborns, and the M-O for the adults and subadults)
nMOP <- sum(subset(samp_pop, samp_pop$ageclass=="N")$mother %in% subset(samp_pop, samp_pop$sex=="F" & samp_pop$ageclass %in% c("A","Y"))$ID) +
sum(subset(samp_pop, samp_pop$ageclass %in% c("A","Y"))$mother %in% subset(samp_pop, samp_pop$sex=="F" & samp_pop$ageclass %in% c("A","Y"))$ID)

if (nFOP==0 & nMOP==0) {
return(rbind(NA, NA, NA, NA)) 
} else {

#----# Muller et al. estimator
N_males <- nAm*nJ/nFOP
N_females <- nAf*nJ/nMOP
N_adults <- N_males + N_females
N_juveniles <- N_adults * ((nJ)/(nAm+nAf))
# N_juveniles <- N_females*((nJ)/(nAf))
# N_juveniles <- N_males*((nJ)/(nAm))
N_total <- nJ * ((nAf/nMOP)+(nAm/nFOP)) * (1+(nJ/(nAm+nAf)))

# Put everything in a DF 
estim <- data.frame(estimates = c(round(N_males), round(N_females), round(N_adults), round(N_juveniles), round(N_total)))
rownames(estim) <- c("N_males","N_females","N_adults","N_juveniles","N_total")

#----# True pop numbers
N_malesT <- nrow(subset(true_pop, true_pop$sex=="M" & true_pop$ageclass %in% c("A","Y")))
N_femalesT <- nrow(subset(true_pop, true_pop$sex=="F" & true_pop$ageclass %in% c("A","Y")))
N_adultsT <- nrow(subset(true_pop, true_pop$ageclass %in% c("A","Y")))
N_juvenilesT <- nrow(subset(true_pop, true_pop$ageclass=="N"))
N_totalT <- nrow(true_pop)

# Put everything in a new DF 
estim2 <- data.frame(estim, True_pop = c(N_malesT, N_femalesT, N_adultsT, N_juvenilesT, N_totalT))
estim2

return(estim2)
}
}



############################################################
# Inference from single occasion capture experiments using genetic markers - Chathurika and Hettiarachchige 2017
# Moment estimator
############################################################

# The package gsl is needed for this estimator
# install.packages("gsl", repos="https://cloud.r-project.org") 
library("gsl")

moment_O <- function(samp_pop, true_pop) {

#--------------------------------------------------------------------------------#
#----# Number of sampled individuals
n <- nrow(samp_pop)

#--------------------------------------------------------------------------------#
#----# Number of identified mothers (I count the mother only once even if it had several offspring during different breeding seasons)
xID <- unique(samp_pop[samp_pop$ID %in% samp_pop$mother,]$ID)

# I count the number of mothers sampled
x <- length(xID)

if (x==0) {
return(c(NA, NA)) 
} else {

#--------------------------------------------------------------------------------#
#----# Number of identified mother-daughter pairs (= number of identified daughters)
dID <- c()

for (a in unique(samp_pop$age)) {

# I subset each generation
temp <- subset(samp_pop, samp_pop$age==a)

# I search for mother in the older generations, if a parent is parent of several offspring of the same age, I count it only once
dID <- c(dID, temp$mother[temp$mother %in% samp_pop[samp_pop$age %in% unique(samp_pop$age)[unique(samp_pop$age)>a],]$ID])

}

# I count the number of identified mother-daughter pairs
d <- length(dID)

if (d==0) {
return(c(NA, NA)) 
} else {


#--------------------------------------------------------------------------------#
#----# Compute the moment estimator (not mine, function provided in the sup mat of the paper)
if(d==x){mu.p.1 <- d/x+0.01} else {mu.p.1 <- d/x}
t.1 <- mu.p.1/exp(mu.p.1)
mu.p.hat <- mu.p.1+lambert_W0(-t.1) 
if(n*(1-exp(-mu.p.hat))>2*x){mu.hat <- n*(1-exp(-mu.p.hat))/x-1} else {mu.hat <- 1}     #estimated mu
p.hat <- mu.p.hat/mu.hat      		#p.hat
N.hat <- mu.hat*n*(n-1)/(d*(1+mu.hat)^2) 


list(N.hat=N.hat,mu.hat=mu.hat,p.hat=p.hat)

#----# Put the estimation and the true pop size 
moment <- data.frame(N_mother_estim=round(N.hat), N_mother_true_pop = nrow(true_pop[true_pop$sex=="F" & true_pop$ageclass!="N",]))
return(moment)

}
}
}