R/CalculationFunctions.R

In USEPA/useeior: USEEIO R modeling software

# Model Calculations

#' Calculate total emissions/resources (LCI) and total impacts (LCIA) for an EEIO model
#' for a given perspective and demand vector.
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param perspective Perspective of the model, can be "DIRECT" or "FINAL". "DIRECT" perspective
#' aligns results with the sectors in which they are produced, while "FINAL" perspective aligns
#' results with the sectors consumed by the final user.
#' @param demand A demand vector, can be name of a built-in model demand vector, e.g. "Production" or "Consumption",
#' or an actual demand vector with names as one or more model sectors and
#' numeric values in USD with the same dollar year as model.
#' @param location, str optional location code for demand vector, required for two-region models
#' @param use_domestic_requirements A logical value: if TRUE, use domestic demand and L_d matrix;
#' if FALSE, use complete demand and L matrix.
#' @param household_emissions, bool, if TRUE, include calculation of emissions from households
#' @param show_RoW, bool, if TRUE, include rows for commodities in RoW, e.g. `111CA/RoW` in result objects.
#' Only valid currently for models with ExternalImportFactors.
#' @export
#' @return A list with LCI and LCIA results (in data.frame format) of the EEIO model.
calculateEEIOModel <- function(model, perspective, demand = "Production", location = NULL,
                               use_domestic_requirements = FALSE, household_emissions = FALSE, show_RoW = FALSE) {
  if (!is.null(model$specs$ExternalImportFactors) && model$specs$ExternalImportFactors) {
    result <- calculateResultsWithExternalFactors(model, perspective, demand, location = location,
                                                  use_domestic_requirements = use_domestic_requirements,
                                                  household_emissions = household_emissions, show_RoW = show_RoW)
  } else {
    # Standard model results calculation
    result <- calculateStandardResults(model, perspective, demand, use_domestic_requirements, location, household_emissions)
  }
  
  logging::loginfo("Result calculation complete.")
  return(result)
}

#' Prepare demand vector for EEIO model results calculations
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param demand A demand vector, can be name of a built-in model demand vector, e.g. "Production" or "Consumption",
#' or an actual demand vector with names as one or more model sectors and
#' numeric values in USD with the same dollar year as model.
#' @param location, str optional location code for demand vector, required for two-region models
#' @param use_domestic_requirements A logical value: if TRUE, use domestic demand and L_d matrix;
#' if FALSE, use complete demand and L matrix.
prepareDemandVectorForStandardResults <- function(model, demand = "Production",
                                                  location = NULL, use_domestic_requirements = FALSE) {
  if (is.character(demand)) {
    #assume this is a model build-in demand 
    #try to load the model vector
    meta <- model$DemandVectors$meta
    if (demand %in% meta$Type) {
      demand_name <- ifelse(use_domestic_requirements,
                            paste0("Domestic", demand),
                            paste0("Complete", demand))
      # Get vector name (ID) from the meta table
      if (is.null(location)) {
        id <- meta[which(meta$Name==demand_name),"ID"]
        if (length(id)>1) {
          stop("Unique demand vector not found, consider passing location")
        }
      } else {
        id <- meta[which(meta$Name==demand_name &
                         meta$Location==location),"ID"]
      }
      d <- model$DemandVectors$vectors[[id]]
    } else {
      stop(paste0("'", demand, "' is not a valid demand vector name in model."))
    }
  } else {
    # Assume this is a user-defined demand vector
    #! Need to check that the given demand 
    if (isDemandVectorValid(demand, model$L)) {
      d <- formatDemandVector(demand, model$L)
    } else {
      stop("Format of the demand vector is invalid. Cannot calculate result.")
    }
  }
  # Convert demand vector into a matrix
  f <- as.matrix(d)  
  
  return(f)
  
}


#' Prepare demand vector for EEIO model results calculations
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param demand A demand vector, can be name of a built-in model demand vector, e.g. "Production" or "Consumption",
#' @param location, str optional location code for demand vector, required for two-region models
prepareDemandVectorForImportResults <- function(model, demand = "Production", location = NULL) {
  if (is.character(demand)) {
    # assume this is a built-in demand 
    if(is.null(location)) {
      location <- "US"
    }
    # Calculate import demand vector y_m. 
    if(demand == "Production") {
      # This option left in for validation purposes.
      logging::loginfo(paste0("Warning: Production demand vector not recommended ",
                              "for estimating results for models with external ",
                              "Import Factors."))
      y_m <- prepareImportProductionDemand(model, location = location)
    } else if(demand == "Consumption") {
      y_m <- prepareImportConsumptionDemand(model, location = location)
    }
  } else {
    # Assume this is a user-defined demand vector
    if (isDemandVectorValid(demand, model$L)) {
      y_m <- formatDemandVector(demand, model$L)
    } else {
      stop("Format of the demand vector is invalid. Cannot calculate result.")
    }
  }
  return(as.matrix(y_m))

}

