GetCatDiet <- function(aCateg,
                       allInp)
  {
  
  Diet <- allInp$Diet
  Category <- allInp$Category
  
  diets <- lapply(allInp$mthnames, function(m) Diet[[Category[[aCateg]][[m]]]])
  names(diets) <- allInp$mthnames
  
  return(diets)
}

#sets the age dynamics of a cohort in a monthly calendar basis
getCohort <- function(maxAge, 
                      birthMth, 
                      growCatNames = c("x"), 
                      allInp,
                      catnames)
  {
  
  Category <- allInp$Category
  n <- maxAge + 1                        # nrows of the matrix to allow for a given maximum age
  age <- seq(0, maxAge, 1)               # age of the category (months) 
  mth <- ((birthMth + age -1) %% 12) + 1 # month of the year
  mthdays <- allInp$DinMth[mth]
  yr <- ((birthMth + age -1) %/% 12) + 1 # age in years
  catName <- vector(mode = "character", length = length(age))
  
  growingAnimals <- c("SucklingCalves", "WeanedCalves", "Heifer", "Heifer", "Heifer")
  
  catName <- SetAgeAttr(age, list(weaning  = allInp$wAge, 
                                 calving1 = length(age), 
                                 calving2 = length(age)), 
                       growingAnimals, allInp, Num = T)
  
  catSex <- unlist(unname(sapply(catName, function(catnm) Category[[catnm]][["Sex"]])))
  
  diets <- unlist(unname(multiApply(FUN = function(catName, m) Category[[catName]][[m]], catName, allInp$mthnames[mth])))
  
  return(cbind.data.frame(age = age, mth = mth, mthdays = mthdays, yr = yr,
                          sex = catSex, diet = diets))
}

getAttr <- function(cohdata, 
                    diets,
                    wAge,
                    attribute)
  {
  
  out  <- cohdata$age
  nage <- which(out < wAge)
  out[out < wAge]  <- sapply(cohdata$mth[out < wAge], function(n) as.numeric(diets$SucklingDiet[[n]][[attribute]]))
  out[out >= wAge] <- sapply(cohdata$mth[out >= wAge], function(n) as.numeric(diets$GrowDiet[[n]][[attribute]]))
  
  return(out)
}

getDietAttrMths <- function(diets,
                            attribute,
                            mthnames)
{
  out <- sapply(allInp$mthnames, function(mn) as.numeric(diets[[mn]][attribute]))
  return(out)
}

SetAgeAttr <- function(age, aAges, attr, allInp, Num = T){
  
  v <- vector("character", length(age))
  v[age < aAges$weaning] <- attr[1]
  v[(age >= aAges$weaning) & (age <= 12)] <- attr[2]
  v[(age > 12) & (age < aAges$calving1)] <- attr[3]
  DryAge <- aAges$calving2 - allInp$PreCalvingDryMths
  v[(age >= aAges$calving1) & (age < DryAge)] <- attr[4]
  v[age >= DryAge] <- attr[5]
  
  return(v)
}

getAges <- function(cumLWG, HfCalvIntv, allInp) {
  calving1 <- which(cumLWG > allInp$Calving1 - allInp$Birth)[1]
  preg <- calving1 - allInp$pregLen
  calving2 <- calving1 + HfCalvIntv
  
  return(list(preg = preg, calving1 = calving1, calving2 = calving2))
}

getPhysStates <- function(pregAge,
                          dryLen,
                          HfpregAge,
                          fstCalvAge,
                          HfCalvIntv,
                          n,
                          allInp)
  {
  # month of pregnancy. NA is used for non-pregnant
  pregLen <- allInp$pregLen
  pMth <- seq(1,pregLen,1)
  
  pregMth <- c(rep(NA, pregAge), pMth, rep(0, HfCalvIntv - pregLen), pMth)
  
  # TODO: use % lactating cows 
  lactLen <- HfCalvIntv - dryLen 
  
  # month of lactation. NA is used for non-lactating
  lMth <- seq(1,lactLen,1) 
  if (n > fstCalvAge){
    lactMth <- c(rep(NA, fstCalvAge), lMth, rep(NA,n-fstCalvAge-lactLen))
  }else{
    lactMth <- NA
  }
  
  return(cbind.data.frame(pregMth = pregMth, lactMth = lactMth))
}  

getPopDyn <- function(initPop,
                      CatIdx,
                      n,
                      allInp)
  {
  Category <- allInp$Category
  
  # get survival rates  
  Morty <- sapply(CatIdx, function(idx) Category[[idx]][["Morty"]])
  Morty <- Morty/100
  Survv <- (1 - as.numeric(Morty)/12)
  # get population dynamics
  Pop <- numeric(n)
  Pop[1] <- initPop
  for (i in 1:(length(Survv)-1)){
    Pop[i+1] <- Pop[i] * Survv[i]
  }
    
  return(cbind.data.frame(Pop = Pop, Morty = Morty, Survv = Survv))
}

GrwFcohMth <- function(cohGEEData,
                       allInp){

  sumList <- list()
  
  varNames = c("FI", "CH4", "N2O")
  
  for(v in varNames){
    dfRes <- data.frame()
    
    for(mn in allInp$mthnames){
      DF <- cohGEEData[[mn]]
      dfList <- split(DF[, v], f = DF$CatStr)
      dfRes  <- rbind(dfRes, as.data.frame(lapply(dfList, mean)))
    }
    
    sumList[[v]] <- as.data.frame(sapply(colnames(dfRes), function(cn) mean(dfRes[,cn])), nm = v)
  }
  
  resDF <- do.call(cbind, sumList)
  
  return(resDF)
}