.packageName <- "AIM"
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/aimfuns-supp.R"
backfit.lm.main=function(x, y, cutp, mincut=0)
{n=length(y)
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}


 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x[,i])}


 muhat=mean(y)
 sigma=sum((y-muhat)^2)/(n-1)

 y.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 
 score.mat=t((t(x)<cutp)*1)
 score.tot=apply(score.mat,1,sum)       

 zscore=sum((y-muhat)*score.tot)/sqrt(sigma*(sum(score.tot^2)-(sum(score.tot))^2/n))
     
 count=1
 loop=1
 while(loop==1)
      {
       score.mat=t((t(x)<cutp)*1)
       score.tot=apply(score.mat,1,sum)       
   
       score.pool=muhat.pool=matrix(0, n, p)
       for(i in 1:p)
          {idx=order.index[,i] 
           score.pool[,i]=(score.tot-score.mat[,i])[idx]
           muhat.pool[,i]=muhat[idx]}
       
       score.stat=t(t(apply(y.pool-muhat,2,cumsum))+apply((y.pool-muhat)*score.pool,2,sum))
       v.stat=sigma*(t(t(2*apply(score.pool,2,cumsum)+1:n)+apply(score.pool^2,2,sum))-(t(matrix(rep(apply(score.pool,2,sum),n),p,n))+1:n)^2/n)

       for(i in 1:p)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p)
          for(i in 1:p)
             {
              test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }
       loop=0
       if(max(mtest)>zscore[count]+1e-8)
         {loop=1
          i=(1:p)[mtest==max(mtest)][1]
          i0=max(effect.range[test[effect.range,i]==mtest[i]])+1
          if(i0>n)
            mcut=Inf
          if(i0<=n)
            mcut=x[order.index[,i],i][i0]
          cutp[i]=mcut
          zscore=c(zscore, max(mtest)) 
          count=count+1           
         }
        
       }
       
       
  return(list(cutp=cutp, zscore=zscore))
 }

############################################################################################################################################
############################################################################################################################################
############################################################################################################################################

backfit.lm.interaction=function(x, trt, y, cutp, mincut=0)
{n=length(y)
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}


 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x[,i])}

 fit=lm(y~trt)
 muhat=fit$fitt
 sigma=sum((y-muhat)^2)/(n-2)
 sigma.sum=sum(solve(t(cbind(1,trt))%*%cbind(1,trt))*sigma)

 y.pool=trt.pool=muhat.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     trt.pool[,i]=trt[order.index[,i]]
     muhat.pool[,i]=muhat[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 
 score.mat=t((t(x)<cutp)*1)
 score.tot=apply(score.mat,1,sum)       
 zscore=sum((y-muhat)*score.tot*trt)/sqrt(sigma*sum(score.tot^2*trt)-(sum(score.tot*trt))^2*sigma.sum)
     
 count=1
 loop=1
 while(loop==1)
      {
       score.mat=t((t(x)<cutp)*1)
       score.tot=apply(score.mat,1,sum)       
   
       score.pool=matrix(0, n, p)
       for(i in 1:p)
          {idx=order.index[,i] 
           score.pool[,i]=(score.tot-score.mat[,i])[idx]}

       score.stat=t(t(apply((y.pool-muhat.pool)*trt.pool,2,cumsum))+apply((y.pool-muhat.pool)*score.pool*trt.pool,2,sum))
       v11.stat=sigma*t(t(apply((2*score.pool+1)*trt.pool,2,cumsum))+apply(score.pool^2*trt.pool,2,sum))
       v22.stat=t((apply(score.pool*trt.pool,2,sum)+t(apply(trt.pool,2,cumsum)))^2)*sigma.sum
       v.stat=(v11.stat-v22.stat)


       for(i in 1:p)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p)
          for(i in 1:p)
             {
              test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }
       loop=0
       if(max(mtest)>zscore[count]+1e-8)
         {loop=1
          i=(1:p)[mtest==max(mtest)][1]
          i0=max(effect.range[test[effect.range,i]==mtest[i]])+1
          if(i0>n)
            mcut=Inf
          if(i0<=n)
            mcut=x[order.index[,i],i][i0]
          cutp[i]=mcut
          zscore=c(zscore, max(mtest)) 
          count=count+1           
         }
        
       }
       
       
  return(list(cutp=cutp, zscore=zscore))
 }

############################################################################################################################################
############################################################################################################################################
############################################################################################################################################

backfit.logistic.main=function(x, y, cutp, mincut=0, wt=wt)
{n=length(y)
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}


 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x[,i])}

 wt.pool=y.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     wt.pool[,i]=wt[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 phat=sum(y*wt)/sum(wt)
 sigma=as.numeric(summary(glm(y~1, family="binomial", weight=wt))$cov.un*(phat*(1-phat))^2*n)

 score.mat=t((t(x)<cutp)*1)
 score.tot=apply(score.mat,1,sum)       

 zscore=sum((y-phat)*score.tot)/sqrt(phat*(1-phat)*(sum(score.tot^2)-(sum(score.tot))^2/n))
     
 count=1
 loop=1
 while(loop==1)
      {
       score.mat=t((t(x)<cutp)*1)
       score.tot=apply(score.mat,1,sum)       
   
       score.pool=matrix(0, n, p)
       for(i in 1:p)
          {idx=order.index[,i] 
           score.pool[,i]=(score.tot-score.mat[,i])[idx]}

       score.stat=t(t(apply((y.pool-phat)*wt.pool,2,cumsum))+apply((y.pool-phat)*(score.pool*wt.pool),2,sum))
       v.stat=sigma*(t(t(2*apply(score.pool*wt.pool^2,2,cumsum)+apply(wt.pool^2,2,cumsum))+apply(score.pool^2*wt.pool^2,2,sum))-(t(matrix(rep(apply(score.pool*wt.pool,2,sum),n),p,n))+apply(wt.pool,2,cumsum))^2/n)

       for(i in 1:p)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p)
          for(i in 1:p)
             {
              test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }
       loop=0
       if(max(mtest)>zscore[count]+1e-8)
         {loop=1
          i=(1:p)[mtest==max(mtest)][1]
          i0=max(effect.range[test[effect.range,i]==mtest[i]])+1
          if(i0>n)
            mcut=Inf
          if(i0<=n)
            mcut=x[order.index[,i],i][i0]
          cutp[i]=mcut
          zscore=c(zscore, max(mtest)) 
          count=count+1           
         }
        
       }
       
       
  return(list(cutp=cutp, zscore=zscore))
 }

############################################################################################################################################
############################################################################################################################################
############################################################################################################################################

backfit.logistic.interaction=function(x, trt, y, cutp, mincut=0, wt=wt)
{n=length(y)
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}


 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x[,i])}

 wt.pool=trt.pool=y.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     wt.pool[,i]=wt[order.index[,i]]
     trt.pool[,i]=trt[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 fit=glm(y~trt, family="binomial", weight=wt)
 phat=fit$fitt
 sigma=sum(summary(fit)$cov.scale)

 score.mat=t((t(x)<cutp)*1)
 score.tot=apply(score.mat,1,sum)       

 zscore=sum((y-phat)*trt*score.tot)/sqrt(sum(score.tot^2*trt*phat*(1-phat))-(sum(score.tot*trt*phat*(1-phat)))^2*sigma)
     
 count=1
 loop=1
 while(loop==1)
      {score.mat=t((t(x)<cutp)*1)
       score.tot=apply(score.mat,1,sum)       
   
       score.pool=phat.pool=matrix(0, n, p)
       for(i in 1:p)
          {idx=order.index[,i] 
           score.pool[,i]=(score.tot-score.mat[,i])[idx]
           phat.pool[,i]=phat[idx]}

       score.stat=t(t(apply(trt.pool*wt.pool*(y.pool-phat.pool),2,cumsum))+apply((y.pool-phat.pool)*trt.pool*score.pool*wt.pool,2,sum))
       v11.stat=t(t(2*apply(score.pool*trt.pool*wt.pool^2*phat.pool*(1-phat.pool),2,cumsum)+apply(trt.pool*wt.pool^2*phat.pool*(1-phat.pool),2,cumsum))+apply(score.pool^2*wt.pool^2*trt.pool*phat.pool*(1-phat.pool),2,sum))
       v22.stat=(t(t(apply(trt.pool*wt.pool*phat.pool*(1-phat.pool),2,cumsum))+apply(score.pool*trt.pool*wt.pool*phat.pool*(1-phat.pool),2,sum)))^2*sigma
       v.stat=(v11.stat-v22.stat)

       for(i in 1:p)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p)
          for(i in 1:p)
             {
              test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }
       loop=0
       if(max(mtest)>zscore[count]+1e-8)
         {loop=1
          i=(1:p)[mtest==max(mtest)][1]
          i0=max(effect.range[test[effect.range,i]==mtest[i]])+1
          if(i0>n)
            mcut=Inf
          if(i0<=n)
            mcut=x[order.index[,i],i][i0]
          cutp[i]=mcut
          zscore=c(zscore, max(mtest)) 
          count=count+1           
         }
        
       }
       
       
  return(list(cutp=cutp, zscore=zscore))
 }



