R/toolAggregate.R

In madrat: May All Data be Reproducible and Transparent (MADRaT) *

#' toolAggregate
#'
#' (Dis-)aggregates a magclass object from one resolution to another based on a
#' relation matrix or mapping
#'
#' Basically toolAggregate is doing nothing more than a normal matrix
#' multiplication which is taking into account the 3 dimensional structure of
#' MAgPIE objects. So, you can provide any kind of relation matrix you would
#' like. However, for easier usability it is also possible to provide weights
#' for a weighted (dis-)aggregation as a MAgPIE object. In this case rel must
#' be a 1-0-matrix or a mapping between both resolutions. The weight
#' needs to be provided in the higher spatial aggregation, meaning for
#' aggregation the spatial resolution of your input data and in the case of
#' disaggregation the spatial resolution of your output data. The temporal and
#' data dimension must be either identical to the resolution of the data set
#' that should be (dis-)aggregated or 1. If the temporal and/or data dimension
#' is 1 this means that the same transformation matrix is applied for all years
#' and/or all data columns. In the case that a column should be just summed up
#' instead of being calculated as a weighted average you either do not provide
#' any weight (then all columns are just summed up) or your set this specific
#' weighting column to NA and mixed_aggregation to TRUE.
#'
#' @param x magclass object that should be (dis-)aggregated
#' @param rel relation matrix, mapping or file containing a mapping in a format
#' supported by \code{\link{toolGetMapping}} (currently csv, rds or rda).
#' A mapping object consists of any number of columns, where one column contains
#' all the elements in x. These elements are mapped to the corresponding values
#' in another column, as described below (see parameters 'from' and 'to').
#' It is possible to not set \code{rel} as long as \code{to} is set and \code{dim}
#' is chosen appropriately. In that case the relation mapping is extracted from
#' the dimnames of the corresponding dimension, e.g. if your data contains a
#' spatial subdimension "country" you can aggregate to countries via
#' \code{toolAggregate(x, to = "country", dim = 1)}.
#' @param weight magclass object containing weights to be used for a weighted
#' (dis)aggregation. The provided weight does not need to be normalized, any
#' number >= 0 is allowed. See "details" section below.
#' @param from Name of source column to be used in rel if it is a
#' mapping (if not set the first column matching the data will be used).
#' @param to Name of the target column to be used in rel if it is a
#' mapping (if not set the column following column \code{from} will be used
#' If column \code{from} is the last column, the column before \code{from is
#' used}). If data should be aggregated based on more than one column these
#' columns can be specified via "+", e.g. "region+global" if the data should
#' be aggregated to column regional as well as column global.
#' In the to column, empty values (in the form of "") will not be included in
#' the aggregated data (but used as normal aggregation targets).
#' If \code{rel} is missing \code{to} refers to the aggregation target dimension name.
#' @param dim Specifying the dimension of the magclass object that should be
#' (dis-)aggregated. Either specified as an integer
#' (1=spatial,2=temporal,3=data) or if you want to specify a sub dimension
#' specified by name of that dimension or position within the given dimension
#' (e.g. 3.2 means the 2nd data dimension, 3.8 means the 8th data dimension).
#' @param wdim Specifying the according weight dimension as chosen with dim
#' for the aggregation object. If set to NULL the function will try to
#' automatically detect the dimension.
#' @param partrel If set to TRUE allows that the relation matrix does contain
#' less entries than x and vice versa. These values without relation are lost
#' in the output.
#' @param negative_weight Describes how a negative weight should be treated. "allow"
#' means that it just should be accepted (dangerous), "warn" returns a warning and
#' "stop" will throw an error in case of negative values
#' @param mixed_aggregation boolean which allows for mixed aggregation (weighted
#' mean mixed with summations). If set to TRUE weight columns filled with NA
#' will lead to summation.
#' @param verbosity Verbosity level of messages coming from the function: -1 = error,
#' 0 = warning, 1 = note, 2 = additional information, >2 = no message
#' @param zeroWeight Describes how a weight sum of 0 for a category/aggregation target should be treated.
#' "allow" accepts it and returns 0 (dangerous), "setNA" returns NA, "warn" throws a warning, "stop" throws an error,
#' "fix" will set zero weights to 10^-30 only where needed, see \code{\link{toolFixWeight}}.
#' @return the aggregated data in magclass format
#'
#' @author Jan Philipp Dietrich, Ulrich Kreidenweis, Pascal Sauer
#' @seealso \code{\link{calcOutput}}
#' @family aggregation
#' @examples
#'
#' # create example mapping
#' p <- magclass::maxample("pop")
#' mapping <- data.frame(from = magclass::getItems(p, dim = 1.1),
#'                       region = rep(c("REG1", "REG2"), 5),
#'                       global = "GLO")
#' print(mapping)
#'
#' # run aggregation
#' toolAggregate(p, mapping)
#' # weighted aggregation
#' toolAggregate(p, mapping, weight = p)
#' # combined aggregation across two columns
#' toolAggregate(p, mapping, to = "region+global")
#'
#' @export
toolAggregate <- function(x,
                          rel, weight = NULL, from = NULL, to = NULL, dim = 1, wdim = NULL, partrel = FALSE,
                          negative_weight = "warn", mixed_aggregation = FALSE, # nolint: object_name_linter.
                          verbosity = 1, zeroWeight = "warn") {

  if (!is.magpie(x)) {
    stop("Input is not a MAgPIE object, x has to be a MAgPIE object!")
  }

  xComment <- getComment(x)

  # create missing rel information from dimension names if argument "to" is set
  if (missing(rel) && !is.null(to)) {
    if (round(dim) != dim) {
      stop("Subdimensions in dim not supported if relation mapping is missing!")
    }
    rel <- data.frame(c(dimnames(x)[dim], getItems(x, dim = dim, split = TRUE, full = TRUE)))
    if (dim == 1 && is.null(rel$global)) {
      rel$global <- "GLO"  # add global column
    }
    if (is.null(rel$all)) {
      rel$all <- "all"
    }
  }

  isMatrix <- (any(class(rel) %in% c("dgCMatrix", "dsCMatrix", "dsyMatrix"))
               || "Matrix" %in% attr(class(rel), "package"))

  if (!is.numeric(rel) && !isMatrix) {
    if (length(to) == 1 && grepl("+", to, fixed = TRUE)) {
      tmprel <- NULL
      toSplit <- strsplit(to, "+", fixed = TRUE)[[1]]
      for (t in toSplit) {
        tmp <- toolGetAggregationMatrix(rel, from = from, to = t, items = getItems(x, dim = dim), partrel = partrel)
        tmprel <- rbind(tmprel, tmp)
      }
      rel <- tmprel
    } else {
      rel <- toolGetAggregationMatrix(rel, from = from, to = to, items = getItems(x, dim = dim), partrel = partrel)
      if (is.null(to)) {
        to  <- names(dimnames(rel))[1]
      }
    }
  }

  # translate dim to dim code
  dim <- dimCode(dim, x, missing = "stop")


  # allow the aggregation, even if the relation matrix is missing some entries from the initial dataset
  if (partrel) {
    datnames <-  getItems(x, dim)

    common <- intersect(datnames, colnames(rel))
    x <- x[common, dim = floor(dim)]

    # datanames not in relnames
    noagg <- datnames[!datnames %in% colnames(rel)]
    if (length(noagg) > 0) {
      vcat(verbosity, paste0(noagg, collapse = ", "), " not mapped in aggregation!")
    }
    rel <- rel[, common, drop = FALSE]
    rel <- rel[Matrix::rowSums(rel) > 0, , drop = FALSE]
  }

  if (!is.null(weight)) {
    result <- toolAggregateWeighted(x = x, rel = rel, weight = weight, from = from, to = to, dim = dim,
                                    wdim = wdim, partrel = partrel, negativeWeight = negative_weight,
                                    mixedAggregation = mixed_aggregation, zeroWeight = zeroWeight,
                                    xComment = xComment)
  } else {
    result <- toolAggregateUnweighted(x = x, rel = rel, to = to, dim = dim, xComment = xComment)
  }
  return(removeEmptyResultItems(result, dim = dim))
}

