CVbandwidth<-function(fname, lbandw=0.1, ubandw=1.2, nfold=10,out=F){
dat<-read.table(fname)
LC<-c(dat$V1)
RC<-c(dat$V2)
N<-length(LC)

bandw.list<-seq(lbandw, ubandw, by=0.1)
LKLHD<-rep(0,length(bandw.list))

#CV
for(v in 1:length(bandw.list)){
NCV<-N-nfold
CVtime<-N/nfold
  for (i in 1:CVtime){
    noind<-seq((i-1)*nfold+1, (i-1)*nfold+nfold, by=1)
    lc<-c(0, LC[-noind])
    rc<-c(0, RC[-noind])
    storage.mode(lc)<-"double"
    storage.mode(rc)<-"double"
    storage.mode(NCV)<-"integer"
    res<-.C("IcmNPMLE", x=lc, P=rc, n=NCV)
    x<-rep(res$x[2:(res$n+1)], res$P[2:(res$n+1)]*N+0.5)
    ks<-ksmooth(x, ker="parzen", bandwidth=bandw.list[v], 
                n.points=101, range.x=c(0,2))
    for (j in noind)
      LKLHD[v]<-LKLHD[v] + lklhd(ks, LC[j], RC[j])
   }
if (out==F)
   cat("LKLHD[",v,"]=", LKLHD[v],"  ")
}
if (out==T){
   maxlkl<-max(LKLHD)
   v<-bandw.list[LKLHD==maxlkl]
   sink("KS1/CVLKbandw1.1", append=T)
   cat(v, "\n")
   sink()
    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)
    x0<-seq(0,2,length=101)
    s0<-1
    sink("ICM/surv1.1", append=T)
    i<-2
    for (j in 1:101){
      if (x0[j]<res$x[i]) cat(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
              }
            cat(s0," ")
           }
           else cat(s0," ")
      }
    cat("\n")
    sink()
    x<-rep(res$x[2:(res$n+1)], res$P[2:(res$n+1)]*N+0.5)
   ks<-ksmooth(x, ker="parzen", bandwidth=v, 
                n.points=101, range.x=c(0,2))
         surv<-rep(1, length(ks$y));
         for (i in 2:length(ks$y))
             surv[i]<-surv[i-1]-(ks$y[i-1]+ks$y[i])*(ks$x[i]-ks$x[i-1])/2
   sink("KS1/CVLKsurv1.1", append=T)
   cat(surv, "\n")
   sink()
   }

#plot(ks)

}

lklhd2<-function(res, lc, rc){
ilc<-irc<-lkl<-0
for (i in 1:length(res$x))
  if (lc<=res$x[i]) {
     ilc<-i;
     break;
     }
if (ilc==0) ilc<-length(res$x)
for (i in 1:length(res$x))
  if (rc<=res$x[i]) {
     irc<-i;
     break;
     }
if (irc==0) irc<-length(res$x)
lkl<-(sum(res$y[ilc:irc])+sum(res$y[(ilc-1):(irc-1)]))*(res$x[irc]-res$x[ilc]);
if (lkl>1) cat("lkl=",lkl,"\n")
if (lkl<1e-100) return(-1000000)
else return(log(lkl))
}

lklhd<-function(res, lc, rc){
lkl<-0
lkl<-sum(res$y[(res$x>=lc)&(res$x<=rc)])*(res$x[3]-res$x[2]);
      # +sum(res$y[(ilc-1):(irc-1)]))*
if (lkl>1) {
   lkl<-1 
   # cat("lkl=",lkl,"\n")
   }
if (lkl<=0) return(-100000000)
else return(log(lkl))
}

DrawKEsurv<-function(fname, bandw=1, what="s"){
dat<-read.table(fname)
lc<-c(0,dat$V2)
rc<-c(0,dat$V3)
storage.mode(lc)<-"double"
storage.mode(rc)<-"double"
N<-length(dat$V2)
storage.mode(N)<-"integer"
res<-.C("IcmNPMLE", x=lc, P=rc, n=N)
x<-rep(res$x[2:res$n], res$P[2:res$n]*N+0.5)
  ks0<-ksmooth(x, ker="parzen", bandwidth=bandw, n.points=101, range.x=c(0,2))
if (what=="d") lines(ks0)
else if (what=="s") {
         # draw NPMLE
         res$x[1]<-0
         res$P[1]<-1
         for (i in 2:(res$n+1))
           res$P[i]<-res$P[i-1]-res$P[i]
         lines(res$x[1:(res$n+1)], res$P[1:(res$n+1)])
         surv<-rep(1, length(ks0$y));
         for (i in 2:length(ks0$y))
             surv[i]<-surv[i-1]-(ks0$y[i-1]+ks0$y[i])*(ks0$x[i]-ks0$x[i-1])/2
         lines(ks0$x, surv, lty=2)
         }

}

runCVbandw<-function(from=1, to=200){
for (i in from:to)
  CVbandwidth(paste("g",i,sep=""), out=T)
}