function [totalDist] = tsp2(mapFile)
% Traveling Salesperson Problem (TSP)
% mapFile is the name (string) of a jpeg file showing a map.
% totalDist is the tour distance of the TSP based on user-selected cities.
% This version calculates a distance matrix.

% Show map, calibrate distance, and select cities
close all
m= imread(mapFile);
image(m)
axis equal off tight
set(gcf,'position',[50 50 800 600])
% Calibrate distance
  title('Click on the two ends of the map scale', 'fontsize', 16)
  [x,y]= ginput(2);
  mapUnit= input('Enter the unit of measurement: ', 's');
  scaleDistance= input('Enter the distance: ');
  mapScale= scaleDistance/abs(x(2)-x(1));
% Prompt user to select cities to visit
  prompt= 'Select cities to visit, beginning with the starting point.  ';
  title([prompt 'Press <Enter>.'], 'fontsize', 16)
  [x,y]= ginput;

n= length(x);              % Number of cities
D= allPairDistances(x,y);  % Distance matrix for all cities

% Determine the "shortest" path to visit all selected cities
path= 1;          % Path begins with city#1
unvisited= 2:n;   % Cities not yet visited
totalDist= 0;     % Total distance traveled so far
hold on
plot(x(1),y(1),'md','markersize',10,'MarkerFaceColor','m')
title(sprintf('Distance traveled = %.1f', totalDist), 'fontsize', 18)

for k= 2:n  % Decide on the kth city of the tour
    pause(.5)
    % Find nearest city from city last visited, path(k-1)
    [city, distance]= nearest(path(k-1),unvisited,D(path(k-1),:));
    unvisited= removePt(city,unvisited);
    totalDist= totalDist + distance;
    path= [path city];
    plot(x(path(k)),y(path(k)),'k^',...
         [x(path(k)) x(path(k-1))], [y(path(k)) y(path(k-1))], ...
                                               'k-','Linewidth',2)
    title(sprintf('Distance traveled = %.0f %s', ...
          totalDist*mapScale, mapUnit), 'fontsize', 18)
end

% Go back to the starting city
plot([x(path(n)) x(1)], [y(path(n)) y(1)], 'k-','Linewidth',2)
totalDist= totalDist + D(path(n),1);

% Scale final distance to mapScale
totalDist= totalDist*mapScale;
title(sprintf('Distance traveled = %.0f %s', ...
          totalDist, mapUnit), 'fontsize', 18)


%%%%%%%%
function d = allPairDistances(x,y)
% x,y are vectors of the same length and represent x- and y-coordinates.
% d is a matrix of pairwise distances.
n= length(x);
d= zeros(n,n);
% Calculate only upper triangle since d(r,c) is same as d(c,r)
for r= 1:n;
    for c= r+1:n;
        d(r,c)= calDist(x(r),y(r), x(c),y(c));
    end
end 
d= d+d  % Fill in lower triangle (mirror image of upper triangle)


%%%%%%%%
function [c, minDist] = nearest(h,p,distVec)
% Find the point that is nearest to point h.
% [p(1) p(2) ... p(n)] are the indices of the candidate points.
% distVec is a vector of distances such that 
%   the distance between point h and point p(k) is distVec(p(k)).
% c is the index of the point nearest to point h.
% minDist is the distance between points h and c.

minDist= distVec(p(1));
c= p(1);
for k= 2:length(p)
    if distVec(p(k))<minDist
        minDist= distVec(p(k));
        c= p(k);
    end
end

    
%%%%%%%%
function d = calDist(x1,y1, x2,y2)
% d is the distance between (x1,y1) and (x2,y2)
d= sqrt( (x1-x2)^2 + (y1-y2)^2 );
   
    
%%%%%%%%
function w = removePt(val, vec)
% w is vector vec the value val removed
w= [];
for k= 1:length(vec)
    if vec(k)~=val
        w= [w vec(k)];
    end
end
% Note:
% The whole function could be a single statement using vectorized code:
%    w= vec(vec~=val);