#' Calculate total emissions/resources (LCI) and total impacts (LCIA) for an EEIO model that has external import factors
#' for a given demand vector.
#' Note that for this calculation, perspective is always FINAL
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param perspective Perspective of the model, can be "DIRECT" or "FINAL". "DIRECT" perspective
#' aligns results with the sectors in which they are produced, while "FINAL" perspective aligns
#' results with the sectors consumed by the final user.
#' @param demand A demand vector, can be name of a built-in model demand vector, e.g. "Production" or "Consumption"
#' @param location, str optional location code for demand vector, required for two-region models
#' @param use_domestic_requirements bool, if TRUE, return only domestic portion of results
#' @param household_emissions, bool, if TRUE, include calculation of emissions from households
#' @param show_RoW, bool, if TRUE, include rows for commodities in RoW, e.g. `111CA/RoW` in result objects.
#' @return A list with LCI and LCIA results (in data.frame format) of the EEIO model.
calculateResultsWithExternalFactors <- function(model, perspective = "FINAL", demand = "Consumption", location = NULL,
                                                use_domestic_requirements = FALSE, household_emissions = FALSE,
                                                show_RoW = FALSE) {
  result <- list()
  y_d <- prepareDemandVectorForStandardResults(model, demand, location = location, use_domestic_requirements = TRUE)
  y_m <- prepareDemandVectorForImportResults(model, demand, location = location)

  if(show_RoW) {
    if(model$specs$IODataSource=="stateior") {
      sector_count <- nrow(y_d)/2
      row_names <- c(colnames(model$M_m),
                     gsub("/.*", "/RoW", colnames(model$M_m[, 1:sector_count])))
    } else {
      row_names <- c(colnames(model$M_m),
                     gsub("/.*", "/RoW", colnames(model$M_m)))
    }
  } else {
    row_names <- colnames(model$M_m)
  }

  ## Description of result components apply to both FINAL and DIRECT perspectives
  # see equations 4-7 and 9-12 in EPA 600/R-24/116
  # G_d - Domestic emissions from domestic production
  # G_mi - Emissions from imported goods consumed as intermediate products by domestic industries
  # G_mf - Emissions from imported goods consumed as final products

  if(perspective == "FINAL") {
    # Calculate Final Perspective LCI (a matrix with total impacts in form of sector x flows)
    logging::loginfo("Calculating Final Perspective LCI and LCIA with external import factors...")
    subscript <- "l"
    G_d <- model$B %*% model$L_d %*% diag(as.vector(y_d))
    G_mi <- model$M_m %*% model$A_m %*% model$L_d %*% diag(as.vector(y_d))

  } else { # Calculate direct perspective results.
    # Calculate Direct Perspective LCI (a matrix with total impacts in form of sector x flows)
    logging::loginfo("Calculating Direct + Imported Perspective LCI and LCIA with external import factors...")
    subscript <- "r"
    s <- getScalingVector(model$L_d, y_d)
    G_d <- t(calculateDirectPerspectiveLCI(model$B, s))
    G_mi <- t(calculateDirectPerspectiveLCI(model$M_m, (model$A_m %*% model$L_d %*% y_d)))
  }
  G_mf <- model$M_m %*% diag(as.vector(y_m))

  if (use_domestic_requirements) {
    # zero out the import results
    G_mi[] <- 0
    G_mf[] <- 0
  }
    
  if(show_RoW) {
    if(model$specs$IODataSource=="stateior") {
      # collapse third term for SoI and RoUS
      G_mf <- G_mf[, 1:sector_count] + G_mf[, (sector_count+1):(sector_count*2)]
      
      if(perspective == "DIRECT") {
        # collapse second and third term for SoI and RoUS
        G_mi <- G_mi[, 1:sector_count] + G_mi[, (sector_count+1):(sector_count*2)]
      }
    }
    if(perspective == "DIRECT") {
      LCI <- cbind(G_d, G_mi + G_mf) # Term 2 and Term 3 are assigned to RoW
    } else {
      LCI <- cbind(G_d + G_mi, G_mf) # Term 3 is assigned to RoW
    }
  } else {
    LCI <- G_d + G_mi + G_mf # All three terms combined and regions do not change 
  }
    
  # Calculate LCIA (matrix with direct impacts in form of sector x impacts)
  LCIA <- model$C %*% LCI
  LCI <- t(LCI)
  LCIA <- t(LCIA)
    
  colnames(LCI) <- rownames(model$M_m)
  rownames(LCI) <- row_names
  colnames(LCIA) <- rownames(model$D)
  rownames(LCIA) <- row_names
    
  # Add household emissions to results if applicable
  if(household_emissions) {
    hh <- calculateHouseholdEmissions(model, f=(y_d + y_m), location, characterized=FALSE, show_RoW=show_RoW)
    hh_lcia <- calculateHouseholdEmissions(model, f=(y_d + y_m), location, characterized=TRUE, show_RoW=show_RoW)
    LCI <- rbind(LCI, hh)
    LCIA <- rbind(LCIA, hh_lcia)
  }
  result[[paste0("G_", subscript)]] <- LCI
  result[[paste0("H_", subscript)]] <- LCIA
  return(result)
  
}