toolAggregateWeighted <- function(x, rel, weight, from, to, dim, wdim, partrel,
                                  negativeWeight, mixedAggregation, zeroWeight, xComment) {
  if (!is.magpie(weight)) {
    stop("Weight is not a MAgPIE object, weight has to be a MAgPIE object!")
  }

  # get proper weight dim
  if (is.null(wdim)) {
    wdim <- union(getDim(rownames(rel), weight, fullmatch = TRUE),
                  getDim(colnames(rel), weight, fullmatch = TRUE))
    # wdim must be in same main dimension as dim
    wdim <- wdim[floor(wdim) == floor(dim)]
  }

  if (length(wdim) == 0) {
    stop("Could not detect aggregation dimension in weight (no match)!")
  }
  if (length(wdim) > 1) {
    # if full dimension and subdimension is matched, use only full dimension
    if (any(wdim == floor(dim))) {
      wdim <- floor(dim)
    } else {
      stop("Could not detect aggregation dimension in weight (multiple matches)!")
    }
  }

  if (floor(dim) == dim) {
    wdim <- floor(wdim)
  }

  if (anyNA(weight)) {
    if (!mixedAggregation) {
      stop("Weight contains NAs which is only allowed if mixed_aggregation=TRUE!")
    } else {
      n <- length(getItems(weight, dim = wdim))
      r <- dimSums(is.na(weight), dim = wdim)
      if (!all(r %in% c(0, n))) {
        stop("Weight contains columns with a mix of NAs and numbers which is not allowed!")
      }
    }
  }
  if (nyears(weight) == 1) {
    getYears(weight) <- NULL
  }
  weight <- collapseNames(weight)
  if (negativeWeight != "allow" && any(weight < 0, na.rm = TRUE)) {
    if (negativeWeight == "warn") {
      warning("Negative numbers in weight. Dangerous, was it really intended?")
    } else {
      stop("Negative numbers in weight. Weight should be positive!")
    }
  }

  weightAndWeightSum <- toolZeroWeight(weight, rel, from, to, dim, wdim, partrel, zeroWeight)
  weight <- weightAndWeightSum$weight
  weight2 <- 1 / (weightAndWeightSum$weightSum + 10^-100)

  if (mixedAggregation) {
    weight2[is.na(weight2)] <- 1
    weight[is.na(weight)] <- 1
  }

  if (wdim != floor(wdim)) {
    getSets(weight)[paste0("d", wdim)]  <- getSets(x)[paste0("d", dim)]
    getSets(weight2)[paste0("d", wdim)] <- getSets(x)[paste0("d", dim)]
  }

  if (setequal(getItems(weight, dim = wdim), getItems(x, dim = dim))) {
    out <- toolAggregate(x * weight, rel, from = from, to = to, dim = dim, partrel = partrel) * weight2
  } else {
    out <- toolAggregate(x * weight2, rel, from = from, to = to, dim = dim, partrel = partrel) * weight
  }

  if (zeroWeight == "fix" && !isTRUE(all.equal(sum(out), sum(x)))) {
    warning("total sum is not the same after toolAggregate despite zeroWeight = 'fix', ",
            "please contact a madrat developer.")
  }

  getComment(out) <- c(xComment, paste0("Data aggregated (toolAggregate): ", date()))
  return(out)
}

