block1 <- function( ... ){
  argv = list( ... )
  i = 0
  for(a in argv){
    m = as.matrix(a)
    if(i == 0)
      rmat = m
    else
    { 
      nr = dim(m)[1]
      nc = dim(m)[2]
      aa = cbind(matrix(0,nr,dim(rmat)[2]),m)
      rmat = cbind(rmat,matrix(0,dim(rmat)[1],nc))
      rmat = rbind(rmat,aa)
    }
    i = i+1
  }
  rmat
}

Min.Set1 <- function(GMask1, Scores){
  Grand.Mask <- apply(GMask1, 2, function(x){any(x==1)})
  Grand.Min  <- rep(FALSE, length(Grand.Mask))
  Min.Genes  <- c()
  k = 0
  while (any(Grand.Min != Grand.Mask)){
    k = k + 1
    Sum.ones     <- GMask1 %*% rep(1, ncol(GMask1))
    Temp.Set     <- which(Sum.ones == max(Sum.ones))
    Min.Genes[k] <- Temp.Set[which.min(Scores[Temp.Set])]
    Grand.Min    <- Grand.Min | (GMask1[Min.Genes[k],]==1)
    GMask1[, Grand.Min] = 0
  }
  return(list(Min.Subset = Min.Genes, Covered.Obs = which(Grand.Min)))
}

POS1 <- function(ES, Core1, Y){
  Y <- as.factor(Y)
  levels(Y) = c(1,2)
  # total Core1 length
  t.Length = pmax(Core1[,2], Core1[,4]) - pmin(Core1[,1], Core1[,3])
  # overlap length
  b <- pmin(Core1[,2], Core1[,4])
  a <- pmax(Core1[,1], Core1[,3])
  o.Length = b - a
  o.Length[o.Length < 0] <- 0
  # Number of inlier (non-outlier) observations
  L <- ES
  L[, Y == 1] = (ES[, Y == 1] >= Core1[,1] & ES[, Y == 1] <= Core1[,2])
  L[, Y == 2] = (ES[, Y == 2] >= Core1[,3] & ES[, Y == 2] <= Core1[,4])
  n.inlier = apply(L, 1, sum)
  # Number of overlapping observations
  L[o.Length > 0, ] = L[o.Length > 0, ] & (ES[o.Length > 0,] <= b[o.Length > 0]) & (ES[o.Length > 0,] >= a[o.Length > 0])
  L[o.Length <= 0, ] = FALSE
  n.ovrlap1 <- apply(L[, Y == 1], 1, sum)
  n.ovrlap2 <- apply(L[, Y == 2], 1, sum)
  n.overlap = n.ovrlap1 + n.ovrlap2
  # POS
  POS1 = 4*(o.Length/t.Length)*(n.ovrlap1*n.ovrlap2/n.inlier)
  POS1[POS1 != 0] = POS1[POS1 != 0]/n.overlap[POS1 != 0]
  POS1[is.na(POS1)] = 1   #in Some Proteomic data, all observations may have a specific value for some genes ==> So, t.Length=0 ==> POS=NA
  return(POS1)
}

RDC1 <- function(GMask1, Y){
  Y <- as.factor(Y)
  levels(Y) = c(1,2)
  Relative.DC <- max.col(cbind(rowSums(GMask1[,Y == 1])/sum(Y == 1), rowSums(GMask1[,Y == 2])/sum(Y == 2)))
  return(Relative.DC)
}