#' Calculate total emissions/resources (LCI) and total impacts (LCIA) for an EEIO model
#' that does not have external import factors for a given perspective and demand vector.
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param perspective Perspective of the model, can be "DIRECT" or "FINAL". "DIRECT" perspective
#' aligns results with the sectors in which they are produced, while "FINAL" perspective aligns
#' results with the sectors consumed by the final user.
#' @param demand A demand vector, can be name of a built-in model demand vector, e.g. "Production" or "Consumption",
#' or an actual demand vector with names as one or more model sectors and
#' numeric values in USD with the same dollar year as model.
#' @param use_domestic_requirements A logical value: if TRUE, use domestic demand and L_d matrix;
#' if FALSE, use complete demand and L matrix.
#' @param location, str optional location code for demand vector, required for two-region models
#' @param household_emissions, bool, if TRUE, include calculation of emissions from households
#' @return A list with LCI and LCIA results (in data.frame format) of the EEIO model.
calculateStandardResults <- function(model, perspective, demand, use_domestic_requirements = FALSE,
                                     location = NULL, household_emissions = FALSE) {
  f <- prepareDemandVectorForStandardResults(model, demand, location, use_domestic_requirements)
  # Initialize results list
  result <- list() 
  # Generate Total Requirements (L or L_d) matrix based on whether "use_domestic"
  if (use_domestic_requirements) {
    L <- model$L_d
    M <- model$M_d
    N <- model$N_d
  } else {
    L <- model$L
    M <- model$M
    N <- model$N
  }
  if(household_emissions) {
    hh <- calculateHouseholdEmissions(model, f, location, characterized=FALSE)
    hh_lcia <- calculateHouseholdEmissions(model, f, location, characterized=TRUE)
  }
  # Calculate LCI and LCIA in direct or final perspective
  if (perspective=="DIRECT") {
    # Calculate Direct Perspective LCI (a matrix with direct impacts in form of sector x flows)
    logging::loginfo("Calculating Direct Perspective LCI...")
    s <- getScalingVector(L, f)
    result$G_r <- calculateDirectPerspectiveLCI(model$B, s)
    # Calculate Direct Perspective LCIA (matrix with direct impacts in form of sector x impacts)
    logging::loginfo("Calculating Direct Perspective LCIA...")
    result$H_r <- calculateDirectPerspectiveLCIA(model$D, s)
    if(household_emissions) {
      result$G_r <- rbind(result$G_r, hh)
      result$H_r <- rbind(result$H_r, hh_lcia)
    }
  } else if (perspective=="FINAL") {
    # Calculate Final Perspective LCI (a matrix with total impacts in form of sector x flows)
    logging::loginfo("Calculating Final Perspective LCI...")
    result$G_l <- calculateFinalPerspectiveLCI(M, f)
    # Calculate Final Perspective LCIA (matrix with total impacts in form of sector x impacts)
    logging::loginfo("Calculating Final Perspective LCIA...")
    result$H_l <- calculateFinalPerspectiveLCIA(N, f)
    if(household_emissions) {
      result$G_l <- rbind(result$G_l, hh)
      result$H_l <- rbind(result$H_l, hh_lcia)
    }
  }
  
  return(result)
}


