findnumCox <- function(target,event,time,confound,numcross,totalcut=3,initial_rr=NULL
                       ,initial_cut=NULL,initial_domain=NULL,numgen,gap=0.03){
  libraries1=c("rgenoud","survival","foreach","doParallel","doRNG","xtable")
  lapply(libraries1, library, quietly = TRUE, character.only = TRUE)
  confound <<- confound
  if(is.null(confound)){
    userdata<-na.omit(data.frame(target,event,time))
  }else{
    userdata<-na.omit(data.frame(target,event,time,confound))
    confound<-data.frame(userdata[,4:dim(userdata)[2]])
  }
  attach(userdata)
  N <<- dim(userdata)[1]
  numboot <<- sample(N)
  aictest.apply<-foreach(i = 1:1,.combine = cbind) %do% {
    attach(userdata)
    target=userdata$target
    event=userdata$event
    time=userdata$time
    if(is.null(confound)){
      counfound=NULL
    }else{confound=userdata$confound}
    mt<-replicate(numcross,sample(N))
    result <- cbind(apply(mt, 2,aictest,target=target,event=event,time=time
                          ,confound=confound,totalcut=totalcut,initial_rr=initial_rr
                          ,initial_cut=initial_cut,initial_domain=initial_domain
                          ,numgen=numgen,gap=gap))
  }
  aictest.apply <<- aictest.apply
  
  lik_cut1=NULL;hr1=NULL
  lik_cut2=NULL;hr2=NULL
  lik_cut3=NULL;hr3=NULL
  pvalue1=NULL;pvalue2=NULL;pvalue3=NULL;
  cutpvalue1=NULL;cutpvalue2=NULL;cutpvalue3=NULL;
  #N <- numcross*simulation
  RESULT <- for(i in 1:numcross){
    A=aictest.apply[[i]]
    lik_cut1[i]=aictest.apply[[i]][1,1]
    lik_cut2[i]=aictest.apply[[i]][1,2]
    lik_cut3[i]=aictest.apply[[i]][1,3]
    hr1=rbind(hr1,A[,1]$beta)
    hr2=rbind(hr2,A[,2]$beta)
    hr3=rbind(hr3,A[,3]$beta)
    pvalue1=rbind(pvalue1,A[,1]$overallstatistic.cross)
    pvalue2=rbind(pvalue2,A[,2]$overallstatistic.cross)
    pvalue3=rbind(pvalue3,A[,3]$overallstatistic.cross)
    cutpvalue1=rbind(cutpvalue1,A[,1]$cutpvalue)
    cutpvalue2=rbind(cutpvalue2,A[,2]$cutpvalue)
    cutpvalue3=rbind(cutpvalue3,A[,3]$cutpvalue)
  }
  ##AIC
  aic <- list(aic1=mean((-2)*unlist(lik_cut1)+dim(userdata)[2]-2)
              ,aic2=mean((-2)*unlist(lik_cut2)+dim(userdata)[2]-1)
              ,aic3=mean((-2)*unlist(lik_cut3)+dim(userdata)[2]))
  ##overall p-value
  overall_pvalue <- list(overallpvalue1 = sort(1-pchisq(2*pvalue1,1)),overallpvalue2 =sort(1-pchisq(2*pvalue2,2))
                         ,overallpvalue2 =sort(1-pchisq(2*pvalue3,3)))
  ##HR
  hr <- list(hr1=exp(median(hr1[,1])),hr2=c(exp(median(hr2[,1])),exp(median(hr2[,2])))
             ,hr3=c(exp(median(hr3[,1])),exp(median(hr3[,2])),exp(median(hr3[,3]))))
  ##cut p-value
  cutpvalue <- list(cutpvalue1=cutpvalue1,cutpvalue2=cutpvalue2,cutpvalue3=cutpvalue3)
  aic <<- aic
  HR <<- hr
  Overall_pvalue <<- overall_pvalue
  Cutpvalue <<- cutpvalue
  return(list(aic = aic, HR = HR, Cutpvalue = cutpvalue))
}

