From 14252e98712b8407ed1333bb30373435dbf07712 Mon Sep 17 00:00:00 2001
From: R0slyn <roslyn.henry.08@aberdeen.ac.uk>
Date: Tue, 1 Oct 2019 17:18:49 +0100
Subject: [PATCH] Script tweaks for Frances for generating shock.csv files when
 creating scenarios.

---
 data/shockDistributions.csv | 16 +++++-----
 scripts/createShockFiles.R  | 64 ++++++++++++++++++++-----------------
 2 files changed, 43 insertions(+), 37 deletions(-)

diff --git a/data/shockDistributions.csv b/data/shockDistributions.csv
index 5ebfa56d..bb9d8cdb 100644
--- a/data/shockDistributions.csv
+++ b/data/shockDistributions.csv
@@ -2,22 +2,22 @@ shock,pname,min,max,type
 protectionism,TRADE_BARRIER_MULTIPLIER,0.1,0.2,additive
 protectionism ,TRANSPORT_COST,0.5,1,multiplicative
 protectionism ,TRANSPORT_LOSSES,0.05,0.1,additive
-automation,TRANSPORT_COST,-0.25,-0.5,multiplicative
-automation,TRANSPORT_LOSSES,-0.05,-0.1,additive
+automation,TRANSPORT_COST,-0.5,-0.25,multiplicative
+automation,TRANSPORT_LOSSES,-0.1,-0.05,additive
 automation,TECHNOLOGY_CHANGE_ANNUAL_RATE,0.002,0.004,additive
 automation,MEAT_EFFICIENCY,0.025,0.05,additive
-automation,FERTILISER_COST_PER_T,-0.25,-0.5,multiplicative
-automation,OTHER_INTENSITY_COST,-0.25,-0.5,multiplicative
-automation,IRRIG_COST_MULTIPLIER,-0.25,-0.5,additive
+automation,FERTILISER_COST_PER_T,-0.5,-0.25,multiplicative
+automation,OTHER_INTENSITY_COST,-0.5,-0.25,multiplicative
+automation,IRRIG_COST_MULTIPLIER,-0.5,-0.25,additive
 automation,IRRIG_EFF_MULTIPLIER,0.1,0.2,additive
 diet,RUMINANT_CHANGE_ANNUAL_RATE,0.4,0.6,additive
 diet,MONOGASTRIC_CHANGE_ANNUAL_RATE,0.4,0.6,additive
 cyber,TRANSPORT_COST,0.5,1,multiplicative
 cyber,TRANSPORT_LOSSES,0.1,0.2,additive
-cyber,TECHNOLOGY_CHANGE_ANNUAL_RATE,-0.002,-0.004,additive
-cyber,MEAT_EFFICIENCY,-0.025,-0.05,additive
+cyber,TECHNOLOGY_CHANGE_ANNUAL_RATE,-0.004,-0.002,additive
+cyber,MEAT_EFFICIENCY,-0.05,-0.025,additive
 cyber,FERTILISER_COST_PER_T,0.25,0.5,multiplicative
 cyber,OTHER_INTENSITY_COST,0.25,0.5,multiplicative
 cyber,IRRIG_COST_MULTIPLIER,0.25,0.5,additive
-cyber,IRRIG_EFF_MULTIPLIER,-0.1,-0.2,additive
+cyber,IRRIG_EFF_MULTIPLIER,-0.2,-0.1,additive
 financial,FINANCIAL_SPECULATION_MULTIPLIER,0.5,1,additive
diff --git a/scripts/createShockFiles.R b/scripts/createShockFiles.R
index 37c1bca3..2e106535 100644
--- a/scripts/createShockFiles.R
+++ b/scripts/createShockFiles.R
@@ -9,28 +9,36 @@ 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  
 	
-	paramValues = data.table(Year = seq(2010,2100,1),
-			trendYearValue = c(rep(0,16),seq(0.2,1,0.2),rep(1,70)),
-			TRADE_BARRIER_MULTIPLIER =scenario[,TRADE_BARRIER_MULTIPLIER], 
-			TRANSPORT_COST = scenario[,TRANSPORT_COST],
-			TECHNOLOGY_CHANGE_ANNUAL_RATE = scenario[,TECHNOLOGY_CHANGE_ANNUAL_RATE],
-			MEAT_EFFICIENCY = scenario[,MEAT_EFFICIENCY],
-			FERTILISER_COST_PER_T=scenario[,FERTILISER_COST_PER_T],
-			OTHER_INTENSITY_COST = scenario[,OTHER_INTENSITY_COST],
-			TRANSPORT_LOSSES = 0.05,
-			IRRIG_COST_MULTIPLIER = 1.0, 
-			RUMINANT_CHANGE_ANNUAL_RATE = 0.0, 
-			MONOGASTRIC_CHANGE_ANNUAL_RATE = 0.0,
-			IRRIG_EFF_MULTIPLIER = 1.0,
-			FINANCIAL_SPECULATION_MULTIPLIER = 1.0)
-	
+  ##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), 
-			climate =as.logical(rbinom(91,size=1,prob = shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),climate])), 
 			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,
@@ -45,14 +53,16 @@ createShockFile = function(scenario){
 			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]		
 				}
@@ -63,9 +73,9 @@ createShockFile = function(scenario){
 		}	
 	}
 	
-	
+	##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:6]){
+		for (shockType in names(shocksHappen)[4:5]){
 			if(shocksHappen[Year == year, get(shockType)]){	
 				for(paramName in shockDistFile[shock == shockType]$pname){
 					
@@ -85,7 +95,9 @@ createShockFile = function(scenario){
 		}
 	}
 	
-	#if negative set to zero, or need to change shock distributions so zero can't happen 
+	
+	
+	#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])
@@ -102,14 +114,8 @@ createShockFile = function(scenario){
 
 drawValue <- function(shockType, paramName){
 	
-	randomValue = runif(1, abs(shockDistFile[shock == shockType & pname == paramName]$min), 
-			abs(shockDistFile[shock == shockType & pname == paramName]$max))
-	
-	if (shockDistFile[shock == shockType & pname == paramName]$min <0){
-		randomValue = randomValue*-1
-	}
-	
-	randomValue
+	runif(1, shockDistFile[shock == shockType & pname == paramName]$min, 
+			shockDistFile[shock == shockType & pname == paramName]$max)
 }
 
 plumDirectory = "/exports/csce/eddie/geos/groups/LURG/models/PLUM"
-- 
GitLab