#' Multiply the Leontief inverse L and the demand vector to calculate scaling vector
#' that represents production needed to fulfill the demand.
#' @param L Leontief inverse.
#' @param demand Final demand vector.
#' @return Scaling vector.
#' @references Yang, Yi, Wesley W. Ingwersen, Troy R. Hawkins, Michael Srocka, and David E. Meyer.
#' 2017. “USEEIO: A New and Transparent United States Environmentally-Extended Input-Output Model.”
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 8.
getScalingVector <- function(L, demand) {
  s <- L %*% demand
  return(s)
}

#' The direct perspective LCI aligns flows with sectors consumed by direct use.
#' Multiply the B matrix and the scaling vector s.
#' @param B Marginal impact per unit of the environmental flows.
#' @param s Scaling vector.
#' @return A matrix with direct impacts in form of sector x flows.
#' @references Yang, Yi, Wesley W. Ingwersen, Troy R. Hawkins, Michael Srocka, and David E. Meyer.
#' 2017. “USEEIO: A New and Transparent United States Environmentally-Extended Input-Output Model.”
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 9.
calculateDirectPerspectiveLCI <- function(B, s) {
  G_r <- t(B %*% diag(as.vector(s), nrow(s)))
  rownames(G_r) <- rownames(s)
  return(G_r)
}

#' The final perspective LCI aligns flows with sectors consumed by final users.
#' Multiply the M matrix and the diagonal of demand, y.
#' @param M, a model M matrix, direct + indirect flows per $ output of sector.
#' @param y, a model demand vector
#' @return A matrix with total impacts in form of sectors x flows..
#' @references Yang, Yi, Wesley W. Ingwersen, Troy R. Hawkins, Michael Srocka, and David E. Meyer.
#' 2017. “USEEIO: A New and Transparent United States Environmentally-Extended Input-Output Model.”
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 10.
calculateFinalPerspectiveLCI <- function(M, y) {
  G_l <- t(M %*% diag(as.vector(y)))
  colnames(G_l) <- rownames(M)
  rownames(G_l) <- colnames(M)
  return(G_l)
}

#' The direct perspective LCIA aligns impacts with sectors consumed by direct use.
#' Multiply the D matrix (the product of C matrix and B matrix) and scaling vector s.
#' @param D Direct impact per unit of the environmental flows.
#' @param s Scaling vector.
#' @return A matrix with direct impacts in form of sector x impact categories.
#' @references Yang, Yi, Wesley W. Ingwersen, Troy R. Hawkins, Michael Srocka, and David E. Meyer.
#' 2017. “USEEIO: A New and Transparent United States Environmentally-Extended Input-Output Model.”
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 9.
calculateDirectPerspectiveLCIA <- function(D, s) {
  H_r <- t(D %*% diag(as.vector(s), nrow(s)))
  rownames(H_r) <- rownames(s)
  return(H_r)
}

#' The final perspective LCIA aligns impacts with sectors consumed by final users.
#' Multiply the N matrix and the diagonal of demand, y.
#' @param N, a model N matrix, direct + indirect impact per unit of the environmental flows.
#' @param y, a model demand vector
#' @return A matrix with total impacts in form of sector x impact categories.
#' @references Yang, Yi, Wesley W. Ingwersen, Troy R. Hawkins, Michael Srocka, and David E. Meyer.
#' 2017. “USEEIO: A New and Transparent United States Environmentally-Extended Input-Output Model.”
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 10.
calculateFinalPerspectiveLCIA <- function(N, y) {
  H_l <- t(N %*% diag(as.vector(y)))
  colnames(H_l) <- rownames(N)
  rownames(H_l) <- colnames(N)
  return(H_l)
}

#' Divide/Normalize a sector x flows matrix by the total of respective flow (column sum)
#' @param m A sector x flows matrix.
#' @return A normalized sector x flows matrix.
#' @export
normalizeResultMatrixByTotalImpacts <- function(m) {
  #Use sweep function to prevent error
  m_norm <- sweep(m, MARGIN = 2, FUN = "/", STATS = colSums(m))
  return(m_norm)
}

#' Dot multiplies two vectors to calculate an impact score and the percent contribution each score to the total
#' @param x, numeric vector of length n
#' @param y, numeric vector of length n
#' @return A dataframe sorted by contribution (high-to-low), also showing "x","y","impact" 
calculatePercentContributiontoImpact <- function (x,y) {
  df <- cbind.data.frame(x,y)
  df["impact"] <- df[,"x"]*df[,"y"] 
  df["contribution"] <- df["impact"]/sum(df["impact"])
  df <- df[order(df$contribution,decreasing=TRUE),]
  return(df) 
}

