#########################################################
# --- Regularized ROC with the TGDR for variable selection
# --- The function for estimation after cross validation
# --- with fixed tau; 
# --- the first covariate is the anchor
#     Written by Shuangge Ma (joint with Jian Huang), Dec. 2005
###################################################################

ind<-read.table("indicator.txt")$V1;
z.dat<-read.table("covariate.txt");

sigma.n<-0.50; # --- sigma.n selected using another program
anchor<-1;
step<- 625;    # --- number of iterations estimated from cross validation
increase<-0.01; tau<-1.0;

samplesize<-nrow(z.dat);
num.cov<-ncol(z.dat);

z<-matrix(0, samplesize, num.cov);
for(i in 1:num.cov)
   z[ ,i]<-z.dat[[i]];
z<-z/sigma.n;
aug.cross<-600;
case.z<-z[ind==1, ]; control.z<-z[ind==0, ];

aug.cross.data<-matrix(0, aug.cross, num.cov);
for(i in 1:aug.cross){
   p1<-as.integer(runif(1, min=1, max=(nrow(case.z)+1)));
   p2<-as.integer(runif(1, min=1, max=(nrow(control.z)+1)));
      
   for(j in 1:num.cov){
      aug.cross.data[i, j]<- case.z[p1, j]-control.z[p2, j];
   }
}

beta.est<-rep(0, num.cov);  beta.est[anchor]<- 1.0;
h.v<-rep(0, num.cov); f.v<-rep(0, num.cov); g.v<-rep(0, num.cov);

for(s in 1:step){
      gradient<- rep(0, num.cov);
      for(i in 1:nrow(aug.cross.data)){
         diff.z<-aug.cross.data[i, ];
         s.func<-1/(1+exp(-sum(diff.z*beta.est)));
         gradient<-gradient+s.func*(1-s.func)*diff.z;
      }
         
      gradient[anchor]<- 0;
      
      # -- standardized gradient ---
      if(max(abs(gradient)) >0.000001){
         for(i in 1:length(g.v))
            g.v[i]<-gradient[i]/max(abs(gradient));
      }
      if(max(abs(gradient)) <= 0.000001){
         for(i in 1:length(g.v))
            g.v[i]<-0.0; 
      }
      
      f.v<-ifelse(abs(g.v)>=tau, 1, 0);
 
      for(i in 1:length(g.v))
         h.v[i]<-f.v[i]*g.v[i];
      
      beta.est<-beta.est+increase*h.v;
}

write.table(t(beta.est), file="beta-est.txt", append=F,
            quote=F, row.name=F, col.name=F);