function [S_best, best, err_best, err_S_best, err_S2N] = ACO_S2N(y,S2N,m,na,p,objf)

%% Ant Colony Optimization

%% Data
% features  : y.features
% data      : y.sample
% labels    : y.label

%% Step 1: Define the terms

row         = size(y.sample_tr,1);
const       = 1;
Trail       = const * (ones(row,1));    % contstant
DeltaTrail  = 0 * (ones(row,1));        % initially zero.

Niter       = 500;  % The ants tend to converge around 700 - 800 in some 
                     % iterations and sometimes may need up to 1000.
k           = round(na/3);     
kbest       = 100;
eta         = 1;
kappa       = 1;
rho         = 0.1;     % (1- rho) represents the evaporation of pheromone trial
Snew = zeros(na,m);
   
%% Step 2: If in the first iteration

  S = []; 
  for j = 1:na
      R = randperm(row);
      S = [S;R(1:m)];
  end
St= {};
 for iter=1:Niter
%      iter
       % Step 3    
       for mm     = m-p+1:m
           MaxUSM = [];
           kk     = 1;
           for j=1:na
               U   = [];
               Tg  = [];
               LI  = [];
               USM = [];
                          
               for i=1:row      % Given subset Sj, Choose feature fi that maximizes USM_i^sj
                   U(i) = 1;    % U(i) =1 if i is not an element of Sj, 0 otherwise
                   for jj=1:mm-1
                    if i==Snew(j,jj)
                       U(i)  = 0;
                       Tg(i) = 0;
                       LI(i) = 0;
                    end
                   end

                   if U(i) == 1
                    LI(i) = S2N(i);                                  
                    Tg(i) = Trail(i);
                   end
               end
               
               for i=1:row
                    if U(i) == 0
                       USM(i) = 0;
                    else
                       USM(i) = ((Tg(i)^eta)*(LI(i)^kappa))/sum(((Tg.^eta).*(LI.^kappa)));
                    end
               end
               
               [USMmax,Indexmax] = max(USM);
               
               % PROBABILITY!!!
               
               % (e.g. q=0.8 implies 80% chance to Indexmax and 20% to
               % chance to others according to their USM(i))

               qq = rand;
               if qq <= 0.4
                   sadd = Indexmax;
               else
                  features = USM;   
                  sumfea   = sum(features);
                  rsumfea  = rand*sumfea;
                  aa=0; sadd=1;
                  for i=1:row
                    aa=aa+features(i);
                    if aa <= rsumfea
                       sadd = i+1;
                    end
                  end
               end
               
               Snew(j,mm) = sadd;
           end
           %Snew = [Snew MaxUSM];    
       end
                 
       if iter ~= 1
        S = Snew;
       end
 
%% Step 4: Evaluate the selected subset of each ant using LOO-SVM

        vindex = S';
        err_tr = feval(objf,vindex,y);
        Ssets  = [S err_tr'];           % selected subsets with error as the last column
        S_sort = sortbycol(Ssets,m+1);  % Sort the set based on error.
        S      = S_sort(:,1:m);         % updated S
        
        if S_sort(1,m+1) < kbest
         best  = S_sort(1,:);          % Best features so far
         kbest = S_sort(1,m+1);
        end                          
   
        %% Step 5: Update Pheromone Trail

        %S_best = S_sort(1:k,:); % the best k ants
        S_best      = S_sort(:,:);
        DeltaTrail  = zeros(row,1);
        
        for j = 1:na
            
            if S_best(j,m+1)== 0
                iter = Niter;
            else
                DT            = 1/S_best(j,m+1);
                i             = S_best(j,1:m);
                DeltaTrail(i) = DeltaTrail(i) + DT;
            end
  
        end
  
        Trail = rho*Trail + DeltaTrail; 
        Snew  = zeros(na,m);
       for j = 1:na
             S_best_vec = [];
            for jj=1:k
             S_best_vec = [S_best_vec S_best(jj,1:m)];
            end
            Su_best_vec   = unique(S_best_vec);
            leng          = length(Su_best_vec);
            R             = randperm(leng);
            Snew(j,1:m-p) = Su_best_vec(R(1:m-p));
       end
        
% St =[St; S_best];

% save St3_2 St rho
 end

 vindex     = best(1:m)';
 err_best   = feval(objf,vindex,y);
 vindex     = S_best(1,1:m)';
 err_S_best = feval(objf,vindex,y);
 S2Nsorted  = sortbycold([(1:row)' S2N],2);
 vindex     = S2Nsorted(1:m,1);
 err_S2N    = feval(objf,vindex,y);
    
return