#' Calculate the percent contribution of sectors to direct+indirect impacts by an indicator,
#' using the product of model L matrix (total requirements) and D matrix (direct impacts by indicator).
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param sector, str, index of a model sector for use in the L matrix, e.g. "221100/US".
#' @param indicator, str, index of a model indicator for use in the D matrix, e.g. "Acidification Potential".
#' @param domestic, boolean, sets model to use domestic flow matrix. Default is FALSE.
#' @return A dataframe sorted by contribution (high-to-low), also showing "L", "D", "impact".
## @export
calculateSectorContributiontoImpact <- function (model, sector, indicator, domestic=FALSE) {
  L <- model$L
  if (domestic) {
    L <- model$L_d
  }
  D <- model$D
  df <- calculatePercentContributiontoImpact(L[,sector],D[indicator,])
  # Rename x and y cols
  colnames(df)[colnames(df)==c("x", "y")] <- c("L", "D")
  # Add sector name for easier interpretation of results 
  rownames(df) <- paste(rownames(df),
                        model$Commodities[match(rownames(df),
                                                model$Commodities$Code_Loc), "Name"],
                        sep = " - ")
  return(df)
}

#' Calculate the percent contribution of flows to direct+indirect impacts by an indicator,
#' using model M matrix for total impacts of flows and C matrix for indicator.
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param sector, str, index of a model sector for use in the M matrix, e.g. "221100/US".
#' @param indicator, str, index of a model indicator for use in the C matrix, e.g. "Acidification Potential".
#' @param domestic, boolean, sets model to use domestic flow matrix. Default is FALSE.
#' @return A dataframe sorted by contribution (high-to-low), also showing "M", "C", "impact".
#' @export 
calculateFlowContributiontoImpact <- function (model, sector, indicator, domestic=FALSE) {
  M <- model$M
  C <- model$C
  if (domestic) {
    M <- model$M_d
  }
  df <- calculatePercentContributiontoImpact(M[,sector], C[indicator,])
  # Rename x and y cols
  colnames(df)[colnames(df)==c("x", "y")] <- c("M", "C")
  return(df)
}

#' Aggregate result matrix by rows
#'
#' @param matrix      A result calculation matrix (e.g. `G` or `H`) with sectors as rows
#' @param to_level    The level of BEA code this matrix will be aggregated to
#' @param crosswalk   Sector crosswalk between levels of detail
#' @return An aggregated result calculation matrix with sectors as rows
aggregateResultMatrixbyRow <- function (matrix, to_level, crosswalk) {
  # Determine the columns within MasterCrosswalk that will be used in aggregation
  from_code <- "USEEIO"
  to_code <- paste0("BEA_", to_level)
  # Remove location codes from sectors
  df <- as.data.frame(matrix)
  rownames(df) <- gsub("/.*","",rownames(df))
  matrix <- as.matrix(df)
  # Aggregate by rows
  value_columns <- colnames(matrix)
  df_fromlevel <- merge(matrix, unique(crosswalk[, c(from_code, to_code)]), by.x = 0, by.y = from_code)
  df_fromlevel_agg <- stats::aggregate(df_fromlevel[, value_columns], by = list(df_fromlevel[, to_code]), sum)
  rownames(df_fromlevel_agg) <- df_fromlevel_agg[, 1]
  df_fromlevel_agg[, 1] <- NULL
  matrix_fromlevel_agg <- as.matrix(df_fromlevel_agg)
  return(matrix_fromlevel_agg)
}

#' Aggregate result matrix by rows and columns
#'
#' @param matrix      A matrix of inventory or impact results with sectors as both rows and columns,
#' such as the output from calculateSectorPurchasedbySectorSourcedImpact()
#' @param to_level    The level of BEA code this matrix will be aggregated to
#' @param crosswalk   Sector crosswalk between levels of detail
#' @return An aggregated matrix with sectors as rows and columns
#' @export
aggregateResultMatrix <- function (matrix, to_level, crosswalk) {
  row_agg_matrix <- aggregateResultMatrixbyRow (matrix, to_level, crosswalk)
  col_agg_matrix <- aggregateResultMatrixbyRow (t(row_agg_matrix), to_level, crosswalk)
  agg_matrix <- t(col_agg_matrix)
  return(agg_matrix)
}