############################################################################################################################################
############################################################################################################################################
############################################################################################################################################

backfit.cox.main=function(x, y, delta, cutp, mincut=0)
{n=length(y)
 p=length(x[1,])
 id=order(y)
 y=y[id]
 delta=delta[id]
 x=x[id,]


 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}


 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x[,i])}

 delta.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {delta.pool[,i]=delta[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 
 ntime=length(unique(y))
 time.index=rep(c(1, cumsum(table(y))+1)[-(ntime+1)], table(y))
 num.risk=(n:1)[time.index]
 time.index.plus=rep(cumsum(table(y)), table(y))

 risk.set=matrix(0, n, n)
 for(i in 1:n)
    {risk.set[(n-num.risk[i]+1):n,i]=1}
 
 score.mat=t((t(x)<cutp)*1)
 score.tot=apply(score.mat,1,sum)       
 score.bar=rev(cumsum(rev(score.tot)))/(n:1)

 u0.stat=sum(delta*(score.tot-score.bar[time.index]))
 i0.stat=sum(delta/num.risk*((rev(cumsum(rev(score.tot^2)))-score.bar^2*(n:1)))[time.index])
 zscore=u0.stat/sqrt(i0.stat+1e-8)
 

 u.pool=i1.pool=matrix(0, n, p)
 for(i in 1:p)
    {idx=order.index[,i]
     idy=order(idx)

     u.pool[,i]=(cumsum(delta/num.risk))[time.index.plus][idx]

     risk.set.sub=risk.set[idx,]
     w=apply(risk.set.sub,2,cumsum)
     w=rbind(0, w[-n,])
     w=t(w[idy,])

     i1.pool[,i]=(cumsum(delta*(num.risk-1)/num.risk^2))[time.index.plus][idx]-2*(diag((apply(delta/num.risk^2*w,2,cumsum))[time.index.plus,]))[idx]
     }


 count=1
 loop=1
 while(loop==1)
      {score.mat=t((t(x)<cutp)*1)
       score.tot=apply(score.mat,1,sum)       
   
       score.pool=score.bar.pool=matrix(0, n, p)
       for(i in 1:p)
          {idx=order.index[,i] 
           score.pool[,i]=(score.tot-score.mat[,i])
           score.bar.pool[,i]=rev(cumsum(rev(score.tot-score.mat[,i])))/(n:1)
           }
       
       i.pool=matrix(0, n, p)
       for(i in 1:p)
           {
            idx=order.index[,i]
            i.pool[,i]=2*u.pool[,i]*score.pool[idx,i]-2*(cumsum(delta/num.risk*score.bar.pool[time.index,i]))[time.index.plus][idx]
            } 
       i.pool=i.pool+i1.pool

       u0.stat=apply(delta*(score.pool-score.bar.pool[time.index,]),2,sum)
       i0.stat=apply(delta*(1/(n:1)*((apply((score.pool^2)[n:1,],2,cumsum))[n:1,]-score.bar.pool^2*(n:1)))[time.index,],2,sum)
       
    
       score.stat=t(u0.stat+t(apply(delta.pool-u.pool,2,cumsum)))
       v.stat=t(i0.stat+t(apply(i.pool,2,cumsum)))

       for(i in 1:p)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p)
          for(i in 1:p)
             {
              test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)
             }
          }
       loop=0
       if(max(mtest)>zscore[count]+1e-8)
         {loop=1
          i=(1:p)[mtest==max(mtest)][1]
          i0=max(effect.range[test[effect.range,i]==mtest[i]])+1
          if(i0>n)
            mcut=Inf
          if(i0<=n)
            mcut=x[order.index[,i],i][i0]
          cutp[i]=mcut
          zscore=c(zscore, max(mtest)) 
          count=count+1           
         }        
       }
       
       
  return(list(cutp=cutp, zscore=zscore))
 }



############################################################################################################################################
############################################################################################################################################
############################################################################################################################################

