Merge pull request #111 from johanneskoch94/policyCosts

Add 'policyCosts.R', an output.R comparison script to compute policy costs

scripts/output/comparison/policyCosts.R

+# |  (C) 2006-2019 Potsdam Institute for Climate Impact Research (PIK)
+# |  authors, and contributors see CITATION.cff file. This file is part
+# |  of REMIND and licensed under AGPL-3.0-or-later. Under Section 7 of
+# |  AGPL-3.0, you are granted additional permissions described in the
+# |  REMIND License Exception, version 1.0 (see LICENSE file).
+# |  Contact: remind@pik-potsdam.de
+
+# policyCosts.R
+#
+# This script produces a pdf in which the policy costs of policy-runs with 
+# respect to specified reference runs are displayed.
+# It can be called via the output.R script. If that is the case, the order 
+# in which the runs are selected is important:
+#   1: policy run 1
+#   2: reference run 1
+#   3: policy run 2
+#   4: reference run 2
+# and so on... 
+# 
+
+suppressPackageStartupMessages(library(tidyverse))
+
+# TODO: Loading magclass shouldn't actually be necessary... but running on the 
+# cluster without doing so throws an error...
+suppressPackageStartupMessages(library(magclass)) 
+
+# Function defintions
+get_run_names <- function(filepaths) {
+  # The run names are found between a "/" and the time-stamp of the run.
+  run_names <- str_match(filepaths, pattern = "^.*\\/(.*)(_....-..-..*)$")[,2]
+  
+  # If for some reason there isn't a time-stamp, return the name of the folder
+  if (any(is.na(run_names))) {
+    run_names <- str_match(filepaths, pattern = "^.*\\/(.*)\\/.*$")[,2]
+  }
+  return(run_names)
+}
+
+policy_costs_pdf <- function(policy_costs,  fileName="PolicyCost.pdf") {
+  
+  cat(paste0("A pdf with the name ",crayon::green(fileName)," is being created.\n"))
+  
+  template <-  c("\\documentclass[a4paper,landscape,twocolumn]{article}",
+                 "\\setlength{\\oddsidemargin}{-0.8in}",
+                 "\\setlength{\\evensidemargin}{-0.5in}",
+                 "\\setlength{\\topmargin}{-0.8in}",
+                 "\\setlength{\\parindent}{0in}",
+                 "\\setlength{\\headheight}{0in}",
+                 "\\setlength{\\topskip}{0in}",
+                 "\\setlength{\\headsep}{0in}",
+                 "\\setlength{\\footskip}{0.2in}",
+                 "\\setlength\\textheight{0.95\\paperheight}",
+                 "\\setlength\\textwidth{0.95\\paperwidth}",
+                 "\\setlength{\\parindent}{0in}",
+                 "\\usepackage{float}",
+                 "\\usepackage[bookmarksopenlevel=section,colorlinks=true,linkbordercolor={0.9882353 0.8352941 0.7098039}]{hyperref}",
+                 "\\hypersetup{bookmarks=true,pdfauthor={GES group, PIK}}",
+                 "\\usepackage{graphicx}",
+                 "\\usepackage[strings]{underscore}",
+                 "\\usepackage{Sweave}",
+                 "\\begin{document}",
+                 "<<echo=false>>=",
+                 "options(width=110)",
+                 "@")
+  
+  # Create temporaray folder in which to create the policyCost pdf
+  system("mkdir tmp_policyCost")
+  
+  # Open stream in tmp_folder
+  sw <- lusweave::swopen(fileName, folder = "tmp_policyCost", template = template)
+  
+  # Write title
+  lusweave::swlatex(sw,"\\section{Policy Costs}")
+  
+  # Loop over subsections and create plots
+  sub_section_variables <- grep("Policy Cost", names(policy_costs), value = T)
+  sub_section_titles <- str_match(sub_section_variables, "\\|(.*) \\(|\\|")[,2]
+  
+  my_ggplot <- function(data, y_var) {
+    gg1 <- ggplot(data) +
+      geom_line(aes(x=period, y=!!sym(y_var), color = `Model Output`)) +
+      geom_point(aes(x=period, y=!!sym(y_var), color = `Model Output`)) +
+      geom_vline(aes(xintercept = min(period)), linetype = 2) +
+      facet_wrap("region", ncol = 3, scales = "free_y") +
+      xlab("Year") +
+      theme_bw() +
+      theme(legend.position = "bottom", 
+            legend.direction = "vertical",
+            legend.title.align=0.5,
+            legend.title = element_text(face="bold"),
+            axis.title.x = element_text(face="bold"),
+            axis.title.y = element_text(face="bold"))
+    return(gg1)
+  }
+  
+  for(i in 1:length(sub_section_titles)){
+    lusweave::swlatex(sw, paste0("\\subsection{", sub_section_titles[i], "}"))
+    
+    p <- my_ggplot(filter(policy_costs, region=="GLO"), sub_section_variables[i])
+    lusweave::swfigure(sw,print,p,sw_option="height=8,width=8")
+    
+    p <- my_ggplot(filter(policy_costs, region!="GLO"), sub_section_variables[i])
+    lusweave::swfigure(sw,print,p,sw_option="height=9,width=8")
+  }
+  
+  # Close stream to-pdf
+  lusweave::swclose(sw)
+  
+  # Copy pdf from tmp folder to remind folder and delete tmp folder
+  system(paste0("mv tmp_policyCost/",fileName," ."))
+  system("rm -r tmp_policyCost")
+  
+}
+
+# Check for an object called "source_include". If found, that means, this script
+# is being called from another (output.R most likely), and the input variables
+# are already in the environment. If not found, the input variables are given
+# default values, and made over-writable with command line values.
+if(!exists("source_include")) {
+  # Set default value
+  outputdirs <- c("../../femulator_hq/remind_runs_lab/base_allT_lab_0point7_2020-03-12_09.36.25",
+                  "../../femulator_hq/remind_runs_lab/base_allT_lab_0point95_2020-03-16_15.35.20",
+                  "../../femulator_hq/remind_runs_lab/base_allT_lab_0point7_2020-03-12_09.36.25",
+                  "../../femulator_hq/remind_runs_lab/base_allT_lab_0point8_2020-03-12_09.36.24")   
+  # Make over-writtable from command line
+  lucode::readArgs("outputdirs")
+  
+}
+
+# Check that the input variable "outputdirs" has an even number of entries.
+if (length(outputdirs) %% 2!=0) {
+  cat(paste0(crayon::red("\nERROR: "), "The number of directories is not even!\n"))
+  cat("Remember, the order in which you choose the directories should be:\n")
+  cat("\t1: policy run 1\n\t2: reference run 1\n\t3: policy run 2\n\t4: reference run 2\nand so on...\n")
+  cat(crayon::red("\nStopping execution now.\n\n"))
+  stop("Number of directories is not even!")
+}
+
+# Get gdx paths
+pol_gdxs <- paste0(outputdirs[seq(1,length(outputdirs),2)], "/fulldata.gdx")
+ref_gdxs <- paste0(outputdirs[seq(2,length(outputdirs),2)], "/fulldata.gdx")
+
+# Get run names
+pol_names <- get_run_names(pol_gdxs)
+ref_names <- get_run_names(ref_gdxs)
+
+# Define pol-ref, policyCost pair names
+pc_pairs <- paste0(pol_names, "_w.r.t_",ref_names)
+
+# If scrpit was called from output.R, check with user if the pol-ref pairs are the 
+# ones he wanted. 
+if(exists("source_include")) {
+  cat(crayon::blue("\nPlease confirm the set-up:\n"))
+  cat("From the order with which you selected the directories, the following policy-cost curves will be created:\n")
+  cat(crayon::green(paste0("\t", pc_pairs ,"\n")))
+  cat("Is that what you intended?\n")
+  cat(paste0("Type '",crayon::green("y"),"' to continue, '",crayon::red("n"),"' to abort: "))
+  user_input <- get_line()
+  if(!user_input %in% c("y","Y","yes")) {
+    cat(crayon::red("\nShame... \n"))
+    cat("Remember, the order in which you choose the directories should be:\n")
+    cat("\t1: policy run 1\n\t2: reference run 1\n\t3: policy run 2\n\t4: reference run 2\nand so on...\n")
+    cat(crayon::red("\nStopping execution now.\n\n"))
+    stop("Wrong set up.")
+  } else cat(crayon::green("Great!\n"))
+}
+
+# Tell the user what's going on
+cat(crayon::blue("\nPolicy cost computations:\n"))
+
+# Get Policy costs for every policy-reference pair
+tmp_policy_costs <- mapply(remind::reportPolicyCosts, pol_gdxs, ref_gdxs, SIMPLIFY = FALSE) %>% 
+  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)
+
+# Tell the user what's going on
+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\n"))