toolZeroWeight <- function(weight, rel, from, to, dim, wdim, partrel, zeroWeight) {
  weightSum <- toolAggregate(weight, rel, from = from, to = to, dim = wdim, partrel = partrel, verbosity = 10)
  if (zeroWeight != "allow" && any(weightSum == 0, na.rm = TRUE)) {
    if (zeroWeight == "warn") {
      warning("Weight sum is 0, so cannot normalize and will return 0 for some ",
              "aggregation targets. This changes the total sum of the magpie object! ",
              'If this is really intended set zeroWeight = "allow", or "setNA" to return NA.')
    } else if (zeroWeight == "setNA") {
      weightSum[weightSum == 0] <- NA
    } else if (zeroWeight == "fix") {
      weight <- toolFixWeight(weight, toolMapFromRel(rel, from, to), dim)
      weightSum <- toolAggregate(weight, rel, from = from, to = to, dim = wdim, partrel = partrel, verbosity = 10)
      if (any(weightSum == 0, na.rm = TRUE)) {
        warning("Critical: toolAggregate zeroWeight = fix did not fix the weight! Please contact a madrat developer.")
      }
    } else {
      stop("Weight sum is 0, so cannot normalize. This changes the total sum of the magpie object!")
    }
  }
  return(list(weight = weight, weightSum = weightSum))
}

