FarmDynR.R

# FarmDynR.R ----
#
# Copyright (c) 2022 Hugo Scherer - Wageningen Economic Research
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.


## Modes ----
#' Retrieve mode of vector
#'
#' This function returns the mode of a vector. If the vector contains a character or factor, the most common character/factor is returned. Numbers written as characters will be compatible with non-character numbers (i.e. doubles/numeric), but the function returns a character.
#'
#' @param x vector from which to retrieve the mode from.
#' @return same class as 'x'.
#' @examples
#' Modes(x = c(1, 1, 3, 0, 2, 4, 2, 1, 5, 2, 1))
#' Modes(x = c('a','b', 'c', 'a', 'c', 'a'))
#' Modes(x = c('a', 2, 'x', 7895, 1, '2', 't', 2, 1))
#' @seealso tabulate()
#' @export Modes
Modes <- function(x) { # Function found on StackOverflow made by Ken Williams and expanded by digEmAll
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  ux[tab == max(tab)] # Takes the highest incidence value
}


## gdxbinwider ----
#' Join BIN data together, make joined dataset wider, and group by a mapping
#'
#' @description
#' This function has been conceived with the Dutch FADN in mind, please use `fadn2fd()` for EU FADN data.
#' 
#' The `gdxbinwider()` function takes in a GDX file with BIN data as parameters p_farmData_NL and p_farmData2GUI, and a mapping as a set.
#' Then the data is widened, and the output is a tibble.
#'
#' @param filename Name of the GDX file with BIN data and mappings.
#' @param BINDir Directory where the FADN data is located.
#' @param gdxmap Name of the set in the GDX file that contains the mapping (e.g. Regs2BINID)
#' @param mapping Column name of the characteristic/variable to be grouped by (e.g. "Regions" or "Regs")
#' @return A tibble `tbl_df`
#' @seealso NULL
#' @export gdxbinwider
gdxbinwider <- function(filename, BINDir, gdxmap, mapping){

  if ('BINDir' %in% ls(envir = .GlobalEnv) & missing(BINDir)) { # Checks if BINDir is in Global Environment and uses it
    BINDir <- get('BINDir', envir = .GlobalEnv)
  } else {
    BINDir
  }

fd2guicolnames <- c("all_binid", "item1", "item2", "value")

fdnlcolnames <- c( "year", "all_binid","item1", "value")

fd2gui <- (gdxrrw::rgdx.param(file.path(BINDir, filename), "p_farmData2GUI", names=fd2guicolnames, compress=FALSE, ts=FALSE,
                                   squeeze=FALSE, useDomInfo = TRUE, check.names = TRUE))

fdnl <- (gdxrrw::rgdx.param(file.path(BINDir, filename), "p_farmData_NL", names=fdnlcolnames, compress=FALSE, ts=FALSE,
                                 squeeze=FALSE, useDomInfo = TRUE, check.names = TRUE))


gdxmapping <- (gdxrrw::rgdx.set(file.path(BINDir,filename), gdxmap, compress=FALSE, ts=FALSE,
                                     useDomInfo = TRUE, check.names = TRUE))

# Pick entries for non-numerical globals (1.6e+303 is code for numerical globals).

dummies <- fd2gui[which(fd2gui$all_binid=="dummy"),]

nonnumdummies <- dummies[dummies$value != 1.6e+303,]

# Separation of dummy from fd2gui

fd2gui <- fd2gui[which(fd2gui$all_binid!="dummy"),]

# Widening of the data for better and easier viewing and handling

fd2gui <- fd2gui %>%
  tidyr::pivot_wider(id_cols = "all_binid", names_from = c("item1", "item2"), values_from = "value", names_sep = "%")


fdnl <- fdnl %>%
  tidyr::pivot_wider(id_cols = c("all_binid", "year"), names_from = c("item1"), values_from = "value", names_sep = "%")

# Since subsetting with !duplicated() is done following the data order, reorder the data to start from youngest to oldest
# Result: all_binids are preserved for youngest years (2020, 2019...) and removed from oldest years (2013, 2014...)

fdnl <- fdnl[order(-as.numeric(fdnl$year)),]

fdnl <- fdnl[!duplicated(fdnl$all_binid),]

# Note: Some of the Weights appear in 2019 (n-1), but not 2020 (n = earliest data year), leading to NAs later on. Make sure at least the latest year has weights.

# Keep list of variables that are NOT numeric, BUT binary or other variables with other meanings.

tokeep <- c(paste("global", nonnumdummies$item2, sep = "%"), "global%soilTypeFirm", "global%derogatie")

# Finding index of global to keep only NON-numerical values based on "tokeep"

keepmatch <- match(tokeep, colnames(fd2gui))

# Making these factors
fd2gui[keepmatch] <- lapply(fd2gui[keepmatch], as.factor)


# Joining all data together

fdnl2gui <- dplyr::left_join(fd2gui, fdnl[,c('all_binid', 'Weight')], by='all_binid')

fdnl2gui[] <- lapply(fdnl2gui, function(x) if(is.factor(x)) as.factor(x) else x)

map2gui <- dplyr::right_join(fdnl2gui, gdxmapping, by='all_binid')

map2gui <- map2gui %>% dplyr::group_by(dplyr::across(dplyr::all_of({{ mapping }})))

return(map2gui)
}


