automatically check project summations (AR6, NAVIGATE), bugfixes ar6climate, fixOnRef, policyCosts

CHANGELOG.md

@@ -20,6 +20,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **37_industry** add process-based steel model as alternative to CES-tree branch
 - **47_regipol** add support for delaying quantity targets and improving regional emission tax convergence
 - **core** change of preference parameters and associated computation of interest rates/mark ups 	
+- **scripts** add script to automatically check project summations from piamInterfaces
+    [[#1587](https://github.com/remindmodel/remind/pull/1587)]
 
 ### fixed
 - fixed weights of energy carriers in `pm_IndstCO2Captured`

DESCRIPTION

@@ -51,7 +51,7 @@ Imports:
     nleqslv,
     optparse,
     piamenv (>= 0.4.0),
-    piamInterfaces (>= 0.12.18),
+    piamInterfaces (>= 0.13.14),
     plotly,
     purrr,
     quitte (>= 0.3128.4),

config/default.cfg

@@ -114,9 +114,13 @@ cfg$pythonEnabled <- "off"  # def <- off
 # Should REMIND output be shown in console? (0:no, 2:logfile, 3:yes)
 cfg$logoption <- 2
 
-# Just list the name of the output scripts that should be used by output.R
-# At the moment there are several R-scripts located in scripts/output/
-cfg$output <- c("reporting","reportCEScalib","rds_report","fixOnRef") #"ar6Climate","emulator"
+# List the name of output scripts from ./scripts/output/ to be used by output.R
+# You can adjust that via an 'output' column in your scenario_config.
+# To overwrite, enter 'reporting,whatever'. To add, enter 'cfg$output,additionalreporting'.
+cfg$output <- c("reporting","reportCEScalib","rds_report","fixOnRef","checkProjectSummations")
+
+# automatically fix remaining differences between run and path_ref
+cfg$fixOnRefAuto <- FALSE
 
 # Set the format for the results folder, type string :date: in order to use the current time stamp in the folder name (e.g. "results:date:") use :title: to use the current title in the folder name
 cfg$results_folder <- "output/:title::date:"

modules/21_tax/on/input/files

@@ -6,4 +6,3 @@ f21_tau_pe_sub.cs4r
 f21_max_pe_sub.cs4r
 f21_tax_convergence.cs4r
 p21_tau_xpres_tax.cs4r
-f21_vehiclesSubsidies.cs4r

output.R

@@ -92,10 +92,14 @@ choose_slurmConfig_output <- function(output) {
   if (!isSlurmAvailable())
     return("direct")
 
-  # Modify slurm options for ar6 reporting, since we want to run MAGICC in parallel and we'll need a lot of memory
-  if ("ar6Climate" %in% output) slurm_options <- paste(slurm_options[1:3], "--tasks-per-node=12 --mem=32000")
-  # reporting.R, in particular remind2::convGDX2MIF, requires at least --mem=8000 of memory
-  if ("reporting" %in% output) slurm_options <- grep("--mem=[0-9]*[0-9]{3}", slurm_options, value = TRUE)
+  # Modify slurm options for reporting options that run in parallel (MAGICC) or need more memory
+  if ("ar6Climate" %in% output) {
+    slurm_options <- paste(slurm_options[1:3], "--tasks-per-node=12 --mem=32000")
+  } else if ("nashAnalysis" %in% output) {
+    slurm_options <- paste(slurm_options[1:3], "--mem=32000")
+  } else if ("reporting" %in% output) {
+    slurm_options <- grep("--mem=[0-9]*[0-9]{3}", slurm_options, value = TRUE)
+  }
 
   if (length(slurm_options) == 1) {
     return(slurm_options[[1]])
@@ -275,12 +279,13 @@ if (comp %in% c("comparison", "export")) {
     } else {
       message("\nStarting output generation for ", outputdir, "\n")
       name <- paste0(output, ".R")
+      scriptsfound <- file.exists(paste0("scripts/output/single/", name))
       if ("--test" %in% flags) {
         message("Test mode, not executing scripts/output/single/", paste(name, collapse = ", "))
-      } else if (all(file.exists(paste0("scripts/output/single/", name)))) {
+      } else {
         if (slurmConfig == "direct") {
           # execute output script directly (without sending it to slurm)
-          for (n in name) {
+          for (n in name[scriptsfound]) {
             message("Executing ", n)
             tmp.env <- new.env()
             tmp.error <- try(sys.source(paste0("scripts/output/single/", n), envir = tmp.env))
@@ -302,9 +307,10 @@ if (comp %in% c("comparison", "export")) {
         }
         # finished
         message("\nFinished ", ifelse(slurmConfig == "direct", "", "starting job for "), "output generation for ", outputdir, "!\n")
-      } else {
-        warning("Skipping ", outputdir, " because some output script selected could not be found ",
-                "in scripts/output/single: ", name[! name %in% dir("scripts/output/single")])
+      }
+      if (any(! scriptsfound)) {
+        warning("Skipping those output script selected that could not be found in scripts/output/single: ",
+                name[! scriptsfound])
       }
     }
 

scripts/output/comparison/policyCosts.R

@@ -297,21 +297,21 @@ if (!"3" %in% special_requests) {
 
   tmp_policy_costs <- tmp_policy_costs_magpie %>%
     lapply(quitte::as.quitte) %>%
-    lapply(select, region, period, data, value)
+    lapply(select, "region", "period", "variable", "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"))
+        left_join(policy_costs, tmp_policy_costs[[i]], by = c("region", "period", "variable"))
     }
   }
   # and do some pivotting
   policy_costs <- policy_costs %>%
     pivot_longer(cols = matches(".*w\\.r\\.t.*"), names_to = "Model Output") %>%
-    pivot_wider(names_from = data)
+    pivot_wider(names_from = "variable")
 
   # By default, plots are only created until 2100
   if (!"4" %in% special_requests) {

scripts/output/single/ar6Climate.R

@@ -224,6 +224,8 @@ climateAssessmentData <- read.quitte(climateAssessmentOutput) %>%
   interpolate_missing_periods(usePeriods, expand.values = FALSE) %>%
   write.mif(remindReportingFile, append = TRUE)
 
+deletePlus(remind_reporting_file, writemif = TRUE)
+
 logmsg <- paste0(
   date(), " postprocessing done! Results appended to REMIND mif '", remindReportingFile, "'\n",
   "ar6Climate.R finished\n"

scripts/output/single/checkProjectSummations.R

@@ -0,0 +1,77 @@
+library(piamInterfaces)
+library(quitte)
+suppressPackageStartupMessages(library(tidyverse))
+
+if(! exists("source_include")) {
+  # Define arguments that can be read from command line
+  outputdir <- "."
+  lucode2::readArgs("outputdir")
+}
+
+scen <- lucode2::getScenNames(outputdir)
+mif  <- file.path(outputdir, paste0("REMIND_generic_", scen, ".mif"))
+mifdata <- as.quitte(mif)
+
+stopmessage <- NULL
+
+absDiff <- 0.00001
+relDiff <- 0.01
+
+# to be skipped for regional aggregation as they are no extensive variables
+varGrep <- paste0("^Tech|CES Price|^Price|^Internal|[Pp]er[- ][Cc]apita|per-GDP|Specific|Interest Rate|",
+                  "Intensity|Productivity|Average Extraction Costs|^PVP|Other Fossil Adjusted|Projected|[Ss]hare")
+unitList <- c("%", "Percent", "percent", "% pa", "1", "share", "USD/capita", "index", "kcal/cap/day",
+             "cm/capita", "kcal/capita/day", "unitless", "kcal/kcal", "m3/ha", "tC/tC", "tC/ha", "years",
+             "share of total land", "tDM/capita/yr", "US$05 PPP/cap/yr", "t DM/ha/yr", "US$2010/kW", "US$2010/kW/yr")
+
+# emi variables where bunkers are added only to the World level
+gases <- c("BC", "CO", "CO2", "Kyoto Gases", "NOx", "OC", "Sulfur", "VOC")
+vars <- c("", "|Energy", "|Energy Demand|Transportation", "|Energy and Industrial Processes",
+          "|Energy|Demand", "|Energy|Demand|Transportation")
+gasvars <- expand.grid(gases, vars, stringsAsFactors = FALSE)
+bunkervars <- unique(sort(paste0("Emissions|", gasvars$Var1, gasvars$Var2)))
+
+
+# failing <- mif %>%
+#   checkSummations(dataDumpFile = NULL, outputDirectory = NULL,  summationsFile = "extractVariableGroups",
+#                   absDiff = 5e-7, relDiff = 1e-8) %>%
+#   filter(abs(diff) >= 5e-7, abs(reldiff) >= 1e-8) %>%
+#   df_variation() %>%
+#   droplevels()
+# if (nrow(failing) > 0) stopmessage <- c(stopmessage, "extractVariableGroups")
+
+for (template in c("AR6", "NAVIGATE")) {
+  message("\n### Check project summations for ", template)
+  d <- generateIIASASubmission(mifdata, outputDirectory = NULL, logFile = NULL,
+                               mapping = template, checkSummation = FALSE)
+  failvars <- d %>%
+    checkSummations(template = template, summationsFile = template, logFile = NULL, dataDumpFile = NULL,
+                    absDiff = absDiff, relDiff = relDiff) %>%
+    filter(abs(diff) >= absDiff, abs(reldiff) >= relDiff) %>%
+    df_variation() %>%
+    droplevels()
+
+  csregi <- d %>%
+    filter(! .data$unit %in% unitList, ! grepl(varGrep, .data$variable)) %>%
+    checkSummationsRegional() %>%
+    rename(World = "total") %>%
+    droplevels()
+  checkyear <- 2050
+  failregi <- csregi %>%
+    filter(abs(.data$reldiff) > 0.5, abs(.data$diff) > 0.00015, period == checkyear) %>%
+    filter(! .data$variable %in% bunkervars) %>%
+    select(-"model", -"scenario")
+  if (nrow(failregi) > 0) {
+    message("For those ", template, " variables, the sum of regional values does not match the World value in 2050:")
+    failregi %>% piamInterfaces::niceround() %>% print(n = 1000)
+    print(paste0(failregi$variable, collapse = ", "))
+  } else {
+    message("Regional summation checks are fine.")
+  }
+
+  if (nrow(failvars) > 0 || nrow(failregi) > 0) stopmessage <- c(stopmessage, template)
+}
+
+if (length(stopmessage) > 0) {
+  stop("Failing summation checks for ", paste(stopmessage, collapse = ", "), ", see above.")
+}

scripts/output/single/fixOnRef.R

@@ -37,7 +37,7 @@ findRefMif <- function(outputdir, envi) {
   return(refmif)
 }
 
-fixMAGICC <- function(d, dref, startyear) {
+fixMAGICC <- function(d, dref, startyear, scenario) {
   magiccgrep <- "^Forcing|^Temperature|^Concentration"
   message("Fixing MAGICC6 data before ", startyear)
   dnew <-
@@ -46,7 +46,9 @@ fixMAGICC <- function(d, dref, startyear) {
              .data$period < startyear),
       filter(d, ! grepl(magiccgrep, .data$variable) |
              .data$period >= startyear)
-    )
+    ) %>%
+    mutate(scenario = factor(scenario)) %>%
+    droplevels()
   return(dnew)
 }
 
@@ -77,21 +79,28 @@ fixOnMif <- function(outputdir) {
   refname <- basename(dirname(refmif))
   d <- quitte::as.quitte(mifs)
   dref <- quitte::as.quitte(refmif)
-  d <- fixMAGICC(d, dref, startyear)
+  d <- fixMAGICC(d, dref, startyear, title)
   failfile <- file.path(outputdir, "log_fixOnRef.csv")
-  fixeddata <- piamInterfaces::fixOnRef(d, dref, ret = "fixed", startyear = startyear, failfile = failfile)
+  fixeddata <- piamInterfaces::fixOnRef(d, dref, ret = "TRUE_or_fixed", startyear = startyear, failfile = failfile)
 
-  if (exists("flags") && isTRUE("--interactive" %in% flags)) {
-    message("\nDo you want to fix that by overwriting ", title, " mif with reference run ", refname, " for t < ", startyear, "?\nType: y/N")
+  update <- paste0("MAGICC data. ", if (! isTRUE(fixeddata)) "Run output.R -> single -> fixOnRef to fix the rest.")
+  if (! isTRUE(fixeddata) && isTRUE(envi$cfg$fixOnRefAuto)) {
+    d <- fixeddata
+    update <- "data from reference run because cfg$fixOnRefAuto=TRUE."
+  } else if (! isTRUE(fixeddata) && exists("flags") && isTRUE("--interactive" %in% flags)) {
+    message("\nDo you want to fix that by overwriting ", title, " mif with reference run ",
+            refname, " for t < ", startyear, "?\nType: y/N")
     if (tolower(gms::getLine()) %in% c("y", "yes")) {
-      message("Updating ", mifs[[1]])
-      tmpfile <- paste0(mifs[[1]], "fixOnMif")
-      quitte::write.mif(fixeddata, tmpfile)
-      file.rename(tmpfile, mifs[[1]])
-      remind2::deletePlus(mifs[[1]], writemif = TRUE)
-      message("Keep in mind to update the runs that use this as `path_gdx_ref` as well.")
+      d <- fixeddata
+      update <- "data from reference run."
     }
   }
+  message("Updating ", mifs[[1]], " with ", update)
+  tmpfile <- paste0(mifs[[1]], "fixOnMif")
+  quitte::write.mif(d, tmpfile)
+  file.rename(tmpfile, mifs[[1]])
+  remind2::deletePlus(mifs[[1]], writemif = TRUE)
+  message("Keep in mind to update the runs that use this as `path_gdx_ref` as well.")
   return(NULL)
 }