toolAggregateUnweighted <- function(x, rel, to, dim, xComment) {
  # convert rel for better performance
  rel <- Matrix::Matrix(rel)

  # make sure that rel covers a whole dimension (not only a subdimension)
  # expand data if necessary
  # set dim to main dimension afterwards
  if (round(dim) != dim) {
    rel <- toolExpandRel(rel, x, dim)
    dim <- round(floor(dim))
  }

  if (dim(x)[dim] != dim(rel)[2]) {
    if (dim(x)[dim] != dim(rel)[1]) {
      stop("Relation matrix has in both dimensions a different number of entries (",
           dim(rel)[1], ", ", dim(rel)[2], ") than x has cells (", dim(x)[dim], ")!")
    }
    rel <- Matrix::t(rel)
  } else if (dim(x)[dim] == dim(rel)[1] && !setequal(colnames(rel), getItems(x, dim))) {
    rel <- Matrix::t(rel)
  }

  # reorder MAgPIE object based on column names of relation matrix if available
  if (!is.null(colnames(rel))) {
    x <- x[colnames(rel), dim = dim]
  }

  # Aggregate data
  if (anyNA(x) || any(is.infinite(x))) {
    out <- apply(x, which(1:3 != dim), toolMatrixMultiplication, rel)
    if (length(dim(out)) == 2) {
      out <- array(out, dim = c(1, dim(out)), dimnames = c("", dimnames(out)))
    }
  } else {
    withr::local_options(matprod = "blas")
    notdim <- setdiff(1:3, dim)
    out <- rel %*% as.array(wrap(x, list(dim, notdim)))
    out <- array(out, dim = c(dim(rel)[1], dim(x)[notdim]))
    dimnames(out)[2:3] <- dimnames(x)[notdim]
  }

  # Write dimnames of aggregated dimension
  if (!is.null(rownames(rel))) {
    regOut <- rownames(rel)
  } else if (dim == 1) {
    regionList <- as.factor(getItems(x, dim = 1.1, full = TRUE))
    # Compute region vector for outputs after aggregation via sending
    # factor values through the relation matrix
    regOut <- factor(as.vector(round(rel %*% as.numeric(regionList) /
                                       (rel %*% rep(1, dim(rel)[2])))))
    levels(regOut) <- levels(regionList)
  } else {
    stop("Missing dimnames for aggregated dimension")
  }

  if (!any(grepl(".", regOut, fixed = TRUE)) && anyDuplicated(regOut)) {
    regOut <- paste(regOut, seq_len(dim(out)[1]), sep = ".")
  }

  dimnames(out)[[1]] <- regOut

  if (dim == 2) {
    out <- wrap(out, map = list(2, 1, 3))
  } else if (dim == 3) {
    out <- wrap(out, map = list(2, 3, 1))
  }

  sets <- getSets(x, fulldim = FALSE)
  # update set name if number of sub-dimensions reduced to 1
  if (ndim(out, dim = dim) == 1  && ndim(x, dim = dim) > 1) {
    sets[dim] <- ifelse(!is.null(to), to, NA)
  }
  getSets(out, fulldim = FALSE) <- sets

  getComment(out) <- c(xComment, paste0("Data aggregated (toolAggregate): ", date()))
  out <- as.magpie(out, spatial = 1, temporal = 2)
  return(out)
}