## gdxreshape ----
#' Reshape from wide to long and save to GDX
#'
#' @description
#' `gdxreshape()` formats the data to be saved in GDX into long format. It is imported from the gdxrrw package with a few improvements for performance and usability, since there is a risk of it being removed from the gdxrrw package in the future.
#' We would like to thank the R GAMS team for this useful function.
#'
#' @param inDF wide dataframe.
#' @param symDim Dimensions of symbol
#' @param symName Symbol name
#' @param tName "Time"
#' @param gdxName Name of gdx file
#' @param setsToo if sets too
#' @param order order of data
#' @param setNames name of sets
#'
#' @return A tibble `tbl_df`.
#' @seealso
#' \itemize{
#' \item{\code{\link[gdxrrw]{wgdx}}}{Write R data to GDX}
#' \item{\code{\link[gdxrrw]{wgdx.lst}}}{Write multiple symbols to GDX}
#' \item{\code{\link[gdxrrw]{wgdx.reshape}}}{Write multiple symbols to GDX}
#' \item{\code{\link[tidyr]{pivot_longer}}}{Make dataframes longer}
#' }
#' @export gdxreshape
gdxreshape <- function (inDF, symDim, symName=NULL, tName="time",
                        gdxName=NULL, setsToo=TRUE, order=NULL,
                        setNames=NULL) {
  # Function based on gdxrrw::wgdx.reshape of the gdxrrw package, modified for performance and usability.

  nCols <- ncol(inDF)
  timeIdx <- symDim                     # default index position for time aggregate
  if (is.null(order)) {
    idCols <- 1:(symDim-1)

    inDF[idCols] <- lapply(inDF[idCols], as.factor)

    outDF <- (tidyr::pivot_longer(inDF, cols=-dplyr::all_of(idCols)))
  }
  else if ((! is.vector(order)) || (symDim != length(order))) {
    stop ("specified order must be a vector of length symDim")
  }
  else {
    timeIdx <- -1
    if (is.character(order)) {
      stop ("order must be numeric for now")
    }
    else if (! is.numeric(order)) {
      stop ("optional order vector must be numeric or character")
    }

    idCols <- 1:(symDim-1)                                   # for k in 1:symDim
    if (any(duplicated(order))) {
      stop ('duplicate entry in order vector: nonsense')
    }

    if ((symDim-1) != sum(order>0)) {
      stop ('order vector must specify symDim-1 ID columns')
    }
    if (all(order>0)) {
      stop ('order vector must have a non-positive entry to specify the "time" index')
    }

    timeIdx <- match(0, order)

    oo <- c(idCols,(1:nCols)[-idCols])
    df2 <- inDF[oo]
    idCols <- 1:(symDim-1)

    df2[idCols] <- lapply(df2[idCols], factor)

    if (symDim == timeIdx) {     # no need to re-order after reshaping
      outDF <- tidyr::pivot_longer(df2, cols=-dplyr::all_of(idCols))
    }
    else {
      df3 <- tidyr::pivot_longer(df2, cols=-dplyr::all_of(idCols))
      oo <- vector(mode="integer",length=symDim+1)

      oo[1:(timeIdx-1)] = 1:(timeIdx-1)

      oo[timeIdx] = symDim

      oo[(timeIdx+1):symDim] = (timeIdx+1):symDim-1

      oo[symDim+1] = symDim+1
      outDF <- (df3[oo])
    }
  }
  outDF$name <- as.factor(outDF$name)
  if (is.null(symName)) {
    symName <- attr(inDF, "symName", exact=TRUE)
  }
  if (! is.character(symName)) {
    stop ("symName must be a string")
  }
  attr(outDF,"symName") <- symName
  symText <- attr(inDF, "ts", exact=TRUE)
  if (is.character(symText)) {
    attr(outDF,"ts") <- symText
  }
  if (is.character(tName)) {
    names(outDF)[timeIdx] <- tName
  }
  else {
    names(outDF)[timeIdx] <- 'time'
  }
  names(outDF)[symDim+1] <- "value"
  # str(outDF)
  if (setsToo) {
    ## write index sets first, then symName
    outLst <- list()

    length(outLst) <- symDim + 1
    setNamesLen <- 0
    if (! is.null(setNames)) {
      if (! is.character(setNames)) {
        stop ("setNames must be a string or string vector")
      }
      else if (! is.vector(setNames)) {
        stop ("setNames must be a string or string vector")
      }
      else {
        ## guard against zero-length vector
        setNamesLen <- length(setNames)
      }
    }
    kk <- 0

    for (i in 1:symDim) {
      lvls <- levels(as.factor(outDF[[i]]))
      outLst[[i]] <- list(name=names(outDF)[[i]], type='set', uels=c(list(lvls)))
      if (setNamesLen > 0) {            # tack on the next set text
        kk <- kk + 1
        outLst[[i]]$ts <- setNames[[kk]]
        if (kk >= setNamesLen)
          kk <- 0
      }
    }
    outLst[[symDim+1]] <- (outDF)

    if (is.character(gdxName)) {
      gdxrrw::wgdx.lst(gdxName,outLst)
    }
    else {
      return(outLst)
    }
  }
  else {
    if (is.character(gdxName)) {
      gdxrrw::wgdx.lst(gdxName,outDF)
    }
    else {
      return(list(outDF))
    }
  }
} # gdxreshape