#' Calculates sector x sector impacts from a given demand vector.
#' @param y a model demand vector or user-specified demand vector
#' that must have the same dimension with the model demand vector
#' @param model A complete EEIO model: a list with USEEIO model components and attributes
#' @param indicator str, index of a model indicator, e.g. "Acidification Potential"
#' @return A matrix of impacts in the form of sector purchased x sector sourced,
#' negative values should be interpreted as "reduced impacts".
#' @export
calculateSectorPurchasedbySectorSourcedImpact <- function (y, model, indicator) {
  L <- model$L
  total_req <- L %*% diag(y)
  colnames(total_req) <- rownames(total_req)
  impacts <- t(diag(model$D[indicator,]) %*% total_req)
  colnames(impacts) <- rownames(impacts)
  return(impacts)
}

#' Calculate sector margin impacts in the form of M and N Matrix
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @return A list with M_margin and N_margin
#' @export
calculateMarginSectorImpacts <- function(model) {
  if (model$specs$IODataSource == "stateior") {
    stop("Margins not available for two-region models")
  }
  # Calculation fractions of producer price for each margin
  MarginCoefficients <- as.matrix(model$Margins[, c("Transportation", "Wholesale", "Retail")] / 
                                    model$Margins[, c("ProducersValue")])
  rownames(MarginCoefficients) <- model$Margins$SectorCode
  MarginCoefficients[is.na(MarginCoefficients)] <- 0
  
  # Create margin_allocation matrix to allocate fractions by margin sector
  # In the matrix, rows are three margin types and columns are margin sectors
  margin_allocation <- matrix(nrow = 3, ncol = nrow(model$MarginSectors), 0)
  rownames(margin_allocation) <- colnames(MarginCoefficients)
  colnames(margin_allocation) <- model$MarginSectors$Code
  # Assign allocation factors to margin sectors based on total Commodity output
  output_ratio <- calculateOutputRatio(model, output_type="Commodity")
  for (i in rownames(margin_allocation)) {
    codes <- model$MarginSectors[model$MarginSectors$Name==i, "Code"]
    margin_allocation[i, codes] <- output_ratio[match(codes, output_ratio$SectorCode), "toSectorRatio"]
  }
  # Multiply fractions by allocation matrix to get a fraction per margin sector for each commodity
  margins_by_sector <- MarginCoefficients %*% margin_allocation
  
  # Put margins_by_sector into a matrix in the form of A
  A_margin <- model$A
  # Make sure sector ordering is the same
  A_margin[,] <- 0 
  for (i in model$MarginSectors$Code_Loc) {
    A_margin[i, ] <- margins_by_sector[gsub("/.*", "", colnames(A_margin)),
                                       gsub("/.*", "", i)]
  }
  # Multiply M and N by margins_by_sector to derive M_margin and N_margin
  ls <- list("M_margin" = model$M %*% A_margin,
             "N_margin" = model$N %*% A_margin)
  logging::loginfo("Model margin impacts calculated.")
  return(ls)
}


#' Calculate household emissions from B_h
#' @param model A complete EEIO model: a list with USEEIO model components and attributes.
#' @param f A demand vector with names as one or more model sectors and
#' numeric values in USD with the same dollar year as model.
#' @param location, str optional location code for demand vector, required for two-region models
#' @param characterized, bool, TRUE to characterize using C matrix, FALSE to show LCI
#' @param show_RoW, bool, if TRUE, include rows for commodities in RoW, e.g. `111CA/RoW` in result objects.
#' Only valid currently for models with ExternalImportFactors.
#' @return A result vector with rows for final demand sector(s)
calculateHouseholdEmissions <- function(model, f, location, characterized=FALSE, show_RoW=FALSE) {
  if(!"B_h" %in% names(model)) {
    logging::logwarn("Household emissions not found in this model")
    return(NULL)
  }
  if(length(model$specs$ModelRegionAcronyms) == 1) {
    # Set location as NULL for single region model
    location <- NULL
  }
  codes <- model$FinalDemandMeta[model$FinalDemandMeta$Group%in%c("Household"), "Code_Loc"]
  if (!is.null(location)) {
    other_code <- codes[!grepl(location, codes)]
    codes <- codes[grepl(location, codes)]
  }

  if(characterized) {
    hh = t(model$C %*% as.matrix(model$B_h[, codes])) * colSums(f)
  } else {
    hh = t(as.matrix(model$B_h[, codes])) * colSums(f)
  }
  rownames(hh) <- codes

  # Create a matrix of 0 values for potential addition to household emissions matrix
  mat <- matrix(0, nrow=nrow(hh), ncol=ncol(hh))
  
  if(!is.null(location)) {
    # add in 0 values for RoUS
    rownames(mat) <- other_code
    hh <- rbind(hh, mat)
  }
  if(show_RoW) {
    # add in 0 values for RoW
    rownames(mat) <- gsub("/.*", "/RoW", codes)
    hh <- rbind(hh, mat)
  }
  return(hh)
}


