# Forming the Gaussian kernel matrix
KGaussian<-function(gamma,A,tildeA){
   nrA<-nrow(A); nrAt<-nrow(tildeA)
   A1<-rowSums(A^2)%*%matrix(1,nrow=1,ncol=nrAt)
   A2<-rowSums(tildeA^2)%*%matrix(1,nrow=1,ncol=nrA)
   K<-A1+t(A2)-2*A%*%(t(tildeA))
   K<-exp(-K*gamma)
   K}
 
# The RLS-SVR
zerossvr<-function(A,y,w0,b0,C){
   rx<-nrow(A); cx<-ncol(A)
   e<-matrix(1,nrow=rx,ncol=1)
   flag<-1
   T<-diag(1,cx+1)
   q<-t(A)%*%e
   t1<-cbind(t(A)%*%A,q)
   t2<-cbind(t(q),t(e)%*%e)
   temp<-rbind(t1,t2)
   hessian<-T+C*temp
   while (flag>10^(-4)){
      h<-C*(y-A%*%w0-e%*%b0)
      gradient<-rbind(w0-t(A)%*%h,b0-t(e)%*%h)
      norm<-max(svd(gradient)$d)
      if (norm/sqrt(rx)>10^(-5)){
         z<-solve(hessian)%*%(-gradient)
         w0<-w0+z[1:cx,]; b0<-b0+z[cx+1,]
         flag<-max(svd(z)$d)}
      else {flag<-0}}
   w<-w0;b<-b0
   list(w=w,b=b)
  }

# gene selection
var_select<-function(gamma,C,A,y,select_time=2,final_size=10,intermedia_size=final_size*10){
   rx<-nrow(A); cx<-ncol(A)
   std<-sqrt(rowSums((A-rowMeans(A))^2)/cx)
   stdx<-(A-rowMeans(A))/std
   IND<-c(1:cx)
   kr<-KGaussian(gamma,stdx,stdx)
   rx<-nrow(kr); cx<-ncol(kr)
   w0<-matrix(0,nrow=cx,1); b0<-0
   ww<-zerossvr(kr,y,w0,b0,C)$w
   wsums<-abs(t(ww)%*%stdx)
   IX1<-order(wsums,decreasing=TRUE)
   B1<-wsums[IX1]; IND=IND[IX1]
   if (select_time==1){
      A1<-stdx[,IX[1:final_size]]
      kr1<-KGaussian(gamma,A1,A1)
      ww<-zerossvr(kr1,y,w0,b0,C)$w
      wsums<-abs(t(ww)%*%A1)
      IX<-order(wsums,decreasing=TRUE)
      B<-wsums[IX]; IND<-IND[IX]
   }
   if (select_time==2){
      A1<-stdx[,IX1[1:intermedia_size]]
      kr1<-KGaussian(gamma,A1,A1)
      ww<-zerossvr(kr1,y,w0,b0,C)$w
      wsums<-abs(t(ww)%*%A1)
      IX2<-order(wsums,decreasing=TRUE)
      B2<-wsums[IX2]; IND<-IND[IX2]
      A2<-A1[,IX2[1:final_size]]
      kr2<-KGaussian(gamma,A2,A2)
      ww<-zerossvr(kr2,y,w0,b0,C)$w
      wsums<-abs(t(ww)%*%A2)
      IX<-order(wsums,decreasing=TRUE)
      B=wsums[IX]; IND=IND[IX]
      }
   list(B=B,IND=IND)
}