function tsp

%%% Using a Genetic Algorithm to Find an Approximate Solution to the Traveling Salesman Problem.

%%%By Evante Garza-Licudine for EPS 109, "Computer Simulations in Earth and
%%%Planetary Science." Taught by Burkhard Militzer, fall 2009 semester at
%%%UC Berkeley.

%%%%%

clear
clf
clc

global n pathtotal rm rc generation

%Parameters
n=40; %Mostly for size of population, or number of cities.
pathtotal=15; %Number of paths you want in the initial population minus one (so if you want 2 paths, set this to 1).
rm=0.05; %Set the mutation rate as a probability
rc=0.85; %Set the crossover rate as a probability
locx=zeros(1,n);
locy=locx;
generation=0;
fraem=0;

%Create n randomly-placed cities on an nxn grid.
for city=1:n
    locx(city)=rand*n;
    locy(city)=rand*n;
end

%determine initial population of paths
initialpaths

newpath=fpath;

%the main loop
while generation<=10000;
    fpath=calcorder(fpath); %order the current population by fitness
    if mod(generation,200)==0
        plot(fpath(1,1:n),fpath(2,1:n),'o-',fpath(1,1),fpath(2,1),'y*') %plot the most fit path
        h1=figure(2);
        xlim([0 n])
        ylim([0 n])
        fraem=fraem+1;
        saveas(h1,strcat('evante_garza_',num2str(999+fraem),'.jpg'));
        fprintf('Frame %4d of %4d\n',fraem,10000/200);
        hold=calcfit(fpath(1,:),fpath(2,:));
        
        generation;
        hold;
    end
    
    for q=1:pathtotal %calculate reproduction, mutation
        
        if rand<rc
            fpath(2*q+1:2*(q+1),:)=calccross(fpath(2*q-1,:),fpath(2*q,:),newpath(2*q+1,:),newpath(2*(q+1),:));
        end
        
    end
    
    for q=1:pathtotal+1
        
        if rand<rm
            newpath(2*q-1:2*q,:)=calcmut(newpath(2*q-1,:),newpath(2*q,:));
        end
        
    end
    
    generation=generation+1;
    
end

    function plotpath %plots every path in matrix fpath
        for k=1:pathtotal+1
            subplot(pathtotal+1,1,k)
            plot(fpath(2*k-1,1:n),fpath(2*k,1:n),'-o',fpath(1,1),fpath(2,1),'y*') %if you plot to 1:n+1, the path will be complete
            xlim([0 n])
            ylim([0 n])
        end
    end

    function initialpaths %creates an initial population of paths
        for pathnum=0:pathtotal %number of paths
    
            drawsetx=locx;
            drawsety=locy;
            pathx(1)=drawsetx(1);
            pathy(1)=drawsety(1);
    
            for in=1:n %number of cities in each path
        
                A=ceil(rand*(n-in))+1;
        
                pathx(in+1)=drawsetx(A); %Choose a random place to permute
                storeA=drawsetx(A); %remember Ath value
                storeN=drawsetx(n+1-in); %remember last value
                drawsetx(n)=storeA; %switch last value with Ath value
                drawsetx(A)=storeN; %switch Ath value with last value
                drawsetx=drawsetx(1:n-in); %discard the used Ath value

                pathy(in+1)=drawsety(A); %Choose a random place to permute
                storeA=drawsety(A); %remember Ath value
                storeN=drawsety(n+1-in); %remember last value
                drawsety(n)=storeA; %switch last value with Ath value
                drawsety(A)=storeN; %switch Ath value with last value
                drawsety=drawsety(1:n-in); %discard the used Ath value

            end

            fpath(2*pathnum+1,:)=pathx; %store x-values in one giant matrix
            fpath(2*(pathnum+1),:)=pathy; %store y-values in one giant matrix

        end
    end

    function mutation=calcmut(x,y) %mutates a path
        xprime=x(2:n);
        yprime=y(2:n);
        
        a=ceil(rand*(n-1));
        b=ceil(rand*(n-1));
        while b==a
            b=ceil(rand*(n-1));
        end
        
        storeXa=xprime(a);
        storeYa=yprime(a);
        storeXb=xprime(b);
        storeYb=yprime(b);
        
        xprime(a)=storeXb;
        yprime(a)=storeYb;
        xprime(b)=storeXa;
        yprime(b)=storeYa;
        
        mutation=[x(1) xprime x(n+1);y(1) yprime y(n+1)];
    end

    function crossover=calccross(xm,ym,xf,yf)%performs a chromosomal crossover between two parents
        xmprime=xm(2:n);
        ymprime=ym(2:n);
        xfprime=xf(2:n);
        yfprime=yf(2:n);
        
        a=ceil(rand*(n-3)); %curiously, for a close to n, child is always a clone
        
        xchild=xmprime(1,1:a);
        ychild=ymprime(1,1:a);
        
        complementx=[];
        complementy=[];
        
        if a<n-1 %if a=n-1, then offspring will just be a maternal clone
            truth=zeros(a,n-1);
            k=0;
            for i=1:n-1 %loop over father's alleles
                for j=1:a %loop over offpsring's current alleles
                    truth(j,i)=xfprime(i)==xchild(j);
                end
                if max(truth(:,i))==0
                    k=k+1;
                    complementx(k)=xfprime(i);
                    complementy(k)=yfprime(i);
                end
            end
            
            xchild=[xchild complementx];
            ychild=[ychild complementy];
        end
        
        crossover=[xm(1) xchild xm(n+1);ym(1) ychild ym(n+1)];
    end

    function fitnessmatrix=calcorder(matrix) %order a matrix by its fitness in decreasing order
        M=size(matrix);
        My=M(1);
        fitnesses=zeros(1,My/2);
        finalmatrix=zeros(My,M(2));
        
        for loop=1:My/2
            fitnesses(loop)=calcfit(matrix(2*loop-1,:),matrix(2*loop,:));
        end
        
        sortedpaths=sortvector(fitnesses);
        S=size(sortedpaths);
        Sx=S(2);

        for pool=1:Sx
            for j=1:My/2
                if sortedpaths(pool)==calcfit(matrix(2*j-1,:),matrix(2*j,:))
                    finalmatrix(2*pool-1,:)=matrix(2*j-1,:);
                    finalmatrix(2*pool,:)=matrix(2*j,:);
                end
            end
        end
        fitnessmatrix=finalmatrix;
    end

    function fitness=calcfit(x,y) %calculates fitness of a path
        sum=0;
        sumnew=0;
        p=size(x);
        p=p(2);
        for Q=1:p-1
            sumnew=norm([x(Q) y(Q)]-[x(Q+1) y(Q+1)]);
            sum=sum+sumnew;
        end
        fitness=sum;
    end

    function sorted=sortvector(a)

    A=size(a);
    N=A(2);
    k=0;

        while N>=1
            for p=1:N
                if a(p)==min(a)
                    k=k+1;
                    best=a(p);
                    last=a(N);
                    a(p)=last;
                    a(N)=best;
                    S(k)=best;
                    a=a(1:N-1);
                    break
                end
            end
            N=N-1;
        end

    sorted=S;
    end


end