--- do_experiments.m
% Usage of heuristics in console:
% [p,l,d]=shortest_path('d198.txt');
% [p2,l2,d2]=insertion_heuristics('d198.txt',@farthest_insertion);
% [p2,l2,d2]=insertion_heuristics('d198.txt',@arbitrary_insertion);
 
function results=do_experiments()
data_files={'d198.txt';'d1291.txt';'rat575.txt'};
number_of_experiments=10;
results=[];
for idx_data=1:length(data_files)
    for i=1:number_of_experiments
        data=char(data_files(idx_data));
        [path,total_length(1),dist]=shortest_path(data);
        [path,total_length(2)]=opt2(path,dist);
        [path,total_length(3),dist]=insertion_heuristics(data,@farthest_insertion);
        [path,total_length(4)]=opt2(path,dist);
        [path,total_length(5),dist]=insertion_heuristics(data,@arbitrary_insertion);
        [path,total_length(6)]=opt2(path,dist);
        results(i+(idx_data-1)*number_of_experiments,:)=total_length;
    end
end
 
--- shortest_path.m
% calculates shortest path using nearest neighborhood
function [path,total_length,distances] = shortest_path(data_file,idx_initial_city)
distances=distance(read_location_data(data_file));
n=length(distances);
if nargin < 2
    idx_initial_city=ceil(rand(1)*n);
end
path=[];
outside=(1:n);
[path,outside]=add_to_path(path,outside,outside(idx_initial_city));
for i = 1:n-1
    entering_city=nearest_neighbor(distances,path(end),outside);
    [path,outside]=add_to_path(path,outside,entering_city);
end
total_length=total_length_of_cycle(distances,path);
 
function [path,outside] = add_to_path(path,outside,city)
path(end+1)=city;
outside(outside==city)=[];
 
--- distance.m
function d = distance(a)
n=length(a);
d=zeros(n,n);
for i = 1:n
    for j = i+1:n
        d(i,j)=norm(a(i,:)-a(j,:));
        d(j,i)=d(i,j);
    end;
end;
 
--- read_location_data.m
function positions = read_location_data(data_file)
load(data_file);
filename=regexp(data_file,'\w*(?=\.)','match');
positions=eval(filename{1});
 
--- nearest_neighbor.m
function nearest_city = nearest_neighbor(distances,last_city,outside)
    min_distance = inf;
    for j=1:length(outside)
        neighbors_distance = distances(outside(j),last_city);
        if neighbors_distance < min_distance
            min_distance = neighbors_distance;
            nearest_city = outside(j);
        end
    end
end
 
--- total_length_of_cycle.m
function total_length = total_length_of_cycle(distances,path)
total_length = 0;
for i = 1:length(path)-1
    total_length = total_length + distances(path(i),path(i+1));
end
total_length = total_length + distances(path(end),path(1));
 
--- insertion_heuristics.m
% generic insertion algorithm
function [path,total_length,distances] = insertion_heuristics(data_file,insertion_rule,idx_initial_city)
distances=distance(read_location_data(data_file));
n=length(distances);
if nargin < 3
    idx_initial_city=ceil(rand(1)*n);
end
path=[];
outside=(1:n);
[path,outside]=add_to_path(path,outside,outside(idx_initial_city),1);
for i = 1:n-1
    entering_city=insertion_rule(path,outside,distances);
    entry_position=find_entry_position(distances,entering_city,path);
    [path,outside]=add_to_path(path,outside,entering_city,entry_position);
end
total_length=total_length_of_cycle(distances,path);
 
function [path,outside] = add_to_path(path,outside,entering_city,entry_position)
old_path=path;
outside(outside==entering_city)=[];
path(entry_position)=entering_city;
for i=entry_position:length(old_path)
    path(i+1)=old_path(i);
end
 
function entry_position=find_entry_position(distances,entering_city,path)
min_increase_in_length = inf;
global path_length
path_length = length(path);
for i=1:path_length
    before = distances(path(i),path(r(i+1)));
    after = distances(entering_city,path(i))+distances(entering_city,path(r(i+1)));
    increase_in_length = after-before;
    if increase_in_length < min_increase_in_length
        min_increase_in_length = increase_in_length;
        entry_position = r(i+1);
    end
end
 
% if index is greater than the path length, turn the index for one round
function result = r(index)
global path_length
if index > path_length
    result = index - path_length;
else
    result = index;
end
 
--- farthest_insertion.m
function entering_city=farthest_insertion(path,outside,distances)
max_distance=0;
for i=1:length(path)
    for j=1:length(outside)
        if distances(path(i),outside(j)) > max_distance
            max_distance=distances(path(i),outside(j));
            entering_city=outside(j);
        end
    end
end
 
--- arbitrary_insertion.m
function entering_city = arbitrary_insertion(path,outside,distances)
entering_city=outside(ceil(rand(1)*length(outside)));
 
--- opt2.m
function [path,total_length] = opt2(path,distances)
global n 
n = length(path);
for i=1:n
    for j=1:n-3
        if change_in_path_length(path,i,j,distances) < 0 
            path=change_path(path,i,j);
        end
    end
end
total_length = total_length_of_cycle(distances,path);
 
function result = change_in_path_length(path,i,j,distances)
before=distances(path(r(i)),path(r(i+1)))+distances(path(r(i+1+j)),path(r(i+2+j)));
after=distances(path(r(i)),path(r(i+1+j)))+distances(path(r(i+1)),path(r(i+2+j)));
result=after-before;
 
function path = change_path(path,i,j)
old_path=path;
% exchange edge cities
path(r(i+1))=old_path(r(i+1+j));
path(r(i+1+j))=old_path(r(i+1));
% change direction of intermediate path 
for k=1:j-1
    path(r(i+1+k))=old_path(r(i+1+j-k));
end
 
% if index is greater than the path length, turn index one round off
function result = r(index)
global n
if index > n
    result = index - n;
else
    result = index;
end