function varargout = tsp_ga(varargin)
%TSP_GA  Finds a (near) optimal solution to the Traveling Salesman Problem (TSP)
%   by setting up a Genetic Algorithm (GA) to search for the shortest
%   path (least distance needed to travel to each city exactly once)
%
%   TSP_GA(NUM_CITIES) where NUM_CITIES is an integer representing the number
%   of cities there are (default = 50)
%
%   TSP_GA(CITIES) where CITIES is an Nx2 matrix representing the X,Y
%   coordinates of user specified cities
%
%   TSP_GA(..., OPTIONS) or TSP_GA(OPTIONS) where OPTIONS include one or
%   more of the following in any order:
%     '-NOPLOT'     turns off the plot showing the progress of the GA
%     '-RESULTS'    turns on the plot showing the final results
%   as well as the following parameter pairs:
%     'POPSIZE', VAL  sets the number of citizens in the GA population
%               VAL should be a positive integer (divisible by 4)
%                 -- default = 100
%     'MRATE', VAL  sets the mutation rate for the GA
%               VAL should be a float between 0 and 1, inclusive
%                 -- default = 0.8
%     'NUMITER', VAL  sets the number of iterations (generations) for the GA
%               VAL should be a positive integer
%                 -- default = 500
%
%   Example:
%     % Solves the TSP for 20 random cities using a population size of 60, 
%     % a 75% mutation rate, and 250 GA iterations
%     tsp_ga(20, 'popsize', 60, 'mrate', 0.75, 'numiter', 250);
%
%   Example:
%     % Solves the TSP for 30 random cities without the progress plot
%     [sorted_cities, best_route, distance] = tsp_ga(30, '-noplot');
%
%   Example:
%     % Solves the TSP for 40 random cities using 1000 GA iterations and
%     % plots the results
%     cities = 10*rand(40, 2);
%     [sorted_cities] = tsp_ga(cities, 'numiter', 1000, '-results');
%
%   NOTE: It is possible for TSP_GA to continue where it left off on a
%   previous set of cities by using the sorted city output matrix as an
%   input, as in the following example:
%     cities = 10*rand(60, 2);
%     sorted_cities = tsp_ga(cities, 'numiter', 100);
%     figure; plot(sorted_cities(:, 1), sorted_cities(:, 2), '.-')
%     sorted_cities2 = tsp_ga(sorted_cities);
%     figure; plot(sorted_cities2(:, 1), sorted_cities2(:, 2), '.-')
% 
% Author: Joseph Kirk
% Email: jdkirk630 at gmail dot com
% Release: 1.2
% Release Date: 5/22/07

error(nargchk(0, 9, nargin));
num_cities = 50; num_dims = 2;
cities = 10*rand(num_cities, num_dims);
pop_size = 100; num_iter = 500; mutate_rate = 0.8;
show_progress = 1; show_results = 0;

% Process Inputs
cities_flag = 0; option_flag = 0;
for var = varargin
    if option_flag
        if ~isfloat(var{1}), error(['Invalid value for option ' upper(option)]); end
        switch option
            case 'popsize', pop_size = 4*ceil(real(var{1}(1))/4); option_flag = 0;
            case 'mrate', mutate_rate = min(abs(real(var{1}(1))), 1); option_flag = 0;
            case 'numiter', num_iter = round(real(var{1}(1))); option_flag = 0;
            otherwise, error(['Invalid option ' upper(option)])
        end
    elseif ischar(var{1})
        switch lower(var{1})
            case '-noplot', show_progress = 0;
            case '-results', show_results = 1;
            otherwise, option = lower(var{1}); option_flag = 1;
        end
    elseif isfloat(var{1})
        if cities_flag, error('CITIES or NUM_CITIES may be specified, but not both'); end
        if length(var{1}) == 1
            num_cities = round(real(var{1}));
            if num_cities < 2, error('NUM_CITIES must be an integer greater than 1'); end
            cities = 10*rand(num_cities, num_dims); cities_flag = 1;
        else
            cities = real(var{1});
            [num_cities, num_dims] = size(cities); cities_flag = 1;
            if or(num_cities < 2, num_dims ~= 2)
                error('CITIES must be an Nx2 matrix of floats, with N > 1')
            end
        end
    else
        error('Invalid input argument.')
    end
end

% Construct the Distance Matrix
mat3d1 = reshape(cities, 1, num_cities, num_dims);
mat3d2 = reshape(cities, num_cities, 1, num_dims);
dist_matx = sqrt(sum((mat3d1(ones(num_cities, 1), :, :) - mat3d2(:, ones(num_cities, 1), :)).^2, 3));

% Plot Cities and Distance Matrix in a Figure
if show_progress
    pfig = figure;
    subplot(2, 2, 1)
    plot(cities(:, 1), cities(:, 2), 'b.')
    if num_cities < 75
        for c = 1:num_cities
            text(cities(c, 1), cities(c, 2), [' ' num2str(c)], 'Color', 'k', 'FontWeight', 'b')
        end
    end
    title([num2str(num_cities) ' Cities'])
    subplot(2, 2, 2)
    imagesc(dist_matx)
    title('Distance Matrix')
    colormap(flipud(gray))
end

% Initialize Population
pop = zeros(pop_size, num_cities);
pop(1, :) = (1:num_cities);
for k = 2:pop_size
    pop(k, :) = randperm(num_cities);
end

