# Runs a cohort of animals born at a month given by birth month

getCohortMtx <- function(initPop,          # initial population  (float)
                         birthMth = 1,     # month of birth of the cohort (mth: integer 1-12)
                         allInp)
  {

  Category   <- allInp$Category
  wAge       <- allInp$wAge           # weaning age (mth: integer)
  Birth_LW   <- allInp$Birth          # birth liveweight (kg: float)  
  pregLen    <- allInp$pregDuration   # duration of pregnancy (mths: integer)
  HfCalvIntv <- allInp$HfCalvIntv     # interval between first and second calving (mth: integer)  
  maxSize    <- allInp$maxSize        # maximum size of the age vector
  
  if(allInp$BullsPerCow == 0){
    Category$Bull <- NULL
  }
  
  catnames <- names(Category)
  
  # Generates a dataframe with basic cohort data  
  coh <- getCohort(maxSize, birthMth, growCatNames = c("SucklingCalves", "WeanedCalves", "Heifer"), allInp = allInp, catnames = catnames)
  
  # Get the daily liveweight gain (kg/d) for each month from cohort birth
  DLWG <- as.numeric(unlist(sapply(coh$diet, function(dietname) allInp$Diet[[dietname]]["LWG.Ref"])))
  
  # Calculate cumulative liveweight gain (kg) in order to determine 
  # the ages at conception, first and second calving
  tempcumLWG    <- cumsum(DLWG*allInp$meanDiM) 
  
  lAges         <- getAges(tempcumLWG, allInp$HfCalvIntv, allInp)
  lAges$weaning <- wAge
  
  if(lAges$calving1 < 22){
    stop("First calving age is less than 22 months. Please verify growing categories' weight gain.")
  }
  
  nTrim   <- lAges$calving1     # nrows of the dataframe until immediately before the second calving
  ageBool <- (coh$age < nTrim)  # boolean vector of the valid ages
  coh     <- coh[ageBool,]      # resize the cohort dataframe  (remove extra rows)
  DLWG    <- DLWG[ageBool]      # resize the DLWG vector
  
  # set indexes of the Category list for the categories of growing animals
  growingAnimals <- c("SucklingCalves", "WeanedCalves", "Heifer", "LactCow", "DryCow")
  CatIdx <- SetAgeAttr(coh$age, lAges, growingAnimals, allInp)
  
  # get digestibility (Dig) and crude protein (CP) of diets depending on age
  Dig <- as.numeric(unlist(sapply(coh$diet, function(dietname) allInp$Diet[[dietname]]["Dig"]))) 
  CP  <- as.numeric(unlist(sapply(coh$diet, function(dietname) allInp$Diet[[dietname]]["CP"]))) 
  
  # get physiological states (mth of pregnancy and mth of lactation)
  PhyState <- getPhysStates(pregAge = lAges$preg, dryLen = allInp$PreCalvingDryMths, HfpregAge =  lAges$preg,
                            fstCalvAge = lAges$calving1, HfCalvIntv, n = nTrim, allInp = allInp)

  # calculate LW from DLWG and birth weight
  LW <- allInp$Birth + c(0, cumsum(DLWG[-length(DLWG)]*allInp$meanDiM ))
  LW <- pmin(LW, allInp$ECRefLW)
  
  # calculate number of animals
  PopDyn <- getPopDyn(initPop, CatIdx, n = nTrim, allInp = allInp)
  
  coh <- cbind.data.frame(coh, 
                          CatStr   = names(Category[CatIdx]), 
                          CategSex = coh$sex, 
                          LW       = LW, 
                          LWG      = DLWG, 
                          PhyState[1:nTrim,], 
                          PopDyn, 
                          Dig = Dig, 
                          CP  = CP)
  return(coh)
}

CalcHerdStruct <- function(GrwFcohGEEData,
                           allInp){
  
  # taxa de aptidão das fêmeas para substituir as vacas
  femSurvRate  <- mean(sapply(GrwFcohGEEData, function(l) l[["Pop"]][length(l[["Pop"]])]))/100
  # média dos 3 primeiros meses, para todos os meses / 100
  suckSurvRate <- mean(sapply(GrwFcohGEEData, function(l) l[["Pop"]][1:allInp$wAge]))/100
  calvingHeifers <- allInp$CowsReplRate * allInp$BreedingCows              # Number of pregnant heifers needed per year
  weanedCalves   <- calvingHeifers / (allInp$femAptitudeRate * femSurvRate)       # Number of 1-year-old heifers needed to have the number of pregnant heifers
  sucklingCalves <- weanedCalves * (12 / allInp$CowsCalvIntv) * (allInp$wAge / 12) * suckSurvRate
  dryCows        <- allInp$BreedingCows * allInp$dryCowsRate  
  lactCows       <- allInp$BreedingCows * allInp$lactCowsRate  
  
  return(list(
    Struct = c(Bull           = 0, 
               DryCow         = dryCows,
               Heifer         = calvingHeifers,
               LactCow        = lactCows, 
               SucklingCalves = sucklingCalves,
               WeanedCalves   = weanedCalves),   
    Stats = c("NA" = NA,
              "NA" = NA,
              "NA" = NA,
              "NA" = NA,
              "NA" = NA,
              "NA" = NA,
              "NA" = NA,
              FPCM      = allInp$FPCM,
              milkProd  = allInp$milkProd,
              pastArea  = ifelse(is.null(allInp$pastArea), 0, allInp$pastArea),
              supplArea = ifelse(is.null(allInp$supplArea), 0, allInp$supplArea))
    ))
}