library(survival)
survival.prob<-function(time,fired,time.threshold)
{
   # If no data, then returen 0 so err on the side of making everyone fire.
   if (length(time)==0) 
   {
     return(0)
   }
   S<-survfit(Surv(time,fired),type="kaplan-meier")
   ntimes<-length(S$time)
   if ((ntimes==1)||(min(S$time)>=time.threshold))
   {
      return(0)
   }
   if (max(S$time)<=time.threshold)
   {
      return(S$surv[length(S$surv)])
   }
   index1<-max(which(S$time<time.threshold)) 
   index2<-index1+1
   p1<-S$surv[index1]
   p2<-S$surv[index2]
   lambda<-(time.threshold-S$time[index1])/(S$time[index2]-S$time[index1])
   p<-p1+lambda*(p2-p1)
   return(p)
}

build.quadrant.classifier<-function(x,y,implanted,days.to.firing,tx.list,ty.list,sens,time.threshold)
{
   lx<-length(x)
   fired<-rep(0,lx)
   time<-rep(0,lx)
   for (i in seq(1,lx))
   {
      if (is.na(days.to.firing[i]))
      {
         fired[i]<-0
         time[i]<-implanted[i]
      } else 
      {
         if (days.to.firing[i]<time.threshold)
         {
          fired[i]<-1
            time[i]<-days.to.firing[i]
         }else
         
         {
            fired[i]<-0
            time[i]<-implanted[i]
         }
      }
   }
   d.temp<-data.frame(x=x,y=y,time=time,fired=fired)
   d.temp<-d.temp[complete.cases(d.temp),]
   best.L<-list(spec=0)
   for (tx in tx.list)
   {
      for (ty in ty.list)
      {
        data<-list()
        data[[1]]<-d.temp[(d.temp$x<tx)&(d.temp$y<ty),]
        data[[2]]<-d.temp[(d.temp$x<tx)&(d.temp$y>ty),]
        data[[3]]<-d.temp[(d.temp$x>tx)&(d.temp$y<ty),]
        data[[4]]<-d.temp[(d.temp$x>tx)&(d.temp$y>ty),]
        p.cells<-c( 
            dim(data[[1]])[1],
            dim(data[[2]])[1],
            dim(data[[3]])[1],
            dim(data[[4]])[1])
        p.cells<-p.cells/sum(p.cells)
        p.fire<-c() 
        p.notfire<-c()
        for (i in seq(1,4)) 
        {
             time<-data[[i]]$time
             fired<-data[[i]]$fired
             p<-survival.prob(time,fired,time.threshold)
             p.notfire<-c(p.notfire,p)  
             p.fire<-c(p.fire,1-p)
             
        }
        p<-p.fire*p.cells
        p<-p/sum(p)
        q<-p.notfire*p.cells
        q<-q/sum(q)
        lr<-p/q
        perm<-order(lr, decreasing=TRUE)
        fprob<-rep(0,4)
        psum<-0
        i<-1
        while(psum+p[perm[i]]<sens) 
        {
           # If a cell is included, assign fprob for that cell 1.
           fprob[perm[i]]<-1
           # Add one more into psum
           psum<-psum+p[perm[i]]
           i<-i+1 
        }
        fprob[perm[i]]<-(sens-psum)/p[perm[i]]
        spec<-sum((1-fprob)*q)
        if (spec>best.L$spec)
            {    
               best.L<-list(fprob=fprob,spec=spec,tx=tx,ty=ty,p=p,q=q)
            }
        }
    }
    best.L
}