Merge pull request #201 from johanneskoch94/policy_costs

  Add special request features

scripts/output/comparison/policyCosts.R

@@ -20,7 +20,8 @@
 
 suppressPackageStartupMessages(library(tidyverse))
 
-# Function defintions
+###########################################################################
+# ###### START FUNCTION DEFINITONS ########################################
 rm_timestamp <- function(strings,
                          name_timestamp_seperator = "_",
                          timestamp_format = "%Y-%m-%d_%H.%M.%S") {
@@ -36,6 +37,7 @@ rm_timestamp <- function(strings,
   return(my_strings_wo_timeStamp)
 }
 
+
 policy_costs_pdf <- function(policy_costs, 
                              fileName="PolicyCost.pdf") {
   
@@ -113,13 +115,14 @@ policy_costs_pdf <- function(policy_costs,
   
 }
 
+
 write_new_reporting <- function(mif_path, 
                                 scen_name, 
                                 new_polCost_data) {
   
   new_mif_path <- paste0(substr(mif_path,1,nchar(mif_path)-4),"_adjustedPolicyCosts.mif")
   
-  cat(paste0("A mif file with the name ",crayon::green(paste0("REMIND_generic_",scen_name,"_adjustedPolicyCosts.mif"))," is being created in the ",scen_name," outputfolder.\n"))
+  cat(paste0("A mif file with the name ",crayon::green(paste0("REMIND_generic_",scen_name,"_adjustedPolicyCosts.mif"))," is being created in the ",scen_name," output folder.\n"))
   
   my_data <- magclass::read.report(mif_path) 
   my_variables <- grep("Policy Cost", magclass::getNames(my_data[[1]][[1]]), value = TRUE, invert = T)
@@ -134,8 +137,78 @@ write_new_reporting <- function(mif_path,
   
   magclass::write.report(my_data, file = new_mif_path, ndigit = 7, skipempty = FALSE)
   
-  #magclass::write.report2(my_data, file=new_mif_path, ndigit=7, skipempty = FALSE)
+  return(new_mif_path)
+}
+
+
+report_transfers <- function(pol_mif, ref_mif) {
+  
+  # Read in reporting files
+  pol_run <- magclass::read.report(pol_mif,as.list = F)
+  ref_run <- magclass::read.report(ref_mif,as.list = F)
+  
+  # Get model and scenario names
+  md <- magclass::getItems(pol_run,3.2)
+  sc <- magclass::getItems(pol_run,3.1)
+  
+  # Tell the user what's going on
+  cat(paste0("Adding ",crayon::green("transfers")," to ",paste0("REMIND_generic_",sc,"_adjustedPolicyCosts.mif"),".\n"))
+  
+
+  # Get gdploss 
+  gdploss <- pol_run[,,"Policy Cost|GDP Loss (billion US$2005/yr)"]
+  # Add rel gdploss (not in percent)
+  gdploss_rel <- magclass::setNames(pol_run[,,"Policy Cost|GDP Loss|Relative to Reference GDP (percent)"]/100,
+                                    "Policy Cost|GDP Loss|Relative to Reference GDP")
+  # Get gdp
+  gdp_ref <- ref_run[,,"GDP|MER (billion US$2005/yr)"]
+  gdp_policy <- pol_run[,,"GDP|MER (billion US$2005/yr)"]
+  
+  # Calculate difference to global rel gdploss  
+  delta_gdploss <- gdploss_rel[,,] - gdploss_rel["GLO",,]
+  # Calculate transfer required to equalize rel gdploss across regions
+  delta_transfer <- magclass::setNames(delta_gdploss * gdp_ref, 
+                                       "Policy Cost|Transfers (billion US$2005/yr)")
+  delta_transfer_rel <- 100*magclass::setNames(delta_transfer/gdp_ref, 
+                                               "Policy Cost|Transfers|Relative to Reference GDP (percent)")
+  
+
+  # Calculate new gdp variables
+  gdp_withtransfers <- magclass::setNames(gdp_policy + delta_transfer,
+                                          "GDP|MER|w/ transfers (billion US$2005/yr)")
+  gdploss_withtransfers <- magclass::setNames(gdp_ref - gdp_withtransfers, 
+                                              "Policy Cost|GDP Loss|w/ transfers (billion US$2005/yr)")
+  gdploss_withtransfers_rel <- 100*magclass::setNames(gdploss_withtransfers/gdp_ref,
+                                                      "Policy Cost|GDP Loss|w/ transfers|Relative to Reference GDP (percent)")
+  
+  # Correct sets
+  magclass::getSets(delta_transfer, fulldim = F)[3] <- "variable"
+  magclass::getSets(delta_transfer_rel, fulldim = F)[3] <- "variable"
+  magclass::getSets(gdp_withtransfers, fulldim = F)[3] <- "variable"
+  magclass::getSets(gdploss_withtransfers, fulldim = F)[3] <- "variable"
+  magclass::getSets(gdploss_withtransfers_rel, fulldim = F)[3] <- "variable"
+
+  # Bind together
+  my_transfers <- NULL
+  my_transfers <- magclass::mbind(my_transfers, delta_transfer) %>% 
+    magclass::mbind(delta_transfer_rel) %>% 
+    magclass::mbind(gdp_withtransfers) %>% 
+    magclass::mbind(gdploss_withtransfers) %>% 
+    magclass::mbind(gdploss_withtransfers_rel)
+  
+  pol_run <- magclass::read.report(pol_mif)
+  pol_run <- magclass::mbind(pol_run[[1]][[1]][,,], my_transfers) %>% 
+    magclass::add_dimension(dim=3.1,add = "model",nm = md) %>% 
+    magclass::add_dimension(dim=3.1,add = "scenario",nm = sc)
+  
+  magclass::write.report(pol_run, file = pol_mif, ndigit = 7, skipempty = FALSE)
+  
+  return(my_transfers)
 }
+# ###### END FUNCTION DEFINITONS ########################################
+###########################################################################
+
+
 
 
 # Check for an object called "source_include". If found, that means, this script
@@ -144,13 +217,13 @@ write_new_reporting <- function(mif_path,
 # default values, and made over-writable with the command line.
 if(!exists("source_include")) {
   # Set default value
-  outputdirs <- c("../../../output/default_2020-03-03_09.38.01",
-                  "../../../output/default_slim_2020-03-03_11.37.51",
-                  "../../../output/default_slim_2020-03-03_11.37.51",
-                  "../../../output/default_2020-03-03_09.38.01")   
+  outputdirs <- c("base_noEffChange_2020-03-09_17.16.28/",
+                  "base_allT_lab_1point25_2020-03-27_16.12.35/",
+                  "base_allT_lab_1point25_2020-03-27_16.12.35/",
+                  "base_noEffChange_2020-03-09_17.16.28/")   
+  special_requests <- c("2")
   # Make over-writtable from command line
-  lucode::readArgs("outputdirs")
-  
+  lucode::readArgs("outputdirs","special_requests")
 }
 
 
@@ -192,10 +265,22 @@ while (!happy_with_input) {
     cat(crayon::green(paste0("\t", pc_pairs ,"\n")))
     cat("Is that what you intended?\n")
     cat(paste0("Type '",crayon::green("y"),"' to continue, '",crayon::blue("r"),"' to reselect output directories, '",crayon::red("n"),"' to abort: "))
+    
     user_input <- get_line()
+    
     if(user_input %in% c("y","Y","yes")) {
       happy_with_input <- TRUE
       cat(crayon::green("Great!\n"))
+      
+      # Get special requests from user
+      cat(crayon::blue("\nDo you have any special requests?\n"))
+      cat("1: Skip creation of adjustedPolicyCost reporting\n")
+      cat("2: Add transfers to adjustedPolicyCost reporting\n")
+      cat("3: Skip plot creation\n")
+      cat("4: Plot until 2150 in pdf\n")
+      cat("Type the number (or numbers seperated by a comma) to choose the special requests, or nothing to continue without any: ")
+      special_requests <- get_line() %>% str_split(",",simplify = T) %>% as.vector()
+      
     } else if(user_input %in% c("r","R","reselect")) {
       if (exists("choose_folder")) {
         cat("Remember, the order in which you choose the directories should be:\n")
@@ -211,45 +296,71 @@ while (!happy_with_input) {
       cat(crayon::red("\nStopping execution now.\n\n"))
       stop("I can't figure this **** out. I give up. ")
     }
+  } else {
+    happy_with_input <- TRUE
   }
 }
 
 
-# Tell the user what is going on
-cat(crayon::blue("\nPolicy cost computations:\n"))
-
 # Get Policy costs for every policy-reference pair
+cat(crayon::blue("\nComputing Policy costs:\n"))
 tmp_policy_costs_magpie <- mapply(remind::reportPolicyCosts, pol_gdxs, ref_gdxs, SIMPLIFY = FALSE) 
+cat(crayon::green("Done!\n"))
 
-tmp_policy_costs <- tmp_policy_costs_magpie %>% 
-  lapply(quitte::as.quitte) %>% 
-  lapply(select, region, period, data, value)
-
-# Combine results in single tibble, with names like "Pol_w.r.t_Ref"
-policy_costs <- rename(tmp_policy_costs[[1]], !!sym(paste0(pol_names[1], "_w.r.t_",ref_names[1])):=value)
-if (length(tmp_policy_costs)>1){
-  for (i in 2:length(tmp_policy_costs)) {
-    policy_costs <- tmp_policy_costs[[i]] %>% 
-      rename(!!sym(paste0(pol_names[i], "_w.r.t_",ref_names[i])):=value) %>% 
-      left_join(policy_costs, tmp_policy_costs[[i]], by=c("region", "period", "data"))
-  }
+
+# Create "adjustedPolicyCost" reporting file
+if (!"1" %in% special_requests) {
+  cat(crayon::blue("\nCreating new reportings:\n"))
+  pol_mifs <- paste0(dirname(pol_gdxs), "/REMIND_generic_", pol_names, ".mif")
+  new_reporting_files <- mapply(write_new_reporting, pol_mifs, pol_names, tmp_policy_costs_magpie)
+  cat(crayon::green("Done!\n"))
 }
-# and do some pivotting
-policy_costs <- policy_costs %>% 
-  pivot_longer(cols = matches(".*w\\.r\\.t.*"), names_to = "Model Output") %>% 
-  pivot_wider(names_from = data)
 
-cat(crayon::green("Done!\n"))
 
-# Create Pdf
-cat(crayon::blue("\nPdf creation:\n"))
-time_stamp <- format(Sys.time(), "_%Y-%m-%d_%H.%M.%S")
-policy_costs_pdf(policy_costs, fileName = paste0("PolicyCost",time_stamp,".pdf"))
-cat(crayon::green("Done!\n"))
+# Add transfer variables to "adjustedPolicyCost" reporting file
+if ("2" %in% special_requests && !"1" %in% special_requests) {
+  cat(crayon::blue("\nComputing transfers:\n"))
+  ref_mifs <- paste0(dirname(ref_gdxs), "/REMIND_generic_", ref_names, ".mif")
+  transfer_info <- mapply(report_transfers, new_reporting_files, ref_mifs, SIMPLIFY = F)
+  cat(crayon::green("Done!\n"))
+}
 
 
-# Create new reporting file
-cat(crayon::blue("\nMif creation:\n"))
-pol_mifs <- paste0(dirname(pol_gdxs), "/REMIND_generic_", pol_names, ".mif")
-return_check <- mapply(write_new_reporting, pol_mifs, pol_names, tmp_policy_costs_magpie)
-cat(crayon::green("Done!\n\n"))
+# Create Pdf
+if (!"3" %in% special_requests) {
+  cat(crayon::blue("\nCreating plots:\n"))
+  
+  # Add transfers, if they exist
+  if (exists("transfer_info")) {
+    tmp_policy_costs_magpie <- mapply(magclass::mbind, tmp_policy_costs_magpie, transfer_info, SIMPLIFY = F)
+  }
+  
+  tmp_policy_costs <- tmp_policy_costs_magpie %>% 
+    lapply(quitte::as.quitte) %>% 
+    lapply(select, region, period, data, value)
+  
+  # Combine results in single tibble, with names like "Pol_w.r.t_Ref"
+  policy_costs <- rename(tmp_policy_costs[[1]], !!sym(paste0(pol_names[1], "_w.r.t_",ref_names[1])):=value)
+  if (length(tmp_policy_costs)>1){
+    for (i in 2:length(tmp_policy_costs)) {
+      policy_costs <- tmp_policy_costs[[i]] %>% 
+        rename(!!sym(paste0(pol_names[i], "_w.r.t_",ref_names[i])):=value) %>% 
+        left_join(policy_costs, tmp_policy_costs[[i]], by=c("region", "period", "data"))
+    }
+  }
+  # and do some pivotting
+  policy_costs <- policy_costs %>% 
+    pivot_longer(cols = matches(".*w\\.r\\.t.*"), names_to = "Model Output") %>% 
+    pivot_wider(names_from = data)
+  
+  # By default, plots are only created until 2100
+  if (!"4" %in% special_requests) {
+    policy_costs <- policy_costs %>% filter(period<=2100)
+  }
+  
+  time_stamp <- format(Sys.time(), "_%Y-%m-%d_%H.%M.%S")
+  policy_costs_pdf(policy_costs, fileName = paste0("PolicyCost",time_stamp,".pdf"))
+  cat(crayon::green("Done!\n"))
+}
+
+cat("\n")
\ No newline at end of file