% cities2
x = [0 2 6 7  15   12 14 9.5 7.5 0.5];
y = [1 3 5 2.5 -0.5 3.5 10 7.5 9 10];

evaporation = 0.1;
tau0 = 0.1;
alpha = 1;
beta = 5;
ro = 0.5;
m = 10; % number of ants
n = m; % number of towns equal to the number of ants
Tmax = 200; % number of tours (cycles)

mostestShortest = inf; % best path
% Pheromone matrix with tau0 on every path
T(1:n, 1:n) = tau0;
T = T - diag(diag(T));
D = zeros(n, n); % distances between towns
L = zeros(m, 1); % length of each ants path

% distance between towns
for i=1:n
    for j=1:n
        D(i,j) = pdist([x(i),y(i); x(j),y(j)]);
    end
end

for i=1:Tmax
    %towns each ant has visited
    C = zeros(m, n);
    %choose a random initial town for each ant
    for k=1:m
        C(k, 1) = randi(n);
    end
    for a=1:m
        for p=2:n
            town = C(a,(find(C(a,1:n), 1, 'last'))); % find in which town the ant is right now, which is the last non-zero value
            remTowns = 1:size(T, 1);  % vector of towns the ant can still go to
            remTowns = setdiff(remTowns, C(a,1:n)); % difference of those two vectors
            A = zeros(n); % decision table
            %Calculate sum (denominator)
            denom = 0;
            for j=1:size(remTowns, 2)
                denom = denom + (T(town, remTowns(j))^alpha)*(1/D(town,remTowns(j))^beta);
            end
            % find the decision vector A for this ant at this moment
        	% assign a value to every available town
            for j=1:size(remTowns, 2)
                A(remTowns(j)) = ((T(town,remTowns(j))^alpha)*(1/D(town, remTowns(j))^beta))/denom;
            end
            %probability for choosing each path
            Aprob = A/sum(A);
            % choose one path
            C(a,p)= randsample( 1:n, 1, true, Aprob );
        end
        %Calculates the length of this tour
        for z=1:n-1
            L(a) = L(a) + D(C(a,z), C(a,z+1));
            if(z==n-1)
                L(a) = L(a) + D(C(a,z), C(a,1));
            end
        end
    end
    %Finds the optimal tour in this cycle
    [shortestLength, minI] = min(L);
    %check if the optimal tour in this cycle is also the best overall
    if shortestLength < mostestShortest
        mostestShortest = shortestLength;
        bestPath = C(minI, :);
        bestPath(n+1) = bestPath(1);
    end
    tauD = zeros(n);

    %for every ant
    for a=1:m
        %for every arc in the path
        for p=1:n-1
            %how much pheromone is added to this arc
            tauD(C(a,p), C(a,p+1)) = tauD(C(a,p), C(a,p+1)) + 1/L(a);
        end
        %ant returns to the initial town
        tauD(C(a,n), C(a,1)) = tauD(C(a,n), C(a,1)) + 1/L(a);
    end

    %update the pheromones
    for p=1:n
        for k=1:n
            T(p, k) = (1-evaporation)*T(p,k) + tauD(p, k);
        end
    end
end

bestPath
mostestShortest

%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(x(bestPath), y(bestPath), 'c');