backfit.cox.interaction=function(x, trt, y, delta, cutp, mincut=0)
{n=length(y)
 p=length(x[1,])
 id=order(y)
 y=y[id]
 delta=delta[id]
 trt=trt[id]
 x=x[id,]

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}

 ntime=length(unique(y))
 time.index=rep(c(1, cumsum(table(y))+1)[-(ntime+1)], table(y))
 num.risk=(n:1)[time.index]
 time.index.plus=rep(cumsum(table(y)), table(y))


 fit=coxph(Surv(y, delta)~trt, method="breslow")
 beta=fit$coef
 sigma0=fit$var
 risk=exp(beta*trt)

 S0=(rev(cumsum(rev(risk))))[time.index]
 Sr=(rev(cumsum(rev(risk*trt))))[time.index]


 order.index=risk.pool=delta.pool=trt.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {idx=order(x[,i])
     order.index[,i]=idx
     risk.pool[,i]=risk[idx]
     delta.pool[,i]=delta[idx]
     trt.pool[,i]=trt[idx]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))
     }

 ds0=(cumsum(delta/S0))[time.index.plus]
 drs0=(cumsum(delta*Sr/S0^2))[time.index.plus]
 dss0=(cumsum(delta/S0^2))[time.index.plus]

 ds0.pool=drs0.pool=dss0.pool=matrix(0, n, p)
 for(i in 1:p)
    {idx=order(x[,i])
     ds0.pool[,i]=ds0[idx]
     dss0.pool[,i]=dss0[idx]
     drs0.pool[,i]=drs0[idx]
     }

 score.mat=t((t(x)<cutp)*1)
 score.tot=apply(score.mat,1,sum)       

 Swwrr0=(rev(cumsum(rev(risk*trt^2*score.tot^2))))[time.index]
 Swrr0=(rev(cumsum(rev(risk*trt^2*score.tot))))[time.index]
 Swr0=(rev(cumsum(rev(risk*trt*score.tot))))[time.index]

 i1=sum(delta*Swwrr0/S0)
 i2=sum(delta*Swr0^2/S0^2)
 i3=sum(delta*Swrr0/S0)
 i4=sum(delta*Swr0*Sr/S0^2)
 v.stat=(i1-i2)-(i3-i4)^2*sigma0[1,1]
 score.stat=sum(delta*score.tot*trt)-sum(delta*Swr0/S0)
 zscore=score.stat/sqrt(v.stat+1e-8)
 
 
 subject=matrix(0, n, n)
 for(i in 1:n)
    {subject[(n-num.risk[i]+1):n,i]=1}
 
 i2.add=apply(trt.pool^2*risk.pool^2*dss0.pool,2,cumsum)
 i3.pool=apply(trt.pool^2*risk.pool*ds0.pool,2,cumsum)
 i4.pool=apply(trt.pool*risk.pool*drs0.pool,2,cumsum)

 w1.pool=matrix(0, n, p)
 for(i in 1:p)
    {idx=order.index[,i]
     subject.id=subject[idx,]
     temp1=apply(risk[idx]*trt[idx]*subject.id,2,cumsum)
     temp1=rbind(0, temp1[-n,])    
     w1.pool[,i]=diag((apply(t(temp1)*delta/S0^2,2,cumsum))[time.index.plus,][idx,])
     }

 score.stat.pool=apply(delta.pool*trt.pool,2,cumsum)-apply(trt.pool*risk.pool*ds0.pool,2,cumsum)

 count=1
 loop=1
 while(loop==1)
       {score.mat=t((t(x)<cutp)*1)
        score.tot=apply(score.mat,1,sum)       
   
             
        score.pool=score.pool0=matrix(0, n, p)
        for(i in 1:p)
           {idx=order.index[,i]
            score.pool[,i]=(score.tot-score.mat[,i])[idx] 
            score.pool0[,i]=(score.tot-score.mat[,i])}


        Swwrr0=((apply((risk*trt^2*score.pool0^2)[n:1,],2,cumsum))[n:1,])[time.index,]
        Swrr0=((apply((risk*trt^2*score.pool0)[n:1,],2,cumsum))[n:1,])[time.index,]
        Swr0=((apply((risk*trt*score.pool0)[n:1,],2,cumsum))[n:1,])[time.index,]


        w.pool=matrix(0, n, p)
        for(i in 1:p)
           {idx=order.index[,i]  
            temp2=(cumsum((rev(cumsum(rev(risk*trt*score.pool0[,i]))))[time.index]*delta/S0^2))[time.index.plus]
            w.pool[,i]=w1.pool[,i]+temp2[idx]
            }

        i10=apply(delta*Swwrr0/S0,2,sum)
        i20=apply(delta*Swr0^2/S0^2,2,sum)
        i30=apply(delta*Swrr0/S0,2,sum)
        i40=apply(delta*Swr0*Sr/S0^2,2,sum)   

        i1.pool=apply((2*score.pool+1)*trt.pool^2*risk.pool*ds0.pool,2,cumsum)
        i2.pool=apply(2*trt.pool*risk.pool*w.pool,2,cumsum)+i2.add
       
        i1=t(t(i1.pool)+i10)
        i2=t(t(i2.pool)+i20)
        i3=t(t(i3.pool)+i30)
        i4=t(t(i4.pool)+i40)

        v.stat=(i1-i2)-(i3-i4)^2*sigma0[1,1]

        score.stat0=apply(delta.pool*score.pool*trt.pool,2,sum)-apply(delta*Swr0/S0,2,sum)
        score.stat=t(t(score.stat.pool)+score.stat0)

        for(i in 1:p)
            {idx=x.replicate[,i]
             v.stat[,i]=v.stat[idx ,i]   
             score.stat[,i]=score.stat[idx, i]
            } 

         test=score.stat/sqrt(v.stat+1e-8)

         if(mincut==0)
           {mtest=apply(test,2,max)}

         if(mincut>0)
           {mtest=rep(0, p)
            for(i in 1:p)
               {
                test.post=test[max(effect.range)+1,i]
                test0=test[effect.range,i]
                test0=test0[test0!=test.post]
                mtest[i]=max(test0)
             }
            }
         loop=0
         if(max(mtest)>zscore[count]+1e-8)
           {loop=1
            i=(1:p)[mtest==max(mtest)][1]
            i0=max(effect.range[test[effect.range,i]==mtest[i]])+1
            if(i0>n)
               mcut=Inf
            if(i0<=n)
              mcut=x[order.index[,i],i][i0]
            cutp[i]=mcut
            zscore=c(zscore, max(mtest)) 
            count=count+1           
            }        
          }
       
       
     return(list(cutp=cutp, zscore=zscore))
    }
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cox.interaction.R"
cox.interaction=function(x, trt, y, delta, nsteps=8, mincut=0.1, backfit=F, maxnumcut=1, dirp=0)
{
n=length(y)
id=order(y)
y=y[id]
x=x[id,]
delta=delta[id]
trt=trt[id]

x0=x
p.true=length(x0[1,])

x=cbind(x0, -x0)
marker.x=rep(1:p.true, 2)
direction.x=c(rep(-1, p.true), rep(1, p.true))

if(sum(abs(dirp))>0)
    {index1=(1:p.true)[dirp==-1 | dirp==0]
     index2=(1:p.true)[dirp==1 | dirp==0]
     x=cbind(x0[,index1], x0[,index2])
     marker.x=c((1:p.true)[index1], (1:p.true)[index2])
     direction.x=c(rep(-1, length(index1)), rep(1, length(index2)))
     }                   
p=length(x[1,])
 
if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
if(mincut==0)
   {effect.range=1:n}

score.current=rep(0,n)
x.out=x

ntime=length(unique(y))
time.index=rep(c(1, cumsum(table(y))+1)[-(ntime+1)], table(y))
num.risk=(n:1)[time.index]
time.index.plus=rep(cumsum(table(y)), table(y))


id.in=NULL
id.out=1:p.true
p.out=p
zvalue=NULL
cut.value=NULL
direction=NULL
imax=NULL
res=as.list(NULL)
num.cut=rep(0, p.true)

fit=coxph(Surv(y, delta)~trt, method="breslow")
beta=fit$coef
sigma0=fit$var
risk=exp(beta*trt)

S0=(rev(cumsum(rev(risk))))[time.index]
Sr=(rev(cumsum(rev(risk*trt))))[time.index]


order.index=risk.pool=delta.pool=trt.pool=x.replicate=matrix(0, n, p)
for(i in 1:p)
   {idx=order(x[,i])
    order.index[,i]=idx
    risk.pool[,i]=risk[idx]
    delta.pool[,i]=delta[idx]
    trt.pool[,i]=trt[idx]
    x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))
    }

ds0=(cumsum(delta/S0))[time.index.plus]
drs0=(cumsum(delta*Sr/S0^2))[time.index.plus]
dss0=(cumsum(delta/S0^2))[time.index.plus]

ds0.pool=drs0.pool=dss0.pool=matrix(0, n, p)
for(i in 1:p)
   {idx=order(x[,i])
    ds0.pool[,i]=ds0[idx]
    dss0.pool[,i]=dss0[idx]
    drs0.pool[,i]=drs0[idx]
    }


subject=matrix(0, n, n)
for(i in 1:n)
    {subject[(n-num.risk[i]+1):n,i]=1}


i2.diff.pool=apply(trt.pool^2*risk.pool^2*dss0.pool,2,cumsum)
i3.diff.pool=apply(trt.pool^2*risk.pool*ds0.pool,2,cumsum)
i4.diff.pool=apply(trt.pool*risk.pool*drs0.pool,2,cumsum)
score.diff.pool=apply(delta.pool*trt.pool,2,cumsum)-apply(trt.pool*risk.pool*ds0.pool,2,cumsum)

w1.pool=matrix(0, n, p.out)
for(i in 1:p.out)
   {idx=order.index[,i]
    subject.id=subject[idx,]
    temp1=apply(risk[idx]*trt[idx]*subject.id,2,cumsum)
    temp1=rbind(0, temp1[-n,])
    w1.pool[,i]=diag((apply(t(temp1)*delta/S0^2,2,cumsum))[time.index.plus,][idx,])
    }
     
