### Function to estimate annual fecundity by nesting attempts

setwd("d:/Amazon Drive/Documents/Academia/Collaboration/Andre/Data Andre/")

library(tcltk, pos=15)
library(aplpack, pos=15)
library(pacman)
pacman::p_load(Rcmdr, nlme, mgcv, chron, cowplot, gplots, Hmisc, doBy, tidyverse, MASS, devtools, reshape2, dplyr, tidyr, purrr) 

FEC <- function(DATA, XX=DATA[,1], name=DATA[,2], NPairs=DATA[,3], BSL=DATA[,4], CS=DATA[,5], I=DATA[,6],N=DATA[,7],PS=DATA[,8],PF=DATA[,9],DSR=DATA[,10],NSR=DATA[,11]) {
  SPDATA <- data.frame(XX, name, NPairs, BSL, CS, I, N, PS, PF, DSR, NSR)
  Fecundity <- function(XX, name, NPairs, BSL, CS, I, N, PS, PF, DSR, NSR) {
    Pair <- rep(NA,NPairs)
    Fecundity <- rep(NA,NPairs)
    output <- data.frame()
    
    for (r in 1:XX) {
      
      for (t in 1:NPairs){
        
        di <- 0 #day counter during incubation 
        dn <- 0 #day counter during nestling period
        j <- 0 #day counter for BSL
        att <- 0 #counter for number of nesting attempts
        
        repeat{
          x <- runif(1)
          s <- 0
          
          
          if(di <= I) {
            
            j <- j + 1
            
            if (x <= DSR){
              di <- di + 1
              i <- di + dn
              
            }
            if (x > DSR){         
              j <- j + PF - 1
              di <- di + 1
              i <- di + dn 
              Pair <- t
              Success <- s
              Nday <- i
              att <- att + 1
              Attempt <- att
              Fecundity <- CS * s
              result <- data.frame(r,name,Pair,BSL,CS,I,N,PS,PF,DSR,NSR,Attempt,Success,Nday,Fecundity,row.names=NULL)  
              output <- rbind(output,result)
              di <- 0
              dn <- 0
              i <- di+dn
              next        
            }
          }
          
          if(di > I) {
            
            j <- j + 1
            
            if (x <= NSR){
              dn <- dn + 1
              i <- di - 1 + dn
              if(dn >= N) {
                s <- 1
                Pair <- t
                Success <- s
                Nday <- i
                att <- att + 1
                Attempt <- att
                Fecundity <- CS * s
                result <- data.frame(r,name,Pair,BSL,CS,I,N,PS,PF,DSR,NSR,Attempt,Success,Nday,Fecundity,row.names=NULL)  
                output <- rbind(output,result)
                j <- j + PS - 1
                di <- 0
                dn <- 0
                i <- di + dn
              }
            }
            if (x > NSR){
              dn  <-  dn + 1
              i <- di - 1 + dn
              Pair <- t
              Success <- s
              Nday <- i
              att <- att + 1
              Attempt <- att
              Fecundity <- CS * s
              result <- data.frame(r,name,Pair,BSL,CS,I,N,PS,PF,DSR,NSR,Attempt,Success,Nday,Fecundity,row.names=NULL)  
              output <- rbind(output,result)
              j <- j + PF - 1
              di <- 0
              dn <- 0
              i <- di + dn
              next
            }
          }
          if (j >= BSL && i >= 3) {next}
          if (j >= BSL && i <= 2) {break}
          
        }
        #cat("\n", "pair  - ", t,"\n")
        
      }
      cat("\n", "iter  - ", r,"\n")
    }
    return(output)
  }
  final <- list()
fulltable <- data.frame()

  for (sp in 1:nrow(DATA)) {  
    final[[sp]] <- Fecundity(SPDATA[sp,1],SPDATA[sp,2],SPDATA[sp,3],SPDATA[sp,4],SPDATA[sp,5],SPDATA[sp,6],SPDATA[sp,7],SPDATA[sp,8],SPDATA[sp,9],SPDATA[sp,10],SPDATA[sp,11])
    fulltable <- rbind(fulltable,final[[sp]])
    cat("\n", "Model  - ", sp,"\n")
  }
  write.table(fulltable, file=paste(SPDATA[1,2],"thamout.txt", row.names = FALSE, col.names = TRUE, sep=","))
  return(fulltable)
}