display_rate = 20;
if num_iter < 50, display_rate = 2; end
fitness = zeros(1, pop_size);
best_fitness = zeros(1, num_iter);
for iter = 1:num_iter
    for p = 1:pop_size
        d = dist_matx(pop(p, 1), pop(p, num_cities));
        for city = 2:num_cities
            d = d + dist_matx(pop(p, city-1), pop(p, city));
        end
        fitness(p) = d;
    end
    [best_fitness(iter) index] = min(fitness);
    best_route = pop(index, :);

    % Plots
    if ~mod(iter, display_rate) && show_progress
        figure(pfig)
        subplot(2, 2, 3)
        route = cities([best_route best_route(1)], :);
        plot(route(:, 1), route(:, 2)', 'b.-')
        title(['Best GA Route (dist = ' num2str(best_fitness(iter)) ')'])
        subplot(2, 2, 4)
        plot(best_fitness(1:iter), 'r', 'LineWidth', 2)
        axis([1 max(2, iter) 0 max(best_fitness)*1.1])
    end

    % Genetic Algorithm Search
    pop = genetic_algorithm(pop, fitness, mutate_rate);
end

if show_progress
    figure(pfig)
    subplot(2, 2, 3)
    route = cities([best_route best_route(1)], :);
    plot(route(:, 1), route(:, 2)', 'b.-')
    title(['Best GA Route (dist = ' num2str(best_fitness(iter)) ')'])
    subplot(2, 2, 4)
    plot(best_fitness(1:iter), 'r', 'LineWidth', 2)
    title('Best Fitness')
    xlabel('Generation')
    ylabel('Distance')
    axis([1 max(2, iter) 0 max(best_fitness)*1.1])
end

if show_results
    figure(2)
    imagesc(dist_matx)
    title('Distance Matrix')
    colormap(flipud(gray))
    figure(3)
    plot(best_fitness(1:iter), 'r', 'LineWidth', 2)
    title('Best Fitness')
    xlabel('Generation')
    ylabel('Distance')
    axis([1 max(2, iter) 0 max(best_fitness)*1.1])
    figure(4)
    route = cities([best_route best_route(1)], :);
    plot(route(:, 1), route(:, 2)', 'b.-')
    for c = 1:num_cities
        text(cities(c, 1), cities(c, 2), [' ' num2str(c)], 'Color', 'k', 'FontWeight', 'b')
    end
    title(['Best GA Route (dist = ' num2str(best_fitness(iter)) ')'])
end

[not_used indx] = min(best_route);
best_ga_route = [best_route(indx:num_cities) best_route(1:indx-1)];
if best_ga_route(2) > best_ga_route(num_cities)
    best_ga_route(2:num_cities) = fliplr(best_ga_route(2:num_cities));
end
varargout{1} = cities(best_ga_route, :);
varargout{2} = best_ga_route;
varargout{3} = best_fitness(iter);

% --- subfunction: genetic algorithm
function new_pop = genetic_algorithm(pop, fitness, mutate_rate)
[p, n] = size(pop);

% Tournament Selection - Round One
new_pop = zeros(p, n);
ts_r1 = randperm(p);
winners_r1 = zeros(p/2, n);
tmp_fitness = zeros(1, p/2);
for k = 2:2:p
    if fitness(ts_r1(k-1)) > fitness(ts_r1(k))
        winners_r1(k/2, :) = pop(ts_r1(k), :);
        tmp_fitness(k/2) = fitness(ts_r1(k));
    else
        winners_r1(k/2, :) = pop(ts_r1(k-1), :);
        tmp_fitness(k/2) = fitness(ts_r1(k-1));
    end
end

% Tournament Selection - Round Two
ts_r2 = randperm(p/2);
winners = zeros(p/4, n);
for k = 2:2:p/2
    if tmp_fitness(ts_r2(k-1)) > tmp_fitness(ts_r2(k))
        winners(k/2, :) = winners_r1(ts_r2(k), :);
    else
        winners(k/2, :) = winners_r1(ts_r2(k-1), :);
    end
end
new_pop(1:p/4, :) = winners;
new_pop(p/2+1:3*p/4, :) = winners;

% Crossover
crossover = randperm(p/2);
children = zeros(p/4, n);
for k = 2:2:p/2
    parent1 = winners_r1(crossover(k-1), :);
    child = winners_r1(crossover(k), :);
    ndx = ceil(n*sort(rand(1, 2)));
    while ndx(1) == ndx(2)
        ndx = ceil(n*sort(rand(1, 2)));
    end
    tmp = parent1(ndx(1):ndx(2));
    for kk = 1:length(tmp)
        child(child == tmp(kk)) = 0;
    end
    child = [child(1:ndx(1)) tmp child(ndx(1)+1:n)];
    child(child == 0) = [];
    children(k/2, :) = child;
end
new_pop(p/4+1:p/2, :) = children;
new_pop(3*p/4+1:p, :) = children;

% Mutate
mutate = randperm(p/2);
num_mutate = round(mutate_rate*p/2);
for k = 1:num_mutate
    ndx = ceil(n*sort(rand(1, 2)));
    while ndx(1) == ndx(2)
        ndx = ceil(n*sort(rand(1, 2)));
    end
    if rand < 0.75 % swap segment between two cities
        new_pop(p/2+mutate(k), ndx(1):ndx(2)) = ...
            fliplr(new_pop(p/2+mutate(k), ndx(1):ndx(2)));
    else % swap two cities
        new_pop(p/2+mutate(k), [ndx(1) ndx(2)]) = ...
            new_pop(p/2+mutate(k), [ndx(2) ndx(1)]);
    end
end

