##########################################################################################################################################
##########################################################################################################################################
##########################################################################################################################################
# ------ S1 Script to simulate pedigrees
##########################################################################################################################################
##########################################################################################################################################
##########################################################################################################################################

# The function is used to simulate populations. The output is a list, 
# each element of this list is a list of two elements: 
# (i) "pedigrees" which is a the population for a given generation and contains information about the individuals (age, sex, age-class, ID and parents' ID),
# and (ii) "adult_mortality" which is the density-depend adult mortality values


##########################################################################################################################################
#----# Variables
##########################################################################################################################################

# Npopstart = Starting pop size (used only to create the initial population)

# Npoplong = Long term pop size 

# sexratiostart = Starting sex-ratio (used only to create the initial population)

# sexratiolong = Long term sex-ratio

# agestart = Starting pop age ("random" (but only adults) or pre-defined number)

# sexmatumale = Males sexual maturity age

# sexmatufemale = FeMales sexual maturity age

# fecV = Average number of offsprings per female

# fecvarV = Variation of the average number of offsprings per female

# repromaleperc = Percentage of sexualy mature males mating

# reprofemaleperc = Percentage of sexualy mature females mating

# multimaleperc = Percentage of sexualy mature males mating with > 1 partner

# multifemaleperc = Percentage of sexualy mature females mating with > 1 partner

# femalepartner = Average partners number for the polyandrous females

# lifeexpV = Max life expectancy

# mortA/mortY/mortN = Age class mortality (Adult, Yearling or Newborn)

# thresh = Mortality threshold (to bound the density-dependent adult mortality rate, otherwise it can be too low or too high and the pop explodes or collapses

# ngen = Number of generation simulated


##########################################################################################################################################
#----# Function
##########################################################################################################################################

