function out = tsp(loc)
% TSP Traveling salesman problem (TSP) using SA (simulated annealing).
% 	TSP by itself will generate 20 cities within a unit cube and
%	then use SA to slove this problem.
%
%	TSP(LOC) solve the traveling salesman problem with cities'
%	coordinates given by LOC, which is an M by 2 matrix and M is
%	the number of cities.
%
% 	For example:
%
%		loc = rand(50, 2);
%		tsp(loc);

if nargin == 0,
% The following data is from the post by Jennifer Myers (jmyers@nwu.edu)
% to comp.ai.neural-nets. It's obtained from the figure in
% Hopfield & Tank's 1985 paper in Biological Cybernetics
% (Vol 52, pp. 141-152).

	loc = [200,250;200,450;500,300;326,18;200,300;200,400;
           180,350;300,450;300,90;400,60;410,80;400,100;
           350,100;390,110;320,125;360,145;340,160;605,140;
           495,154;550,173;595,178;500,183;475,206;500,220;
           480,220;450,230;530,214;550,204;570,225;568,250;
           516,245;545,255;480,240;480,275;400,250;440,254;
           300,235;345,278;342,228;360,237;430,300;300,305;
           250,290;385,500;300,250;158,275;554,295;200,260;
           280,270;350,260;120,335;300,350;545,325;525,315;
           450,275;375,275;345,300;385,305;520,280;474,307;
           500,335;214,483;410,495;450,478;454,450;454,415;
           480,345;160,380;260,460;380,453;450,345;415,350;
           450,325;470,330;450,380;420,315;415,430;407,330;
           330,290;236,355;245,375;255,320;235,338;240,440;
           360,320;380,320;380,340;315,325;325,340;360,350;
           400,360;270,420;340,450;300,430;340,420;260,390;
           300,370;345,330;270,340;280,290;350,275;360,330;
           390,290;320,370;340,355;325,390;310,400;390,370;
           420,390;410,390;430,400;480,410;360,400;360,380;
           360,370;390,430;380,380];
end
 
NumCity = length(loc);		% Number of cities
distance = zeros(NumCity);	% Initialize a distance matrix
% Fill the distance matrix
for i = 1:NumCity,
	for j = 1:NumCity,
		distance(i, j) = norm(loc(i, :) - loc(j, :));
	end
end

% To generate energy (objective function) from path
%path = randperm(NumCity);
%energy = sum(distance((path-1)*NumCity + [path(2:NumCity) path(1)]));

% Find typical values of dE
count = 125;
all_dE = zeros(count, 1);
for i = 1:count
	path = randperm(NumCity);
	energy = sum(distance((path-1)*NumCity + [path(2:NumCity) path(1)]));
	new_path = path;
	index = round(rand(2,1)*NumCity+.5);
	inversion_index = (min(index):max(index));
	new_path(inversion_index) = fliplr(path(inversion_index));
	all_dE(i) = abs(energy - ...
		sum(sum(diff(loc([new_path new_path(1)],:))'.^2)));
end
dE = max(all_dE);

temp = 10*dE;	% Choose the temperature to be large enough
fprintf('Initial energy = %f\n\n',energy);

% Initial plots
out = [path path(1)];
plot(loc(out(:), 1), loc(out(:), 2),'r.', 'Markersize', 20);
axis square; hold on
h = plot(loc(out(:), 1), loc(out(:), 2)); hold off
 
MaxTrialN = NumCity*100;		% Max. # of trials at a temperature
MaxAcceptN = NumCity*10;		% Max. # of acceptances at a temperature
StopTolerance = 0.005;		% Stopping tolerance
TempRatio = 0.5;		% Temperature decrease ratio
minE = inf;			% Initial value for min. energy
maxE = -1;			% Initial value for max. energy
x_temp = [];
x_energy = [];
x_minmaxE = [];
x_AcceptTrialN = [];
% Major annealing loop
while (maxE - minE)/maxE > StopTolerance,
	minE = inf;
	maxE = 0;
	TrialN = 0;		% Number of trial moves
	AcceptN = 0;		% Number of actual moves 
	while TrialN < MaxTrialN & AcceptN < MaxAcceptN,
		new_path = path;
		index = round(rand(2,1)*NumCity+.5);
		inversion_index = (min(index):max(index));
		new_path(inversion_index) = fliplr(path(inversion_index)); 
		new_energy = sum(distance( ...
			(new_path-1)*NumCity+[new_path(2:NumCity) new_path(1)]));
		if rand < exp((energy - new_energy)/temp),	% accept it!
			energy = new_energy;
			path = new_path;
			minE = min(minE, energy);
			maxE = max(maxE, energy);
			AcceptN = AcceptN + 1;
		end
		TrialN = TrialN + 1;
	end

	% Update plot
	out = [path path(1)];
	set(h, 'xdata', loc(out(:), 1), 'ydata', loc(out(:), 2));
	drawnow;
	% Print information in command window 
	fprintf('temp. = %f\n', temp);
	tmp = sprintf('%d ',path);
	fprintf('path = %s\n', tmp);
	fprintf('energy = %f\n', energy);
	fprintf('[minE maxE] = [%f %f]\n', minE, maxE);
	fprintf('[AcceptN TrialN] = [%d %d]\n\n', AcceptN, TrialN);
	% Lower the temperature
	x_temp=[x_temp;temp];
    x_energy = [x_energy;energy];
    x_minmaxE = [x_minmaxE; [minE maxE]];
    x_AcceptTrialN = [x_AcceptTrialN; [AcceptN TrialN]];
    temp = temp*TempRatio;
    pause(1.2)
end
figure
k=1:length(x_temp);
subplot(2,2,1), plot(x_temp),grid;ylabel('temp');xlabel('iteraciones')
subplot(2,2,2), plot(k,x_energy,k,x_minmaxE,'-.'),grid;ylabel('energy');xlabel('iteraciones')
subplot(2,2,3), plot(x_minmaxE),grid;ylabel('minE - maxE');xlabel('iteraciones')
subplot(2,2,4), plot(x_AcceptTrialN),grid;ylabel('AcceptN - TrialN');xlabel('iteraciones')
% Print sequential numbers in the graphic window
%for i = 1:NumCity,
%	text(loc(path(i),1)+0.01, loc(path(i),2)+0.01, int2str(i), ...
%		'fontsize', 8);
%end
