#  File containing my commonly used R functions.  This one contains goodness-of-fit
#  functions.
#
#  source("/home/waichler/R/functions2.R")
#
# Alternative goodness-of-fit statistics are from
# Legates, D.R., and G.J. McCabe, 1999, "Evaluating the use of 'goodness-of-fit' measures
# in hydrologic and hydroclimatic model validation," Water Resources Research 35(1),
# pp. 233-241.  One of their recommendations is to use a baseline adjusted mean of
# observed values, instead of a grand mean.
# The baseline timestep is logically a seasonal timestep, say a month.  It could be
# continally changing; i.e., a new value for each month in the period of record, but
# the approach I take here is to use the monthly mean (average across all years).  

# Baseline mean computed with respect to months.  
MONTHLY.BASELINE.MEAN <- function(x, month) {
  monthly.baseline.mean <- NULL
  for(i in 1:12) {
    monthly.baseline.mean[i] <- mean(x[ month == i ])
  }
  return(monthly.baseline.mean)
}

# Factor baseline mean.  For each factor (distinct value of flag vector),
# the mean is computed.  A vector the same length as the input is returned.
# Flag could be month, so this is a more general implementation than above.
FACTOR.BASELINE.MEAN <- function(x, flag) {
  ind.valid <- which(!is.na(x))
  unique.baseline.mean <- NULL
  baseline.mean <- NULL
  unique.flags <- unique(flag[ind.valid])
  for(i in 1:length(unique.flags)) {
    unique.baseline.mean[i] <- mean(x[ !is.na(x) & flag == unique.flags[i] ])
  }
  
  for(i in 1:length(x)) {
    baseline.mean[i] <- unique.baseline.mean[unique.flags == flag[i]]
  }
  return(baseline.mean)
}

# E1PRIME()
# Returns first-degree, baseline-adjusted efficiency E1'.
E1PRIME <- function(obs, sim, month) {
  baseline.mean <- MONTHLY.BASELINE.MEAN(obs, month)
  stat <- 1 - (sum(abs(obs - sim)) / sum(abs(obs - baseline.mean[month])))
  return( round(stat, 3) )
}

# E1PRIME.v2()
# Returns first-degree, baseline-adjusted efficiency E1'.  This is different
# from the base E1' function in using a character vector as the basis for grouping
# observations to compute the baseline mean.  Each group (having same flag) has its
# own mean.
E1PRIME.v2 <- function(obs, sim, flag) {
  baseline.mean <- FACTOR.BASELINE.MEAN(obs, flag)
  stat <- 1 - (sum(abs(obs - sim)) / sum(abs(obs - baseline.mean)))
  return( round(stat, 3) )
}

# E1PRIME.v3()
# Returns first-degree, baseline-adjusted efficiency E1'.  This is different
# from the base E1' function by providing the baseline mean as an argument and
# having the same length as the obs, sim vectors.
E1PRIME.v3 <- function(obs, sim, baseline.mean) {
  stat <- 1 - (sum(abs(obs - sim)) / sum(abs(obs - baseline.mean)))
  return( round(stat, 3) )
}

# D1()
# Returns first-degree modified index of agreement d1.
D1 <- function(obs, sim) {
  grand.mean <- mean(obs)
  stat <- 1 - ( sum(abs(obs - sim)) /
      (sum(abs(obs - grand.mean)) + sum(abs(sim - grand.mean))) )
  return( round(stat, 3) )
}

# D1PRIME()
# Returns first-degree baseline-adjusted modified index of agreement d1'.
D1PRIME <- function(obs, sim, month) {
  baseline.mean <- MONTHLY.BASELINE.MEAN(obs, month)
  stat <- 1 - (sum(abs(obs - sim)) /
      (sum(abs(obs - baseline.mean[month])) + sum(abs(sim - baseline.mean[month]))) )
  return( round(stat, 3) )
}

# D1PRIME.v2()
# Returns first-degree baseline-adjusted modified index of agreement d1'.  This is different
# from the base D1' function in using a character vector as the basis for grouping
# observations to compute the baseline mean.  Each group (having same flag) has its
# own mean.
D1PRIME.v2 <- function(obs, sim, flag) {
  baseline.mean <- FACTOR.BASELINE.MEAN(obs, flag)
  stat <- 1 - (sum(abs(obs - sim)) /
      (sum(abs(obs - baseline.mean)) + sum(abs(sim - baseline.mean))) )
  return( round(stat, 3) )
}

# D1PRIME.v3()
# Returns first-degree baseline-adjusted modified index of agreement d1'.  This is different
# from the base D1' function by providing the baseline mean as an argument and
# having the same length as the obs, sim vectors.
D1PRIME.v3 <- function(obs, sim, baseline.mean) {
  stat <- 1 - (sum(abs(obs - sim)) /
      (sum(abs(obs - baseline.mean)) + sum(abs(sim - baseline.mean))) )
  return( round(stat, 3) )
}

# R2()
# Returns square of Pearson's product-moment correlation coefficient
# (R-squared).
R2 <- function(obs, sim) {
  stat <- (cor(sim, obs))^2
  return( round(stat, 3) )
}

# E1()
# Returns first-degree efficiency E1.
E1 <- function(obs, sim) {
  stat <- 1 - (sum( abs(obs - sim) ) / sum( abs(obs - mean(obs)) ))
  return( round(stat, 3) )
}

# E2()
#
# Returns Nash and Sutcliffe efficiency E.
E2 <- function(obs, sim) {
  stat <- 1 - (sum((obs - sim)^2) / sum((obs - mean(obs))^2))
  return( round(stat, 3) )
}

# E2prime()
#
# Uses baseline mean for comparison.
E2PRIME <- function(obs, sim, baseline.mean) {
  stat <- 1 - (sum((obs - sim)^2) / sum((obs - baseline.mean)^2))
  return( round(stat, 3) )
}

# ME()
# Returns mean error.
ME <- function(obs, sim) {
  stat <- mean(sim - obs)
  return( round(stat, 4) )
}

# MAE()
# Returns mean absolute error.
MAE <- function(obs, sim) {
  stat <- mean( abs(sim - obs) )
  return( round(stat, 4) )
}

# RMSE()
# Returns root mean square error.
RMSE <- function(obs, sim) {
  mse <- mean( (sim - obs)^2 )
  stat <- sqrt(mse)
  return( round(stat, 4) )
}


# end of functions2.R
