createShockFiles.R

library(data.table)

createShockFile = function(scenario){
  
  EnsembleParamtDt= paramDt[Ensemble==scenario[,Ensemble]]
  if(nrow(EnsembleParamtDt) == 0){
    print("shock ensemble name does not match scenario ensemble name")
  }
  else{
    dt=data.table(yearCol = seq(2010,endYear,1))
    for(param in EnsembleParamtDt[,parameter]){
      
      if(param %in% names(scenario)){
        baseValue = scenario[,get(param)]
        setDT(dt)[, c(param) := rep(baseValue,length(yearCol))]
        
        shockYears=EnsembleParamtDt[parameter==param,yearsOfChange]
        recoveryTime=EnsembleParamtDt[parameter==param,recoverBy]
        
        if(grepl( ":", shockYears, fixed = TRUE)){
          
          startEndYear= as.numeric(tstrsplit(shockYears, ":"))
          startOfChange=startEndYear[1]
          endOfChange=startEndYear[2]
          yearsC = seq(startOfChange,endOfChange,1)
          
          shockValue=calculateShockValue(EnsembleParamtDt[parameter==param],baseValue, 1)
          
          yearlyChanges = (shockValue-baseValue)/length(yearsC)
          
          setDT(dt)[yearCol %in% yearsC, get('param') := yearlyChanges*(1+yearCol-startOfChange)+get(param)]
          
          if(!is.na(recoveryTime)){
            
            recovered = endOfChange+recoveryTime
            yearlyRecovery = (shockValue-baseValue)/recoveryTime
            
            setDT(dt)[yearCol %in% seq(endOfChange+1,recovered,1), get('param') := shockValue - yearlyRecovery*(yearCol-endOfChange)]
          }
          else{ 
            setDT(dt)[yearCol > max(yearsC), get('param') := shockValue]
          }
        }
        else{
          yearsC = as.numeric(tstrsplit(shockYears, ","))
          shockValues=calculateShockValue(EnsembleParamtDt[parameter==param],baseValue, length(yearsC))
          
          setDT(dt)[yearCol %in% yearsC, get('param') := shockValues]
          
          if(!is.na(recoveryTime)){
            
            recovered = yearsC+recoveryTime
            yearlyRecovery = (shockValues-baseValue)/recoveryTime
            index=1
            
            for(recoverN in recovered){
              setDT(dt)[yearCol %in% seq(yearsC[index]+1,recoverN,1), get('param') := shockValues[index] - yearlyRecovery[index]*(yearCol-yearsC[index])]
              index=index+1
            }
          }
        }
        
      }
      else{
        print(paste0("Cannot find parameter ", param, " value in scenario table"))
        
      }
    }
    
    ensDir=file.path(scenario[,Ensemble],scenario[,Scenario])
    
    if (dir.exists(file.path(baseOutputDir,ensDir))) {
      print(paste(ensDir, "exists.  Creating shock file"))
      write.csv(dt, file.path(baseOutputDir,ensDir,"shocks.csv"), row.names=FALSE, quote=FALSE)
    }
    else{
      print(paste(ensDir, "does not exist.  Stopping"))
    }
  }
}

calculateShockValue = function(ensembleRow, baseValue, nValues){
  
  if (ensembleRow[,distribution] == "value") {
    s = rep(ensembleRow[,value],nValues)
  }
  else if(ensembleRow[,distribution] == "unif"){
    s = runif(nValues, ensembleRow[,xmin], ensembleRow[,xmax])
  }
  else{
    print(paste0("unexpected distribution ", dist))
    s='NA'
  }
  
  if(ensembleRow[,type]=='percentage'){
    s=baseValue*(1+s)
  }
  else{
    s=s+baseValue
  }
  
  s
}

#This is the shock parameter and distribution file for the ensemble
paramDt = fread(file.path(plum_data_dir, 'paramTest.csv'))
#This is the ensemble file generated by monte carlo script
scenarioTable = fread("C:\\Users\\rhenry2\\eclipse-workspace\\montecarlo_param_files\\juliette_storyline_final.csv")
#scenarioTable = fread(file.path(plumDirectory,ensembleFile))
baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output" # "~/Downloads"
#baseOutputDir="C:\\Users\\rhenry2\\eclipse-workspace"
ensembleFile = commandArgs(trailingOnly = TRUE)[1]
endYear = commandArgs(trailingOnly = TRUE)[2]

for(i in 1:nrow(scenarioTable)){
  #first check if shocks are intended
  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])
    }
  }
}