#Species data

Asmall_wo2 <- data.frame(XX = 100, name = "Asmall_wo2", Npairs = 30, BSL = 39, CS = 2, I = 11, N = 12, PS = 27, PF = 35, DSR = 0.85, NSR = 0.85)
Bsmall_wo3 <- data.frame(XX = 100, name = "Bsmall_wo3", Npairs = 30, BSL = 39, CS = 3, I = 13, N = 12, PS = 27, PF = 35, DSR = 0.85, NSR = 0.85)
Csmall_w2 <- data.frame(XX = 100, name = "Csmall_w2", Npairs = 30, BSL = 39, CS = 2, I = 11, N = 12, PS = 27, PF = 22, DSR = 0.85, NSR = 0.85)
Dsmall_w3 <- data.frame(XX = 100, name = "Dsmall_w3", Npairs = 30, BSL = 39, CS = 3, I = 13, N = 12, PS = 27, PF = 22, DSR = 0.85, NSR = 0.85)
Emid_wo2 <- data.frame(XX = 100, name = "Emid_wo2", Npairs = 30, BSL = 54, CS = 2, I = 11, N = 12, PS = 27, PF = 36, DSR = 0.88, NSR = 0.88)
Fmid_wo3 <- data.frame(XX = 100, name = "Fmid_wo3", Npairs = 30, BSL = 54, CS = 3, I = 13, N = 12, PS = 27, PF = 36, DSR = 0.88, NSR = 0.88)
Gmid_w2 <- data.frame(XX = 100, name = "Gmid_w2", Npairs = 30, BSL = 54, CS = 2, I = 11, N = 12, PS = 13, PF = 14, DSR = 0.88, NSR = 0.88)
Hmid_w3 <- data.frame(XX = 100, name = "Hmid_w3", Npairs = 30, BSL = 54, CS = 3, I = 13, N = 12, PS = 13, PF = 14, DSR = 0.88, NSR = 0.88)
Ilarge_wo2 <- data.frame(XX = 100, name = "Ilarge_wo2", Npairs = 30, BSL = 69, CS = 2, I = 11, N = 12, PS = 27, PF = 53, DSR = 0.97, NSR = 0.97)
Jlarge_wo3 <- data.frame(XX = 100, name = "Jlarge_wo3", Npairs = 30, BSL = 69, CS = 3, I = 13, N = 12, PS = 27, PF = 53, DSR = 0.97, NSR = 0.97)
Klarge_w2 <- data.frame(XX = 100, name = "Klarge_w2", Npairs = 30, BSL = 69, CS = 2, I = 11, N = 12, PS = 29, PF = 26, DSR = 0.97, NSR = 0.97)
Llarge_w3 <- data.frame(XX = 100, name = "Llarge_w3", Npairs = 30, BSL = 69, CS = 3, I = 13, N = 12, PS = 29, PF = 26, DSR = 0.97, NSR = 0.97)

TABLE <- data.frame()
TABLE <- rbind(TABLE, Asmall_wo2, Bsmall_wo3, Csmall_w2, Dsmall_w3, Emid_wo2, Fmid_wo3, Gmid_w2, Hmid_w3, Ilarge_wo2, Jlarge_wo3, Klarge_w2, Llarge_w3)

TABLE2<- data.frame(TABLE[,1],TABLE[,2], TABLE[,3], TABLE[,4], TABLE[,5], TABLE[,6], TABLE[,7], TABLE[,8], TABLE[,9], TABLE[,10], TABLE[,11])

thcasim <- FEC(TABLE2)
names(thcasim)
Aggdata<- aggregate(Fecundity ~ name + r + CS + Pair, data=thcasim, FUN=sum)
aggdata2 <- aggregate(name + Fecundity, data = Aggdata, FUN = length)