count=1
loop=1
while(loop==1){


Swwrr0=(rev(cumsum(rev(risk*trt^2*score.current^2))))[time.index]
Swrr0=(rev(cumsum(rev(risk*trt^2*score.current))))[time.index]
Swr0=(rev(cumsum(rev(risk*trt*score.current))))[time.index]


score.pool=w.pool=matrix(0, n, p.out)
for(i in 1:p.out)
    {idx=order.index[,i]   
     score.pool[,i]=score.current[idx]
     temp=(cumsum((rev(cumsum(rev(risk*trt*score.current))))[time.index]*delta/S0^2))[time.index.plus]
     w.pool[,i]=w1.pool[,i]+temp[idx]
     }
     

i1=sum(delta*Swwrr0/S0)+apply((2*score.pool+1)*trt.pool^2*risk.pool*ds0.pool,2,cumsum)
i2=sum(delta*Swr0^2/S0^2)+apply(2*trt.pool*risk.pool*w.pool,2,cumsum)+i2.diff.pool
i3=sum(delta*Swrr0/S0)+i3.diff.pool
i4=sum(delta*Swr0*Sr/S0^2)+i4.diff.pool


v.stat=(i1-i2)-(i3-i4)^2*sigma0[1,1]
score.stat=sum(delta*score.current*trt)-sum(delta*Swr0/S0)+score.diff.pool

for(i in 1:p.out)
   {idx=x.replicate[,i]
    v.stat[,i]=v.stat[idx ,i]   
    score.stat[,i]=score.stat[idx, i]
    } 

test=score.stat/sqrt(v.stat+1e-8)

if(mincut==0)
  {mtest=apply(test,2,max)}

if(mincut>0)
  {mtest=rep(0, p.out)
   for(i in 1:p.out)
      {
       test.post=test[max(effect.range)+1,i]
       test0=test[effect.range,i]
       test0=test0[test0!=test.post]
       mtest[i]=max(test0)}
       }


mcut=rep(0, p.out)
i0=rep(0, p.out)
for(i in 1:p.out)
   {i0[i]=max(effect.range[test[effect.range,i]==mtest[i]])+1
    if(i0[i]>n)
       mcut[i]=Inf
    if(i0[i]<=n)
       mcut[i]=x.out[order.index[,i],i][i0[i]]
    }
         

i.sel=which.max(mtest)
score.current=score.current+(x.out[,i.sel]<mcut[i.sel])

marker.sel=marker.x[i.sel]
id.in=c(id.in, marker.sel)
direction=c(direction, direction.x[i.sel])
cut.value=c(cut.value, -mcut[i.sel]*direction.x[i.sel])
num.cut[marker.sel]=num.cut[marker.sel]+1


if(num.cut[marker.sel]==maxnumcut)
  {id.exclude=(1:p.out)[marker.x==marker.sel]

   x.out=x.out[,-id.exclude, drop=F]
   order.index=order.index[, -id.exclude, drop=F]
   trt.pool=trt.pool[, -id.exclude, drop=F]
   risk.pool=risk.pool[, -id.exclude, drop=F]
   ds0.pool=ds0.pool[, -id.exclude, drop=F]
   i2.diff.pool=i2.diff.pool[, -id.exclude, drop=F]
   i3.diff.pool=i3.diff.pool[, -id.exclude, drop=F]
   i4.diff.pool=i4.diff.pool[, -id.exclude, drop=F]
   score.diff.pool=score.diff.pool[, -id.exclude, drop=F]
   w1.pool=w1.pool[, -id.exclude, drop=F]
   x.replicate=x.replicate[, -id.exclude, drop=F]
   
   direction.x=direction.x[-id.exclude]
   marker.x=marker.x[-id.exclude]   
   p.out=p.out-length(id.exclude)                  
   }

       
if(backfit==T)
  {if(count>1)
     {x.adj=x0[, id.in]
      x.adj=-t(t(x.adj)*direction)
      cutp=-cut.value*direction
      fit=backfit.cox.interaction(x.adj, trt, y, delta, cutp, mincut=mincut)
         
      jmax=id.in
      cutp=-fit$cutp*direction
      maxdir=direction
             
      res[[count]]=cbind(jmax, cutp, maxdir)
      zvalue=c(zvalue,  max(fit$zscore)) 
      score.current=apply((t(x.adj)<fit$cutp),2,sum)
      }

    if(count==1)
      {jmax=id.in
       cutp=cut.value
       maxdir=direction
       res[[count]]=cbind(jmax, cutp, maxdir)
       zvalue=c(zvalue,  max(mtest))
       }
    }
if(backfit==F)
  {cutp=cut.value
   maxdir=direction
   jmax=id.in
   zvalue=c(zvalue, max(mtest))
   maxsc=zvalue
   res[[count]]=cbind(jmax, cutp, maxdir, maxsc)
   }


count=count+1
loop=(length(id.in)<nsteps)
}

       
return(list(res=res, maxsc=zvalue))
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cox.main.R"