Sel.Features1 <- function(ES, Y, K = "Min", Verbose = FALSE){
  Y <- as.factor(Y)
  levels(Y) = c(1,2)
  Core1   <- CI.emprical1(ES, Y)
  G.Mask <- GMask1(ES, Core1, Y)
  Pos1    <- POS1(ES, Core1, Y)
  Min.S  <- Min.Set1(G.Mask, Pos1)
  if (K == "Min"){
    MinSub <- cbind(Feature = Min.S$Min.Subset, Pos1 = Pos1[Min.S$Min.Subset])
    return(list(Features = MinSub, Covered.Obs = Min.S$Covered.Obs))
  }
  else{
    R.DC         <- RDC1(G.Mask, Y)
    Net.Features <- cbind(Feature = c(1:length(Pos1)), Pos1 = Pos1, RDC1 = R.DC)[-Min.S$Min.Subset,]
    # dividing feauters into two ordered groups (one for each class)
    Cl1.Features <- Net.Features[Net.Features[,"RDC1"] == 1,][order(Net.Features[Net.Features[,"RDC1"] == 1,][,"Pos1"]), ]
    Cl2.Features <- Net.Features[Net.Features[,"RDC1"] == 2,][order(Net.Features[Net.Features[,"RDC1"] == 2,][,"Pos1"]), ]
    # craete an empty vector to save list of ranked features within
    G.Rank <- vector(mode = "integer", length = nrow(Cl1.Features) + nrow(Cl2.Features))
    # get the start position for class whose features' number is less - 1st index in "G.Rank" is for less Pos measure
    Start <- order(c(Cl1.Features[1,"Pos1"], Cl2.Features[1,"Pos1"]))[which.min(c(nrow(Cl1.Features), nrow(Cl2.Features)))]
    # generate an indices sequence for class whose features' number is less (those indices are used to apply the alternating fashion)
    Squns <- seq(from = Start, by = 2, length.out = min(nrow(Cl1.Features), nrow(Cl2.Features)))
    
    if(which.min(c(nrow(Cl1.Features), nrow(Cl2.Features))) == 1){
      G.Rank[Squns]  <- Cl1.Features[,"Feature"]
      G.Rank[-Squns] <- Cl2.Features[,"Feature"]
    }
    else{
      G.Rank[Squns]  <- Cl2.Features[,"Feature"]
      G.Rank[-Squns] <- Cl1.Features[,"Feature"]
    }
    
    Final.Rank <- c(Min.S$Min.Subset, G.Rank)[1:K]
    nMin.Genes <- length(Min.S$Min.Subset)
    if (Verbose){
      Status  <- c(rep("Min.Set1", nMin.Genes), R.DC[G.Rank])[1:K]
      Details <- data.frame(Features = Final.Rank, Pos1 = Pos1[Final.Rank], Status = Status)
      return(list(Features = Final.Rank, nMin.Features = nMin.Genes, Measures = Details))
    }
    else{return(list(Features = Final.Rank, nMin.Features = nMin.Genes))}
  }
}

CI.emprical1 <- function(ES, Y){
  Y <- as.factor(Y)
  levels(Y) = c(1,2)
  C1Q1 <- rowQ(ES[,Y == 1], round(sum(Y == 1)/4))
  C1Q3 <- rowQ(ES[,Y == 1], round(sum(Y == 1)*3/4))
  C2Q1 <- rowQ(ES[,Y == 2], round(sum(Y == 2)/4))
  C2Q3 <- rowQ(ES[,Y == 2], round(sum(Y == 2)*3/4))
  MADC1 <- mad(ES[,Y == 1])
  MADC2 <- mad(ES[,Y == 2])
  Q    = do.call(cbind, list(C1Q1, C1Q3, C2Q1, C2Q3))
  ToAdd <- matrix(rep(c(-0.9*MADC1,0.9*MADC1,-0.9*MADC2,0.9*MADC2),nrow(Q)),nrow = nrow(Q),byrow = T)
  Core1  = data.frame(Q +ToAdd)
  Core1
}

GMask1 <- function(ES, Core1, Y){
  Y <- as.factor(Y)
  levels(Y) = c(1,2)
  ES[, Y == 1] = (ES[, Y == 1] >= Core1[,1] & ES[, Y == 1] <= Core1[,2]) & !(ES[, Y == 1] >= Core1[,3] & ES[, Y == 1] <= Core1[,4])
  ES[, Y == 2] = (ES[, Y == 2] >= Core1[,3] & ES[, Y == 2] <= Core1[,4]) & !(ES[, Y == 2] >= Core1[,1] & ES[, Y == 2] <= Core1[,2])
  mask = ES
  mode(mask) = "integer"
  return(mask)
}