createShockFiles.R

# TODO: Add comment
# 
# Author: rhenry2
###############################################################################
library(data.table)

createShockFile = function(scenario){
	
#setting all parameters to default values. If parameters are changed according to ssp then pull value from samp table, if not
#set as default same as in config file or to a value that makes no change  
	
  ##This creates the shocks.csv table that PLUM needs to read in to update the parameter values over time
  ##the parameters below and their order must match the order of the method in PLUM in ModelConfig readInShocksFile()
  ##if the parameter has a baseline value as specified in the scenario table then use that otherwise use the PLUM default value
  ##as in ModelConfig 
  paramValues = data.table(Year = seq(2010,2100,1),
                           trendYearValue = c(rep(0,scenarioTable[i,SHOCK_START_YEAR]-2010),
                                              seq(0,1,1/scenarioTable[i,SHOCK_END_YEAR-SHOCK_START_YEAR]),
                                              rep(1,2100-scenarioTable[i,SHOCK_END_YEAR])),
                           TRADE_BARRIER_MULTIPLIER =scenario[, if(exists("TRADE_BARRIER_MULTIPLIER")) TRADE_BARRIER_MULTIPLIER else 1.0],
                           TRANSPORT_COST = scenario[, if(exists("TRANSPORT_COST")) TRANSPORT_COST else 0.1],
                           TECHNOLOGY_CHANGE_ANNUAL_RATE = scenario[, if(exists("TECHNOLOGY_CHANGE_ANNUAL_RATE")) TECHNOLOGY_CHANGE_ANNUAL_RATE else 0.002],
                           MEAT_EFFICIENCY = scenario[, if(exists("MEAT_EFFICIENCY")) MEAT_EFFICIENCY else 1.0],
                           FERTILISER_COST_PER_T = scenario[, if(exists("FERTILISER_COST_PER_T")) FERTILISER_COST_PER_T else 1.4],
                           OTHER_INTENSITY_COST = scenario[, if(exists("OTHER_INTENSITY_COST")) OTHER_INTENSITY_COST else 0.8],
                           TRANSPORT_LOSSES = scenario[, if(exists("TRANSPORT_LOSSES")) TRANSPORT_LOSSES else 0.05],
                           IRRIG_COST_MULTIPLIER = scenario[, if(exists("IRRIG_COST_MULTIPLIER")) IRRIG_COST_MULTIPLIER else 1.0],
                           RUMINANT_CHANGE_ANNUAL_RATE = scenario[, if(exists("RUMINANT_CHANGE_ANNUAL_RATE")) RUMINANT_CHANGE_ANNUAL_RATE else 0.0],
                           MONOGASTRIC_CHANGE_ANNUAL_RATE = scenario[, if(exists("MONOGASTRIC_CHANGE_ANNUAL_RATE")) MONOGASTRIC_CHANGE_ANNUAL_RATE else 0.0],
                           IRRIG_EFF_MULTIPLIER = scenario[, if(exists("IRRIG_EFF_MULTIPLIER")) IRRIG_EFF_MULTIPLIER else 1.0],
                           FINANCIAL_SPECULATION_MULTIPLIER = scenario[, if(exists("FINANCIAL_SPECULATION_MULTIPLIER")) FINANCIAL_SPECULATION_MULTIPLIER else 1.0]
  )
  
	paramValues[, names(paramValues) := lapply(.SD, as.numeric)]
	
	##This decides if and when a cyber shock or financial shock will happen according to the shockProbabilities.csv table and the SSP
	shocksHappen = data.table(SSP_Scenario = scenario[,SSP_SCENARIO], Scenario = scenario[,Scenario], Year = seq(2010,2100,1), 
			cyber=as.logical(rbinom(91,size=1,prob = shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),cyber])), 
			financial =as.logical(rbinom(91,size=1,prob = shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),financial])))
	
	##this creates the table for the shocksRecords.csv file. so contains only the shock value, hence where no shock the value is 0
	shockValues = data.table(Year = rep(seq(2010,2100,1),6),shock = rep(c("protectionism","automation","diet","cyber","climate","financial"), each=91), 
			TRADE_BARRIER_MULTIPLIER = 0, 
			TRANSPORT_COST = 0,
			TECHNOLOGY_CHANGE_ANNUAL_RATE = 0,
			MEAT_EFFICIENCY = 0,
			FERTILISER_COST_PER_T=0,
			OTHER_INTENSITY_COST = 0,
			TRANSPORT_LOSSES = 0,
			IRRIG_COST_MULTIPLIER = 0, 
			RUMINANT_CHANGE_ANNUAL_RATE = 0, 
			MONOGASTRIC_CHANGE_ANNUAL_RATE = 0,
			IRRIG_EFF_MULTIPLIER = 0,
			FINANCIAL_SPECULATION_MULTIPLIER = 0)
	
	##add trend values to shocks.csv file by drawing a random value from the distributions file for each parameter for the type of trend
	for (trendType in c("protectionism", "automation", "diet")){
		if(as.logical(rbinom(1,size=1,prob =shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),get(trendType)]))){
			for(paramName in shockDistFile[shock == trendType]$pname){
				#draw a random parameter value from the distribution file
				randomValue = drawValue(trendType, paramName)
				#put in records table
				shockValues[shock == trendType, c(paramName) := randomValue * paramValues[,trendYearValue]]
				##add the value to the shocks.csv table and adjust depending on whether multiplicative or additive, also adjust value
				##according to year of implementation 
				if(shockDistFile[shock == trendType & pname == paramName]$type == 'additive'){
					paramValues[,c(paramName) := get(paramName) + randomValue * trendYearValue]		
				}
				else{
					paramValues[,c(paramName) := get(paramName) * (1 + randomValue*trendYearValue)]
				}
			}		
		}	
	}
	
	##do the same as above for individual year shocks, don't adjust by implementation year as cyber shocks and financial shocks are single year occurances
	for(year in paramValues[,Year]){
		for (shockType in names(shocksHappen)[4:5]){
			if(shocksHappen[Year == year, get(shockType)]){	
				for(paramName in shockDistFile[shock == shockType]$pname){
					
					randomValue = drawValue(shockType, paramName)
					#put in records table
					shockValues[shock == shockType & Year == year, c(paramName) := randomValue]
					
					#shock is changed from baseline value not change from previous year 
					if(shockDistFile[shock == shockType & pname == paramName]$type == 'additive'){
						paramValues[Year == year,c(paramName) := get(paramName) + randomValue]
					}
					else{
						paramValues[Year == year,c(paramName) := get(paramName) * (1 + randomValue)]
					}
				}		
			}
		}
	}
	
	
	
	#if negative set to zero
	for(col in names(paramValues)) set(paramValues, i=which(paramValues[[col]]<0), j=col, value=0)
	
	ensDir=file.path(scenario[,Ensemble],scenario[,Scenario])
	
	if (dir.exists(file.path(baseOutputDir,ensDir))) {
		print(paste(ensDir, "exists.  Creating shock file"))
		write.csv(paramValues, file.path(baseOutputDir,ensDir,"shocks.csv"), row.names=FALSE, quote=FALSE)
		write.csv(shockValues, file.path(baseOutputDir,ensDir,paste0("shocksRecord_",scenario[,Scenario],".csv")), row.names=FALSE, quote=FALSE)
	}
	else{
		print(paste(ensDir, "does not exist.  Stopping"))
	}
}

drawValue <- function(shockType, paramName){
	
	runif(1, shockDistFile[shock == shockType & pname == paramName]$min, 
			shockDistFile[shock == shockType & pname == paramName]$max)
}

plumDirectory = "/exports/csce/eddie/geos/groups/LURG/models/PLUM"
baseDataDir = "/exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/data"
baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output" # "~/Downloads"

shockDistFile = fread(file.path(baseDataDir, "shockDistributions.csv"))
shockProbs = fread(file.path(baseDataDir,"shockProbabilities.csv"))

ensembleFile = commandArgs(trailingOnly = TRUE)[1]

scenarioTable = fread(file.path(plumDirectory,ensembleFile))

for(i in 1:nrow(scenarioTable)){
#first if column exists
	if("SHOCKS_POSSIBLE" %in% colnames(scenarioTable)){
		if(scenarioTable[i,SHOCKS_POSSIBLE]){
			print(paste("Shocks possible in ",scenarioTable[i,Ensemble],scenarioTable[i,Scenario]))
			createShockFile(scenario = scenarioTable[i])
		}
	}
}