function [bestCost,bestTour]=ACO_TSP(a,genMax,alpha,betap,rho,nameFile)
// a = number of ants
// genMax = number of generations
// alfa, beta = parameters used in the computation of the transition probability
// rho = evaporation rate for pheromones
// example: ACO_TSP(50,1000,0.5,0.5,0.75,"eil51.tsp")
    if (argn(2)==6)
        rez=csvRead(nameFile,' ')';
        x=rez(2:3,:);
        n=size(x,2);
    else    
        [x,n]=generare();                   // select coordinates of nodes
    end   
    // best tour initialization 
    for i=1:n
        bestTour(i)=i;
    end
    bestCost=cost(x,bestTour);
    // initialization of the pheromone matrix
    tau=grand(n,n,'unf',0,1); // symmetric matrix
    gen=0;
    costMatrix=computeCostMatrix(n,x);
    disp(costMatrix);
    while gen<genMax
        gen=gen+1;
        for k=1:a // each agent constructs a tour 
            tour=zeros(1,n);
            tour(1)=1;  // all tours start from node 1
            for i=2:n
                tour(i)=selectNext(i,n,tour,tau,alpha,betap,costMatrix);
            end        
            tourCost=cost(x,tour);
            if (tourCost<bestCost)
                bestTour=tour;                
                bestCost=tourCost;
            end
        end
        tau=update(tau,rho,n,bestTour,costMatrix);
//        disp(tau,"tau=");
        disp(bestCost,"cost=");
//        disp(bestTour,"tur=");
    end 
    //plotting the solution
    clf
    plot(x(1,:),x(2,:),'ko');
    for i=1:n-1
         plot([x(1,bestTour(i)),x(1,bestTour(i+1))],[x(2,bestTour(i)),x(2,bestTour(i+1))],'r*-');
    end
    plot([x(1,bestTour(n)),x(1,bestTour(1))],[x(2,bestTour(n)),x(2,bestTour(1))],'r*-');
endfunction

// selection of nodes coordinates
function [x,n]=generare()
   // data input
   clf;
   plot2d(0,0,0,rect=[0,0,1,1]);
   xgrid(3);
   p=locate(-1,1);   
   n=size(p,2); 
   x = p;
endfunction  

// selection of the next node

function node=selectNext(i,n,tour,tau,alpha,betab,costMatrix)
  poz=zeros(1,n);
  for j=1:i-1
      poz(tour(j))=1;
  end
  suma=0
  for j=1:n   
      if(tour(i-1)~=j) 
        suma=suma+(1-poz(j))*(tau(tour(i-1),j)^alpha)/(costMatrix(tour(i-1),j)^betab);
      end  
  end
  prob=zeros(1,n);
  for j=1:n
      if (poz(j)==0) then
          prob(j)=(tau(tour(i-1),j)^alpha/(costMatrix(tour(i-1),j)^betab))/suma;
      end
  end
  node=selectNode(prob);
endfunction        
    
// computation of the tour cost
function ct=cost(x,t)
    n=length(t);
    ct=0;
    for i=1:n-1
        ct=ct+sqrt((x(1,t(i))-x(1,t(i+1)))^2+(x(2,t(i))-x(2,t(i+1)))^2);
    end    
    ct=ct+sqrt((x(1,t(n))-x(1,t(1)))^2+(x(2,t(n))-x(2,t(1)))^2);   
endfunction    

function costMatrix=computeCostMatrix(n,x);
    costMatrix=zeros(n,n);
    for i=1:n
        for j=1:n
            costMatrix(i,j)=sqrt((x(1,i)-x(1,j))^2+(x(2,i)-x(2,j))^2);
        end
    end
endfunction

//  selection node (proportional)
function i=selectNode(prob)
    i=1;
    s=prob(1);
    u=grand(1,1,"unf",0,1);
    while u>s
        i=i+1;
        s=s+prob(i);
    end    
endfunction

// update of the pheromone matrix
function tau=update(tau,rho,n,bestTour,costMatrix)
    for i=1:n
        for j=1:n
            tau(i,j)=(1-rho)*tau(i,j)+updateTerm(i,j,n,bestTour,costMatrix);
        end
    end
endfunction

function termen=updateTerm(i,j,n,bestTour,costMatrix)
    termen=0;
    for k=1:n-1
        if ((i==bestTour(k))&(j==bestTour(k+1))) 
            termen=1.0/costMatrix(i,j);
        end
    end
endfunction
// ACO_TSP(a,genMax,alpha,beta,rho,nameFile)
//[bestCost,bestTour]=ACO_TSP(50,100,0.5,0.5,0.25,"eil51.tsp")
[bestCost,bestTour]=ACO_TSP(20,100,0.5,0.5,0.25,"eil51.tsp")
