tours = 5;
alpha = 1;
beta = 5;
ro = 0.5; %evaporation coeff
%cities 3
x = [0 3 6 7  15   10 16 5 8 1.5];
y = [1 2 1 4.5 -1 2.5 11 6 9 12];
N = 10; %number of cities
%number of ants is equal to number of cities
m = N;
%xy = [3,1; 2,4; 12,2; 7,4.5; 9,9; 3,1.5; 16,11; 11,8; 9,10;2,7];
%distance between cities
D = zeros(N);
for i=1:N
    for j=1:N
        D(i,j) = pdist([x(i),y(i); x(j),y(j)]);
    end
end
%display(D);
tau0 = 1/(min(max(D))); %initial pheromones
tau = zeros(N);
for i=1:N
    for j=1:N
        if(i~=j) %if not the same city
            tau(i,j) = tau0;
        end
    end
end

for t=1:tours
    %creating ants
    ant = zeros(m,3+N); %starting city = 1, current city = 2, distance traveled = 3
    for i=1:m
        ant(i,1) = randi([1 10],1,1); %choosing starting city
        ant(i,2) = ant(i,1); %ant is in the same city
        ant(i,3+ant(i, 1)) = 1; %first city in sequence
    end
    dtau = zeros(N); %increase of pheromone
    for iant=1:m
        for icity=2:N+1
            if(icity<N+1)
                numerator = zeros(N);
                denominator = zeros(N,1);
                A = zeros(N);
                for city=1:N
                    for dest=1:N
                        if(city==dest)
                            numerator(city,dest) = 0; %can't go to city you're in
                        else
                            if(ant(iant,3+dest)>0)
                                numerator(city,dest) = 0; %can't go to city you've been
                            else
                                numerator(city,dest) =(tau(city,dest)^alpha)*((1/D(city,dest))^beta);
                                denominator(city) =denominator(city)+numerator(city,dest);
                            end
                        end
                    end
                end
                for city=1:N
                    for dest=1:N
                        if(numerator(city,dest)==0)
                            A(city,dest)=0;
                        else
                            A(city,dest)=numerator(city,dest)/denominator(city);
                        end
                    end
                end

                p = zeros(N);
                pd = zeros(N,1);
                %probability matrix
                for city=1:N
                    for dest=1:N
                        if(ant(iant,3+dest)==0) %if city not yet visited
                            pd(city)=pd(city)+A(city,dest);
                            %denominator sum for probability
                        end
                    end
                end
                for city=1:N
                    for dest=1:N
                        p(city,dest)=A(city,dest)/pd(city);
                        %probability
                    end
                end
                %choosing destination
                choice = 0;
                prob = 0;

                for dest=1:N
                    if(ant(iant,3+dest)==0) %if city not visited
                        if(choice==0)
                            prob = prob + p(ant(iant,2),dest);
                            %probability of current -> dest
                            if(prob>rand
                                choice = dest; %chosen destination

                                break
                            end
                        end
                    end
                end
                dest=choice;
            else
                dest=ant(iant,1); %no cities left, go home
            end
            ant(iant,3) = ant(iant,3) + D(ant(iant,2),dest);
            %distant increases
            ant(iant,2) = dest; %new current city
            if(icity < 11) %end of journey, return to starting city
                ant(iant, 3+dest) = icity;
            end
        end
        %pheromone trail
        for city=1:N
            for dest=1:N
                if(city~=dest)
                    if(ant(iant, 3+dest)-ant(iant, 3+city)==1)
                        dtau(city, dest) = dtau(city, dest) + (1/ant(iant, 3));
                        dtau(dest, city) = dtau(dest, city) + (1/ant(iant, 3)); %direction doesn't matter
                    end %calc of pheromone increase after total distance for ant is known
                end
            end
        end
    end
    %pheromone update
    for city=1:N
        for dest=1:N
            tau(city, dest) = (1 - ro) * tau(city, dest) + dtau(city, dest);
        end
    end
end

%mrufka miszcz
best = 1;
for iant=1:m
    if(ant(iant, 3)<ant(best, 3))
        best = iant;
    end
end

%result
%ant(best, :)
distance = ant(best, 3);
x1 = zeros(11, 1);
y1 = zeros(11, 1);
sequence = zeros(1, 10);
for city=1:N
    for i=4:13
        if(ant(best, i)==city)
            sequence(city) = i-3;
            x1(city) = x(i-3);
            y1(city) = y(i-3);
        end
    end
end

x1(11) = x1(1); %back to starting city
y1(11) = y1(1);
sequence
distance
%plot
hold on
scatter(x, y, 'm')
a = [1:10]'; %labels
b = num2str(a);
c = cellstr(b);
text(x+0.1, y+0.1, c); %0.1 not to overlap
plot(x1, y1, 'c');