#Functions common to Model Evaluation and Model Diagnostics tools
r.common.functions.version = "4.0"
svn.version.common.functions = "$ProjectRevision: SVN 1942 (2016-01-29) $"
svn.prj.name.common.functions = "$ProjectName: /cerc/Products/Utilities/ModelEvaluationToolkit/Tags/ModelEvaluationToolkit_4_0_0 $"

# Function to read in user defined averaging time and statistic for Model Evaluation tool
# process output averaging and merge with pollutants data frame 
process.output.averaging = function(output.averaging.file, pol)  {
  # Firstly read in the data, checking its format
  output.averaging = read.csv(file=output.averaging.file)
  names(output.averaging) = tolower(names(output.averaging))
  # make sure each column exists
  if (is.null(output.averaging$pollutant) | is.null(output.averaging$output.avg.time.hours) | is.null(output.averaging$output.statistic) | is.null(output.averaging$output.daily.max) ) {
    stop("You have chosen to perform concentration evaluation but you have not included all of the required information in the output averaging times file.")
    quit("no")
  }
  output.averaging$pollutant = tolower(as.character(output.averaging$pollutant))
  output.averaging$output.avg.time.hours = as.numeric(output.averaging$output.avg.time.hours)
  output.averaging$output.statistic = tolower(as.character(output.averaging$output.statistic))
  output.averaging$output.daily.max = tolower(as.character(output.averaging$output.daily.max))
  # check each pollutant in the pollutant data frame (which may be a single pollutant) is included in the user defined averaging file
  for (j in 1:length(pol$pollutant)) {   
    if (length(which(output.averaging$pollutant==pol$pollutant[j]))==0 ) {
      stop(paste("Pollutant",toupper(pol$pollutant[j]),"is missing from the output averaging times file."))
      quit("no")
    }
    if (length(which(output.averaging$pollutant==pol$pollutant[j]))>1 ) {
      stop(paste("Pollutant",toupper(pol$pollutant[j]),"is listed more than once in the output averaging times file."))
      quit("no")
    }
  }
  # now merge with pollutants file
  pol = merge(pol,output.averaging,by="pollutant")
  #Make sure the statistics entered are allowed
  allowed.statistics = c("max","mean","rolling mean")
  for (j in 1:length(pol$output.statistic)) {
    if (length(which(allowed.statistics==pol$output.statistic[j]))==0 ) {
      stop(paste("Averaging statistics for",toupper(pol$pollutant[j]),"must be either max, mean or rolling mean"))
      quit("no")
    }
  }
  # Make sure the user defined averaging time and statistic can be calculated from the modelled averaging time
  #Make sure that the modelled averaging time and statistic is either same as user-defined or 1 hour mean
  not.same = subset(pol, mod.statistic!=output.statistic | mod.avg.time.hours!=output.avg.time.hours)
  not.same.not.1hrmean = subset(not.same, mod.avg.time.hours!=1 | mod.statistic!="mean")
  if (length(not.same.not.1hrmean$pollutant) !=0 ) {
    stop("The modelled data must be a 1-hour mean or the modelled data and the user defined averaging time and statistic must be the same.")
    quit("no")
  }
  #Make sure that if the daily max is to be calculated, then the averaging time is less than or equal to 
  #24 hours, unless statistic is rolling mean.
  for (j in 1:length(pol$pollutant)) {
    if (is.na(pol$output.daily.max[j])) {
      stop(paste("output.daily.max setting for",toupper(pol$pollutant[j]),"is missing"))
      quit("no")      
    }
    if (pol$output.daily.max[j]=="yes" & pol$output.statistic[j] != "rolling mean" & pol$output.avg.time.hours[j]>24) {
      stop(paste("Averaging time for",toupper(pol$pollutant[j]),"must be less than or equal to 24 hours for the daily max of a non-rolling statistic to be calculated"))
      quit("no")
    } else if (pol$output.daily.max[j]!="yes" & pol$output.daily.max[j]!="no") {
      stop(paste("output.daily.max setting for",toupper(pol$pollutant[j]),"must be yes or no"))
    }
  }
  return(pol)
}

