##=========================================================================##
##                                                                         ##
##                 Start function "make.target.table" version 2            ##
##                                 -----------------                       ##
##  Function to make a table suitable for target diagram plots containing  ##
##             uRMSD and nBIAS for selected categories of data             ##
##                                                                         ##
##       Input:   formula: variables ~ .  (variables to be grouped by)     ##
##                 df (dataframe containing data                      ##
##                 val_obs (column with observed values)                   ##
##                 val_mod (column with modelled values)                   ##
##      Output:    df.target (dataframe with uRMSD and nbias)              ##
##   Reference:    Jolliff(2009) J Mar Sys, 76(1-2), 64-82                 ##
##      Author:    willem.stolte@deltares.nl                               ##
##  webaddress:    https://svn.oss.deltares.nl/repos/openearthtools/       ##
##                 trunk/r/applications/Delft3D/waq/target-function.R      ##
##  testscript:    https://svn.oss.deltares.nl/repos/openearthtools/       ##
##                 trunk/r/applications/Delft3D/waq/target-diagram.R       ##
##   copyright:    Deltares                                                ##
##                                                                         ##
##=========================================================================##

make.target.table3 <- function (formulax, df, val_obs, val_mod, logtrans = F) {

#   TESTDATA TO RUN THE FUNCTION AS SCRIPT
#    df = read.csv("stattable.csv")  
#   formulax = ~ substance + location 
#   val_obs = "value.x"
#   val_mod = "value.y"
#    logtrans = T
  
  require(plyr)

#   ## Do transformation
#   if(logtrans) {
#     min_obs <- (min(df$val_obs))
#     min_mod <- (min(val_mod))
#     min_all <- (min(min_obs, min_mod))
#     val_obs <- (log(val_obs) + min_obs + 1)
#     val_mod <- (log(val_mod) + min_mod + 1)
#   }
   if(logtrans)  {
     print("log transformation used") }   else {
     print("no transformation")  }

  ## calculate square differences (SD)
   if(logtrans)  {
     df.summary <- ddply(df, formulax, here(summarise),
                          observed = log(get(val_obs) + 1),
                          modelled = log(get(val_mod) + 1),
                          SD = ((log(get(val_obs) + 1) - mean(log(get(val_obs) + 1))) - (log(get(val_mod) + 1) - mean(log(get(val_mod) + 1))))^2
                          )
   } else  {
     df.summary <- ddply(df, formulax, here(summarise),
                          observed = get(val_obs),
                          modelled = get(val_mod),
                          SD = ((get(val_obs) - mean(get(val_obs))) - (get(val_mod) - mean(get(val_mod))))^2
                          )
      }
     
   
  ## calculate normalized root mean square difference (uRMSD)
  ## and normalized bias (nBIAS)

      df.target <- ddply(df.summary, formulax, summarise,
                         uRMSD = (sqrt(mean(SD))*sign(sd(modelled)-sd(observed)))/sd(observed),
                         nBIAS = (mean(modelled) - mean(observed))/sd(observed)
        )
  
  return(df.target)
}

####################### end function #########################################

# df.stat = read.csv("d:/weeber/Documents/Laptop/OpenEarthTools/Delft3D/waq/stattable.csv")
# str(df)
# 
# make.target.table2(formulax = ~substance + location, df = df.stat, val_obs = "value.x",val_mod = "value.y")