function traveling_salesman()
global dist x y N col
N1 = 8;
Nclusters = 9;
N = N1*Nclusters;
if mod(N1,2) == 0
    N2 = N/2;
else 
    N2 = (N + 1)/2;
end
c1 = winter(N2);
c = [c1; c1(N2 : -1 : 1,:)];
col = c(1 : N,:);


c = [1, 5, 9];
[xc,yc] = meshgrid(c,c);
xc = xc(:);
yc = yc(:);

x = zeros(N,1);
y = zeros(N,1);
dist = zeros(N);
for j = 1 : Nclusters
    ind = [(j - 1)*N1 + 1 : j*N1];
    x(ind) = xc(j) + 2*rand([N1,1]) - 1;
    y(ind) = yc(j) + 2*rand([N1,1]) - 1;
end
ax = x*ones(1,N);
dx = abs(ax - ax');
ay = y*ones(1,N);
dy = abs(ay - ay');
dist = dx + dy; % distance in the Manhattan metric d(a,b) = |xa - xb| + |ya - yb|


% generate initial path
minpath = randperm(N);
Cmin = mycost(minpath);
fprintf('mincost = %d\n',Cmin);
figure(1);
clf;
hold on;
drawpath(minpath);

% simulated annealing
tem = [12, 8, 4, 2, 1, 0.5, 0.1];
ntem = length(tem);

nflip = 1; % flip nflip + 1 subsequence 
Nf = N - nflip; 
nsmax = N*5000;
namax = 10000*N;
for j = 1 : ntem
    ns = 0; % the # of succeesses
    na = 0; % the # of attempts
    path = minpath;
    C = Cmin;
    fprintf('j = %d\t tem = %d\n',j,tem(j));
    fprintf('path:\n');
    reshape(path,Nclusters,N1)
    fprintf('\n');
    while ns < nsmax | na < namax
% pick randomly  two concequetive cities and reverse the order in which they are passed        
%         flip = ceil(Nf*rand);
%         newpath = path;
%         newpath(flip : flip + nflip) = path(flip + nflip : -1 : flip);
% pick randomly two cities and switch them 
        newpath = path;
        site = ceil(N*rand);
        city = newpath(site);
        newsite = ceil(N*rand);
        newpath(site) = newpath(newsite);
        newpath(newsite) = city;
        Cnew = mycost(newpath);
        na = na + 1;
        if Cnew <= C
            path = newpath;
            C = Cnew;
            if C < Cmin
                Cmin = C;
                minpath = path;
                drawpath(minpath);
                fprintf('mincost = %d\t acceptance ratio = %d\n',Cmin,(ns+1)/na);

            end
            ns = ns + 1;
        else
            eta = rand;
            if eta < exp(-(Cnew - C)/tem(j))
                path = newpath;
                C = Cnew;
                ns = ns + 1;
            end
        end
    end
end
fprintf('path:\n');
reshape(path,Nclusters,N1)
fprintf('\n');
     
end
%% 
function f = mycost(path)
global dist N

p = circshift(path,[0,1]);
ind = path + (p - 1)*N;
f = sum(dist(ind));
end

%%
function drawpath(path)
global x y N col

clf;hold on;
p = circshift(path,[0,1]);
for j = 1 : N
    plot([x(path(j)),x(p(j)),x(p(j))],[y(path(j)),y(path(j)),y(p(j))],...
        'LineWidth',2,'color',col(j,:));
end
plot(x,y,'.','color','r','Markersize',30);
drawnow;
end








    