Merge branch 'linkstateiorImportF' into imports_2region

USEPA · Jan 12, 2024 · 7c2ed04 · 7c2ed04
2 parents d89737d + 37d0313
commit 7c2ed04
Show file tree

Hide file tree

Showing 9 changed files with 82 additions and 1,064 deletions.
diff --git a/R/BuildModel.R b/R/BuildModel.R
@@ -137,9 +137,9 @@ constructEEIOMatrices <- function(model, configpaths = NULL) {
     }
     logging::loginfo("Calculating N_d matrix (total environmental impacts per dollar from domestic activity)...")
     model$N_d <- model$C %*% model$M_d
-    if(!is.null(model$M_m)) {
+    if(!is.null(model$Q_t)) {
       logging::loginfo("Calculating N_m matrix (total environmental impacts per dollar from imported activity)...")
-      model$N_m <- model$C %*% model$M_m # I don't think this is correct. See ImportFactorsTests.rmd, "Calculate the N matrix for the standard model using coupled model approach" chunk.
+      model$N_m <- model$C %*% model$Q_t # I don't think this is correct. See ImportFactorsTests.rmd, "Calculate the N matrix for the standard model using coupled model approach" chunk.
     }
   }
 

diff --git a/R/CalculationFunctions.R b/R/CalculationFunctions.R
@@ -122,15 +122,15 @@ calculateResultsWithExternalFactors <- function(model, demand = "Consumption", l
   if (use_domestic_requirements) {
     result$LCI_f <- (model$B %*% model$L_d %*% y_d)
   } else {
-    result$LCI_f <- (model$B %*% model$L_d %*% y_d) + (model$M_m %*% model$A_m %*% model$L_d %*% y_d + model$M_m %*% y_m)
+    result$LCI_f <- (model$B %*% model$L_d %*% y_d) + (model$Q_t %*% model$A_m %*% model$L_d %*% y_d + model$Q_t %*% y_m)
   }
   result$LCIA_f <- model$C %*% result$LCI_f
 
   result$LCI_f <- t(result$LCI_f)
   result$LCIA_f <- t(result$LCIA_f)
 
-  colnames(result$LCI_f) <- rownames(model$M_m)
-  rownames(result$LCI_f) <- colnames(model$M_m)
+  colnames(result$LCI_f) <- rownames(model$Q_t)
+  rownames(result$LCI_f) <- colnames(model$Q_t)
 
   colnames(result$LCIA_f) <- rownames(model$N_m)
   rownames(result$LCIA_f) <- colnames(model$N_m)

diff --git a/R/IOFunctions.R b/R/IOFunctions.R
@@ -266,14 +266,14 @@ buildModelwithImportFactors <- function(model, configpaths = NULL) {
   logging::loginfo("Calculating M_d matrix (total emissions and resource use per dollar from domestic activity)...")
   model$M_d <- model$B %*% model$L_d 
 
-  logging::loginfo("Calculating M_m matrix (total emissions and resource use per dollar from imported activity)...")
-  M_m <- loadExternalImportFactors(model, configpaths)
+  logging::loginfo("Calculating Q_t matrix (total emissions and resource use per dollar from imported activity)...")
+  Q_t <- loadExternalImportFactors(model, configpaths)
 
-  # Fill in flows for M_m not found in Import Factors but that exist in model and align order
-  M_m <- rbind(M_m, model$M_d[setdiff(rownames(model$M_d), rownames(M_m)),])
-  M_m <- M_m[rownames(model$M_d), ]
+  # Fill in flows for Q_t not found in Import Factors but that exist in model and align order
+  Q_t <- rbind(Q_t, model$M_d[setdiff(rownames(model$M_d), rownames(Q_t)),])
+  Q_t <- Q_t[rownames(model$M_d), ]
 
-  model$M_m <- M_m
+  model$Q_t <- Q_t
 
   return(model)
 }

diff --git a/R/LoadExternalImportFactors.R b/R/LoadExternalImportFactors.R
@@ -4,7 +4,7 @@
 #' @param model An EEIO form USEEIO model object with model specs loaded
 #' @param configpaths str vector, paths (including file name) of model configuration file
 #' and optional agg/disagg configuration file(s). If NULL, built-in config files are used.
-#' @return M_m, matrix of import coefficients (flow x sector).
+#' @return Q_t, matrix of import coefficients (flow x sector).
 loadExternalImportFactors <- function(model, configpaths = NULL) {
 
   # Read in file with Import factors
@@ -50,12 +50,12 @@ loadExternalImportFactors <- function(model, configpaths = NULL) {
     IFTable['Location'] <- "US"
   }
 
-  M_m <- standardizeandcastSatelliteTable(IFTable, model)
+  Q_t <- standardizeandcastSatelliteTable(IFTable, model)
   # standardizeandcast prepares df for industries, convert to commodities
-  M_m[, setdiff(model$Commodities$Code_Loc, colnames(M_m))] <- 0
+  Q_t[, setdiff(model$Commodities$Code_Loc, colnames(Q_t))] <- 0
   # Adjust column order to be the same with V_n rownames
-  M_m <- M_m[, model$Commodities$Code_Loc]
-  M_m <- as.matrix(M_m)
+  Q_t <- Q_t[, model$Commodities$Code_Loc]
+  Q_t <- as.matrix(Q_t)
 
-  return(M_m)
+  return(Q_t)
 }
diff --git a/R/ValidateModel.R b/R/ValidateModel.R
@@ -323,13 +323,13 @@ printValidationResults <- function(model) {
   print(paste("Number of flow totals by commodity failing:",q_x_val$N_Fail))
   print(paste("Sectors with flow totals failing:", paste(unique(q_x_val$Failure$rownames), collapse = ", ")))
 
-if (model$specs$CommodityorIndustryType=="Commodity") {
-  print("Validate that commodity output equals to domestic use plus production demand")
-  q_val <- compareCommodityOutputandDomesticUseplusProductionDemand(model, tolerance = 0.01)
-  print(paste("Number of flow totals by commodity passing:",q_val$N_Pass))
-  print(paste("Number of flow totals by commodity failing:",q_val$N_Fail))
-  print(paste("Sectors with flow totals failing:", paste(unique(q_val$Failure$rownames), collapse = ", ")))
-}
+  if (model$specs$CommodityorIndustryType=="Commodity") {
+    print("Validate that commodity output equals to domestic use plus production demand")
+    q_val <- compareCommodityOutputandDomesticUseplusProductionDemand(model, tolerance = 0.01)
+    print(paste("Number of flow totals by commodity passing:",q_val$N_Pass))
+    print(paste("Number of flow totals by commodity failing:",q_val$N_Fail))
+    print(paste("Sectors with flow totals failing:", paste(unique(q_val$Failure$rownames), collapse = ", ")))
+  }
 }
 
 #' Removes hybrid processes form a model object for successful validation
@@ -354,152 +354,6 @@ removeHybridProcesses <- function(model, object) {
   return(object)
 }
 
-#' Validate Leontief matrix (L) of two-region model and final demand against
-#' SoI and RoUS output.
-#' @param model A complete EEIO model: a list with USEEIO model components and attributes
-#' @param state A text value specifying state of interest.
-#' @return A list of validation components and result.
-validate2RegionLagainstOutput <- function(model, state=NULL) {
-  startLogging()
-  #Get model data
-  iolevel <- model$specs$BaseIOLevel
-  ioschema <- model$specs$BaseIOSchema
-  year <- model$specs$IOYear
-  state_abb <- sub(".*-","",model$specs$ModelRegionAcronyms[1]) ## Extract characters after -
-
-  # Define industries and commodities
-  industries <- getVectorOfCodes(ioschema, iolevel, "Industry")
-  commodities <- getVectorOfCodes(ioschema, iolevel, "Commodity")
-  ita_column <- ifelse(iolevel == "Detail", "F05100", "F051")
-  # Define state abbreviation
-  logging::loginfo(paste0("Generating A matrix of ",state_abb," (SoI) Make table ..."))
-  # SoI Make
-  TwoRegionMake <- as.data.frame(model$V) # will need to modify this when this part is moved to a new buildTwoRegionA matrix function
-  SoI_Make <- TwoRegionMake[endsWith(rownames(TwoRegionMake), state_abb),
-                            endsWith(colnames(TwoRegionMake), state_abb)]
-  # SoI commodity output
-  TwoRegionCommodityOutput <- model$q
-  SoI_Commodity_Output <- TwoRegionCommodityOutput[endsWith(names(TwoRegionCommodityOutput),
-                                                            state_abb)]
-  # SoI A matrix
-  SoI_A <- normalizeIOTransactions(SoI_Make, SoI_Commodity_Output)
-  # Check column sums of SoI_A
-  if (all(abs(colSums(SoI_A) - 1) < 1E-3)) {
-    logging::loginfo("FACT CHECK: column sums of A matrix of SoI Make table == 1.")
-  } else {
-    logging::logwarn("Column sums of A matrix of SoI Make table != 1")
-  }
-
-  logging::loginfo("Generating A matrix of RoUS Make table ...")
-  # RoUS Make
-  RoUS_Make <- TwoRegionMake[endsWith(rownames(TwoRegionMake), "RoUS"),
-                             endsWith(colnames(TwoRegionMake), "RoUS")]
-  # RoUS commodity output
-  RoUS_Commodity_Output <- TwoRegionCommodityOutput[endsWith(names(TwoRegionCommodityOutput),
-                                                             "RoUS")]
-
-  # RoUS A matrix
-  RoUS_A <- normalizeIOTransactions(RoUS_Make, RoUS_Commodity_Output)
-  # Check column sums of RoUS_A
-  if (all(abs(colSums(RoUS_A) - 1) < 1E-3)) {
-    logging::loginfo("FACT CHECK: column sums of A matrix of RoUS Make table == 1.")
-  } else {
-    logging::logerror("Column sums of A matrix of RoUS Make table != 1")
-  }
-
-  # Two-region A matrix
-  #                       Commodity                Commodity                  Industry                     Industry
-  #                          SoI                     RoUS                        SoI                         RoUS
-  #               +-----------------------+-------------------------+---------------------------+----------------------------+
-  # Commodity SoI |           0           |            0            |norm(SoI2SoI_Use,  SoI_TIO)|norm(SoI2RoUS_Use, RoUS_TIO)|
-  #               +-----------------------+-------------------------+---------------------------+----------------------------+
-  # Commodity RoUS|           0           |            0            |norm(RoUS2SoI_Use, SoI_TIO)|norm(RoUS2RoUS_Use,RoUS_TIO)|
-  #               +-----------------------+-------------------------+---------------------------+----------------------------+
-  # Industry  SoI |norm(SoI_Make, SoI_TCO)|            0            |             0             |             0              |
-  #               +-----------------------+-------------------------+---------------------------+----------------------------+
-  # Industry  RoUS|           0           |norm(RoUS_Make, RoUS_TCO)|             0             |             0              |
-  #               +-----------------------+-------------------------+---------------------------+----------------------------+
-  #                           Total column sum must equal 1                 Total column sum can but shouldn't equal 1
-
-  logging::loginfo("Generating two-region Domestic Use tables ...") 
-  ls <- model$DomesticUseTransactionswithTrade
-  TwoRegionIndustryOutput <- model$x
-  SoI_Industry_Output <- TwoRegionIndustryOutput[endsWith(names(TwoRegionIndustryOutput),
-                                                          state_abb)]
-  RoUS_Industry_Output <- TwoRegionIndustryOutput[endsWith(names(TwoRegionIndustryOutput),
-                                                           "RoUS")]
-  # If industry/comm output == 0, it's not viable to generate A matrix, hence set it to 1.
-  SoI_Industry_Output[SoI_Industry_Output == 0] <- 1
-
-  logging::loginfo("Generating A matrix of SoI2SoI Domestic Use table ...")
-  SoI2SoI_A <- normalizeIOTransactions(ls[["SoI2SoI"]][, industries],
-                                                 SoI_Industry_Output)
-
-  logging::loginfo("Generating A matrix of RoUS2SoI Domestic Use table ...")
-  RoUS2SoI_A <- normalizeIOTransactions(ls[["RoUS2SoI"]][, industries],
-                                                  SoI_Industry_Output)
-
-  logging::loginfo("Generating A matrix of SoI2RoUS Domestic Use table ...")
-  SoI2RoUS_A <- normalizeIOTransactions(ls[["SoI2RoUS"]][, industries],
-                                                  RoUS_Industry_Output)
-
-  logging::loginfo("Generating A matrix of RoUS2RoUS Domestic Use table ...")
-  RoUS2RoUS_A <- normalizeIOTransactions(ls[["RoUS2RoUS"]][, industries],
-                                                   RoUS_Industry_Output)
-
-  logging::loginfo("Assembling the complete A matrix ...")
-  # Assemble A matrix
-  A_top <- cbind(diag(rep(0, length(commodities)*2)),
-                 cbind(rbind(SoI2SoI_A, RoUS2SoI_A),
-                       rbind(SoI2RoUS_A, RoUS2RoUS_A)))
-  colnames(A_top) <- c(1:ncol(A_top))
-  A_btm <- cbind(as.matrix(Matrix::bdiag(list(as.matrix(SoI_A),
-                                              as.matrix(RoUS_A)))),
-                 diag(rep(0, length(industries)*2)))
-  A <- as.matrix(rbind(A_top, setNames(A_btm, colnames(A_top))))
-  rownames(A) <- paste(c(rep(c(state, "RoUS"), each = length(commodities)),
-                         rep(c(state, "RoUS"), each = length(industries))),
-                       c(rep(commodities, 2), rep(industries, 2)),
-                       c(rep("Commodity", length(commodities)*2),
-                         rep("Industry", length(industries)*2)),
-                       sep = ".")
-  colnames(A) <- rownames(A)
-
-  logging::loginfo("Generating the L matrix ...")
-  # Calculate L matrix
-  I <- diag(nrow(A))
-  L <- solve(I - A, tol = 1E-20)
-  #----------------- END OF BUILD A, L MATRICES
-  logging::loginfo("Calculating y (Final Demand totals) of SoI and RoUS ...")
-  # Calculate Final Demand (y)
-  FD_columns  <- unlist(sapply(list("HouseholdDemand", "InvestmentDemand",
-                                    "ChangeInventories", "Export", "Import",
-                                    "GovernmentDemand"),
-                               getVectorOfCodes, iolevel = "Summary"))
-  ita_column <- ifelse(iolevel == "Detail", "F05100", "F051")
-  SoI2SoI_y   <- rowSums(ls[["SoI2SoI"]][, c(FD_columns, ita_column, "ExportResidual")])
-  SoI2RoUS_y  <- rowSums(ls[["SoI2RoUS"]][, c(FD_columns, ita_column)])
-  RoUS2SoI_y  <- rowSums(ls[["RoUS2SoI"]][, c(FD_columns, ita_column)])
-  RoUS2RoUS_y <- rowSums(ls[["RoUS2RoUS"]][, c(FD_columns, ita_column, "ExportResidual")])
-  y <- c(SoI2SoI_y + SoI2RoUS_y, RoUS2SoI_y + RoUS2RoUS_y, rep(0, length(industries)*2))
-  names(y) <- rownames(L)
-
-  logging::loginfo("Validating L*y == commodity and industry output ...")
-  # Validate L * y == Output
-  # Output = c(SoI_TCO, RoUS_TCO, SoI_TIO, RoUS_TIO)
-  output <- c(SoI_Commodity_Output, RoUS_Commodity_Output,
-              SoI_Industry_Output, RoUS_Industry_Output)
-  validation <- as.data.frame(L %*% y - output)
-  colnames(validation) <- "L*y-output"
-
-  validation$rel_diff <- validation$`L*y-output`/output
-  validation$Ly <- as.numeric(L %*% y)
-  validation$output <- output
-  validation[validation$output == 1, "rel_diff"] <- 0
-
-  logging::loginfo("Validation complete.")
-  return(list(A = A, L = L, y = y, Validation = validation))
-}
 
 #' Compare commodity or industry output calculated from Make and Use tables.
 #' @param model A model list object with model specs and IO tables listed
@@ -568,12 +422,12 @@ validateImportFactorsApproach <- function(model, demand = "Consumption"){
   # Calculate LCI using coupled model approach
   # Revised equation from RW email (2023-11-01):
   # LCI <- (s_d * L_d * Y_d) + (s*L*A_m*L_d*Y_d + s*L*Y_m). I.e., s in RW email is analogous to model$B
-  # For validation, model$M = model$M_m, whereas in normally we'd be using model$M_m instead of model$M
+  # For validation, we use M as a stand-in for import emissions , whereas in normally we'd be using model$Q_t
 
   LCI_dm <- (model$M_d %*% y_d) + (M %*% model$A_m %*% model$L_d %*% y_d + M %*% y_m)
 
   cat("\nTesting that LCI results are equivalent between standard and coupled model approaches (i.e., LCI = LCI_dm) when\n")
-  cat("assuming model$M = model$M_m.\n")
+  cat("assuming model$M = model$Q_t.\n")
   print(all.equal(LCI, LCI_dm))
 
   # Calculate LCIA using standard approach
@@ -592,7 +446,7 @@ validateImportFactorsApproach <- function(model, demand = "Consumption"){
   rownames(LCIA_dm) <- colnames(model$N_m)
 
   cat("\nTesting that LCIA results are equivalent between standard and coupled model approaches (i.e., LCIA = LCIA_dm) when\n")
-  cat("assuming model$M = model$M_m.\n")
+  cat("assuming model$M = model$Q_t.\n")
   all.equal(LCIA_dm, LCIA)
 
 }
diff --git a/R/WriteModel.R b/R/WriteModel.R
@@ -209,7 +209,7 @@ writeModelMetadata <- function(model, dirs) {
   model_desc <- file.path(dirs$data, "models.csv")
   ID <- model$specs$Model
   Name <- model$specs$Model
-  Location <- model$specs$ModelRegionAcronyms
+  Location <- model$specs$ModelRegionAcronyms[1]
   Description <- ""
   #Add in sector schema for model
   if (is.null(model$specs$DisaggregationSpecs)) {
@@ -276,8 +276,8 @@ writeModelMetadata <- function(model, dirs) {
   flows <- flows[order(flows$ID),]
   flows$Index <- c(1:nrow(flows)-1)
   flows <- flows[, fields$flows]
-  checkNamesandOrdering(flows$ID, rownames(model$B),
-                        "flows in flows.csv and rows in B matrix")
+  #checkNamesandOrdering(flows$ID, rownames(model$B),
+  #                      "flows in flows.csv and rows in B matrix")
   utils::write.csv(flows, paste0(dirs$model, "/flows.csv"), na = "",
                    row.names = FALSE, fileEncoding = "UTF-8")
 

diff --git a/inst/extdata/VisualizationEssentials.yml b/inst/extdata/VisualizationEssentials.yml
@@ -94,3 +94,55 @@ Indicators:
     - "#428E55" # green, Eonomic & Social
   - "Jobs Supported":
     - "#428E55"
+
+# State color palette is from https://medialab.github.io/iwanthue/
+StateColorPalette: ["#d93762",
+"#44c96a",
+"#9969e4",
+"#89c13d",
+"#8b45b3",
+"#43a132",
+"#db4bb1",
+"#c5ba34",
+"#4d5ac6",
+"#e09129",
+"#6982ec",
+"#919f35",
+"#d974e1",
+"#3d8339",
+"#e43588",
+"#56bd8a",
+"#ac3e9b",
+"#84b96b",
+"#b13670",
+"#5bcbb9",
+"#d73440",
+"#4ab4dd",
+"#d04818",
+"#406dac",
+"#ee5e42",
+"#3b9987",
+"#af3733",
+"#347c52",
+"#e7699f",
+"#5f7f24",
+"#ba84d6",
+"#c99e3b",
+"#74569e",
+"#bdb76b",
+"#8a99dd",
+"#e47d45",
+"#a05e90",
+"#596222",
+"#d88bc1",
+"#876c1d",
+"#964969",
+"#88894b",
+"#e1646d",
+"#b37e48",
+"#e68c95",
+"#a74f1b",
+"#e6a573",
+"#a64d54",
+"#8c532c",
+"#d07259"]