function [c]=GA_TSP(m,genMax,pm,nameFile)
// m = numar de elemente in populatie (obs:  trebuie sa fie numar par)
// genMax = nr generatii
// pm = probabilitate de mutatie    
// exemplu apel: GA_TSP(50,1000,0.5,"eil51.tsp")
    if argn(2)==4
        rez=csvRead(nameFile,' ')';
        x=rez(2:3,:);
        n=size(x,2);
    else    
        [x,n]=generare(); // coordonate locatii 
    end    
    p=initPopulation(m,n);                  // initializare populatie
    c=evalPopulation(x,p);                  // evaluare populatie
    [tbest,cbest]=bestElement(p,c);         // determinarea celui mai bun element
    gen=0;
    while gen<genMax
        gen=gen+1;
        pnou=[];
        cnou=[];
        disp(gen,"gen="); 
        prob=computeProb(c);                // calcul probabilitati selectie
        for i=1:m/2
            i1=selectParent(prob);            // selectie proportionala
            i2=selectParent(prob);
//        i1=selectParentTournament(c,m);   // selectie de tip turneu
//        i2=selectParentTournament(c,m);
            [f1,f2]=crossover(p,i1,i2);     // incrucisare
//  eliminare incrucisare        
            f1=p(i1,:);
            f2=p(i2,:);
            if grand(1,1,"unf",0,1)<pm      // aplicare mutatie
                f1=mutation(f1);
            end    
            if grand(1,1,"unf",0,1)<pm
                f2=mutation(f2);
            end    
            c1=cost(x,f1);                 // calcul cost urmasi 
            c2=cost(x,f2);
            pnou=[pnou;f1;f2];             
            cnou=[cnou;c1;c2];
        end   
        [p,c]=selectSurvivors(p,c,pnou,cnou);
        [tbest,cbest]=bestElement(p,c);
        disp(cbest,"cost=");
    end 
    //vizualizare traseu aproximat de catre retea
    clf
    plot(x(1,:),x(2,:),'ko');
    for i=1:n-1
         plot([x(1,tbest(i)),x(1,tbest(i+1))],[x(2,tbest(i)),x(2,tbest(i+1))],'r*-');
    end
    plot([x(1,tbest(n)),x(1,tbest(1))],[x(2,tbest(n)),x(2,tbest(1))],'r*-');
endfunction

// generarea coordonatelor punctelor
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  

// initializarea populatiei
function p=initPopulation(m,n)
    p=[];
    for i=1:m
        p=[p; initTour(n)];
    end    
endfunction
    
// initializarea unui traseu
function t=initTour(n)
    for i=1:n
        t(1,i)=i;
    end   
    ninv=grand(1,1,"uin",1,n);
    for i=1:ninv
        i1=grand(1,1,"uin",1,n);
        i2=grand(1,1,"uin",1,n);
        aux=t(1,i1); t(1,i1)=t(1,i2); t(1,i2)=aux;
    end     
endfunction 

// evaluarea populatiei
function c=evalPopulation(x,p)
    c=[];
    for i=1:size(p,1)
        c=[c; cost(x,p(i,:))];
    end
endfunction    

// calcul cost pentru un traseu
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    

// determinarea celui mai bun  element
function [tbest,cbest]=bestElement(p,c)
    tbest=p(1,:);
    cbest=c(1,:);
    for i=2:size(p,1)
        if c(i)<cbest
            cbest=c(i,:);
            tbest=p(i,:);
        end
    end    
endfunction

// calcul probabilitati de selectie
function prob=computeProb(c)
    prob=zeros(size(c,1));
    s=0;
    for i=1:size(c,1)
        s=s+1.0/c(i);
    end
    for i=1:size(c,1)
        prob(i)=(1.0/c(i,:))/s;
    end    
endfunction

// selectie proportionala (parinti)
function i=selectParent(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

// selectie de tip turneu (parinti)
function i=selectParentTournament(c,m)
    i1=grand(1,1,"uin",1,m);
    i2=random(1,1,"uin",1,m);
    if c(i1,:)<c(i2,:)
        i=i1;
    else
        i=i2;
    end    
endfunction

function r=search(t,v)
    i=1;
    if length(t)==0
        r=0;
    else    
        while ((v~=t(i))&(i<length(t)))
            i=i+1;
        end
        if v==t(i)
            r=1;
        else 
            r=0;
        end    
    end    
endfunction

// incrucisare
function [f1,f2]=crossover(p,i1,i2)
    p1=p(i1,:);
    p2=p(i2,:);
    n=length(p1);
    k=grand(1,1,"uin",1,n);
    for i=1:k
        f1(1,i)=p1(i);
        f2(1,i)=p2(i);
    end
    j1=k;
    j2=k;
    for i=1:n
        r1=search(p2(k+1:n),p1(i));
        if r1==1
            j1=j1+1;
            f2(1,j1)=p1(i);
        end
        r2=search(p1(k+1:n),p2(i));
        if r2==1
            j2=j2+1;
            f1(1,j2)=p2(i);
        end
    end    
endfunction

// generarea unui nou traseu; transformare de tip 2-opt
function tnou=mutation(t)
    n=length(t);
    i1=grand(1,1,"uin",1,n); i2=grand(1,1,"uin",1,n);
    if i1>i2
      aux=i1; i1=i2; i2=aux;
    end
    for i=1:i1-1
        tnou(:,i)=t(:,i);
    end
    for i=i1:i2
        tnou(:,i)=t(:,i1+i2-i);
    end
    for i=i2+1:n
        tnou(:,i)=t(:,i);
    end    
endfunction

function [psel,csel]=selectSurvivors(p,c,pnou,cnou)
    psel=[];
    csel=[];
    m=size(p,1);
    for i=1:m
        i1=grand(1,1,"uin",1,2*m);
        if i1<=m
            e1=p(i1,:);
            c1=c(i1,:);
        else    
            e1=pnou(i1-m,:);
            c1=cnou(i1-m,:);
        end    
        i2=grand(1,1,"uin",1,2*m);
        if i2<=m
            e2=p(i2,:);
            c2=c(i2,:);
        else    
            e2=pnou(i2-m,:);
            c2=cnou(i2-m,:);
        end
        if c1<c2
            psel=[psel;e1];
            csel=[csel;c1];
        else
            psel=[psel;e2];
            csel=[csel;c2];
        end    
    end    
endfunction

GA_TSP(50,500,0.5,"eil51.tsp")