#function to average data in Model Evaluation tool to user defined output time and statistic
average.data = function(original.dataframe, pol, j, stns, threshold) {
  # create time.avg string for timeAverage(avg.time=...) option.
  if (pol$output.avg.time.hours[j]==24) {
    time.avg = "day"
  } else {
    time.avg = paste(pol$output.avg.time.hours[j],"hour")
  }
  # applying averaging to subset of data for pollutant and station ID
  if (pol$output.avg.time.hours[j]==pol$mod.avg.time.hours[j] & pol$output.statistic[j]==pol$mod.statistic[j]) {
    data.averaged = original.dataframe
  } else {
    if (pol$output.statistic[j]=="rolling mean") {
      data.averaged.obs = rollingMean(original.dataframe,pollutant="obs",width=pol$output.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right")
      data.averaged.obs = subset(data.averaged.obs,select=c("date","site","rolling"))
      names(data.averaged.obs)[names(data.averaged.obs)=="rolling"]="obs"
      data.averaged.mod = rollingMean(original.dataframe,pollutant="mod",width=pol$output.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right")
      data.averaged.mod = subset(data.averaged.mod,select=c("date","site","rolling"))
      names(data.averaged.mod)[names(data.averaged.mod)=="rolling"]="mod"
      data.averaged = merge(data.averaged.obs,data.averaged.mod)
    } else {  #(pol$output.statistic[j]=="mean" | pol$output.statistic[j]=="max")
      data.averaged = timeAverage(original.dataframe,avg.time=time.avg,statistic=pol$output.statistic[j],data.thresh=threshold,type="site")
    } 
  }
  #Convert to daily maximum if required
  if (pol$output.daily.max[j]=="yes" & time.avg != "day") {
    data.averaged = timeAverage(data.averaged,avg.time="day",statistic="max",data.thresh=threshold,type="site")
  }
  # add station type back in by merge, to mimic all.data format
  names(stns)[names(stns) == "station"] = "site"
  data.averaged$site = as.character(data.averaged$site)
  data.averaged = merge(data.averaged, stns, by="site")
  return(data.averaged)  # because of placement in ModelEvaluation code, actually data averaged by pollutant
}

#Process index file and add to pollutants file
#Inputs index scales file and pollutants, outputs pollutants (merged with index scales)
process.index.scales = function(index.scales,pol) {
  #Read and process the index scales file
  indices = read.csv(file=index.scales)  
  #Make sure the index columns exist and there is at least one index band defined for each pollutant
  if (is.null(indices$pollutant) | is.null(indices$index.units) | is.null(indices$index.statistic) | is.null(indices$index.avg.time.hours) | is.null(indices$index.daily.max) | is.null(indices$i1)) {
    stop("You have chosen to assess forecast indices but you have not included all the required forecast index information in the index scales file.")
    quit("no")
  }
  # convert strings to lower case for consistency within R
  indices$pollutant = tolower(as.character(indices$pollutant))
  indices$index.units = tolower(as.character(indices$index.units))   
  indices$index.statistic=tolower(as.character(indices$index.statistic))
  indices$index.daily.max = tolower(as.character(indices$index.daily.max))
  #Make sure that the units entered are allowed
  for (j in 1:length(indices$pollutant)) {
    if (length(which(allowed.units==indices$index.units[j]))==0) {
      stop(paste("Concentration units for",toupper(indices$pollutant[j]),"must be either ppb, ppm, ug/m3, mg/m3 or g/m3"))
      quit("no")
    } 
  }
  #Make sure the statistics entered are allowed
  for (j in 1:length(indices$pollutant)) {
    if (length(which(allowed.statistics==indices$index.statistic[j]))==0 ) {
      stop(paste("Averaging statistics for",toupper(indices$pollutant[j]),"must be either max, mean or rolling mean"))
      quit("no")
    } 
  }
  #Make sure that the averaging time is numeric
  for (j in 1:length(indices$pollutant)) {
    if (!(is.numeric(indices$index.avg.time.hours[j]))) {
      stop(paste("Averaging time for",toupper(indices$pollutant[j]),"is not numeric."))
      quit("no")
    } 
  }
  # now merge the error-checked indices data with the pollutants data.frame
  pol = merge(pol,indices,by="pollutant")
  
  #Make sure that the modelled averaging time and statistic is either same as index or 1 hour mean
  not.same = subset(pol, mod.statistic!=index.statistic | mod.avg.time.hours!=index.avg.time.hours)
  not.same.not.1hrmean = subset(not.same, mod.avg.time.hours!=1 | mod.statistic!="mean")
  if (length(not.same.not.1hrmean$pollutant) !=0 ) {
    stop("The modelled data must either be a 1 hour mean or the same averaging time and statistic as the index definition.")
    quit("no")
  }
  #Make sure that if the daily max is to be calculated, then the averaging time is less than or equal to 
  #24 hours, unless statistic is rolling mean.
  for (j in 1:length(pol$pollutant)) {
    if (pol$index.daily.max[j]=="yes" & pol$index.statistic[j] != "rolling mean" & pol$index.avg.time.hours[j]>24) {
      stop(paste("Averaging time for",toupper(pol$pollutant[j]),"must be less than or equal to 24 hours for the daily max of a non-rolling statistic to be calculated"))
      quit("no")
    } else if (pol$index.daily.max[j]!="yes" & pol$index.daily.max[j]!="no") {
      stop(paste("index.daily.max setting for",toupper(pol$pollutant[j]),"must be yes or no"))
    }
  }
  return(pol)
}

#Obtain the lower bounds of each index
#Inputs pollutants, outputs lower bounds of each index
find.index.lower.bounds = function(pol) {
  
  #Find out the number of forecast index bands
  #We already know that there's an 'i1' column, so look for 'i2', 'i3' etc
  j=1
  while (length(which(names(pol)==paste("i",j,sep="")))!=0) {
    j=j+1
  }
  number.index.bands = j-1
  index.headers = paste("i",seq(from=1, to=number.index.bands, by=1),sep="")
  lower.bounds.of.index = pol[,which(names(pol)==index.headers[1])]
  if (number.index.bands > 1) {
    for (j in 2:number.index.bands) {
      start = (j-1)*length(pol$pollutant)+1
      finish = j*length(pol$pollutant)
      lower.bounds.of.index[start:finish] = pol[,which(names(pol)==index.headers[j])]
    }
  }
  lower.bounds.of.index = matrix(lower.bounds.of.index,nrow=length(pol$pollutant),ncol=number.index.bands)
  #Check that index lower bounds are given in ascending order
  for (j in 1:length(pol$pollutant)) {
    if (any(lower.bounds.of.index[j,]!=sort(lower.bounds.of.index[j,]))) {
      stop(paste("Index lower bounds are not in ascending order for pollutant",toupper(pol$pollutant[j])))
      quit("no")
    }
  }
  return(lower.bounds.of.index)
}

#Calculate indices of data
#Inputs all data, pollutants, number of pollutants and lower bounds and outputs data with indices applied
find.data.for.indexing = function(data,pol,j,lower.bounds.of.index,threshold) {
  
  data = subset(data, select=c("date","site","obs","mod")) 
  
  #GET OBS AND MOD TO SAME AVERAGING TIME AND STATISTIC AS INDEX
  if (pol$index.avg.time.hours[j]==24) {
    index.avg.time = "day"
  } else {
    index.avg.time = paste(pol$index.avg.time.hours[j],"hour")
  }
  index.statistic = pol$index.statistic[j]
  index.daily.max = pol$index.daily.max[j]
  mod.statistic = pol$mod.statistic[j]
  mod.avg.time = pol$mod.avg.time[j]
  if (index.avg.time==mod.avg.time & index.statistic==mod.statistic) {
    data.averaged = data
  } else {
    if (index.statistic=="rolling mean") {
      data.averaged.obs = rollingMean(data,pollutant="obs",width=pol$index.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right")
      data.averaged.obs = subset(data.averaged.obs,select=c("date","site","rolling"))
      names(data.averaged.obs)[names(data.averaged.obs)=="rolling"]="obs"
      data.averaged.mod = rollingMean(data,pollutant="mod",width=pol$index.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right")
      data.averaged.mod = subset(data.averaged.mod,select=c("date","site","rolling"))
      names(data.averaged.mod)[names(data.averaged.mod)=="rolling"]="mod"
      data.averaged = merge(data.averaged.obs,data.averaged.mod)
    } else {
      data.averaged = timeAverage(data,avg.time=index.avg.time,statistic=index.statistic,data.thresh=threshold,type="site")
    }
  }
  #Convert to daily maximum if required
  if (index.daily.max=="yes" & index.avg.time != "day") {
    data.for.indexing = timeAverage(data.averaged,avg.time="day",statistic="max",data.thresh=threshold,type="site")
  } else {
    data.for.indexing = data.averaged
  }
  #Convert concentrations to index units
  names(data.for.indexing)[which(names(data.for.indexing)=="obs")]="obs.output.units"
  names(data.for.indexing)[which(names(data.for.indexing)=="mod")]="mod.output.units"
  data.for.indexing$obs.index.units = convert.units(data.for.indexing$obs.output.units,pol$output.units[j],pol$index.units[j],pol$conv.ugm3.ppb[j])
  data.for.indexing$mod.index.units = convert.units(data.for.indexing$mod.output.units,pol$output.units[j],pol$index.units[j],pol$conv.ugm3.ppb[j])
  
  #APPLY INDEX BANDS TO CONCENTRATIONS
  data.for.indexing$index.obs = findInterval(data.for.indexing$obs.index.units,lower.bounds.of.index[j,])
  data.for.indexing$index.mod = findInterval(data.for.indexing$mod.index.units,lower.bounds.of.index[j,])
  
  #return data with site names as character strings for future manipulation 
  data.for.indexing$site = as.character(data.for.indexing$site)
  
  return(data.for.indexing)
}
