##############################################################################################
#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 = 200,  #The number of equilidistant distributed combination of (number of drivers - 1) drivers for each separate
    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
    

){
  
  
  
  
  if(way == 3){
    #randomly distributed
    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, equally distributed input
    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 <- 1
      
        k <- match(#H.tc5.lr01$contributions$var[ii], 
                   names(name_order)[ii] , H.tc5.lr01$gbm.call$predictor.names
                   )
        var.name <- H.tc5.lr01$gbm.call$predictor.names[k]
        x.label <- paste(edrp_labx[k]," ",edrp_labx_unit[k], "  (", round(name_order[ii][[1]] * 100,1),"%)", sep="")
        
        #if(H.tc5.lr01$contributions$var[ii] != "SW"){
        if(FALSE){ 
         if(H.tc5.lr01$contributions$var[ii] == "zm_zi"){
          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 = 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)
          }
        } 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 = cexvar*1.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 = cexvar,
                   #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 = cexvar*2.4)
            }
        
      
        
        if(var.name == "zm_zi" & vari == "H_en"){
          
          
#           temp.lo <- loess(colMeans(rpo_zmzi) ~ ide[[var.name]], span = 0.3)
#           
#           if(common_scale){
#             plot(ide[[var.name]],colMeans(rpo_zmzi), type='l',   lwd = cexvar*1.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 = cexvar,
#                  #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 = cexvar*2.4)
#           }
#           
          
#           pred.day.zmzi <- sapply(1:50, function(zz) mean(pred.matrix[ zz, 1:50]))
#           pred.ngt.zmzi <- sapply(1:50, function(zz) mean(pred.matrix[ zz,31:50]))
#           
#           temp.lo <- loess(pred.day.zmzi ~ ide[[var.name]], span = 0.3)
#           lines(ide[[var.name]], fitted(temp.lo), lty = 1 , col="green", lwd = 1.3)
#           
#           temp.lo <- loess(pred.ngt.zmzi ~ ide[[var.name]], span = 0.3)
#           lines(ide[[var.name]], fitted(temp.lo), lty = 1 , col="blue", lwd = 1.3)          

          
          
          
          #randomly distributed for zm_zi
           rd_zmzi <- sapply(1:length(H.tc5.lr01$var.levels), function(x)
            #x <- 1
            runif(25, min = range(H.tc5.lr01$var.levels[x])[1], max = range(H.tc5.lr01$var.levels[x])[2])
          )
          rd_zmzi <- as.data.frame(rd_zmzi)
          dimnames(rd_zmzi)[[2]] <- H.tc5.lr01$gbm.call$predictor.names
          
          
          #response plot output
          rpo_zmzi <- 
            sapply(1:n, function(y) {
              #y <- 1
              sapply(1:50, function(x) {
                #x <- 2
                edi_sub <- rd_zmzi
                edi_sub$zm_zi <- ide$zm_zi[y]
                edi_sub$cos_azi <- ide$cos_azi[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)
          
          #lines(ide[[var.name]], colMeans(rpo_zmzi[11:50,]), lty = 1 , col="black", lwd = 1.3)
          
          temp.lo <- loess(colMeans(rpo_zmzi[1:8,]) ~ ide[[var.name]], span = 0.3)
          lines(ide[[var.name]], fitted(temp.lo), lty = 1 , col="green", lwd = 1.3)
          
          aa <- loess(rpo[[var.name]] ~ ide[[var.name]], span = 0.3)
          lines(ide[[var.name]], 2*fitted(aa) -fitted(temp.lo), lty = 1 , col="blue", lwd = 1.3)
          
#           temp.lo <- loess(colMeans(rpo_zmzi[9:50,]) ~ ide[[var.name]], span = 0.3)
#           lines(ide[[var.name]], fitted(temp.lo), lty = 1 , col="blue", lwd = 1.3)
#           
#           
#           image.plot(x = ide[[var.name]], y = ide$cos_azi[1:50], z = t(rpo_zmzi[1:50,]), zlim =range(rpo_zmzi)
#                      , xlab = "zm/zi", ylab="cos(azi)"#
#                      #, col=PREcol 
#           )
          
          
        }
        
        
        
        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()