aictest<-function(target,event,time,confound=NULL,totalcut=totalcut,initial_rr=NULL,initial_cut=NULL,initial_domain=NULL,numgen=numgen,gap=NULL,numboot){
  confound <<- confound
  if(is.null(confound)){
    userdata<-na.omit(data.frame(target,event,time))
  }else{
    userdata<-na.omit(data.frame(target,event,time,confound))
    confound<-userdata[,4:dim(userdata)[2]]
  }
  target<-userdata$target
  censor<-userdata$event
  Ti<-userdata$time
  N<-dim(userdata)[1]
  aicfunction<-function(numcut,numgen,numboot){
    set1 <- numboot[1:round(N/2)]
    set2 <- numboot[(round(N/2)+1):N] #cross validation
    numcut <<- numcut
    
    generate1 <- maxloglik(target[set1],numcut,Ti[set1],censor[set1],confound[set1],numgen,domain_range(target[set1],initial_domain,numcut,confound[set1]),initial(target[set1],initial_rr,initial_cut,numcut),gap)
    generate2 <- maxloglik(target[set2],numcut,Ti[set2],censor[set2],confound[set2],numgen,domain_range(target[set2],initial_domain,numcut,confound[set2]),initial(target[set2],initial_rr,initial_cut,numcut),gap)
    cut1 <- generate1$par_corr
    cut2 <- generate2$par_corr
    group1 <- cut(target[set1],breaks = c(0,cut2[1:numcut],max(target[set1])),labels = c(0:numcut),quantile=FALSE)
    group2 <- cut(target[set2],breaks = c(0,cut1[1:numcut],max(target[set2])),labels = c(0:numcut),quantile=FALSE)
    if(is.null(confound)){
      fit <- coxph(Surv(time[c(set1,set2)],censor[c(set1,set2)])~factor(c(group1,group2))+strata(rep(c(1,2),c(round(N/2),N-round(N/2)))))
      overallstatistic.cross <- as.numeric(summary(fit)[["logtest"]]["test"])
    }else{
      fit <- coxph(Surv(time[c(set1,set2)],censor[c(set1,set2)])~factor(c(group1,group2))+strata(rep(c(1,2),c(round(N/2),N-round(N/2))))+as.matrix(confound[c(set1,set2)]))
      fit.re <- coxph(Surv(time[c(set1,set2)],censor[c(set1,set2)])~strata(rep(c(1,2),c(round(N/2),N-round(N/2))))+as.matrix(confound[c(set1,set2)]))
      overallstatistic.cross <- fit$loglik[2]-fit.re$loglik[2]
    }
    #Wald Test
    varcov = fit$var
    cutpvalue<-as.numeric(summary(fit)$coefficients[,5])
    beta<-as.numeric(fit$coefficients)
    return(list(aic=(-2*fit$loglik[2]+2*numcut) ,loglik=fit$loglik[2],cutpvalue=cutpvalue,beta=beta,
                cut1=cut1,cut2=cut2,overallstatistic.cross=overallstatistic.cross,
                chiWald=varcov))
  }
  catchwrongaic <- function (numcut,baby,rand) {
    out <- tryCatch(aicfunction(numcut,baby,rand), 
                    error = function(e) mget(c("aic","loglik","pvalue","beta","cut1","cut2","overall","wald"),new.env(),
                                             ifnotfound=as.list(list(c("wrong"),c("wrong"),rep(c("wrong"),numcut),
                                                                     rep(c("wrong"),numcut),
                                                                     rep(c("wrong"),2*numcut),
                                                                     c("wrong"),c("wrong"))))
    )
    return(out)
  }
  return(sapply(c(1:totalcut),aicfunction,numgen=numgen,numboot=numboot))
}

maxloglik<-function(target,numcut,time,censor,confound=NULL,baby,domain,initial,gap=NULL){
  time_global<<-time
  censor_global<<-censor
  target_global<<-target
  nvars<<-2*numcut+dim(data.frame(confound))[2]
  if(is.null(confound)){
    confound_global <<- NULL
  }else{
    confound_global<<-as.matrix(confound)
  }
  numcut_global<<-numcut
  target_max<<-max(target)
  if(is.null(gap)){
    gap_global <<- quantile(sort(diff(sort(na.omit(target)))),probs = 0.5)
  }else{
    gap_global<<-gap
  }
  ccc<-genoud(obj, nvars, max=TRUE, pop.size=100, max.generations=baby, wait.generations=10,
              hard.generation.limit=TRUE, starting.values=initial, MemoryMatrix=TRUE,
              Domains=domain, solution.tolerance=0.001,print.level = 0,
              gr=NULL, boundary.enforcement=2, lexical=FALSE, gradient.check=TRUE)
  ccc$par_corr<-ccc$par #the coefficients of genoud
  ccc$par_corr[1:numcut]<-sort(ccc$par[1:numcut])+seq(0,gap_global*(numcut-1),by=gap_global) #sort cutpoint
  return(ccc)
}

initial<-function(target,initial_rr=NULL,initial_cut=NULL,numcut){
  if(is.null(initial_cut)||is.null(initial_rr)){
    incut <- quantile(target,probs = seq(0,1,1/(numcut+1)))
    initial <- c(incut[2:(numcut+1)],initial_rr)
  }else{
    incut<-array(matrix(c(initial_cut),ncol=3,byrow=TRUE),dim=c(1,3))
    inrr<-array(matrix(c(initial_rr),ncol=3,byrow=TRUE),dim=c(1,3))
    initial <- c(as.numeric(unlist(incut[1,numcut])),as.numeric(unlist(inrr[1,numcut])))
  }
  return(initial)
}

domain_range<-function(target,initial_domain=NULL,numcut,confound=NULL){
  if(is.null(confound)){
    vars<-numcut
  }else{
    vars<-numcut+dim(data.frame(confound))[2]
  }
  if(is.null(initial_domain)){
    cutdomain <- c(quantile(target,probs = 0.1),quantile(target,probs = 0.99))
    rrdomain <- rep(c(-10,10),vars)
    ran<-if(numcut==1){
      c(1,2)
    }else if(numcut==2){
      c(1,2,1,2)
    }else{rep(1:2,times=numcut)}
    indomain <- matrix(c(cutdomain[ran],rrdomain),ncol = 2,byrow = TRUE)
  }else{
    indomain <- initial_domain[[numcut]]
  }
  return(indomain)
}

obj <- function(xx){  
  cutoff <- xx[1:numcut_global] #cutpoint
  cut_design <- cut(target_global,breaks=c(0,sort(cutoff)+seq(0,gap_global*(length(cutoff)-1),by=gap_global),target_max),quantile=FALSE,labels=c(0:numcut_global))
  beta <- xx[(numcut+1):nvars]  #coefficients of parameters
  if(is.null(confound)){
    beta1 <- coxph(Surv(time_global,censor_global)~cut_design, init = beta,iter.max=0)#iteration is zero
  }else{
    beta1 <- coxph(Surv(time_global,censor_global)~cut_design+confound_global, init = beta,iter.max=0)#iteration is zero
  }
  return(beta1$loglik[2])
}