## groupstats ----
#' Generate descriptive statistics and save to GDX
#'
#' @description
#' This function has been conceived with the Dutch FADN in mind, please use `fd_dec()` for EU FADN data.
#' 
#' `groupstats()` returns descriptive statistics per group based on the mapping given. For example, if your mapping
#' is 'regions', this function will give you the weighted mean, median, min, max, number of observations per variable for each
#' region based on the individual farm data. When `writegdx` is `TRUE`, it writes the GDX in the format 'farmStats_(mapping).gdx'
#'
#' @inheritParams gdxbinwider
#' @param cols Which columns to derive descriptive statistics from
#' @param w Column with the Weights for the weighted mean
#' @param writegdx Logical. If `TRUE`, it writes a GDX with the descriptive statistics.
#' @param filtern Logical. If `TRUE`, results will be limited to more than 15 observations per variable for reporting
#'
#' @return A tibble `tbl_df`.
#' @examples
#' BINDir <- "inst/extdata/GAMS"
#' datafile <- 'FarmDynRexampledata.gdx'
#' groupstats('FarmDynRexampledata.gdx',
#'             BINDir="inst/extdata/GAMS/",
#'             gdxmap = 'map2binid',
#'             mapping = 'mapping',
#'             cols = c('a', 'b'),
#'             w='Weight')
#' @seealso
#' \itemize{
#' \item{\code{\link{summary}}}{summary statistics}
#' \item{\code{\link[psych]{describe}}}{Descriptive statistics}
#' \item{\code{\link[gdxrrw]{wgdx}}}{Write R data to GDX}
#' \item{\code{\link[gdxrrw]{wgdx.lst}}}{Write multiple symbols to GDX}
#' \item{\code{\link[tidyr]{pivot_longer}}}{Make dataframes longer}
#' }
#' @export groupstats
#FIXME fix the example data!! Not working
#TODO Include options for different types of aggregation based on different years and if only farms that are available every year should be selected.
groupstats <- function(filename, BINDir, gdxmap, mapping, cols, w, writegdx = TRUE, filtern = FALSE) {

  if ("BINDir" %in% ls(envir = .GlobalEnv) & missing(BINDir)) {
    BINDir <- get("BINDir", envir = .GlobalEnv)
  } else {
    BINDir
  }


  data <- gdxbinwider(filename, BINDir, gdxmap, mapping)
  names(data) <- gsub(x=names(data), pattern = "\\w*%", '')


  groupmap <- data %>%
    dplyr::select(dplyr::all_of(mapping), dplyr::all_of(cols), dplyr::all_of(w))  %>%
    dplyr::summarise(weightedmean = dplyr::across(dplyr::all_of(cols),~ stats::weighted.mean(as.numeric(as.character(.x)), w=.data[[w]],na.rm=TRUE)), # Make a new column named weightedmean where the values are the weighted means of only the numeric columns (otherwise error)
                     min = dplyr::across(dplyr::all_of(cols),~ min(as.numeric(as.character(.x))), na.rm =TRUE), # Same as weightedmean but with min
                     max = dplyr::across(dplyr::all_of(cols),~ max(as.numeric(as.character(.x))), na.rm =TRUE), # Idem
                     median = dplyr::across(dplyr::all_of(cols),~ stats::median(as.numeric(as.character(.x))), na.rm =TRUE), # Idem
                     mode = dplyr::across(dplyr::all_of(cols),~ Modes(as.numeric(as.character(.x)))[1]),
                     n = dplyr::across(dplyr::all_of(cols),~ sum(!is.na(.x))), # Make a column with n of each variable in the group
                     .groups = 'keep' # .keep is to keep the grouped groups as in the original mapping (otherwise it will group with only one group, but the results are the same)
  )

  groupmap <- groupmap %>% tidyr::unpack(cols = c(weightedmean, min, max, median, mode, n), names_sep = "%") %>% # Data cleaning activities and making a simple table
    tidyr::pivot_longer(cols = -dplyr::all_of(mapping), names_to = c('item1', 'variable'), names_sep = "%") %>%
    tidyr::pivot_wider(id_cols = c(dplyr::all_of(mapping), variable), names_from = item1, values_from = value)


  groupmap[] <- lapply(groupmap, function(x) if(is.numeric(x)) round(x, digits = 2) else x) # Rounding the results, comment if not needed
  groupmap[] <- lapply(groupmap, function(x) if(is.factor(x)) as.factor(x) else x) # For some reason, factors include ALL factors  (BIN_IDs, soil type, random things),
  # this results in gdx file with 1000s of unwanted data that serves no use, factor(x) eliminates that.

  if (filtern == TRUE) {
    groupmap <- dplyr::filter(groupmap, n >= 15) # Keep variables that have more than 15 observations (as seen in column 'n', created earlier)
  }
  if (writegdx == TRUE) {
    allmaps <- paste0(mapping[1:length(mapping)], '.gdx')
    gdxfilename <- paste('farmStats', allmaps, sep = '_')

    gdxreshape(gdxName = gdxfilename, as.data.frame(groupmap), sum(length(mapping), 2), symName = 'p_farmStats',
                 tName = 'colsFarmStats', order = c(1:sum(length(mapping),1),0)
    )
    return(groupmap)
  } else
    return(groupmap)
}