#' For a given impact, provided via indicator or elementary flow label, 
#' disaggregate the total impacts per purchase (indicator: N, flow: M) into 
#' direct impacts (indicator: D, flow: B) and upstream, Tier 1 purchase impacts. 
#' Return a long-format df of exchanges, with sector names mapped to sector codes.
#' @param model A complete EEIO model: a list with USEEIO model components and attributes
#' @param impact str, a model indicator (e.g., "Greenhouse Gases") row index of N, 
#'                    or elementary flow (e.g., "Methane/emission/air/kg") index of M
#' @param opt_impact str {'indicator', 'elemflow'}, string code to specify impact type
#' @export
#' @return A data frame of direct and per-tier-1-purchase sector impacts
disaggregateTotalToDirectAndTier1 <- function(model, impact, opt_impact="indicator") {
  mtx_direct <- c("indicator"="D", "elemflow"="B")[opt_impact]
  if (is.na(mtx_direct)) {
    stop(paste0("'",opt_impact,"' is not a valid opt_impact string code"))
  }
  # get direct sector impacts
  df_direct <- tryCatch({  # catches bad `impact` row label
      tibble::enframe(model[[mtx_direct]][impact,])
  }, error=function(e) {
      stop(paste0("'",impact,"' is not a valid ",opt_impact," label"))
  })
  df_direct <- dplyr::rename(df_direct, impact_per_purchase=value, sector_code=name) 
  # assign "Direct" as purchased commodity label for data-vis & stat convenience
  df_direct <- dplyr::mutate(df_direct, purchased_commodity = 'Direct')
  # get total impacts per Tier 1 purchase by sector
  df_total <- calculateTotalImpactbyTier1Purchases(model, impact, opt_impact) 
  df_total <- tibble::as_tibble(df_total, rownames="purchased_commodity_code") 
  df_total <- reshape2::melt(df_total, id.vars="purchased_commodity_code",
                                       variable.name="sector_code",
                                       value.name="impact_per_purchase") 
  # map sector codes to names
  sector_map <- setNames(model$Commodities$Name, model$Commodities$Code_Loc)
  df_total <- dplyr::mutate(df_total, 
    purchased_commodity = dplyr::recode(purchased_commodity_code, !!!sector_map))
  # concat direct + total impacts  
  df_impacts <- dplyr::bind_rows(df_total, df_direct) 
  df_impacts <- dplyr::mutate(df_impacts, 
    sector = dplyr::recode(sector_code, !!!sector_map))
  return(df_impacts)
}

#' Calculate sector x sector total impacts (single indicator or elementary flow) 
#' for Tier 1 purchases. Multiply each row of the sector by sector A matrix by 
#' the scalar elements of a single-impact by sector array (indicator: N, flow: M)
#' @param model A complete EEIO model: a list with USEEIO model components and attributes
#' @param impact str, a model indicator (e.g., "Greenhouse Gases") row index of N, 
#'                    or elementary flow (e.g., "Methane/emission/air/kg") index of M
#' @param opt_impact str {'indicator', 'elemflow'}, string code to specify impact type
#' @return A sector by sector, impact-per-tier-1-purchase matrix.
calculateTotalImpactbyTier1Purchases <- function(model, impact, opt_impact='indicator') {
  mtx_total <- c("indicator"="N", "elemflow"="M")[opt_impact]
  if (is.na(mtx_total)) {
    stop(paste0("'",opt_impact,"' is not a valid opt_impact string code"))
  }
  df_direct <- tryCatch({  # catches bad `impact` row label
    totalImpactPerPurchase <- model[[mtx_total]][impact,] * model$A
    # totalImpactPerPurchase <- model$N[impact,] * model[["A"]]
  }, error=function(e) {
    stop(paste0("'",impact,"' is not a valid ",opt_impact," label"))
  })
  return(totalImpactPerPurchase)
}


