# 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, # dataframe with the output of the getCohortMtx call?
                           milkCowCategName,
                           femProp = 0.5,  # femProp: proportion of females in calves (useful if calving)
                           allInp)
  {

  Category    <- allInp$Category
  RR          <- allInp$CowsReplRate
  CowsCI      <- allInp$CowsCalvIntv
  HfCI        <- allInp$HfCalvIntv
  PreCvDP     <- allInp$PreCalvingDryMths
  Cows        <- allInp$BreedingCows
  CowsLMP     <- allInp$MilkProd
  BullsPerCow <- allInp$BullsPerCow
  CvDist      <- allInp$CvSsnalDist
  CowsMorty   <- as.numeric(Category[[milkCowCategName]][["Morty"]])
  CowsMorty   <- CowsMorty/100
  
  # Number of calves produced by the cohort per year and month of calving of the cohort
  cohCalvingData <- sapply(GrwFcohGEEData, function(DF) DF[nrow(DF),"Pop"]) 
  
  calvesPropMth <- unlist(lapply(GrwFcohGEEData, function(cohBirth) sum(cohBirth[cohBirth$CatStr == "SucklingCalves", "Pop"])/12))/100
  weanedPropMth <- unlist(lapply(GrwFcohGEEData, function(cohBirth) sum(cohBirth[cohBirth$CatStr == "WeanedCalves", "Pop"])/12))/100
  heiferPropMth <- unlist(lapply(GrwFcohGEEData, function(cohBirth) sum(cohBirth[cohBirth$CatStr == "Heifer", "Pop"])/12))/100
  
  PHCv             <- mean(cohCalvingData)/100                # PHCv: proportion of Heifers Calving
  RpH              <- RR * Cows / (1 - CowsMorty/2)           # number of replacement heifers needed 
  retFemCalves     <- RpH / PHCv                              # number of female calves needed
  cvnbMth          <- allInp$CvSsnalDist * retFemCalves
  LL               <- CowsCI - PreCvDP                        # lactation length (mth)
  LLHf             <- HfCI - PreCvDP                          # lactation length of heifers (mth)
  propLact         <- LL/CowsCI                               # proportion of lactating cows (dimensionless)
  propLactHf       <- LLHf/HfCI
  CowCR            <- 12/CowsCI                               # cows calving rate (calves/cow/yr)
  HfCR             <- 12/HfCI                                 # heifers calving rate (calves/hf/yr)
  LLDays           <- LL * allInp$meanDiM                     # lactation length (days)
  nCv              <- Cows * (1-RR)*CowCR + RpH               # annual calvings (calves/yr)  
  annualFemCalves  <- nCv * femProp                           # annual female calves (calves/yr)  
  annualMaleCalves <- nCv * (1- femProp)                      # annual male calves (calves/yr)
  Bulls            <- Cows * BullsPerCow                      # number of bulls
  avgDMPL          <- CowsLMP/LLDays                          # average daily milk production (kg/cow)
  CvPerMth         <- CvDist * nCv                            # number of calvings per month
  LactHeifers      <- RpH * propLactHf
  Cows2ndLacPlus   <- Cows * propLact - RpH  
  DryCows          <- Cows * (1-propLact) + RpH * (1-propLactHf)  # number of dry cows
  
  return(list(
    Struct = c(Bull           = Bulls, 
               DryCow         = DryCows,
               Heifer         = sum(heiferPropMth * cvnbMth),
               LactCow        = Cows2ndLacPlus + LactHeifers, 
               SucklingCalves = sum(calvesPropMth * cvnbMth),
               WeanedCalves   = sum(weanedPropMth * cvnbMth)),                         
    Stats = c(annualFemCalves       = annualFemCalves, 
              retFemCalves          = retFemCalves, 
              femCalvesSurplus      = annualFemCalves - retFemCalves, 
              annualMaleCalves      = annualMaleCalves,
              replaceCalvingHeifers = RpH, 
              cowsLactationLen      = LL, 
              heifersLactationLen   = LLHf,
              FPCM      = allInp$FPCM,
              milkProd  = allInp$milkProd,
              pastArea  = allInp$pastArea,
              supplArea = allInp$supplArea)))
}