toolMatrixMultiplication <- function(y, x) {
  if (any(is.infinite(y))) {
    # Special Inf treatment to prevent that a single Inf in x
    # is setting the full output to NaN (because 0*Inf is NaN)
    # Infs are now treated in a way that anything except 0 times Inf
    # leads to NaN, but 0 times Inf leads to NaN
    for (i in c(-Inf, Inf)) {
      j <- (is.infinite(y) & (y == i))
      x[, j][x[, j] != 0] <- i
      y[j] <- 1
    }
  }
  if (any(is.na(y)) && !all(is.na(y))) {
    # Special NA treatment to prevent that a single NA in x
    # is setting the full output to NA (because 0*NA is NA)
    # NAs are now treated in a way that anything except 0 times NA
    # leads to NA, but 0 times NA leads to 0
    x[, is.na(y)][x[, is.na(y)] != 0] <- NA
    y[is.na(y)] <- 0
  }
  return(as.array(x %*% y))
}

toolGetAggregationMatrix <- function(rel, from = NULL, to = NULL, items = NULL, partrel = FALSE) {
  if ("tbl" %in% class(rel)) {
    rel <- data.frame(rel)
  }
  if (!(is.matrix(rel) || is.data.frame(rel))) {
    if (length(rel) > 1) {
      stop("Malformed relation mapping!")
    }
    if (!file.exists(rel)) {
      stop("Cannot find region mapping file: ", rel, " (working directory ", getwd(), ")")
    }
    rel <- toolGetMapping(rel, where = "local")
  }
  if (is.matrix(rel)) {
    rel <- as.data.frame(rel)
  }

  if (length(rel) < 2) {
    stop("relation mapping has only ", length(rel), " column!")
  }

  if (is.null(from)) {
    if (partrel) {
      from <- as.integer(which(vapply(lapply(rel, intersect, items), length, numeric(1)) > 0))
    } else {
      from <- as.integer(which(vapply(rel, setequal, logical(1), items)))
    }
    if (length(from) == 0) {
      maxMatchColumn <- which.max(vapply(lapply(rel, intersect, items), length, numeric(1)))
      unmappedItems <- setdiff(items, rel[[maxMatchColumn]])
      missingItems <- setdiff(rel[[maxMatchColumn]], items)

      if (length(unmappedItems) > 0) {
        warning(paste("The following items were not found in the mapping:", toString(unmappedItems)))
      }
      if (length(missingItems) > 0) {
        warning(paste("The mapping expected the following items, but they are missing:", toString(missingItems)))
      }
      stop("Complete mapping failed. If you want a partial mapping, call toolAggregate(..., partrel = TRUE).")
    }
    if (length(from) > 1) {
      from <- from[1]
    }
  }
  if (is.null(to)) {
    if (from == length(rel)) {
      to <- from - 1
    } else {
      to <- from + 1
    }
  }

  toItems <- as.character(unique(rel[, to]))
  fromItems <- as.character(unique(rel[, from]))
  m <- Matrix::Matrix(data = 0, nrow = length(toItems), ncol = length(fromItems),
                      dimnames = list(toItems = toItems, fromItems = fromItems))
  m[cbind(match(rel[, to], rownames(m)), match(rel[, from], colnames(m)))] <- 1
  if (is.numeric(to)) {
    to <- dimnames(rel)[[2]][to]
  }
  if (is.numeric(from)) {
    from <- dimnames(rel)[[2]][from]
  }
  names(dimnames(m)) <- c(to, from)
  return(m)
}