## updateFarmData ----
#' Reshape from wide to long and save to GDX
#' Create sample farms
#>>>>>>> 1ebeb0906d519cd31fea60640b1541768ae43ad8
#'
#' @description
#' This function has been conceived with the Dutch FADN in mind, please use `fadn2fd()` for EU FADN data.
#' 
#' `updateFarmData()` creates sample farms by aggregating data based on the weighted mean and the selected mapping for use in FarmDyn.
#' For non-numerical globals, it summarises based on the mode using the `Modes()` function. When `writegdx` is `TRUE`, it writes the GDX in the format 'farmData_(mapping).gdx'.
#'
#' @inheritParams groupstats
#' @param cptcoeffs Logical. When this is set to `TRUE` it calculates farm-specific parameters based on CPT coefficients
#' @param farmchars GDX file containing farm characteristics
#' @param cptcoeffsxl Location of the CPT coefficients excel file
#'
#' @return A tibble `tbl_df`.
#' @examples
#' BINDir <- "inst/extdata/GAMS"
#' datafile <- 'FarmDynRexampledata.gdx'
#' updateFarmData('FarmDynRexampledata.gdx',
#'         &BINDir="inst/extdata/GAMS/",
#'         &gdxmap = 'map2binid',
#'         &mapping = 'mapping',
#'         &w='Weight')
#' @seealso
#' \itemize{
#' \item{\code{\link[FarmDynR]{gdxbinwider}}}{Widens BIN data directly from GDX}
#' \item{\code{\link[FarmDynR]{gdxreshape}}}{Lengthens data and saves to GDX}
#' \item{\code{\link[gdxrrw]{wgdx}}}{Write R data to GDX}
#' \item{\code{\link[gdxrrw]{wgdx.lst}}}{Write multiple symbols to GDX}
#' \item{\code{\link[dplyr]{summarise}}}{Summarises data and aggregates to group}
#' \item{\code{\link[stats]{weighted.mean}}}{Calculates weighted mean}
#' }
#' @export updateFarmData
updateFarmData <- function(filename, BINDir,  gdxmap, mapping, writegdx = TRUE, cptcoeffs = FALSE, farmchars = NULL, cptcoeffsxl = NULL) {


  if ("BINDir" %in% ls(envir = .GlobalEnv) & missing(BINDir)) {
    BINDir <- get("BINDir", envir = .GlobalEnv)
  } else {
    BINDir
  }

  map2gui <- gdxbinwider(filename, BINDir, gdxmap, mapping)

  if (cptcoeffs == TRUE) {
    
    farmchars <- gdxload(farmchars, 'p_farmCharBIN', symbol = 'param', symName= "p_farmCharBIN", names=c('all_binid', 'year', 'char', 'value'), compress=FALSE, ts=FALSE,
                         squeeze=FALSE, useDomInfo = TRUE, check.names = TRUE)

    cptcoeffs <- read_xlsx(cptcoeffsxl,sheet = 'dat')

    farmchars <- farmchars[order(-as.numeric(as.character(farmchars$year))),]

    narabland <- map2gui %>% ungroup() %>%  select(c('all_binid', `global%nArabLand`))

    colnames(narabland)[2] <- 'ArabLand'

    farmchars <- farmchars %>% pivot_wider(id_cols = c(all_binid, year), values_from = value, names_from = char) %>%
      mutate('PropSalesArable' = 1)

    farmchars <- farmchars[!duplicated(farmchars$all_binid),]

    farmchars[is.na(farmchars)] <- 0 #NAs are produced when going from long to wide.

    farmchars <- inner_join(farmchars, narabland, by='all_binid') %>% mutate('logArabLandabs' = log(abs(ArabLand)))

    farmchars[] <- lapply(farmchars,function(x) if(is.factor(x)) factor(x) else x) %>% stats::na.omit()

    colnames(cptcoeffs)[1]<- 'char'

    alphaprod <- farmchars[colnames(farmchars) %in% cptcoeffs$char]*cptcoeffs$b_alpha[match(names(farmchars), cptcoeffs$char)][col(farmchars)] %>% stats::na.omit()

    betaprod <- farmchars[colnames(farmchars) %in% cptcoeffs$char]*cptcoeffs$b_beta[match(names(farmchars), cptcoeffs$char)][col(farmchars)] %>% stats::na.omit()

    gammaprod <- farmchars[colnames(farmchars) %in% cptcoeffs$char]*cptcoeffs$b_gamma[match(names(farmchars), cptcoeffs$char)][col(farmchars)] %>% stats::na.omit()

    deltaprod <- farmchars[colnames(farmchars) %in% cptcoeffs$char]*cptcoeffs$b_delta[match(names(farmchars), cptcoeffs$char)][col(farmchars)] %>% stats::na.omit()

    alpha <- rowSums(alphaprod)+cptcoeffs$b_alpha[7]

    beta <- rowSums(betaprod)+cptcoeffs$b_beta[7]

    delta <- rowSums(deltaprod)+cptcoeffs$b_delta[7]

    gamma <- rowSums(gammaprod)+cptcoeffs$b_gamma[7]

    coeffs <- list(data.frame('TKAlpha' = alpha, 'all_binid' = factor(farmchars$all_binid)),data.frame('TKBeta' = beta,'all_binid' = factor(farmchars$all_binid)),
                   data.frame('TKdelta' = delta, 'all_binid' = factor(farmchars$all_binid)), data.frame('TKgamma' = gamma, 'all_binid' = factor(farmchars$all_binid)))

    map2gui <- Reduce(function(x, y) merge(x, y, by="all_binid"), list(map2gui,coeffs))
    map2gui[] <- lapply(map2gui,function(x) if(is.factor(x)) factor(x) else x)

    map2gui <- map2gui %>%
      rename('global%TKAlpha' = 'TKAlpha',
             'global%TKBeta' = 'TKBeta',
             'global%TKdelta' = 'TKdelta',
             'global%TKgamma' = 'TKgamma'
      )
    map2gui$all_binid.1 <- NULL
    map2gui$all_binid.2 <- NULL
    map2gui$all_binid.3 <- NULL
  }

  map2gui <- map2gui %>% dplyr::select(-all_binid) %>%
    summarise(dplyr::across(everything(),~ if(is.numeric(.)) stats::weighted.mean(., w=.data[['Weight']], na.rm = TRUE) else Modes(.)))
  # In RStudio, ignore the (X) error unmatched bracket, everything is fine and all works.

  map2gui$Weight <- NULL

  map2gui[!colnames(map2gui) %in% mapping] <- lapply(map2gui[!colnames(map2gui) %in% mapping], function(x) if(is.factor(x)) (as.numeric(as.character(x))) else x)

  map2gui <- map2gui %>% tidyr::pivot_longer(cols = where(is.numeric), names_to = c('item1', 'item2'), names_sep = "%") %>%
    tidyr::pivot_wider(id_cols = c(dplyr::all_of(mapping), 'item1'), names_from = 'item2', values_from = 'value')
  map2gui[] <- lapply(map2gui, function(x) if(is.numeric(x)) round(x, digits = 2) else x) # Rounding the results, comment if not needed
  map2gui[is.na(map2gui)] <- 0 # By bringing these values from long to wide, there are undefined values and, therefore, NAs. BUT these are not real NAs, plus it's for modelling purposes on GAMS, which will ignore the 0s, so we can safely assign a 0 to the NAs.


  if (writegdx == TRUE) {
    allmaps <- paste0(mapping[1:length(mapping)], '.gdx')
    gdxfilename <- paste('farmData', allmaps, sep = '_')

    gdxreshape(map2gui, symDim = 3, order = c(1,2,0), gdxName = gdxfilename, symName = 'p_farmData')


    return(map2gui)
  } else
      return(map2gui)
}

