source("NestleInputs.R")
source("NestleUtils.R")
source("EmissionsModels.R")
source("NestleHerdDynamics.R")
source("WasteManagement.R")
source("SoilCTransitionMatrix.R")
source("OutputFix.R")

runGEE <- function(allInp, logger){
  log4r::info(logger, "CALCULATING EMISSIONS...")
  
  resList <- list(ProductionUnit = allInp$ProductionUnit)
  
  # run growing animals cohort emissions in order to get herd structure proportions
  GrwFcohGEEData <- multiApply(FUN = calcCohortEmissions, birthMth = 1:12, initPop = 100,
                               SIMPLIFY = F, MoreArgs = list(allInp = allInp))
  
  if(!is.list(GrwFcohGEEData) && GrwFcohGEEData == 1){
    stop("Error identified: check log file.")
  }
  
  names(GrwFcohGEEData) <- allInp$mthnames
  
  log4r::info(logger, "Growing animals cohort OK.")
  
  aveLW <- numeric()
  for(nm in unique(GrwFcohGEEData[[1]]$CatStr)){
    aveLW[[nm]] <- sum(sapply(GrwFcohGEEData, function(x) sum(x[x[["CatStr"]] %in% nm, "LW"])/nrow(x[x[["CatStr"]] %in% nm,])))/length(GrwFcohGEEData)
  }
  
  allInp[["averageLW"]] <- c(aveLW, "LactCow" = allInp$ECRefLW, "DryCow" = allInp$ECRefLW, "Bull" = 1.2*allInp$ECRefLW)
  allInp[["averageLW"]] <- allInp[["averageLW"]][order(names(allInp[["averageLW"]]))]
  
  # then calc herd structure including mature animals
  HerdStruct <- CalcHerdStruct(GrwFcohGEEData   = GrwFcohGEEData,
                               milkCowCategName = "LactCow",
                               allInp           = allInp)
 
  dfHSt <- herdStatisticsOutPut(stats = HerdStruct$Stats)
  resList <- c(resList, HerdStatistics = list(dfHSt))
  
  allHerdStructure <- herdStructureOutPut(struct = round(HerdStruct$Struct, digits = 1), allInp = allInp)
  resList <- c(resList, HerdStructure = list(allHerdStructure))
  
  log4r::info(logger, "Equilibrium herd structure OK.")
  
  #### Emissions for equilibrium herd
  HerdGHGTotal <- getTotalHerdEmissions(CohFun     = calcCohortEmissions, 
                                        allInp     = allInp,
                                        HerdStruct = HerdStruct$Struct,
                                        HerdStats  = HerdStruct$Stats)
  
  emissionList <- calculateHerdEmissions(HerdGHGTotal     = HerdGHGTotal,
                                         allHerdStructure = allHerdStructure, 
                                         allInp           = allInp,
                                         multSimpl        = c(FI = 1e+3, CH4 = 1, N2O = 1))
 
  emissionDF <- emissionListOutPut(emissions = emissionList)
  resList <- c(resList, emissionDF)
  
  log4r::info(logger, "Equilibrium and current herd emissions OK.")
  
  wasteEmissions <- herdWasteManagement(allInp           = allInp, 
                                        allHerdStructure = allHerdStructure)
  
  resList <- c(resList, 
               list(EquilibriumWasteEmissions = wasteEmissions$Equilibrium),
               list(CurrentWasteEmissions     = wasteEmissions$Current))
  
  
  log4r::info(logger, "Waste management emissions OK.")
  
  soilEmissions <- soilTransition(allInp = allInp, logger = logger)
  
  soilEmissionsDF <- transitionMatrixOutPut(stocks = soilEmissions)
  resList <- c(resList, soilEmissionsDF)
  
  log4r::info(logger, "Soil transition OK.")
  
  emissionStats <- GHGSynth(allHerdStructure = allHerdStructure,
                            emissionList     = emissionList,
                            wasteManage      = wasteEmissions,
                            soilEmissions    = soilEmissionsDF$PlotsStatistics, 
                            allInp           = allInp)
  
  emissionStatsDF <- emissionsStatsOutPut(stats = emissionStats)
  resList <- c(resList, EmissionStatistics = list(emissionStatsDF))
  
  log4r::info(logger, "Emissions statistics OK.")
  
  return(resList)
}