cox.main=function(x, y, delta, nsteps=8, mincut=0.1, backfit=F, maxnumcut=1, dirp=0)
{
 n=length(x[,1])
 id=order(y)
 y=y[id]
 delta=delta[id]
 x=x[id,]

 
 x0=x
 p.true=length(x0[1,])
 
 x=cbind(x0, -x0)
 marker.x=rep(1:p.true, 2)
 direction.x=c(rep(-1, p.true), rep(1, p.true))

 if(sum(abs(dirp))>0)
    {index1=(1:p.true)[dirp==-1 | dirp==0]
     index2=(1:p.true)[dirp==1 | dirp==0]
     x=cbind(x0[,index1], x0[,index2])
     marker.x=c((1:p.true)[index1], (1:p.true)[index2])
     direction.x=c(rep(-1, length(index1)), rep(1, length(index2)))
     }                   
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}

 score.current=rep(0,n)
 x.out=x

 
 order.index=delta.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
      {idx=order(x[,i])
       order.index[,i]=idx
       delta.pool[,i]=delta[idx]
       x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))
       } 
 
 id.in=NULL
 id.out=1:p.true
 p.out=p
 zvalue=NULL
 cut.value=NULL
 direction=NULL
 imax=NULL
 res=as.list(NULL)
 num.cut=rep(0, p.true)

 tby=table(y)
 ntime=length(unique(y))
 time.index=rep(c(1, cumsum(tby)+1)[-(ntime+1)], tby)
 num.risk=(n:1)[time.index]
 time.index.plus=rep(cumsum(tby), tby)


 risk.set=matrix(0, n, n)
 for(i in 1:n)
    {risk.set[(n-num.risk[i]+1):n,i]=1}


 
 i1.pool=i3.pool=u.pool=matrix(0, n, p.out)
 u0=(cumsum(delta/num.risk))[time.index.plus]
 i10=(cumsum(delta*(num.risk-1)/num.risk^2))[time.index.plus]
 for(i in 1:p.out)
    {idx=order.index[,i]
     idy=order(idx)

     u.pool[,i]=u0[idx]

     i1.pool[,i]=i10[idx]
      
     risk.set.sub=risk.set[idx,]
     w=apply(risk.set.sub,2,cumsum)
     w=rbind(0, w[-n,])
     w=t(w[idy,])
     i3.pool[,i]=2*(diag((apply(delta/num.risk^2*w,2,cumsum))[time.index.plus,]))[idx]
     }

 
 count=1
 loop=1
 while(loop==1)
  {score.bar=(rev(cumsum(rev(score.current)))/(n:1))
   u0.stat=sum(delta*(score.current-score.bar[time.index]))
   i0.stat=sum(delta/num.risk*((rev(cumsum(rev(score.current^2)))-score.bar^2*(n:1)))[time.index])

 
   i.pool=matrix(0, n, p.out)
   for(i in 1:p.out)
      {
       idx=order.index[,i]                
       i.pool[,i]=i1.pool[,i]+2*u.pool[,i]*score.current[idx]-2*(cumsum(delta/num.risk*score.bar[time.index]))[time.index.plus][idx]-i3.pool[,i]
       } 


     score.stat=u0.stat+apply(delta.pool-u.pool,2,cumsum)
     v.stat=i0.stat+apply(i.pool,2,cumsum)

     for(i in 1:p.out)
        {idx=x.replicate[,i]
         v.stat[,i]=v.stat[idx ,i]   
         score.stat[,i]=score.stat[idx, i]
         } 

     test=score.stat/sqrt(v.stat+1e-8)

     if(mincut==0)
        {mtest=apply(test,2,max)}

     if(mincut>0)
        {mtest=rep(0, p.out)
         for(i in 1:p.out)
             {test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
         }


     mcut=rep(0, p.out)
     i0=rep(0, p.out)
     for(i in 1:p.out)
       {i0[i]=max(effect.range[test[effect.range,i]==mtest[i]])+1
        if(i0[i]>n)
          mcut[i]=Inf
        if(i0[i]<=n)
          mcut[i]=x.out[order.index[,i],i][i0[i]]
       }
         

     i.sel=(1:p.out)[mtest==max(mtest)][1]
     score.current=score.current+(x.out[,i.sel]<mcut[i.sel])

     marker.sel=marker.x[i.sel]
     id.in=c(id.in, marker.sel)
     direction=c(direction, direction.x[i.sel])
     cut.value=c(cut.value, -mcut[i.sel]*direction.x[i.sel])
     num.cut[marker.sel]=num.cut[marker.sel]+1

     if(num.cut[marker.sel]==maxnumcut)
       {id.exclude=(1:p.out)[marker.x==marker.sel]
   
        x.out=x.out[,-id.exclude,drop=F]
        order.index=order.index[, -id.exclude,drop=F]
        delta.pool=delta.pool[,-id.exclude,drop=F]
        x.replicate=x.replicate[, -id.exclude,drop=F]
        i1.pool=i1.pool[,-id.exclude,drop=F]
        i3.pool=i3.pool[,-id.exclude,drop=F]
        u.pool=u.pool[,-id.exclude,drop=F]
        direction.x=direction.x[-id.exclude]
        marker.x=marker.x[-id.exclude]   
        p.out=p.out-length(id.exclude)                  
        }
          

       
     if(backfit==T)
         {if(count>1)
            {x.adj=x0[, id.in]
             x.adj=-t(t(x.adj)*direction)
             cutp=-cut.value*direction
             fit=backfit.cox.main(x.adj, y, delta, cutp, mincut=mincut)
         
             jmax=id.in
             cutp=-fit$cutp*direction
             maxdir=direction
             
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(fit$zscore)) 
             score.current=apply((t(x.adj)<fit$cutp),2,sum)
             }

          if(count==1)
            {jmax=id.in
             cutp=cut.value
             maxdir=direction
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(mtest))
             }
           }
      if(backfit==F)
         {cutp=cut.value
          maxdir=direction
          jmax=id.in
          zvalue=c(zvalue, max(mtest))
          maxsc=zvalue
          res[[count]]=cbind(jmax, cutp, maxdir, maxsc)
          }

       
       count=count+1
       loop=(length(id.in)<nsteps)
       }

       
 return(list(res=res, maxsc=zvalue))
 }
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cv.cox.interaction.R"
cv.cox.interaction=function(x, trt, y, status, K.cv=5, num.replicate=1, nsteps, mincut=0.1, backfit=F, maxnumcut=1, dirp=0){
x=as.matrix(x)
n=length(y)


sc.tot=matrix(0, nrow=K.cv, ncol=nsteps)
preval.tot=matrix(0, nrow=nsteps, ncol=n)

for(b in 1:num.replicate)
   {g=sample(rep(1:K.cv,(n+K.cv-1)/K.cv))
    g=g[1:n]
    cva=vector("list", K.cv)

    

    sc=matrix(0, nrow=K.cv, ncol=nsteps)
    preval=matrix(0, nrow=nsteps, ncol=n)
    pv.fit=rep(0, nsteps)
    for(i in 1:K.cv)
       {cva[[i]]=cox.interaction(x[g!=i ,], trt[g!=i], y[g!=i], status[g!=i], nsteps=nsteps, mincut=mincut, backfit=backfit, maxnumcut=maxnumcut, dirp=dirp)
        for(ii in 1:nsteps)
           {aa=index.prediction(cva[[i]]$res[[ii]],x[g==i,])
            fit=coxph(Surv(y[g==i],status[g==i])~aa*trt[g==i])
            sc[i,ii]=summary(fit)$coef[3,4]
            preval[ii,g==i]=aa
           }
        }
   

    sc.tot=abs(sc.tot)+sc
    preval.tot=preval.tot+preval
    }

meansc=colMeans(sc.tot)
kmax=which.max(meansc)
pvfit.score=rep(0, nsteps)
for(i in 1:nsteps)
   {fit=coxph(Surv(y,status)~preval.tot[i,]*trt)
    pvfit.score[i]=summary(fit)$coef[3,4]
    }


return(list(kmax=kmax, meanscore=meansc/num.replicate, pvfit.score=pvfit.score, preval=preval.tot/num.replicate))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cv.cox.main.R"
cv.cox.main=function(x, y, status, K.cv=5, num.replicate=1, nsteps, mincut=0.1, backfit=F, maxnumcut=1, dirp=0){
x=as.matrix(x)
n=length(y)

sc.tot=matrix(0, nrow=K.cv, ncol=nsteps)
preval.tot=matrix(0, nrow=nsteps, ncol=n)

for(b in 1:num.replicate)
   {g=sample(rep(1:K.cv,(n+K.cv-1)/K.cv))
    g=g[1:n]
    cva=vector("list",K.cv)

    sc=matrix(0, nrow=K.cv, ncol=nsteps)
    preval=matrix(0, nrow=nsteps, ncol=n)
    for(i in 1:K.cv)
       {cva[[i]]=cox.main(x[g!=i ,],  y[g!=i], status[g!=i],mincut=mincut, nsteps=nsteps, backfit=backfit, maxnumcut=maxnumcut, dirp=dirp)
        for(ii in 1:nsteps)
           {aa=index.prediction(cva[[i]]$res[[ii]],x[g==i,])
            fit=coxph(Surv(y[g==i],status[g==i])~aa)
            sc[i,ii]=sign(fit$coef)*sqrt(fit$scor)
            preval[ii,g==i]=aa
            }
        }
    sc.tot=abs(sc.tot)+sc
    preval.tot=preval.tot+preval
    }

meansc=colMeans(sc.tot)
kmax=which.max(meansc)
pvfit.score=rep(0, nsteps)
for(i in 1:nsteps)
   {fit=coxph(Surv(y,status)~preval.tot[i,])
    pvfit.score[i]=sign(fit$coef)*sqrt(fit$scor)
    }


return(list(kmax=kmax, meanscore=meansc/num.replicate, pvfit.score=pvfit.score, preval=preval.tot/num.replicate))
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cv.lm.interaction.R"
cv.lm.interaction=function(x, trt, y, K.cv=5, num.replicate=1, nsteps, mincut=0.1, backfit=F, maxnumcut=1, dirp=0){
x=as.matrix(x)
n=length(y)


sc.tot=matrix(0, nrow=K.cv, ncol=nsteps)
preval.tot=matrix(0, nrow=nsteps, ncol=n)

for(b in 1:num.replicate)
   {g=sample(rep(1:K.cv,(n+K.cv-1)/K.cv))
    g=g[1:n]
    cva=vector("list", K.cv)

    

    sc=matrix(0, nrow=K.cv, ncol=nsteps)
    preval=matrix(0, nrow=nsteps, ncol=n)
    for(i in 1:K.cv)
       {cva[[i]]=lm.interaction(x[g!=i ,], trt[g!=i],  y[g!=i], nsteps=nsteps, mincut=mincut, backfit=backfit, maxnumcut=maxnumcut, dirp=dirp)
        for(ii in 1:nsteps)
           {aa=index.prediction(cva[[i]]$res[[ii]],x[g==i,])
            fit=lm(y[g==i]~aa*(trt[g==i]))
            sc[i,ii]=summary(fit)$coef[4,3]
            preval[ii,g==i]=aa
           }
        }

    sc.tot=abs(sc.tot)+sc
    preval.tot=preval.tot+preval
    }

meansc=colMeans(sc.tot)
kmax=which.max(meansc)
pvfit.score=rep(0, nsteps)
for(i in 1:nsteps)
   {fit=lm(y~preval.tot[i,]*trt)
    pvfit.score[i]=summary(fit)$coef[4,3]
    }


return(list(kmax=kmax, meanscore=meansc/num.replicate, pvfit.score=pvfit.score, preval=preval.tot/num.replicate))
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cv.lm.main.R"
cv.lm.main=function(x, y, K.cv=5, num.replicate=1, nsteps, mincut=0.1, backfit=F, maxnumcut=1, dirp=0){
x=as.matrix(x)
n=length(y)

sc.tot=matrix(0, nrow=K.cv, ncol=nsteps)
preval.tot=matrix(0, nrow=nsteps, ncol=n)

for(b in 1:num.replicate)
   {g=sample(rep(1:K.cv,(n+K.cv-1)/K.cv))
    g=g[1:n]
    cva=vector("list",K.cv)

    sc=matrix(0, nrow=K.cv, ncol=nsteps)
    preval=matrix(0, nrow=nsteps, ncol=n)
    for(i in 1:K.cv)
       {cva[[i]]=lm.main(x[g!=i ,],  y[g!=i], mincut=mincut, nsteps=nsteps, backfit=backfit, maxnumcut=maxnumcut, dirp=dirp)
        for(ii in 1:nsteps)
           {aa=index.prediction(cva[[i]]$res[[ii]],x[g==i,])
            fit=lm(y[g==i]~aa)
            sc[i,ii]=summary(fit)$coef[2,3]
            preval[ii,g==i]=aa
            }
        }
    sc.tot=abs(sc.tot)+sc
    preval.tot=preval.tot+preval
    }

meansc=colMeans(sc.tot)
kmax=which.max(meansc)
pvfit.score=rep(0, nsteps)
for(i in 1:nsteps)
   {fit=lm(y~preval.tot[i,])
    pvfit.score[i]=summary(fit)$coef[2,3]
    }


return(list(kmax=kmax, meanscore=meansc/num.replicate, pvfit.score=pvfit.score, preval=preval.tot/num.replicate))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cv.logistic.interaction.R"
cv.logistic.interaction=function(x, trt, y, K.cv=5, num.replicate=1, nsteps, mincut=0.1, backfit=F, maxnumcut=1, dirp=0, weight=1){
x=as.matrix(x)
n=length(y)


sc.tot=matrix(0, nrow=K.cv, ncol=nsteps)
preval.tot=matrix(0, nrow=nsteps, ncol=n)

for(b in 1:num.replicate)
   {g=sample(rep(1:K.cv,(n+K.cv-1)/K.cv))
    g=g[1:n]
    cva=vector("list", K.cv)

    

    sc=matrix(0, nrow=K.cv, ncol=nsteps)
    preval=matrix(0, nrow=nsteps, ncol=n)
    for(i in 1:K.cv)
       {cva[[i]]=logistic.interaction(x[g!=i ,],  trt[g!=i], y[g!=i], nsteps=nsteps, mincut=mincut, backfit=backfit, maxnumcut=maxnumcut, dirp=dirp, weight=weight)
        for(ii in 1:nsteps)
           {aa=index.prediction(cva[[i]]$res[[ii]],x[g==i,])
            fit=glm(y[g==i]~aa*(trt[g==i]), family="binomial")
            sc[i,ii]=summary(fit)$coef[4,3]
            preval[ii,g==i]=aa
           }
        }

    sc.tot=abs(sc.tot)+sc
    preval.tot=preval.tot+preval
    }

meansc=colMeans(sc.tot)
kmax=which.max(meansc)
pvfit.score=rep(0, nsteps)
for(i in 1:nsteps)
   {fit=glm(y~preval.tot[i,]*trt, family="binomial")
    pvfit.score[i]=summary(fit)$coef[4,3]
    }


return(list(kmax=kmax, meanscore=meansc/num.replicate, pvfit.score=pvfit.score, preval=preval.tot/num.replicate))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/cv.logistic.main.R"
cv.logistic.main=function(x, y,  K.cv=5, num.replicate=1, nsteps, mincut=0.1, backfit=F, maxnumcut=1, dirp=0, weight=1){
x=as.matrix(x)
n=length(y)

sc.tot=matrix(0, nrow=K.cv, ncol=nsteps)
preval.tot=matrix(0, nrow=nsteps, ncol=n)

for(b in 1:num.replicate)
   {g=sample(rep(1:K.cv,(n+K.cv-1)/K.cv))
    g=g[1:n]
    cva=vector("list",K.cv)

    sc=matrix(0, nrow=K.cv, ncol=nsteps)
    preval=matrix(0, nrow=nsteps, ncol=n)
    for(i in 1:K.cv)
       {cva[[i]]=logistic.main(x[g!=i ,],  y[g!=i], mincut=mincut, nsteps=nsteps, backfit=backfit, maxnumcut=maxnumcut, dirp=dirp, weight=weight)
        for(ii in 1:nsteps)
           {aa=index.prediction(cva[[i]]$res[[ii]],x[g==i,])
            fit=glm(y[g==i]~aa, family="binomial")
            sc[i,ii]=summary(fit)$coef[2,3]
            preval[ii,g==i]=aa
            }
        }
    sc.tot=abs(sc.tot)+sc
    preval.tot=preval.tot+preval
    }

meansc=colMeans(sc.tot)
kmax=which.max(meansc)
pvfit.score=rep(0, nsteps)
for(i in 1:nsteps)
   {fit=glm(y~preval.tot[i,], family="binomial")
    pvfit.score[i]=summary(fit)$coef[2,3]
    }


return(list(kmax=kmax, meanscore=meansc/num.replicate, pvfit.score=pvfit.score, preval=preval.tot/num.replicate))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/index.prediction.R"

index.prediction=function(res, x){
k=nrow(res)
scor=rep(0,nrow(x))
for(i in 1:k){
 jmax=res[i,"jmax"]
 cutp=res[i,"cutp"]
 maxdir=res[i,"maxdir"]
  if(maxdir==1){ temp=1*(x[,jmax]> cutp)}
  if(maxdir==-1){ temp=1*(x[,jmax]< cutp)}
 scor=scor+temp
}
return(scor)
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/lm.interaction.R"

lm.interaction=function(x, trt, y, nsteps=8, backfit=F, mincut=0.1, maxnumcut=1, dirp=0)
{n=length(y)
 x0=x
 p.true=length(x0[1,])
 
 x=cbind(x0, -x0)
 marker.x=rep(1:p.true, 2)
 direction.x=c(rep(-1, p.true), rep(1, p.true))

 if(sum(abs(dirp))>0)
    {index1=(1:p.true)[dirp==-1 | dirp==0]
     index2=(1:p.true)[dirp==1 | dirp==0]
     x=cbind(x0[,index1], x0[,index2])
     marker.x=c((1:p.true)[index1], (1:p.true)[index2])
     direction.x=c(rep(-1, length(index1)), rep(1, length(index2)))
     }                   
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}


 score.current=rep(0,n)
 x.out=x

 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x.out[,i])}

 y.pool=trt.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     trt.pool[,i]=trt[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 
 id.in=NULL
 id.out=1:p.true
 p.out=p
 zvalue=NULL
 cut.value=NULL
 direction=NULL
 imax=NULL
 res=as.list(NULL)
 num.cut=rep(0, p.true)


 fit=lm(y~trt)
 muhat=fit$fitted
 sigma=sum((y-muhat)^2)/(n-2)
 sigma.sum=sum(solve(t(cbind(1,trt))%*%cbind(1,trt))*sigma)
 loop=1
 count=1
 

 while(loop==1)
      {score.pool=muhat.pool=matrix(0, n, p.out)
       for(i in 1:p.out)
          {score.pool[,i]=score.current[order.index[,i]]
           muhat.pool[,i]=muhat[order.index[,i]]}

       score.stat=apply(trt.pool*(y.pool-muhat.pool),2,cumsum)+sum((y-muhat)*score.current*trt)
       v11.stat=sigma*(apply((2*score.pool+1)*trt.pool,2,cumsum)+sum(score.current^2*trt))
       v22.stat=(sum(score.current*trt)+apply(trt.pool,2,cumsum))^2*sigma.sum
       v.stat=(v11.stat-v22.stat)

       for(i in 1:p.out)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p.out)
          for(i in 1:p.out)
             {test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }


       mcut=rep(0, p.out)
       i0=rep(0, p.out)
       for(i in 1:p.out)
          {i0[i]=max(effect.range[test[effect.range,i]==mtest[i]])+1
           if(i0[i]>n)
              mcut[i]=Inf
           if(i0[i]<=n)
              mcut[i]=x.out[order.index[,i],i][i0[i]]
           }
         
       i.sel=(1:p.out)[mtest==max(mtest)][1]
       score.current=score.current+1*(x.out[,i.sel]<mcut[i.sel])     

       marker.sel=marker.x[i.sel]
       id.in=c(id.in, marker.sel)
       direction=c(direction, direction.x[i.sel])
       cut.value=c(cut.value, -mcut[i.sel]*direction.x[i.sel])
       num.cut[marker.sel]=num.cut[marker.sel]+1

       if(num.cut[marker.sel]==maxnumcut)
         {id.exclude=(1:p.out)[marker.x==marker.sel]
   
          x.out=x.out[,-id.exclude,drop=F]
          order.index=order.index[, -id.exclude,drop=F]
          y.pool=y.pool[,-id.exclude,drop=F]
          trt.pool=trt.pool[,-id.exclude,drop=F]
          x.replicate=x.replicate[, -id.exclude,drop=F]

          direction.x=direction.x[-id.exclude]
          marker.x=marker.x[-id.exclude]   
          p.out=p.out-length(id.exclude)        
          }

       if(backfit==T)
         {if(count>1)
            {x.adj=x0[, id.in]
             x.adj=-t(t(x.adj)*direction)
             cutp=-cut.value*direction
             fit=backfit.lm.interaction(x.adj, trt, y, cutp, mincut=mincut)
         
             jmax=id.in
             cutp=-fit$cutp*direction
             maxdir=direction
             
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(fit$zscore)) 
             score.current=apply((t(x.adj)<fit$cutp),2,sum)
             }

          if(count==1)
            {jmax=id.in
             cutp=cut.value
             maxdir=direction
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(mtest))
             }
           }
       if(backfit==F)
         {cutp=cut.value
          maxdir=direction
          jmax=id.in
          zvalue=c(zvalue, max(mtest))
          maxsc=zvalue
          res[[count]]=cbind(jmax, cutp, maxdir, maxsc)
          }

       count=count+1
       loop=(length(id.in)<nsteps)
       }

       
 return(list(res=res, maxsc=zvalue))
 }
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/lm.main.R"
lm.main=function(x, y,  nsteps=8, backfit=F, mincut=0.1, maxnumcut=1, dirp=0)
{n=length(y)
 x0=x
 p.true=length(x0[1,])
 
 x=cbind(x0, -x0)
 marker.x=rep(1:p.true, 2)
 direction.x=c(rep(-1, p.true), rep(1, p.true))

 if(sum(abs(dirp))>0)
    {index1=(1:p.true)[dirp==-1 | dirp==0]
     index2=(1:p.true)[dirp==1 | dirp==0]
     x=cbind(x0[,index1], x0[,index2])
     marker.x=c((1:p.true)[index1], (1:p.true)[index2])
     direction.x=c(rep(-1, length(index1)), rep(1, length(index2)))
     }                   
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}

 score.current=rep(0,n)
 num.cut=rep(0, p.true)

 x.out=x

 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x.out[,i])}

 y.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 
 id.in=NULL
 zvalue=NULL
 cut.value=NULL
 direction=NULL
 res=as.list(NULL)
 p.out=p

 muhat=mean(y)
 sigma=sum((y-muhat)^2)/(n-1)
 loop=1
 count=1
 

 while(loop==1)
      {score.pool=matrix(0, n, p.out)
       for(i in 1:p.out)
          {score.pool[,i]=score.current[order.index[,i]]}

       score.stat=apply(y.pool-muhat,2,cumsum)+sum((y-muhat)*score.current)
       v.stat=sigma*(2*apply(score.pool,2,cumsum)+1:n+sum(score.current^2)-(sum(score.current)+1:n)^2/n)

       for(i in 1:p.out)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p.out)
          for(i in 1:p.out)
             {test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }


       mcut=rep(0, p.out)
       i0=rep(0, p.out)
       for(i in 1:p.out)
          {i0[i]=max(effect.range[test[effect.range,i]==mtest[i]])+1
           if(i0[i]>n)
              mcut[i]=Inf
           if(i0[i]<=n)
              mcut[i]=x.out[order.index[,i],i][i0[i]]
           }
         
       i.sel=(1:p.out)[mtest==max(mtest)][1]
       score.current=score.current+1*(x.out[,i.sel]<mcut[i.sel])     

       marker.sel=marker.x[i.sel]
       id.in=c(id.in, marker.sel)
       direction=c(direction, direction.x[i.sel])
       cut.value=c(cut.value, -mcut[i.sel]*direction.x[i.sel])
       num.cut[marker.sel]=num.cut[marker.sel]+1

       if(num.cut[marker.sel]==maxnumcut)
         {id.exclude=(1:p.out)[marker.x==marker.sel]
   
          x.out=x.out[,-id.exclude,drop=F]
          order.index=order.index[, -id.exclude,drop=F]
          y.pool=y.pool[,-id.exclude,drop=F]
          x.replicate=x.replicate[, -id.exclude,drop=F]
          direction.x=direction.x[-id.exclude]
          marker.x=marker.x[-id.exclude]   
          p.out=p.out-length(id.exclude)        
          }

       if(backfit==T)
         {if(count>1)
            {x.adj=x0[, id.in]
             x.adj=-t(t(x.adj)*direction)
             cutp=-cut.value*direction
             fit=backfit.lm.main( x.adj, y,cutp, mincut=mincut)
         
             jmax=id.in
             cutp=-fit$cutp*direction
             maxdir=direction
             
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue, max(fit$zscore)) 
             score.current=apply((t(x.adj)<fit$cutp),2,sum)
             }

          if(count==1)
            {jmax=id.in
             cutp=cut.value
             maxdir=direction
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(mtest))
             }
           }
       if(backfit==F)
         {cutp=cut.value
          maxdir=direction
          jmax=id.in
          zvalue=c(zvalue, max(mtest))
          maxsc=zvalue
          res[[count]]=cbind(jmax, cutp, maxdir, maxsc)
          }

       count=count+1
       loop=(length(id.in)<nsteps)
       }

       
 return(list(res=res, maxsc=zvalue))
 }
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/logistic.interaction.R"