## str_firstlast ----
#' @name str_line_replace
#' @rdname str_line_replace
#'
#' @title Replace first or last line in strings
#'
#' @description
#' These functions serve to change the first or last line of strings which match a specific pattern (regex).
#' `str_firstLine_replace()` replaces the first line that matches the pattern.
#' `str_lastLine_replace()` replaces the last line that matches the pattern
#' They are useful, for example, when reading a text file with many lines and you want to preserve the lines of that text file.
#' When `which='all'`, it is a wrapper for `stringr::str_replace()`.
#'
#' @param str String with pattern to make replacement
#' @param pattern Regular expression to replace
#' @param replacement What to replace the pattern with
#' @param which which one? first, last, all or the poles (first AND last)
#'
#' @return string
#' @examples
#'
#' somelines <- c('AAAAA', 'textytext', 'BBBBB', \\n
#'               'AAAAA', 'writingwriting', 'AAAAA', \\n
#'               'etc', 'etc', 'BBBBB')
#'
#' str_firstLine_replace(somelines, 'AAAAA', 'changedfirstline')
#'
#' str_lastLine_replace(somelines, 'AAAAA', 'changedlastline')
#'
#' str_line_replace(somelines, 'AAAAA', 'changedpoles', which='poles')
#'
#' str_line_replace(somelines, 'AAAAA', 'changedall', which='all')
#'
#' @seealso
#' [stringr::str_replace()]
#' @export
str_line_replace <- function(str, pattern, replacement, which=c('first', 'last', 'poles', 'all')) {

  if (which=='first') {
    str_firstLine_replace(str, pattern, replacement)
  }
  if (which=='last') {
    str_lastLine_replace(str, pattern, replacement)
  }
  if (which=='all') {
    return(stringr::str_replace(str, pattern, replacement))
  }
  if (which=='poles') {
    str_lastLine_replace(str_firstLine_replace(str, pattern, replacement), pattern, replacement)
  }
}

#' @rdname str_line_replace
#' @export str_firstLine_replace
str_firstLine_replace <- function(str, pattern, replacement) {
  str[grepl(pattern=pattern, x=str)][1] <- replacement
  return(str)
}
#' @rdname str_line_replace
#' @export str_lastLine_replace
str_lastLine_replace <- function(str, pattern, replacement) {
  str[grepl(pattern=pattern, x=str)][length(str[grepl(pattern=pattern, x=str)])] <- replacement
  return(str)
}

## writeBatch ----
#' Write batch file for batch file execution mode in FarmDyn
#'
#' @description
#' This function writes the batch file for you. It directly takes the necessary information from runInc.gms in FarmDyn, so the GUI
#' settings remain the same as you have set them.
#'
#' @inheritParams runFarmDynfromBatch
#' @inheritParams gdxbinwider
#' @param farmIds Individual farm Identifiers. This is usually the aggregation (mapping) of your p_farmData file. Whereas here the mapping is just the name of the aggregation, farmIds is a vector of the names of the different farm samples made. For example, if mapping = NUTS0, then famrIds would be AT, BE, BG... etc.
#'
#' @return Writes batch file necessary to run FarmDyn
#' @examples NULL
#'
#' @seealso
#' \code{\link[FarmDynR]{runFarmDynfromBatch}}
#'
#' @export writeBatch