pedisim <- function(Npopstart, Npoplong, 
					sexratiostart, sexratiolong, 
					agestart, 
					sexmatumale, sexmatufemale, 
					fec, fecvar, 
					repromaleperc, reprofemaleperc, 
					multimaleperc, multifemaleperc, femalepartner, 
					lifeexp, 
					mortA, mortY, mortN, 
					ngen, thresh,
					silent) {

##########################################################################################################################################
#----# Initialization
##########################################################################################################################################

# Define the number of males and females as a function of the initial sex-ratio and the initial population size
sexIni <- c(rep("M", Npopstart*sexratiostart), rep("F", Npopstart*(1-sexratiostart)))

# Define the age of starting inidividuals (I added the function invisible to not print the output)
invisible(ifelse(agestart=="random", ageIni <- sample(3:lifeexp, Npopstart, replace = T), 
ifelse(agestart<=lifeexp, ageIni <- rep(agestart, Npopstart), print("!!!cannot set an initial age older than the life expectancy parameter!!!"))))

# Define the age class 
ageclassIni <- ageIni
ageclassIni[ageIni<=1] <- "N"
ageclassIni[ageIni>1 & ageclassIni<=2] <- "Y"
ageclassIni[ageIni>2] <- "A"

# Create the starting population, add mother and father columns needed in the next steps
popIni <- cbind.data.frame(generation=0, ID=1:Npopstart, sex=sexIni, age=ageIni, ageclass=ageclassIni, mother=NA, father=NA)

# Create a list to store each generation after breeding
newpopl <- list()
newpopl[[1]] <- popIni

# Create a list to store each generation after breeding AND mortality to be used to estimate pop size
truepopl <- list()

##########################################################################################################################################
#----# Long-term simulations
##########################################################################################################################################


#--------------------------------------------------------------------------------#
#----# Start the repeat loop
#--------------------------------------------------------------------------------#
mortal <- c()
fecund <- c()

repeat {

#--------------------------------------------------------------------------------#
#----# The new generation beggin with the last generation as potential parents
#--------------------------------------------------------------------------------#
newpop <- newpopl[[length(newpopl)]]

# Keep the ID max, we'll need it to define the ID of the new individuals
IDmax <- max(newpop$ID)

#--------------------------------------------------------------------------------#
#----# Age, adults, yearlings and newborns death
#--------------------------------------------------------------------------------#
# Randomly kill the individuals as a function of the class mortality rate (I assume individuals died before mating, could change that if needed)
Atodie <- newpop[newpop$ageclass=="A",]$ID 
Ytodie <- newpop[newpop$ageclass=="Y",]$ID
Ntodie <- newpop[newpop$ageclass=="N",]$ID

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Density-dependant adult mortality
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

#I use the first 3 generations as a burn-in 
if (length(truepopl)+1 > 3) { 

# Add a threshold to bound the mortality rate, otherwise it can be too low or too high and the pop explodes or collapses
mortA2 <- mortA*(nrow(truepopl[[length(truepopl)]])/Npoplong)^c(fec)
ifelse(mortA2>c(mortA+thresh), mortA2 <- c(mortA+thresh), ifelse(mortA2<c(mortA-thresh), mortA2 <-c(mortA-thresh), mortA2 <- mortA2))

} else {  mortA2 <- mortA }

# Store the mortality values of each generation
mortal <- c(mortal, mortA2)

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@


Atodie <- sample(Atodie, length(Atodie)*mortA2)
Ytodie <- sample(Ytodie, length(Ytodie)*mortY)
Ntodie <- sample(Ntodie, length(Ntodie)*mortN)

newpop <- newpop[!(newpop$ID %in% c(Atodie, Ytodie, Ntodie)),]

# Add 1 year to the individuals and change the age class accordingly, add 1 year to the generation
newpop$age <- newpop$age+1
try(newpop[newpop$age>1 & newpop$age<=2,]$ageclass <- "Y", silent=TRUE)
newpop[newpop$age>2,]$ageclass <- "A"
newpop$generation <- newpop$generation+1

# Filter too old individuals (I assume individuals died before mating)
newpop <- newpop[newpop$age<lifeexp,]

# Store the population
truepopl[[length(truepopl)+1]] <- newpop

# Display the generation and the number of individuals in the pop on the console during the loop
if (silent==FALSE) {
cat("\n", "\n", "####", "generation", max(newpop$generation), "####", sep=" " , "\n")
cat("--->", nrow(newpop), "individuals", sep=" ", "\n")
flush.console()
}

#--------------------------------------------------------------------------------#
#----# Reproduction
#--------------------------------------------------------------------------------#

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Split the newpop by mating system
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Keep only sexualy mature individuals
malerepro <- newpop[newpop$sex=="M" & newpop$age>=sexmatumale,]
femalerepro <- newpop[newpop$sex=="F" & newpop$age>=sexmatufemale,]

# Randomly choose sexualy mature individuals that are going to reproduce
malerepro <- malerepro[sample(1:nrow(malerepro), nrow(malerepro)*repromaleperc),]
femalerepro <- femalerepro[sample(1:nrow(femalerepro), nrow(femalerepro)*reprofemaleperc),]

# Identify males and females with multiple mates (I assume this characteristic randomly changes every year, 
# i.e., an individual can be polygynous one year and monogamous the next year, could change that if needed)
ifelse(multimaleperc>0, malemulti <- malerepro[sample(1:nrow(malerepro), nrow(malerepro)*multimaleperc),], malemulti <- c("NA"))
ifelse(multifemaleperc>0, femalemulti <- femalerepro[sample(1:nrow(femalerepro), nrow(femalerepro)*multifemaleperc),], femalemulti <- c("NA"))

# Identify monogamous males and females (I assume that all the mature individuals are breeding, could add a parameter to change that)
ifelse(multimaleperc==1, malemonog <- c(NA), ifelse(multimaleperc!=1 & multimaleperc>0, malemonog <- malerepro[!(malerepro$ID %in% malemulti$ID),], malemonog <- malerepro))
ifelse(multifemaleperc==1, femalemonog <- c(NA), ifelse(multifemaleperc!=1 & multifemaleperc>0, femalemonog <- femalerepro[!(femalerepro$ID %in% femalemulti$ID),], femalemonog <- femalerepro))

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Monogamous females
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

if (multifemaleperc!=1) {
# Define the number of offspring for each female (sometimes the number can be <0, so I replace negative values by 0)
femalemonogoffnb <- round(rnorm(nrow(femalemonog), mean=fec, sd=fecvar))
femalemonogoffnb[femalemonogoffnb<0] <- 0
femalemonogoffspring <- rep(femalemonog$ID, femalemonogoffnb)

# Create a new DF for the offspring
femalemonogoffspring <- cbind.data.frame(generation=unique(newpop$generation), ID=NA, sex=NA, age=0, ageclass="N", mother=femalemonogoffspring, father=NA)

# Randomly attribute the sex to the offsrping
for (i in levels(as.factor(femalemonogoffspring$mother))) {

	temp <- femalemonogoffspring[femalemonogoffspring$mother==i,]
	tempsex <- c(rep("M", nrow(temp)*sexratiolong), rep("F", nrow(temp)*(1-sexratiolong)))
	ifelse(length(tempsex)<nrow(temp), tempsex <- c(tempsex, sample(c("M", "F"),1)), tempsex <- c(tempsex))

	femalemonogoffspring[femalemonogoffspring$mother==i,]$sex <- sample(tempsex, length(tempsex))

}

# Define the potential male partners
ifelse(multimaleperc==1, availablemales <- c(malemulti$ID), 
ifelse(multimaleperc==0, availablemales <- c(malemonog$ID), availablemales <- c(malemulti$ID, malemonog$ID)))

# Create a table with the couples (random association) 
# If there are not enough available sexually mature males for all the monogamous females, 
# limit the number of breeding monogamous females to the number of available males
if (length(unique(femalemonogoffspring$mother))>length(availablemales)) { 
femalemonogoffspring <- femalemonogoffspring[femalemonogoffspring$mother %in% sample(unique(femalemonogoffspring$mother), length(availablemales)),]
}

femalemonogpartner <- cbind.data.frame(mother=unique(femalemonogoffspring$mother), father=sample(availablemales, length(unique(femalemonogoffspring$mother))))

# Attribute the father to the offsrping
for (i in levels(as.factor(femalemonogoffspring$mother))) {

	femalemonogoffspring[femalemonogoffspring$mother==i,]$father <- femalemonogpartner[femalemonogpartner$mother==i,]$father

}
}

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Polygamous females
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

if (multifemaleperc!=0) {

# Define the number of offspring for each female (sometimes the number can be <0, so I replace negative values by 0)
femalemultioffnb <- round(rnorm(nrow(femalemulti), mean=fec, sd=fecvar))
femalemultioffnb[femalemultioffnb<0] <- 0

# repeat the mothers ID for each child
femalepolyffspring <- rep(femalemulti$ID, femalemultioffnb)

# Define the number of partners for each female (Poisson distribution +1 to avoid having female with offsprings but no partner)
femalepolyffspringpartner <- rep(rpois(nrow(femalemulti), femalepartner)+1, femalemultioffnb)

# Create a new DF for the offspring (I've added a new column with the number of male mates, don't forget to delete it later)
femalepolyffspring <- cbind.data.frame(generation=unique(newpop$generation), ID=NA, sex=NA, age=0, ageclass="N", mother=femalepolyffspring, father=NA, npartner=femalepolyffspringpartner)

# Attribute the sex to the offsrping and check that the number of partners <= to the number of offsprings, othewise, replace the number of partners by the number of offsprings
for (i in levels(as.factor(femalepolyffspring$mother))) {

	temp <- femalepolyffspring[femalepolyffspring$mother==i,]
	tempsex <- c(rep("M", nrow(temp)*sexratiolong), rep("F", nrow(temp)*(1-sexratiolong)))
	ifelse(length(tempsex)<nrow(temp), tempsex <- c(tempsex, sample(c("M", "F"),1)), tempsex <- c(tempsex))

	femalepolyffspring[femalepolyffspring$mother==i,]$sex <- sample(tempsex, length(tempsex))

	if (unique(temp$npartner)>nrow(temp)) {femalepolyffspring[femalepolyffspring$mother==i,]$npartner <- nrow(temp)} 

} 

# Define the potential male partners at the beggining of the loop, 
#differentiate monogamous from polygamous males and remove monogamous males that already reproduced with a monogamous female (step above)
ifelse(multimaleperc!=1, availablemonomales <- malemonog$ID[!c(malemonog$ID %in% unique(femalemonogoffspring$father))], availablemonomales <- c("NA"))
ifelse(multimaleperc!=0, availablemultimales <- malemulti$ID[!c(malemulti$ID %in% unique(femalemonogoffspring$father))], availablemultimales <- c("NA"))
ifelse(multimaleperc==1, availableallmales <- c(availablemultimales),
ifelse(multimaleperc==0, availableallmales <- c(availablemonomales), availableallmales <- c(availablemultimales, availablemonomales)))

# Create a table with the couples 
for (i in levels(as.factor(femalepolyffspring$mother))) {

	temp <- femalepolyffspring[femalepolyffspring$mother==i,]

# Randomly select the males for a given female and attribute them an equal number of her offsprings 
# I've added a "if else" because sometimes, there is no more males available, so I just remove the offspring
# I've added an "ifelse" cause sometimes there is not enought different males, so I just make all the males breed  		
if (length(availableallmales)!=0) {
 	ifelse(unique(temp$npartner)>length(availableallmales), 
	tempmales <- availableallmales,
	tempmales <- rep(sample(availableallmales, unique(temp$npartner)), nrow(temp)/unique(temp$npartner)))

# As it might not be possible to attribute exactly the same number of offsprings to each males, 
# add a step to randomly attribute the remaining offsprings to some of the males already selected    	
		while (length(tempmales)<nrow(temp))
		{
		tempmales <- c(tempmales, sample(tempmales, 1))
		}

	femalepolyffspring[femalepolyffspring$mother==i,]$father <- tempmales

# Remove the monogamous males if they already reproduced
ifelse(multimaleperc!=1, availablemonomales <- availablemonomales[!c(availablemonomales %in% unique(femalepolyffspring$father))], availablemonomales <- c())
ifelse(multimaleperc==1, availableallmales <- c(availablemultimales),
ifelse(multimaleperc==0, availableallmales <- c(availablemonomales), availableallmales <- c(availablemultimales, availablemonomales)))

} else {  femalepolyffspring <- femalepolyffspring[!c(femalepolyffspring$mother %in% c(i)),]

}

}

}

#--------------------------------------------------------------------------------#
# Merge mono- and polygamous females
#--------------------------------------------------------------------------------#

# Merge the old generations, the offsprings of monogamous females and the offsrping of polygamous males (don't forget to remove the last column with the male partners number)
ifelse(multifemaleperc!=0, 
ifelse(multifemaleperc!=1, newoffspring <- rbind.data.frame(femalemonogoffspring, femalepolyffspring[,-8]), 
newoffspring <- rbind.data.frame(femalepolyffspring[,-8])), newoffspring <- rbind.data.frame(femalemonogoffspring))

newoffspring$ID <- c(IDmax+1):c(IDmax+nrow(newoffspring))

newpop <- rbind.data.frame(newpop, newoffspring)

# Concatenate the new generation in the list
newpopl[[length(newpopl)+1]] <- newpop

	
#--------------------------------------------------------------------------------#
# Stop the loop when enough generations
#--------------------------------------------------------------------------------#
if(unique(newoffspring$gen) == ngen) {
      break
   }
}

# Output of the function is a list with three elements: the pedigrees, the mortality and the fecundity values for all the generations
return(list(pedigrees=truepopl, adult_mortality=mortal))

}