#setwd('/pub1/data/mg_projects/projects/web_script/R/')
library(optparse)
option_list <- list(
  make_option(c("-i", "--infile"), type = "character", default = '/pub1/data/mg_projects/projects/web_script/tool_runing/52219422b02ea81046dd5ff10e20b480/input.json',
              action = "store", help = "Input a exp file path!"
  ),
  make_option(c("-o", "--outfile"), type = "character", default = '/pub1/data/mg_projects/projects/web_script/web_file_catche/runing/dae42f7681eca28355c6ba9105a3b5ff',
              action = "store", help = "Input a outfolder path!"
  )
)
logs=c()
#stx=rbind()
tryCatch({
  Args <- commandArgs()
  opt = parse_args(OptionParser(option_list = option_list, usage = "Data press"))
  #logs=c(logs,paste0('geting data:',paste0(paste0(names(opt),'=',opt),collapse = ',')))
  logs=c(logs,paste0('run nomogram.R-',basename(opt$outfile)))
  #library("rjson")
  library(jsonlite)
  data<-jsonlite::stream_in(file(opt$infile),pagesize = 100)
  

library(rms)
getFunList=function(){
  all_fun=list()
  s1=function(x) return(sur1(cut.time[1],x))
  s2=function(x) return(sur1(cut.time[2],x))
  s3=function(x) return(sur1(cut.time[3],x))
  s4=function(x) return(sur1(cut.time[4],x))
  s5=function(x) return(sur1(cut.time[5],x))
  s6=function(x) return(sur1(cut.time[6],x))
  s7=function(x) return(sur1(cut.time[7],x))
  s8=function(x) return(sur1(cut.time[8],x))
  s9=function(x) return(sur1(cut.time[9],x))
  s10=function(x) return(sur1(cut.time[10],x))
  all_fun2=list(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10)
  return(all_fun2[1:length(cut.time)])
}

getNomogramData=function(nom){
  all_point_data1=rbind()
  all_point_data2=rbind()
  all_point_data=rbind()
  for(i in 1:length(nom)){
    if(names(nom)[i]=='total.points'){
      #all_point_data=rbind(all_point_data
      #                     ,c('TOL',)
    }else if(names(nom)[i]=='lp'){
  #    all_point_data=rbind(all_point_data
   #                        ,cbind('LP',nom[[i]][[1]],nom[[i]][[2]]))
    }else if(substr(names(nom)[i],1,1)=='C'){
      all_point_data1=rbind(all_point_data1
                           ,cbind(names(nom)[i]
                                  ,nom[[i]][[1]],nom[[i]][[3]]))
    }else if(substr(names(nom)[i],1,1)=='P'){
      all_point_data2=rbind(all_point_data2
                           ,cbind(names(nom)[i],nom[[i]][[1]],nom[[i]][[2]]))
    }
  }
  return(list(Col=all_point_data1,Row=all_point_data2,LP=cbind(nom$lp$x,nom$lp$x.real),TP=c( min(nom$total.points$x),max(nom$total.points$x))))
}

set.seed(17)

#oldfolder=data$oldFolder
#if(is.null(oldfolder)|oldfolder==''){
  cut.time=as.numeric(unlist(data$cutTime))
  #inputData=lung[,c(-3,-2)]
  tm=as.numeric(unlist(data$times))
  ev=as.numeric(unlist(data$events))
  vals=unlist(data$values)
  valType=as.numeric(unlist(data$colType))
  inputData=(matrix(vals,nrow=length(tm),byrow = F))
  
  #class(inputData[,3])='numeric'
  #inputData[,3]=as.numeric(inputData[,3])
  inputData=as.data.frame(inputData,stringsAsFactors =F)
  for(i in 1:ncol(inputData)){
    if(valType[i]==1){
      inputData[,i]=as.numeric(inputData[,i])
    }else{
      inputData[which(inputData[,i]==''),i]=NA
    }
  }
  logs=c(logs,paste0('获得矩阵，行：',nrow(inputData),',列：',ncol(inputData)))
  
  colnames(inputData)=paste0('C',1:ncol(inputData))
old_names=colnames(inputData)
colnames(inputData)=paste0('C',1:ncol(inputData))
fum=as.formula(paste0('Surv(time,status)~',paste0(colnames(inputData),collapse = '+')))
ndat=cbind(time=as.numeric(tm),status=as.numeric(ev),inputData)
dd <- datadist(ndat)
options(datadist="dd")
logs=c(logs,'开始进行多因素生存回归！')
md1=cph(fum,data=ndat,x=T,y=T,surv=T)
md2=survival::coxph(fum,data=ndat,x=T,y=T)
ps=summary(md2)[[7]][,5]

if(!md1$fail){
  logs=c(logs,'多因素生存回归成功！开始进行nomogram计算')
sur1=Survival(md1)
nom=nomogram(md1,fun=getFunList()
             ,fun.at = seq(0.05,0.95,by = 0.05)
             ,funlabel = paste0('P',1:length(cut.time))
)
#plot(nom)
logs=c(logs,'nomogram计算完成')
nom.data=getNomogramData(nom)
logs=c(logs,'获得nomogram数据')
old_names[match(nom.data$Col[,1],colnames(inputData))]
n.res.dat=cbind(old_names[match(nom.data$Col[,1],colnames(inputData))],
      nom.data$Col[,2:3])
nom.data$Row[,1]=gsub('P','',nom.data$Row[,1])

result<-rbind(cbind(0,t(nom.data$TP),''),
        cbind(1,nom.data$LP,''),
        cbind(2,n.res.dat),
        cbind(3,nom.data$Row),
        cbind(4,1:length(ps),ps,''))
#plot(nom)
logs=c(logs,'保存nomogram数据')
save(md1,cut.time,file = paste0(opt$outfile,'/stat.RData'))
write.table(result,file = paste0(opt$outfile,'/stat.mtx')
            ,row.names = F,col.names = F,quote = F,sep = '\t')
#dim(ndat)
logs=c(logs,'保存风险得分数据')
predict(md1)->riskscore
cindx=survcomp::concordance.index(x=riskscore[which(!is.na(riskscore))], surv.time=ndat$time[which(!is.na(riskscore))]
                            , surv.event=ndat$status[which(!is.na(riskscore))],method="noether")
logs=c(logs,paste0('C-index:',paste0(c(cindx$c.index,cindx$lower,cindx$upper,cindx$p.value),collapse = ',')))

riskscore=cbind(as.character(ndat$time),as.character(ndat$status),riskscore)
write.table(riskscore,file = paste0(opt$outfile,'/riskscore.txt')
            ,row.names = F,col.names = F,quote = F,sep = '\t')
#}else{
#  load(file = paste0(oldfolder,'/stat.RData'))
logs=c(logs,'开始计算校准曲线')
  get_calibrate=function(cox_result,u_time){
    set.seed(123)
    all_cal=c()
    points=c(2,3,4,5,6,7)
    for(m1 in floor(sum(cox_result$n)/points)){
      tryCatch({
        cal <- calibrate(cox_result, u=u_time, cmethod='KM', m=m1, B=1000)
        #plot(cal)
        se <- cal[, "std.err"]
        ciupper <- function(surv, d) ifelse(surv == 0, 0, pmin(1,surv * exp(d)))
        cilower <- function(surv, d) ifelse(surv == 0, 0, surv * exp(-d))
        cal1 <- cal[, "KM"]
        lw=cilower(cal1,1.959964 * se)
        up=ciupper(cal1,1.959964 * se)
        pred <- cal[, "mean.predicted"]
        cal.corrected <- cal[, "KM.corrected"]
        cal.dat=cbind(m1,pred,cal.corrected,cal1,lw,up)
        all_cal=rbind(all_cal,cal.dat)
      }, error = function(e) {
      }, finally = {})
    }
    return(all_cal)
  }
  all.calib=rbind()
  for(ut in cut.time){
    al=get_calibrate(md1,ut)
    if(!is.null(al)){
      al=al[which(!is.na(al[,3])),]
      all.calib=rbind(all.calib,cbind(ut,al))
    }
  }
  logs=c(logs,'获得校准曲线数据')
  write.table(all.calib,file = paste0(opt$outfile,'/calibrate.mtx')
              ,row.names = F,col.names = F,quote = F,sep = '\t')
#}
  logs=c(logs,'succ')
#get_calibrate(md1,365)
}else{
  capture.output(str(cph(fum,data=ndat,x=T,y=T,surv=T)))->putTxt
  logs=c(logs,'在多因素回归是计算失败，数据存在过拟合！')
  logs=c(logs,paste0('ERCPH:',putTxt))
}
},error = function(e) {
  print(conditionMessage(e))
  logs=c(logs,paste0('error:',conditionMessage(e)))
}, finally = {
  write.table(logs,file = paste0(opt$outfile,'/run.log'),quote = F,row.names = T,col.names = T,sep = '\t')
})