writeBatch <- function(FarmDynDir, mapping, farmIds) {

  if ("FarmDynDir" %in% ls(envir = .GlobalEnv) & missing(FarmDynDir)) {
    FarmDynDir <- get("FarmDynDir", envir = .GlobalEnv)
  } else {
    FarmDynDir
  }

  readLines(
    file.path(FarmDynDir, 'gams', 'incgen', 'runInc.gms')
  )[((which(readLines(file.path('/FARMDYNTRUNK', 'gams', 'incgen', 'runInc.gms'))=='* Setting for executing the task in batch file mode'))-1):which(readLines(file.path('/FARMDYNTRUNK', 'gams', 'incgen', 'runInc.gms'))=='* end batch execution file')] %>%
    str_replace(pattern = 'Scenario description = \\w*', replacement = paste0('Scenario description = ', mapping)) %>%
    str_firstLine_replace(pattern='Farm sample file = \\w*', replacement = paste0('  Farm sample file = ', 'farmData_', mapping)) %>%
    str_lastLine_replace(pattern='Farm sample file = \\w*', replacement = paste0('  Farm sample file = ', 'farmData_', mapping, '\n  macro = ', paste(farmIds, collapse = '\\'))) %>%
    str_replace('farmIds = \\w*', replacement = paste0('farmIds = ', farmIds[length(farmIds)])) %>%
    str_replace('execute=Gamsrun', replacement = '  startparallel
    FOR farmidinloop = %allfarms%
      farmIds  = %farmidinloop%
      execute = Gamsrun
    ENDFOR
  collectparallel
*  execute=Gamsrun
') %>%
    writeLines(con = paste(mapping, 'batch.txt', sep = '_', collapse = ''))
}


## runFarmDynfromBatch ----
#' Execute FarmDyn
#'
#' @description
#' `runFarmDynfromBatch()` makes it possible to run FarmDyn from R using the batch file
#'
#' @param FarmDynDir Directory where FarmDyn is located
#' @param IniFile Name of the IniFile
#' @param XMLFile Name of the XML file
#' @param BATCHDir Directory where the .batch file is located
#' @param BATCHFile Name of the .batch file.
#'
#' @return Executes FarmDyn from R
#' @examples NULL
#' @seealso NULL
#' @export runFarmDynfromBatch

runFarmDynfromBatch <- function(FarmDynDir, IniFile, XMLFile, BATCHDir, BATCHFile) {
#
  # make sub directories
  GUIDir <- paste(FarmDynDir,"GUI",sep="/")
  BATCHFilePath <- paste(BATCHDir, BATCHFile, sep = "\\")

  # General JAVA command
  javacmdstrg <- r"(java -Xmx1G -Xverify:none -XX:+UseParallelGC -XX:PermSize=20M -XX:MaxNewSize=32M -XX:NewSize=32M -Djava.library.path=jars -classpath jars\gig.jar de.capri.ggig.BatchExecution)"

  # append specific files to JAVA command
  javacmdparac <- paste(javacmdstrg,IniFile,XMLFile,BATCHFilePath,sep = " ")

  # create bat file
  runbat   = paste0(GUIDir,"/runfarmdyn.bat")
  if (file.exists(runbat)) x=file.remove(runbat)

  b = substr(runbat,1,2)
  b = c(b,paste('cd',gsub("/", "\\\\",GUIDir)))
  b = c(b,c("SET PATH=%PATH%;./jars"))
  b = c(b,javacmdparac)
  writeLines(b,runbat)
  rm(b)

  # execute farmdyn in batch mode
  system(runbat)


}

## gdxload ----
#' Load from GDX
#'
#' @description
#' `gdxload()` is a wrapper for gdxrrw::rgdx. It fixes the issue of having all UELs as factor levels for variables for which they do not belong.
#'
#' @param filename Name of the GDX file and its location
#' @param symbol The symbol type to be loaded (set, parameter, scalar)
#' @param symName The name of the symbol to be loaded
#' @param names Column names of the symbol
#' @param ...  Arguments to be passed to gdxrrw::rgdx(...)
#'
#' @return A dataframe
#'
#' @seealso
#' *gdxrrw::rgdx.set()
#' *gdxrrw::rgdx.param()
#' *gdxrrw::rgdx.scalar()
#'
#' @export gdxload

gdxload <- function(filename, symbol=c('set', 'param', 'scalar'), symName, names = NULL, ...) {
  if (symbol=='set') {
    symb <- gdxrrw::rgdx.set(filename, symName, names)
  }
  if (symbol=='param') {
    symb <- gdxrrw::rgdx.param(filename, symName, names)
  }
  if (symbol=='scalar') {
    symb <- gdxrrw::rgdx.scalar(filename, symName)
  }

  symb[] <- lapply(symb, function(x) if(is.factor(x)) factor(x) else x)
  return(symb)
}

## rgdx.var ----
#' Improved function to load variables from gdx files (from Renger in GAMSWorld forum)
#' @param varname Name of the variable to load
#' @param scen_name Name of the scenarios to load
#' @param res_fold Path to the folder where the gdx files are
#' @return A data.table with the variable
#' @export rgdx.var
#' @references https://forum.gamsworld.org/viewtopic.php?t=9966
rgdx.var <- function(varname, scen_name, res_fold) {
  varname <- rgdx(
    paste(
      res_fold, scen_name,
      sep = "/"
    ),
    requestList = list(name = varname)
  )
  var.data <- data.frame(varname$val)
  var.dim <- length(varname$uels)
  domains <- varname$domains
  for (j in (1:(var.dim))) {
    if (domains[j] == "*") {
      domains[j] <- paste("X", j, sep = "")
    }
  }
  for (i in 1:var.dim) {
    dim <- data.frame(varname$uels[[i]])
    dim$id <- seq_along(dim[, 1])
    index <- varname$domains[i]
    colnames(dim) <- c(index, "id")
    var.colname <- paste("X", i, sep = "")
    var.data <- merge(dim, var.data, by.x = "id", by.y = var.colname)
    var.data <- var.data[, -which(colnames(var.data) == "id")]
  }
  var.data <- var.data[, c(var.dim:1, var.dim + 1)]
  colnames(var.data)[var.dim + 1] <- c("value")
  colnames(var.data)[var.dim] <- "field"
  attributes(var.data)$domains <- varname$domains
  attributes(var.data)$type <- "variable"
  attributes(var.data)$symName <- varname$name
  attributes(var.data)$description <- varname$description
  return(var.data)
}

## vars_dump ----
#' Load variables from dump gdx files and return a data.table
#' @inheritParams rgdx.var
#' @return A data.table with the variable
#' @export vars_dump
vars_dump <- function(res_fold, varname, scen_name) {
  files_in_res <- list.files(res_fold, pattern = "gdx", full.names = FALSE)
  files <- lapply(scen_name, function(x) {
    files_in_res[grepl(paste0("dump_", x, "_", "[[:alnum:]]{4}"), files_in_res)]
  })
  scen_call <- list()
  var_call <- list()
  for (y in seq_along(scen_name)) {
    scen_call[[y]] <- lapply(files[[y]], function(x) {
      rgdx.var(
        varname = varname,
        scen_name = x,
        res_fold = res_fold
      ) %>%
        select(-c("t")) %>%
        mutate(
          farmIds = gsub("^\\w+_", "", x) %>% gsub(".gdx", "", .)
        )
    })
    # Add farm ids (4 last letters and numbers of the file name) to the data.table
    var_call[[y]] <-
      var_call[[y]] <- rbindlist(scen_call[[y]])
    var_call[[y]]$sims <- scen_name[y]
  }
  return(var_call)
}

## scen_analysis ----
#' Data analysis function for repeated tasks
#' This loads p_res from different scenarios
#' @inheritParams rgdx.var
#' @return A data.table with the parameters
#' @export scen_analysis
scen_analysis <- function(res_fold, scen_name) {
  # Load the sims for all the files in files_in_resNUTS2 using rgdx.param and load p_res
  # Then, select the columns "farmIds", "scen", "resItem1", "resItem2", "value" from the data.table
  scen_call <- list()
  for (i in seq_along(scen_name)) {
    scen_call[[i]] <- gdxrrw::param(
      paste(
        res_fold, paste0("res_", scen_name[i], "_until_2010.gdx"),
        sep = "/"
      ), "p_res",
      names = c(
        "sims", "farmIds", "resItem1", "resItem2", "resItem3", "time", "value"
      ), compress = FALSE, ts = FALSE,
      squeeze = FALSE, useDomInfo = TRUE, check.names = TRUE
    )
    scen_call[[i]]$sims <- scen_name[i]
  }
  scen_call <- rbindlist(scen_call)

  scen_call <- scen_call[scen_call$time == "mean", ]

  return(scen_call)
}

## scen_diff ----
#' Calculates the relative difference of the reference and the scenarios in a new column
#' @param data A dataframe
#' @param vars_to_diff A vector with the names of the variables to calculate the difference
#' @return A dataframe with the new columns
#' @export scen_diff
# Function for diffs
scen_diff <- function(data, vars_to_diff) {
  for (i in seq_along(vars_to_diff)) {
    data <- data %>%
      dplyr::group_by(farmIds) %>%
      dplyr::mutate(
        !!paste0(vars_to_diff[i], "_diff") := (get(vars_to_diff[i]) - first(get(vars_to_diff[i]))) / first(get(vars_to_diff[i]))
      )
  }
  return(data)
}

## abs_diff ----
#' Calculates the absolute difference of the reference and the scenarios in a new column
#' @inheritParams scen_diff
#' @return A dataframe with the new columns
#' @export abs_diff
abs_diff <- function(data, vars_to_diff) {
  for (i in seq_along(vars_to_diff)) {
    data <- data %>%
      dplyr::group_by(farmIds) %>%
      dplyr::mutate(
        !!paste0(vars_to_diff[i], "_absolute_diff") := (get(vars_to_diff[i]) - first(get(vars_to_diff[i])))
      )
  }
  return(data)
}

# Function in case I need to use the dumps instead of res
#' `load_dumps()` loads the p_res from the dump gdx files
#' @inheritParams scen_analysis
#' @return A data.table with the parameters
#' @export load_dumps
load_dumps <- function(res_fold, scen_name) {
  files_in_res <- list.files(res_fold, pattern = "gdx", full.names = FALSE)
  files <- lapply(scen_name, function(x) {
    files_in_res[grepl(paste0("dump_", x, "_", "[[:alnum:]]{4}"), files_in_res)]
  })
  scen_call <- list()
  for (i in seq_along(scen_name)) {
    scen_call[[i]] <- lapply(files[[i]], function(x) {
      gdxrrw::param(
        paste(
          res_fold, x,
          sep = "/"
        ), "p_res",
        names = c(
          "sims", "farmIds", "resItem1", "resItem2", "resItem3", "time", "value"
        ), compress = FALSE, ts = FALSE,
        squeeze = FALSE, useDomInfo = TRUE, check.names = TRUE
      )
    })
    scen_call[[i]] <- rbindlist(scen_call[[i]])
    scen_call[[i]]$sims <- scen_name[i]
  }
  return(scen_call)
}

# Function to load any dump parameter
## load_dump_par ----
#' `load_dump_par()` loads any parameter from the dump gdx files
#' @inheritParams scen_analysis
#' @param param Name of the parameter to load
#' @param names Names of the columns to load
#' @return A data.table with the parameters
#' @export load_dump_par
load_dump_par <- function(res_fold, scen_name, param, names = NULL) {
  files_in_res <- list.files(res_fold, pattern = "gdx", full.names = FALSE)
  files <- lapply(scen_name, function(x) {
    files_in_res[grepl(paste0("dump_", x, "_", "[[:alnum:]]{4}"), files_in_res)]
  })
  scen_call <- list()
  for (i in seq_along(scen_name)) {
    scen_call[[i]] <- lapply(files[[i]], function(x) {
      gdxrrw::param(
        paste(
          res_fold, x,
          sep = "/"
        ), param,
        names = names, compress = FALSE, ts = FALSE,
        squeeze = FALSE, useDomInfo = TRUE, check.names = TRUE
      ) %>%
        mutate(
          farmIds = gsub("^\\w+_", "", x) %>% gsub(".gdx", "", .)
        )
    })
    scen_call[[i]] <- rbindlist(scen_call[[i]])
    scen_call[[i]]$sims <- scen_name[i]
  }
  return(scen_call)
}

# Function for other things that rgdx is too stubborn to take out of the dumps
## load_dump_scalar ----
#' `load_dump_scalar()` loads any scalar from the dump gdx files
#' @inheritParams scen_analysis
#' @param scalar Name of the scalar to load
#' @return A data.table with the scalar
#' @export load_dump_scalar
#' @examples
#' load_dump_scalar(res_fold = res_fold, scen_name = scen_name, scalar = "p_Nmin")
load_dump_scalar <- function(res_fold, scen_name, scalar) {
  files_in_res <- list.files(res_fold, pattern = "gdx", full.names = FALSE)
  files <- lapply(scen_name, function(x) {
    files_in_res[grepl(paste0("dump_", x, "_", "[[:alnum:]]{4}"), files_in_res)]
  })
  scen_call <- list()
  for (i in seq_along(scen_name)) {
    scen_call[[i]] <- lapply(files[[i]], function(x) {
      rgdx(
        paste(
          res_fold, x,
          sep = "/"
        ), list(name = scalar)
      )$val
    })
    # scen_call[[i]] <- rbindlist(t(scen_call[[i]]))
    names(scen_call[[i]]) <- gsub("^\\w+_", "", files[[i]]) %>% gsub(".gdx", "", .)
    scen_call[[i]]$sims <- scen_name[i]
    # scen_call[[i]]$farmIds <- gsub("^\\w+_", "", files[[i]]) %>% gsub(".gdx", "", .)
  }
  scen_call <- t(scen_call)
  scen_call <- rbindlist(scen_call)
  scen_call <- pivot_longer(scen_call, cols = c(-sims), names_to = "farmIds", values_to = scalar)
  return(scen_call)
}

# Function for other things that rgdx is too stubborn to take out of the dumps
## load_dump_scalar ----
#' `load_dump_marg()` loads any marginal from the dump gdx files
#' @inheritParams scen_analysis
#' @param marginal Name of the marginal to load
#' @return A data.table with the marginal
#' @export load_dump_marg
load_dump_marg <- function(res_fold, scen_name, marginal) {
  files_in_res <- list.files(res_fold, pattern = "gdx", full.names = FALSE)
  files <- lapply(scen_name, function(x) {
    files_in_res[grepl(paste0("dump_", x, "_", "[[:alnum:]]{4}"), files_in_res)]
  })
  scen_call <- list()
  for (i in seq_along(scen_name)) {
    scen_call[[i]] <- lapply(files[[i]], function(x) {
      rgdx(
        paste(
          res_fold, x,
          sep = "/"
        ), list(name = marginal, field = "m")
      )$val[, 4]
    })
    # scen_call[[i]] <- rbindlist(t(scen_call[[i]]))
    names(scen_call[[i]]) <- gsub("^\\w+_", "", files[[i]]) %>% gsub(".gdx", "", .)
    scen_call[[i]]$sims <- scen_name[i]
    # scen_call[[i]]$farmIds <- gsub("^\\w+_", "", files[[i]]) %>% gsub(".gdx", "", .)
  }
  scen_call <- t(scen_call)
  scen_call <- rbindlist(scen_call)
  scen_call <- pivot_longer(scen_call, cols = c(-sims), names_to = "farmIds", values_to = marginal)
  return(scen_call)
}


# Function for descriptive statistics WITH exclusion of farms with less than 15 farms
## fd_desc ----
#' `fd_desc()` calculates the descriptive statistics of the farm data
#' @param farm_data A dataframe with the p_farmData
#' @param type Type of farm data to analyse (Dairy or Arable farms)
#' @param csv Logical. If TRUE, it saves the results as a csv
#' @param dir If csv is TRUE, where to save the csv
#' @return A dataframe with the descriptive statistics
#' @export fd_desc
fd_desc <- function(farm_data, type = c("arable", "dairy"), csv = FALSE, dir = NULL) {
  ifelse(!"NUTS0" %in% colnames(farm_data), farm_data$NUTS0 <- substr(farm_data$farmIds, 1, 2), farm_data$NUTS0 <- farm_data$NUTS0)
  if (type == "arable") {
    descstats <- farm_data %>%
      group_by(NUTS0) %>%
      summarise(
        `Land [ha]` = stats::weighted.mean(`global%nArabLand`, `misc%weights`, na.rm = TRUE),
        SummerCere = stats::weighted.mean(`misc%SummerCere`, `misc%weights`, na.rm = TRUE),
        Winterbarley = stats::weighted.mean(`misc%Winterbarley`, `misc%weights`, na.rm = TRUE),
        WinterWheat = stats::weighted.mean(`misc%WinterWheat`, `misc%weights`, na.rm = TRUE),
        MaizCorn = stats::weighted.mean(`misc%MaizCorn`, `misc%weights`, na.rm = TRUE),
        MaizSil = stats::weighted.mean(`misc%MaizSil`, `misc%weights`, na.rm = TRUE),
        #        `Annual Work Units` = median(`global%Aks`, `misc%weights`, na.rm = TRUE),
        `Farm Net Value Added [EUR]` = stats::weighted.mean(`misc%net cashflow`, `misc%weights`, na.rm = TRUE),
        `Median FNVA [EUR]` = stats::median(`misc%net cashflow`, na.rm = TRUE),
        `Annual Work Units` = stats::weighted.mean(`global%Aks`, `misc%weights`, na.rm = TRUE),
        `FNVA per AWU` = `Farm Net Value Added [EUR]` / `Annual Work Units`,
        N_use = stats::weighted.mean(`misc%N_use`, `misc%weights`, na.rm = TRUE),
        n = sum(`misc%nFarms`)
      )
    descstats <- descstats %>%