##############################################################################################
#Created by: Ke Xu (xuke2012abroad@gmail.com)
#Last change by: Ke Xu (2015-05-12)
#Limited abroadreement: All rights are with the author. Use, modification, porting and distribution of all or part of the algorithms, as well as publication of results derived with the algorithms only with written consent of the author.
#Disclaimer: The author is not liable for any errors in the algorithms of for any erroneous information resulting from use of the algorithms.
##############################################################################################

#equidistance response sensitivity plot function

# data_dir <- "/sap/kxu/kxu_NEON_140912"
# png(paste(data_dir, "/tower_ERF/out/Phase2_", in_ver, "_", out_ver, "/", analysis,  "/", vari, "gbm_inte_resp_way3_Tv.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(
    H.tc5.lr01 = H.tc5.lr01,  #output from gbm.step.multic()
    way = 3,
    n = 50,   #The number of separates in the x axis for plotting
    m = 1000,  #The number of equilidistant distributed combination of (number of drivers - 1) drivers for each separate
    edrp_labx = edrp_labx,
    edrp_labx_unit = edrp_labx_unit,
    title_list = title_list,
    smooth =TRUE,#
    common_scale = TRUE,
    rug = TRUE,                   # plot a rug of deciles
    rug.side = 3,                  # which axis for rug plot? default (3) is top; bottom=1
    rug.lwd = 1,                   # line width for rug plots
    rug.tick = 0.03,              # tick length for rug plots
    train.data = train.data
    

){
  
  
  
  
  if(way == 3){
    #equally distributed input
    edi <- sapply(1:length(H.tc5.lr01$var.levels), function(x)
    #x <- 1
    runif(m, min = range(H.tc5.lr01$var.levels[x])[1], max = range(H.tc5.lr01$var.levels[x])[2])
    )
    edi <- as.data.frame(edi)
    dimnames(edi)[[2]] <- H.tc5.lr01$gbm.call$predictor.names
    
    
    
    #x axis:
    ide <- sapply(1:length(H.tc5.lr01$var.levels), function(x)
    #x <- 1
    seq(from=range(H.tc5.lr01$var.levels[x])[1], to=range(H.tc5.lr01$var.levels[x])[2], length.out=n)
    )
    ide <- as.data.frame(ide)
    dimnames(ide)[[2]] <- H.tc5.lr01$gbm.call$predictor.names
    
    
    
    #response plot output
    rpo <- sapply(1:length(H.tc5.lr01$var.levels), function(x)
    #x <- 1
    sapply(1:n, function(y)
    #y <- 1
    { edi_sub <- edi
        edi_sub[,x] <- ide[y,x]
        
        mean(predict.gbm(
        object=H.tc5.lr01,
        newdata=as.data.frame(edi_sub),
        n.trees=H.tc5.lr01$gbm.call$best.trees,
        type="response"
        ))
    }
    ))
    
    rpo <- as.data.frame(rpo)
    dimnames(rpo)[[2]] <- dimnames(edi)[[2]]
    
    
    
    
    for (ii in 1:length(dimnames(rpo)[[2]])){
        #ii <- 2
      
        k <- match(H.tc5.lr01$contributions$var[ii] , H.tc5.lr01$gbm.call$predictor.names
                   )
        var.name <- H.tc5.lr01$gbm.call$predictor.names[k]
        x.label <- var.name
         # paste(edrp_labx[k]," ",edrp_labx_unit[k], "  (", round(H.tc5.lr01$contributions[ii,2],1),"%)", sep="")
       
        #if(H.tc5.lr01$contributions$var[ii] != "SW"){
        if(FALSE){
       
          #limx <- c(0, 0.1)
          temp.lo <- loess(rpo[[var.name]] ~ ide[[var.name]], span = 0.3)
          
          if(common_scale){
            plot(ide[[var.name]],rpo[[var.name]], type='l',   lwd = 2,
                 ylim = range(rpo, fitted(temp.lo)),
                 xlim = limx,
                 xlab = x.label, ylab = paste("Fitted ",vari,sep=""))
          } else {
            plot(ide[[var.name]],rpo[[var.name]], type='l',   lwd = 2,
                 #ylim = range(rpo, fitted(temp.lo)),
                 xlab = x.label, ylab = paste("Fitted ",vari,sep=""))
          }
          #smooth <- TRUE
          if (smooth & is.vector(ide[[var.name]])) {
            
            lines(ide[[var.name]],fitted(temp.lo), lty = 1, col = 2, lwd = 4)
          }
          
          
          #predict day and night zm/zi
          
         
         else {}}
        
    
        temp.lo <- loess(rpo[[var.name]] ~ ide[[var.name]], span = 0.3)
        
        if(common_scale){
            plot(ide[[var.name]],rpo[[var.name]], type='l',  # lwd = 2,
            ylim = range(rpo, fitted(temp.lo)),
            xlab = x.label, ylab = title_list[which(vari == varis)])
        } else {
          plot(ide[[var.name]],rpo[[var.name]], type='l',  # lwd = 2,
               #ylim = range(rpo, fitted(temp.lo)),
               xlab = x.label, ylab = title_list[which(vari == varis)])
        }
        #smooth <- TRUE
        if (smooth & is.vector(ide[[var.name]])) {
            
            lines(ide[[var.name]],fitted(temp.lo), lty = 1, col = 2, lwd = 4)
        }
        
        
#         if(H.tc5.lr01$contributions$var[ii] == "zm_zi"){
#           pred.day.zmzi <- sapply(1:50, function(zz) mean(pred.matrix[zz,1:30]))
#           pred.ngt.zmzi <- sapply(1:50, function(zz) mean(pred.matrix[zz,31:50]))
#           
#           lines(ide, lty = 3 , col="black", lwd = 1.6)
#         }
        
        
        dataframe.name <- H.tc5.lr01$gbm.call$dataframe
        data <- eval(parse(text = dataframe.name))
        
        
        if (rug & is.vector(data[,H.tc5.lr01$gbm.call$gbm.x[k]])
            ) {
          rug(quantile(data[,H.tc5.lr01$gbm.call$gbm.x[k]], probs = seq(0, 1, 0.1), na.rm = TRUE), side = rug.side, lwd = rug.lwd, ticksize = rug.tick)
          
          rug(quantile(seq(from=range(H.tc5.lr01$var.levels[k])[1], to=range(H.tc5.lr01$var.levels[k])[2], length.out=1000), probs = seq(0, 1, 0.1), na.rm = TRUE), side = 1, lwd = rug.lwd, ticksize = rug.tick)
        }
        #}
       # }
    }
  }
  
  if (way == 1){
    #
    n <- 200
    
    edi <- sapply(1:length(H.tc5.lr01$var.levels), function(x) 
      #x <- 1
      runif(n, min = range(H.tc5.lr01$var.levels[x])[1], max = range(H.tc5.lr01$var.levels[x])[2])
                  # seq(from=range(H.tc5.lr01$var.levels[x])[1], to=range(H.tc5.lr01$var.levels[x])[2], length.out=n)
    )
    
    edi <- as.data.frame(edi)
    
    dimnames(edi)[[2]] <- H.tc5.lr01$gbm.call$predictor.names
    
    #       for (cc in dimnames(edi)[[2]]){
    #         #cc <- "par"
    #         edi[[cc]] <- edi[[cc]][order(edi[[cc]])]
    #       }
    
    
    
    
    #response plot output
    rpo <- sapply(1:length(H.tc5.lr01$var.levels), function(x) 
      #x <- 1
      sapply(1:n, function(y) 
        #y <- 1
      { edi_sub <- edi
        edi_sub[,x] <- edi[y,x]
        
        mean(predict.gbm(
          object=H.tc5.lr01, 
          newdata=as.data.frame(edi_sub),
          n.trees=H.tc5.lr01$gbm.call$best.trees,
          type="response"
        ))
      }
      ))
    
    rpo <- as.data.frame(rpo)
    dimnames(rpo)[[2]] <- dimnames(edi)[[2]]
    
    png(paste(data_dir, "/tower_ERF/out/Phase2_", in_ver, "_", out_ver, "/", vari,  "/", driv_set, vari, "gbm_inte_resp_20.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")
    
    #   for (ii in 1:length(dimnames(rpo)[[2]])){
    #     #ii <- 1
    #     k <- match(H.tc5.lr01$contributions$var[ii],H.tc5.lr01$gbm.call$predictor.names)
    #     var.name <- H.tc5.lr01$gbm.call$predictor.names[k]
    #     x.label <- paste(var.name,"  (",round(H.tc5.lr01$contributions[ii,2],1),"%)",sep="")
    #     
    #     plot(edi[[var.name]],rpo[[var.name]], type='l',
    #          xlab = x.label, ylab = paste("Fitted ",vari,sep=""))
    #   }
    
    for (ii in 1:length(dimnames(rpo)[[2]])){
      #ii <- 1
      k <- match(H.tc5.lr01$contributions$var[ii],H.tc5.lr01$gbm.call$predictor.names)
      var.name <- H.tc5.lr01$gbm.call$predictor.names[k]
      x.label <- paste(var.name,"  (",round(H.tc5.lr01$contributions[ii,2],1),"%)",sep="")
      
      plot(edi[[var.name]][order(edi[[var.name]])],rpo[[var.name]][order(edi[[var.name]])], type='l',
           xlab = x.label, ylab = paste("Fitted ",vari,sep=""))
    }
    
    dev.off()
  }
  
  if(way == 2){
    #
    n <- 15
    
    edi <- sapply(1:length(H.tc5.lr01$var.levels), function(x) 
      #x <- 1
      #runif(n, min = range(H.tc5.lr01$var.levels[x])[1], max = range(H.tc5.lr01$var.levels[x])[2])
      seq(from=range(H.tc5.lr01$var.levels[x])[1], to=range(H.tc5.lr01$var.levels[x])[2], length.out=n)
    )
    
    edi <- as.data.frame(edi)
    
    dimnames(edi)[[2]] <- H.tc5.lr01$gbm.call$predictor.names
    
    edi_sub <- sapply(1:n^(length(H.tc5.lr01$var.levels)), function(i) 
      #i <- 1
      # j <- 1
      as.matrix(c(sapply(1:length(H.tc5.lr01$var.levels), function(j) 
        #edi[floor((i-1)%%(n^(length(H.tc5.lr01$var.levels)))/(n^(length(H.tc5.lr01$var.levels) - j)) + 1), j],
        edi[floor((i-1)%%(n^(length(H.tc5.lr01$var.levels) - j + 1))/(n^(length(H.tc5.lr01$var.levels) - j)) + 1), j]
                         #edi[floor((i-1)%%(n^(length(H.tc5.lr01$var.levels) - j + 1))/(n^(length(H.tc5.lr01$var.levels) - j)) + 1), j],
                         #edi[floor((i-1)%%(n^(length(H.tc5.lr01$var.levels) - 3))/(n^(length(H.tc5.lr01$var.levels) - j)) + 1), j],
                         #edi[floor((i-1)%%(n^(length(H.tc5.lr01$var.levels) - 4))/(n^(length(H.tc5.lr01$var.levels) - j)) + 1), j]
      )
      )
      )
    )
    
    edi_sub <- as.data.frame(t(edi_sub))
    
    dimnames(edi_sub)[[2]] <- H.tc5.lr01$gbm.call$predictor.names
    
    rp <- predict.gbm(
      object=H.tc5.lr01, 
      newdata=as.data.frame(edi_sub),
      n.trees=H.tc5.lr01$gbm.call$best.trees,
      type="response"
    )
    
    edi_sub$response <- rp
    
    rpo <- sapply(1:length(H.tc5.lr01$var.levels), function(x) 
      #x <- 1
      sapply((1:n), function(y) 
        #y <- 1
        mean(edi_sub[which(edi_sub[,x] == edi[y,x]),"response"])
      )
    )
    
    rpo <- as.data.frame(rpo)
    dimnames(rpo)[[2]] <- dimnames(edi)[[2]]
    
    png(paste(data_dir, "/tower_ERF/out/Phase2_", in_ver, "_", out_ver, "/", vari,  "/", driv_set, vari, "gbm_inte_resp_15.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")
    
    for (ii in 1:length(dimnames(rpo)[[2]])){
      #ii <- 1
      k <- match(H.tc5.lr01$contributions$var[ii],H.tc5.lr01$gbm.call$predictor.names)
      var.name <- H.tc5.lr01$gbm.call$predictor.names[k]
      x.label <- paste(var.name,"  (",round(H.tc5.lr01$contributions[ii,2],1),"%)",sep="")
      
      plot(edi[[var.name]],rpo[[var.name]], type='l',
           xlab = x.label, ylab = paste("Fitted ",vari,sep=""))
    }
    
    
    
    dev.off()
  }
  
  
}

# ersp(H.tc5.lr01)
# 
# dev.off()