##############################################################################################
#' @title Definition function: 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{erf} 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 def.ersp(erf = 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
#   Stefan Metzger (2016-06-02)
#     moved library(gbm) call to DESCRIPTION file import(gbm)
##############################################################################################


def.ersp <- function(
  erf,
  laby,
  numSpr,
  numCom,
  smth,
  rug
  
){  
  
  data <- eval(parse(text = erf$gbm.call$dataframe))
  
  
  #####equally distributed input########
  InpEd <- as.data.frame(sapply(1:length(erf$var.levels), function(i)
    #i <- 1
    runif(numCom, min = range(erf$var.levels[i])[1], max = range(erf$var.levels[i])[2])
  ))
  
  dimnames(InpEd)[[2]] <- erf$gbm.call$predictor.names
  
  
  
  #########prepare x axis##############
  idep <- as.data.frame(sapply(1:length(erf$var.levels), function(i)
    #i <- 1
    seq(from=range(erf$var.levels[i])[1], to=range(erf$var.levels[i])[2], length.out=numSpr)
  ))
  
  dimnames(idep)[[2]] <- erf$gbm.call$predictor.names
  
  
  
  #####prepare y axis: response plot output######
  depe <- as.data.frame(sapply(1:length(erf$var.levels), function(i)
    #i <- 1
    sapply(1:numSpr, function(j)
      #j <- 1
    { tmpInpEd <- InpEd
      tmpInpEd[,i] <- idep[j,i]
      
      mean(gbm::predict.gbm(
        object=erf,
        newdata=as.data.frame(tmpInpEd),
        n.trees=erf$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(erf$contributions$var[i], erf$gbm.call$predictor.names)
    var <- erf$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(erf$contributions[i,2],1),"%)",sep=""), 
         ylab = laby)
    
    if (rug & is.vector(data[,erf$gbm.call$gbm.x[idx]])) {
      rug(quantile(data[,erf$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) 
  
}