Project 1 - Traveling salesman problem

clc
clear all
citypoints=20;
x=randi(100,citypoints,1);
y=randi(100,citypoints,1);
cities = [x y];
Dmatrix =  toeplitz(randi([0 9],20,1));; %Costs are rounded to integers and symmetric
[NumOfCities,dummy]=size(cities);
Primes=primes(NumOfCities*10);
%     D=diag(inf*ones(1,NumOfCities)); %Assign an inifinite cost to traveling from a city to itself
%     for i=1:NumOfCities %Assign the distances between pairs of cities
%         for j=i+1:NumOfCities
%             D(i,j)=norm(cities(i,:)-cities(j,:));
%              D(j,i)=D(i,j);
%     end
%     end
    D=Dmatrix;

NumOfDataSets=1;
for i=2:NumOfCities
    NumOfDataSets=NumOfDataSets+nchoosek(NumOfCities,i)
end
Data(NumOfDataSets).S=[];
Data(NumOfDataSets).l=0;
Data(NumOfDataSets).cost=inf;
Data(NumOfDataSets).pre=[];
Data(NumOfDataSets).m=[];
LookUpTable(NumOfDataSets)=0;
%Define a data structure that holds the following pieces of data we need
%for later. This data structure uses the same notation used in the paper 
%by Held and Karp (1962):
% S - the set of cities in the tour.
% l - the last city visited in the set S. 
% cost - the cost of a tour, which includes all city in S and ends at l.
%In addition, the following data items are used in the dataset for reducing
%runtime:
% Pre - the index of predecessor dataset, i.e. the one with Spre=S-{l}
% m - the city in S-{l} that yielded the lowest cost C(Spre,m)+D(m,l).
% This index will facilitate the generation of the optimal tour without
% further calculations.
Data(1).S=[1];
Data(1).l=1;
Data(1).cost=0;
Data(1).Pre=[];
Data(1).m=[];
for s=2:NumOfCities
    Data(s).S=[Data(1).S,s];
    Data(s).l=s;
    Data(s).cost=D(s,1);
    Data(s).Pre=1;
    Data(s).m=1;
    LUT=calcLUT(Data(s).S,s,Primes);
    LookUpTable(s)=LUT;
end
IndexStartPrevStep=2; %index into Data that marks the beginning of the previous step
IndexLastStep=NumOfCities; %index into Data that marks the end of the previous step
CurrentData=IndexLastStep; %index into Data that marks the current dataset
%This is the main Dynamic Programming loop
for s=3:NumOfCities
    %generate possible sets with s-1 cities out of the possible N-1 cities
    %(not including city 1)
    TempSets=nchoosek(2:NumOfCities,s-1);
    NumOfSets=size(TempSets);
    for j=1:NumOfSets(1) %iterate through all the sets
        for k=1:NumOfSets(2) %iterate through all the elements in each set
            SminuskSet=[1,TempSets(j,1:k-1),TempSets(j,k+1:NumOfSets(2))]; %this is the set S-{k}               
            candidatecost(2:length(SminuskSet))=inf;
            indices=[];
            for mm=2:length(SminuskSet) %choose a city in S-{k} that will be last
                LUV=calcLUT(SminuskSet,SminuskSet(mm),Primes);
                index=find(LUV==LookUpTable(IndexStartPrevStep:IndexLastStep));  
                index=index+IndexStartPrevStep-1;
                if index==0
                    candidatecost(mm)=inf;                    
                else
                    candidatecost(mm)=Data(index).cost+D(SminuskSet(mm),TempSets(j,k));
                    indices(mm)=index;
                end                
            end
            [mincost,indexcost]=min(candidatecost(2:end));
            CurrentData=CurrentData+1;
            Data(CurrentData).S=[1,TempSets(j,:)];
            Data(CurrentData).l=TempSets(j,k);
            Data(CurrentData).cost=mincost;
            Data(CurrentData).Pre=indices(indexcost+1);
            Data(CurrentData).m=SminuskSet(indexcost+1);
            LookUpTable(CurrentData)=calcLUT(Data(CurrentData).S,TempSets(j,k),Primes);
        end
    end
    IndexStartPrevStep=IndexLastStep+1;
    IndexLastStep=CurrentData;    
end
mm=0;
%Now add the distance back from the last city to city 1
for i=IndexStartPrevStep:IndexLastStep
    mm=mm+1;
    candidatecost(mm)=Data(i).cost+D(Data(i).l,1);
end
%Find the one that minimizes the total distance
[mincost,indexcost]=min(candidatecost);
Temp=Data(IndexStartPrevStep+indexcost-1);
%Generate the optimal tour by traversing back from the last city to its
%predecessors
OptimalTour=1;
while ~isempty(Temp.Pre)
    OptimalTour=[OptimalTour,Temp.l];
    Temp=Data(Temp.Pre);
end
OptimalTour=[OptimalTour,1];
function LUT=calcLUT(vec,last,Primes);
j=length(vec);
LUT=Primes(last);
for i=2:j
    LUT=LUT*Primes(vec(i));
end
end