logistic.interaction=function(x, trt, y, nsteps=8, mincut=0.1, backfit=F, maxnumcut=1, dirp=0, weight=1)
# interaction AIM for logistic
#  y- n vector of outcomes (0-1)
# trt= 0-1 treatment variable- should be centered before input into function
# x- n by p matrix of predcitors
# nsteps= # number terms to add 
# mincut- min  cutting prop (at either end)  
# backfit- should split points be re-adjusted after every term is added?
# maxnumcut- max # of splits per feature
# dirp- directions of splits allowed- a p-vector of 0,1, or -1.
#  1 means split to the left; -1 to the right; 0 means either. dirp=0
#   means all splits are unrestricted
#  weight= obs wt given to second class

# Output is a list res, with each component summarizing the model
# after a term has been added. columns are
#  jmax- predictor#,  cutp- cutpoint,  maxdir- direction of split (1 means le,
# -1 means >) and    maxsc- likelihood score achieved
{n=length(y)
 x0=x
 p.true=length(x0[1,])

 if(length(weight)==1)
   {wt=y*(weight-1)+1}
 if(length(weight)==n)
   {wt=weight} 

 x=cbind(x0, -x0)
 marker.x=rep(1:p.true, 2)
 direction.x=c(rep(-1, p.true), rep(1, p.true))

 if(sum(abs(dirp))>0)
    {index1=(1:p.true)[dirp==-1 | dirp==0]
     index2=(1:p.true)[dirp==1 | dirp==0]
     x=cbind(x0[,index1], -x0[,index2])
     marker.x=c((1:p.true)[index1], (1:p.true)[index2])
     direction.x=c(rep(-1, length(index1)), rep(1, length(index2)))
     }                   
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}

 score.current=rep(0,n)
 num.cut=rep(0, p.true)

 order.index=matrix(0, n, p)
 for(i in 1:p.true)
    {order.index[,i]=order(x[,i])
     order.index[,i+p.true]=rev(order.index[,i])}

 wt.pool=trt.pool=y.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     trt.pool[,i]=trt[order.index[,i]]
     wt.pool[,i]=wt[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 x.out=x 
 id.in=NULL
 p.out=p
 zvalue=NULL
 cut.value=NULL
 direction=NULL
 res=as.list(NULL)

 fit=glm(y~trt, family="binomial", weight=wt)
 phat=fit$fitt
 sigma=sum(summary(fit)$cov.scale)


 count=1
 loop=1
 while(loop==1)
      {score.pool=phat.pool=matrix(0, n, p.out)
       for(i in 1:p.out)
          {score.pool[,i]=score.current[order.index[,i]]
           phat.pool[,i]=phat[order.index[,i]]}

       score.stat=apply(trt.pool*wt.pool*(y.pool-phat.pool),2,cumsum)+sum((y-phat)*trt*score.current*wt)
       v11.stat=2*apply(score.pool*trt.pool*wt.pool^2*phat.pool*(1-phat.pool),2,cumsum)+apply(trt.pool*wt.pool^2*phat.pool*(1-phat.pool),2,cumsum)+sum((score.current*wt)^2*trt*phat*(1-phat))
       v22.stat=(sum(score.current*trt*wt*phat*(1-phat))+apply(trt.pool*wt.pool*phat.pool*(1-phat.pool),2,cumsum))^2*sigma
       v.stat=(v11.stat-v22.stat)

       
       for(i in 1:p.out)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p.out)
          for(i in 1:p.out)
             {
              test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }


       mcut=rep(0, p.out)
       i0=rep(0, p.out)
       for(i in 1:p.out)
          {i0[i]=max(effect.range[test[effect.range,i]==mtest[i]])+1
           if(i0[i]>n)
              mcut[i]=Inf
           if(i0[i]<=n)
              mcut[i]=x.out[order.index[,i],i][i0[i]]
           }
         
       i.sel=(1:p.out)[mtest==max(mtest)][1]
       score.current=score.current+(x.out[,i.sel]<mcut[i.sel])


       marker.sel=marker.x[i.sel]
       id.in=c(id.in, marker.sel)
       direction=c(direction, direction.x[i.sel])
       cut.value=c(cut.value, -mcut[i.sel]*direction.x[i.sel])
       num.cut[marker.sel]=num.cut[marker.sel]+1

       if(num.cut[marker.sel]==maxnumcut)
         {id.exclude=(1:p.out)[marker.x==marker.sel]

          x.out=x.out[,-id.exclude,drop=F]
          order.index=order.index[, -id.exclude,drop=F]
          y.pool=y.pool[,-id.exclude,drop=F]
          wt.pool=wt.pool[,-id.exclude, drop=F]
          trt.pool=trt.pool[,-id.exclude,drop=F]
          x.replicate=x.replicate[, -id.exclude,drop=F]          

          direction.x=direction.x[-id.exclude]
          marker.x=marker.x[-id.exclude]
          p.out=p.out-length(id.exclude)
         }

       
       if(backfit==T)
         {if(count>1)
            {x.adj=x0[, id.in]
             x.adj=-t(t(x.adj)*direction)
             cutp=-cut.value*direction
             fit=backfit.logistic.interaction(x.adj, trt, y, cutp, mincut=mincut, wt=wt)
         
             jmax=id.in
             cutp=-fit$cutp*direction
             maxdir=direction
             
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(fit$zscore)) 
             score.current=apply((t(x.adj)<fit$cutp),2,sum)
             }

          if(count==1)
            {jmax=id.in
             cutp=cut.value
             maxdir=direction
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(mtest))
             }
           }
       if(backfit==F)
         {cutp=cut.value
          maxdir=direction
          jmax=id.in
          zvalue=c(zvalue, max(mtest))
          maxsc=zvalue
          res[[count]]=cbind(jmax, cutp, maxdir, maxsc)
          }
       
       count=count+1
       loop=(length(id.in)<nsteps)
       }

       
 return(list(res=res, maxsc=zvalue))
 }

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpuIiijG/R.INSTALL10d63af1/AIM/R/logistic.main.R"
logistic.main=function(x, y, nsteps=8, mincut=.1, backfit=F, maxnumcut=1, dirp=0, weight=1)
# main effects AIM for logistic
# y- n vector of outcomes (0-1)
# x- n by p matrix of predictors
# nsteps= # number terms to add 
# mincut- min  cutting prop (at either end)  
# backfit- should split points be re-adjusted after every term is added?
# maxnumcut- max # of splits per feature
# dirp- directions of splits allowed- a p-vector of 0,1, or -1.
#  1 means split to the left; -1 to the right; 0 means either. dirp=0
#   means all splits are unrestricted
#  weight= obs wt given to second class

