##############################################################################################
#' @title Equidistant response sensitivity plots

#' @author
#' Ke Xu \email{xuke2012abroad@gmail.com} \cr
#' Andrei Serafimovich \email{andrei.serafimovich@gfz-potsdam.de} \cr
#' Stefan Metzger \email{eddy4R.info@gmail.com}

#' @description 
#' Function defintion. Create equidistant response sensitivity plots for machine learning results to display the response curve of each driver, when each driver is equally distributed within its range. 

#' @param \code{ERFs} The enivronmental response function object reported from either gbm.step() or gbm.step.multic() 
#' @param \code{laby} Type string. The label of y axis shown in the plots, e.g. "fitted H [W m-2]".
#' @param \code{numSpr} Type numeric. The number of equidistant separates in the x axis for plotting, e.g. 50.
#' @param \code{numCom} Number of combination of remaining equidistant distributed drivers for each separate, e.g. 200.
#' @param \code{smth} flag for response function smoothing, e.g. TRUE.
#' @param \code{rug}  flag for plotting a rug of deciles, e.g. TRUE.

#' @return Currently none

#' @references Currently none
#' license: Terms of use of the NEON FIU algorithm repository dated 2015-01-16

#' @keywords environmental response function, plot, equidistant response

#' @examples ersp(ERFs = H.tc5.lr01, laby = "fitted H [W m-2]", numSpr = 50, numCom = 200, smth = TRUE, rug = TRUE)

#' @seealso Currently none

#' @export

# changelog and author contributions / copyrights
#   Ke Xu (TBD)
#     original creation
#   Andrei Serafimovich (2015-07-23)
#     added rug-ticks for uniform distributed samples
#   Stefan Metzger (2015-11-29)
#     added Roxygen2 tags
#   Ke Xu (2016-06-01)
#     apply eddy4R code-style convention
##############################################################################################

#this needs to go into the function call:
data_dir <- "/sap/kxu/kxu_NEON_140912"
png(paste(data_dir, "/tower_ERF/out/Phase2-", in_ver, "-", out_ver, "/", vari,  "/", vari, "_github.png",sep=""), width = 1000, height = 1000, units = "px",
pointsize = 12, bg = "white")
cexvar=3; par(las=0, lwd=cexvar, cex.main=cexvar, cex.axis=cexvar, cex.lab=cexvar, family="times",
pty="s", mfrow=c(3,3), mar=c(6,7,4,2), mgp=c(5,2,0))#, bg = "thistle")

ersp <- function(
  ERFs, #= H.tc5.lr01,  #output from gbm.step.multic()
  laby = "fitted H [W m-2]",
  numSpr = 50,   
  numCom = 200,  
  smth =TRUE,   
  rug = TRUE                    
 
){  
  
  library(gbm)
  data <- eval(parse(text = ERFs$gbm.call$dataframe))
  
  
  #####equally distributed input########
  InpEd <- as.data.frame(sapply(1:length(ERFs$var.levels), function(i)
    #i <- 1
    runif(numCom, min = range(ERFs$var.levels[i])[1], max = range(ERFs$var.levels[i])[2])
  ))
  
  dimnames(InpEd)[[2]] <- ERFs$gbm.call$predictor.names
  
  
  
  #########prepare x axis##############
  idep <- as.data.frame(sapply(1:length(ERFs$var.levels), function(i)
    #i <- 1
    seq(from=range(ERFs$var.levels[i])[1], to=range(ERFs$var.levels[i])[2], length.out=numSpr)
  ))
  
  dimnames(idep)[[2]] <- ERFs$gbm.call$predictor.names
  
  
  
  #####prepare y axis: response plot output######
  depe <- as.data.frame(sapply(1:length(ERFs$var.levels), function(i)
    #i <- 1
    sapply(1:numSpr, function(j)
      #j <- 1
    { tmpInpEd <- InpEd
      tmpInpEd[,i] <- idep[j,i]
      
      mean(predict.gbm(
        object=ERFs,
        newdata=as.data.frame(tmpInpEd),
        n.trees=ERFs$gbm.call$best.trees,
        type="response"
      ))
    }
    )))
  
  dimnames(depe)[[2]] <- dimnames(InpEd)[[2]]
  
  
  
  #####plotting########
  for (i in 1:length(dimnames(depe)[[2]])){
    #i <- 1
    idx <- match(ERFs$contributions$var[i], ERFs$gbm.call$predictor.names)
    var <- ERFs$gbm.call$predictor.names[idx]
    
    #Local Polynomial Regression Fitting
    LPRF <- loess(depe[[var]] ~ idep[[var]], span = 0.3)
    
    plot(idep[[var]],depe[[var]], type='l', lwd = 2,
         ylim = range(depe, fitted(LPRF)),
         xlab = paste(var,"  (",round(ERFs$contributions[i,2],1),"%)",sep=""), 
         ylab = laby)
    
    if (rug & is.vector(data[,ERFs$gbm.call$gbm.x[idx]])) {
      rug(quantile(data[,ERFs$gbm.call$gbm.x[idx]], probs = seq(0, 1, 0.1), na.rm = TRUE), side = 3, lwd = 1, ticksize = 0.03)
      
      rug(quantile(idep[[var]], probs = seq(0, 1, 0.1), na.rm = TRUE), side = 1, lwd = 1, ticksize = 0.03)
    }
    
  
    if (smth & is.vector(idep[[var]])) {      
      lines(idep[[var]],fitted(LPRF), col = 2, lty = 1, lwd = 4)
    }
    
  }
  
  title("Equidistant response sensitivity plots", outer=TRUE) 
  
}

#this needs to go into the function call:
#ersp(H.tc5.lr01)
#dev.off()