#' Function to get backward or forward economic linkages from a USEEIO model,
#' Provides results based on input demand, sorted by total (direct+indirect)
#' Backward linkages use A and L matrices
#' Forward linkages use Ghosh counterparts to A and L
#' @param model An EEIO model object with model specs, IO tables, satellite tables,
#'    and indicators loaded
#' @param demand A demand vector, can be name of a built-in model demand vector,
#'    e.g. "Production" or "Consumption",
#'    or an actual demand vector with names as one or more model sectors and
#'    numeric values in USD with the same dollar year as model.
#' @param type "backward" linkages use A and L matrices "forward" linkages use Ghosh
#'    counterparts to A and L
#' @param use_domestic_requirements A logical value: if TRUE, use domestic demand
#'    and L_d matrix; if FALSE, use complete demand and L matrix.
#' @param location str optional location code for demand vector, required for
#'    two-region models
#' @param cutoff numeric value, shows the cutoff value for sorted results.
#'    Values smaller than cutoff are not shown
#' @return dataframe of direct, indirect, and total economic linkages
#' @export
getSectorLinkages <- function(model, demand, type = "backward",
                              location = NULL, use_domestic_requirements = FALSE,
                              cutoff = 0.01) {
  f <- prepareDemandVectorForStandardResults(model, demand, location,
                                             use_domestic_requirements)
  if(type=="backward") {
    if (use_domestic_requirements) {
      tier1 <- model$A_d %*%f 
      all <- model$L_d %*%f  
    } else {
      tier1 <- model$A%*%f 
      all <- model$L%*%f   
    }
  } else if(type=="forward") {
    Ghosh_B <- calculateGhoshB(model, use_domestic_requirements)
    Ghosh_G <- calculateGhoshG(Ghosh_B)
    tier1 <- Ghosh_B%*%f 
    all <- Ghosh_G%*%f  
  }
  tier2andbeyond <- all-tier1
  
  # Remove the additional "1" from the Leontief or Ghosh to avoid double counting indirect 
  tier2andbeyond <- tier2andbeyond-f
  total <- tier1 + tier2andbeyond
  
  linkages <- data.frame(direct=tier1, indirect=tier2andbeyond, total=total)
  linkages$Name <- model$Commodities[match(rownames(linkages),
                                           model$Commodities$Code_Loc),"Name"]
  # Filter by share > 1%
  linkages <- linkages[linkages$total > cutoff, ]
  # Sort them
  linkages <- linkages[order(linkages$total, decreasing = TRUE), ]
  return(linkages)
}

#' Function to get supply side Ghosh matrix B. B can be calculated based on commodity-by-commodity or 
#' industry-by-industry direct requirement matrix: (Miller and Blair)
#' for industry-by-industry requirements: B = inverse(qhat) * A(industry-by-industry)*qhat,
#' where q is the total commodity output
#' for commodity-by-commodity requirements: B = inverse(xhat) * A(industry-by-industry)*xhat,
#' where x is the total inudstry output
#' @param model An EEIO model object with model specs, IO tables, satellite tables, and indicators loaded
#' @param use_domestic_requirements A logical value: if TRUE, use domestic demand and L_d matrix;
#'    if FALSE, use complete demand and L matrix.
#' @return B, Ghosh B matrix
calculateGhoshB <- function(model, use_domestic_requirements = FALSE) {
  q <- model$q # total commodity output
  x <- model$x # total industry output
  if(use_domestic_requirements) {
    A <- model$A_d
  } else {
    A <- model$A    
  }
  if(model$specs$CommodityorIndustryType == "Commodity") {
    B <- solve(diag(q)) %*% A %*% diag(q)
    row.names(B) <- model$Commodities$Code_Loc
    colnames(B) <- model$Commodities$Code_Loc
  } else if(model$specs$CommodityorIndustryType == "Industry") {
    B <- solve(diag(x)) %*% A %*% diag(x)
    row.names(B) <- model$Industries$Code_Loc
    colnames(B) <- model$Industries$Code_Loc
  }

  return(B)
}

#' Function to calculate inverse of Ghosh matrix B, returns the inverse G
#' @param GhoshB Ghosh matrix B
#' @return G, Ghosh matrix G
calculateGhoshG <- function(GhoshB) {
  G <- useeior:::calculateLeontiefInverse(GhoshB)
  return(G)
}