# Output is a list res, with each component summarizing the model
# after a term has been added. columns are
#  jmax- predictor#,  cutp- cutpoint,  maxdir- direction of split (1 means le,
#-1 means >) and    maxsc- likelihood score achieved

{n=length(y)
 x0=x
 p.true=length(x0[1,])

 if(length(weight)==1)
   {wt=y*(weight-1)+1}
 if(length(weight)==n)
   {wt=weight} 

 x=cbind(x0, -x0)
 marker.x=rep(1:p.true, 2)
 direction.x=c(rep(-1, p.true), rep(1, p.true))

 if(sum(abs(dirp))>0)
    {index1=(1:p.true)[dirp==-1 | dirp==0]
     index2=(1:p.true)[dirp==1 | dirp==0]
     x=cbind(x0[,index1], -x0[,index2])
     marker.x=c((1:p.true)[index1], (1:p.true)[index2])
     direction.x=c(rep(-1, length(index1)), rep(1, length(index2)))
     }                   
 p=length(x[1,])
 

 if(mincut>0)
   {effect.range=ceiling(n*mincut):floor(n*(1-mincut))}
 if(mincut==0)
   {effect.range=1:n}

 score.current=rep(0,n)
 num.cut=rep(0, p.true)

 x.out=x
 id.in=NULL
 p.out=p
 zvalue=NULL
 cut.value=NULL
 direction=NULL
 res=as.list(NULL)
 
 order.index=matrix(0, n, p)
 for(i in 1:p)
    {order.index[,i]=order(x.out[,i])}

 wt.pool=y.pool=x.replicate=matrix(0, n, p)
 for(i in 1:p)
    {y.pool[,i]=y[order.index[,i]]
     wt.pool[,i]=wt[order.index[,i]]
     x.replicate[,i]=rep((1:n)[cumsum(table(x[,i]))], table(x[,i]))}

 phat=sum(y*wt)/sum(wt)
 sigma=as.numeric(summary(glm(y~1, family="binomial", weight=wt))$cov.un*(phat*(1-phat))^2*n)

 loop=1
 count=1
 
 while(loop==1)
      {score.pool=matrix(0, n, p.out)
       for(i in 1:p.out)
          {score.pool[,i]=score.current[order.index[,i]]}

       score.stat=apply((y.pool-phat)*wt.pool,2,cumsum)+sum((y-phat)*(score.current*wt))
       v.stat=sigma*(2*apply(score.pool*wt.pool^2,2,cumsum)+apply(wt.pool^2,2,cumsum)+sum((score.current*wt)^2)-(sum(score.current*wt)+apply(wt.pool,2,cumsum))^2/n)


       for(i in 1:p.out)
          {idx=x.replicate[,i]
           v.stat[,i]=v.stat[idx ,i]   
           score.stat[,i]=score.stat[idx, i]
          } 

       test=score.stat/sqrt(v.stat+1e-8)

       if(mincut==0)
         {mtest=apply(test,2,max)}

       if(mincut>0)
         {mtest=rep(0, p.out)
          for(i in 1:p.out)
             {test.post=test[max(effect.range)+1,i]
              test0=test[effect.range,i]
              test0=test0[test0!=test.post]
              mtest[i]=max(test0)}
          }


       mcut=rep(0, p.out)
       i0=rep(0, p.out)
       for(i in 1:p.out)
          {i0[i]=max(effect.range[test[effect.range,i]==mtest[i]])+1
           if(i0[i]>n)
              mcut[i]=Inf
           if(i0[i]<=n)
              mcut[i]=x.out[order.index[,i],i][i0[i]]
           }
         
       i.sel=(1:p.out)[mtest==max(mtest)][1]
       score.current=score.current+1*(x.out[,i.sel]<mcut[i.sel])     

       
       marker.sel=marker.x[i.sel]
       id.in=c(id.in, marker.sel)
       direction=c(direction, direction.x[i.sel])
       cut.value=c(cut.value, -mcut[i.sel]*direction.x[i.sel])
       num.cut[marker.sel]=num.cut[marker.sel]+1

       if(num.cut[marker.sel]==maxnumcut)
         {id.exclude=(1:p.out)[marker.x==marker.sel]
   
          x.out=x.out[,-id.exclude,drop=F]
          order.index=order.index[, -id.exclude,drop=F]
          y.pool=y.pool[,-id.exclude,drop=F]
          wt.pool=wt.pool[,-id.exclude, drop=F]
          x.replicate=x.replicate[, -id.exclude,drop=F]
          direction.x=direction.x[-id.exclude]
          marker.x=marker.x[-id.exclude]   
          p.out=p.out-length(id.exclude)        
          }
          


       if(backfit==T)
         {if(count>1)
            {x.adj=x0[, id.in]
             x.adj=-t(t(x.adj)*direction)
             cutp=-cut.value*direction
             fit=backfit.logistic.main(x.adj, y, cutp, mincut=mincut, wt=wt)
         
             jmax=id.in
             cutp=-fit$cutp*direction
             maxdir=direction
             
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(fit$zscore)) 
             score.current=apply((t(x.adj)<fit$cutp),2,sum)
             }

          if(count==1)
            {jmax=id.in
             cutp=cut.value
             maxdir=direction
             res[[count]]=cbind(jmax, cutp, maxdir)
             zvalue=c(zvalue,  max(mtest))
             }
           }
       if(backfit==F)
         {cutp=cut.value
          maxdir=direction
          jmax=id.in
          zvalue=c(zvalue, max(mtest))
          maxsc=zvalue
          res[[count]]=cbind(jmax, cutp, maxdir, maxsc)
          }

       
       count=count+1
       loop=(length(id.in)<nsteps)
       }

       
 return(list(res=res, maxsc=zvalue))
 }

