# Comparing two samples using the NPMLEs of their distribution functions
# Author: Wei Pan, weip@biostat.umn.edu
# INPUT: 1) datafile: name of the file containing the data
#        2) B: number of bootstrap replications
#        3) range.x: range of the survival time to be considered
NPMLE2test<-function(datafile, B=500, 
                         range.x=c(0,2), debug=F){
x0<-seq(0,range.x[2],length=101*range.x[2]/2)
  dat<-read.table(datafile)
  n<-length(dat$V2)
  nseq<-1:length(dat$V2)
  sam1<-nseq[dat$V3==0]
  sam2<-nseq[dat$V3==1]
  LC1<-dat$V1[sam1]
  RC1<-dat$V2[sam1]
  LC2<-dat$V1[sam2]
  RC2<-dat$V2[sam2]
  n1<-length(LC1)
  n2<-length(LC2)

  y1<-NPMLEvalues(LC1, RC1, n1, x0)
  y2<-NPMLEvalues(LC2, RC2, n2, x0)
  KSstat0<-max(abs(y1-y2))
  Dstat0<-sum(y1-y2)*(x0[2]-x0[1])
if (debug==T){
#   cat("stat0=",stat0,"\n")
   plot(x0, y1, type="l")
   lines(x0, y2, lty=2)
   }

  #bootstrap
  b<-1
  KSngt<-KSnlt<-0
  Dngt<-Dnlt<-0
  for (b in 1:B){
    subsam1<-sample(1:n, n1, replace=T)
    subsam2<-sample(1:n, n2, replace=T)
    lc1<-dat$V1[subsam1]
    rc1<-dat$V2[subsam1]
    lc2<-dat$V1[subsam2]
    rc2<-dat$V2[subsam2]
    y1<-NPMLEvalues(lc1, rc1, n1, x0)
    y2<-NPMLEvalues(lc2, rc2, n2, x0)
    KSstat<-max(abs(y1-y2))
    Dstat<-sum(y1-y2)*(x0[2]-x0[1])
if (debug==T){
#cat("  ", stat, " ");
plot(x0, y1, type="l")
lines(x0, y2, lty=2)
   }
    if (KSstat0<KSstat) KSnlt<-KSnlt+1
    else if (KSstat0>KSstat) KSngt<-KSngt+1
    if (Dstat0<Dstat) Dnlt<-Dnlt+1
    else if (Dstat0>Dstat) Dngt<-Dngt+1
    }

if (debug==F){
  sink(paste("NPMLE/2TestB",B,sep=""),
        append=T)
  cat(KSnlt," ", KSngt," ", KSstat0,
      " ", Dnlt," ",Dngt," ", Dstat0,"\n")
  sink()
  }
}

NPMLEvalues<-function(LC, RC, N, x0){
    lc<-c(0, LC)
    rc<-c(0, RC)
    storage.mode(lc)<-"double"
    storage.mode(rc)<-"double"
    storage.mode(N)<-"integer"
    res<-.C("IcmNPMLE", x=lc, P=rc, n=N)
    y0<-rep(0, length(x0))
    s0<-1
    i<-2
    for (j in 1:length(x0)){
      if (x0[j]<res$x[i]) y0[j]<-s0
      else if (i<res$n+1) {
            while (x0[j]>=res$x[i] && i<res$n+1) {
              s0<-s0-res$P[i]
              i<-i+1
              }
            y0[j]<-s0
           }
           else y0[j]<-s0
      }
    return(y0)
}