#' toolExpandRel
#'
#' Expands aggregation matrix to full items for matrices that only
#' contain subdim items in from/to.
#'
#' Works by first generating all possible from/to items, then sub-setting
#' to only the ones that are actually valid, i.e. from-items that appear in the
#' original data, and to-items that only contain sub-items that also occur in
#' the from-items, so an aggregated item should only contain sub-items that
#' were also present in all from-items it was aggregated from.
#' @param rel The aggregation matrix with subdim items
#' @param x The data to be aggregated already sub-set in case of a partial relation
#' @param dim The (sub-)dim to be aggregated
#' @return An expanded aggregation matrix
#' @keywords internal
toolExpandRel <- function(rel, x, dim) {
  names <- getItems(x, round(floor(dim)))

  # Expand rel matrix to full dimension if rel is only provided for a subdimension
  if (is.null(colnames(rel)) || is.null(rownames(rel))) {
    stop("colnames and/or rownames missing in relation matrix!")
  }

  if (round(dim) == dim || suppressWarnings(all(colnames(rel) == names))) {
    return(rel)
  }

  subdim <- round((dim - floor(dim)) * 10)
  fullItems <- strsplit(getItems(x, floor(dim)), ".", fixed = TRUE)
  onlyNames <- getItems(x, dim)

  if (!setequal(colnames(rel), onlyNames)) {
    if (!setequal(rownames(rel), onlyNames)) {
      stop("The provided mapping contains entries which could not be found in the data: ",
           paste(setdiff(colnames(rel), onlyNames), collapse = ", "))
    }
    rel <- Matrix::t(rel)
  }

  # Build mapping from subdimension values to their aggregated values
  subdimFrom <- colnames(rel)
  subdimToList <- lapply(subdimFrom, function(from) {
    rownames(rel)[rel[, from] == 1]
  })
  names(subdimToList) <- subdimFrom

  # Generate all valid target items
  toItems <- unique(unlist(lapply(fullItems, function(fullItem) {
    toNames <- subdimToList[[fullItem[subdim]]]
    lapply(toNames, function(toName) {
      fullItem[subdim] <- toName
      paste0(fullItem, collapse = ".")
    })
  })))

  # Order toItems to preserve the order of rownames in rel
  # Extract the subdimension part from toItems and order by their position in rownames(rel)
  toItemsSubdim <- lapply(strsplit(toItems, ".", fixed = TRUE), function(x) x[subdim])
  toItems <- toItems[order(match(toItemsSubdim, rownames(rel)))]

  #
  # Actual Expansion to all possible items
  #
  templates <- unique(lapply(fullItems, function(fullItem) {
    fullItem[subdim] <- ""
    fullItem
  }))
  .expandTemplates <- function(fill) {
    unlist(lapply(templates, function(template) {
      lapply(fill, function(item) {
        template[subdim] <- item
        paste0(template, collapse = ".")
      })
    }))
  }
  cnames <- .expandTemplates(colnames(rel))
  rnames <- .expandTemplates(rownames(rel))

  newRel <- Matrix::Matrix(0, nrow = length(rnames), ncol = length(cnames), dimnames = list(rnames, cnames))

  for (i in seq_along(templates)) {
    newRel[seq_len(nrow(rel)) + (i - 1) * nrow(rel), seq_len(ncol(rel)) + (i - 1) * ncol(rel)] <- rel
  }

  # Only return the rows and columns that contain valid
  # target or source items
  return(newRel[toItems, names, drop = FALSE])
}

toolMapFromRel <- function(rel, from, to) {
  if (is.data.frame(rel)) {
    map <- rel[, c(from, to)]
  } else {
    map <- do.call(rbind, lapply(colnames(rel), function(i) {
      return(data.frame(coarse = i, fine = names(which(rel[, i] == 1))))
    }))
    colnames(map) <- c(from, to)
  }
  return(map)
}

removeEmptyResultItems <- function(m, dim) {
  dim <- floor(dim)
  if (any("" %in% unlist(getItems(m, dim = dim)))) {
    vcat(1, "Aggregation target included \"\". Those items were removed from aggregation result.")

    if (any("" %in% unlist(getItems(m, dim = dim)))) {
      m <- m["", dim = dim, invert = TRUE]
    }

    return(m)
  } else {
    return(m)
  }
}