% gademo.m -- Genetic algorithm and other random search algorithms
% demonatration solving the travelling salesperson 
% N: number of cities. Cities are labeled from 1 to N
% taking advantage of rotation invariance, we start every tour at city #1
% we did not take advantage of symmetricity yet
% implemented algorithms:
% 1. exhaustive search 
% 2. genetic algorithm
% 3. greedy search (Lin and Kernighan)
% 4. random search
% 5. simulated annealing
% call rshift.m, permugen.m
%
% (C) 2001 by Yu Hen Hu
% created: 12/1/2001
% call m-functions: permugen.m, rshift.m, randomize.m, mydist.m
% modified: 12/9/2005

clear all, close all
N=input('# of cities (default = 5, no more than 20) = ');
if isempty(N), N=5; end
if N>20, N=20; end  % not to take too long to run

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (1) generate a TSA problem for N cities
pos=rand(N,2);  % cities located within unit square
citydistance=mydist(pos,pos);
figure(1),clf
subplot(231),plot(pos(:,1),pos(:,2),'o')
axis([0 1 0 1]), axis square, hold on
for i=1:N
   text(pos(i,1)+.03,pos(i,2),int2str(i));
end
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% exhaustive search
disp('********** EXHAUSTIVE SEARCH METHOD'); 

costmin=N; % initialize with an upper bound of cost
if N <=10,
   tmp=permugen([2:N]);
   alltour=[ones(size(tmp,1),1) tmp ];  % each row is a tour
   for i=1:factorial(N-1),
      cost(i)=sum(diag(citydistance(alltour(i,:)',rshift(alltour(i,:))')));
   end
   [costmin,idx]=min(cost);
   tourmin=alltour(idx,:);
   disp(['cost function evaluation: ' int2str(factorial(N-1)) ' times!'])
   disp(['minmum trip length = ' num2str(costmin)])
   disp('optimum tour = ')
   disp(num2str(tourmin))
   
   
   subplot(232),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on
   trip=pos(tourmin',:); trip=[trip; trip(1,:)];
   plot(trip(:,1),trip(:,2),'-')
   title(['Exhaustive search, cost = ' num2str(costmin)])
   hold off
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Genetic search
disp('********** GENETIC ALGORITHM');
ngen=input('# of generations to evolve (default = 1000) = ');
if isempty(ngen), ngen=1000; end
ngpool=input(['# of chromosoms in the gene pool  (default = ' int2str(N) ') = ']); % size of gene pool
if isempty(ngpool), ngpool = N; end
gpool=zeros(ngpool,N); % gene pool
for i=1:ngpool, % intialize gene pool
   gpool(i,:)=[1 randomize([2:N]')'];
end

costmin=N, tourmin=zeros(1,N); cost=zeros(1,ngpool); 
for i=1:ngen, 
   % step 1. evaluate fitness of current chromosoms
   % for TSP problem, it is the trip length. The small the better
   for i=1:ngpool,
      cost(i)=sum(diag(citydistance(gpool(i,:)',rshift(gpool(i,:))')));
   end
   % record current best solution
   [costmin,idx]=min(cost);
   tourmin=gpool(idx,:);
   
   % step 2. cross breeding and mutation
   % since each chromosom is an integer vector, cross-breeding 
   % can not be
   % done by simply combining binary vectors.
   % Instead, we let the off-spring to keep the common genes 
   % of parents, and
   % randomly shuffle genes when parents disagree
   % for simplicity, we only cross-breed the two best solutions 
   
   [csort,ridx]=sort(cost); % sort from small to big.
   mother=gpool(ridx(1),:); father=gpool(ridx(2),:); % parents
   % find index of genes that are the same in father and mother
   sameidx=[father==mother]; % include index #1, a mask vector
   diffidx=find(sameidx==0); % idx of different elements
   if length(diffidx) <=4, % father and mother too close!
      % do mutation
      boy=[1 randomize([2:N]')'];
      girl=[1 randomize([2:N]')'];
   else
      boy=father.*sameidx; boy(diffidx)=randomize(father(diffidx)')';
      girl=mother.*sameidx; girl(diffidx)=randomize(mother(diffidx)')';
   end
   gpool=[gpool(ridx(1:ngpool-2),:);boy;girl]; % replace the worst two
end

disp(['cost function evaluation: ' int2str(ngen*2+ngpool) ' times!'])
disp(['minmum trip length = ' num2str(costmin)])
disp('optimum tour = ')
disp(num2str(tourmin))

subplot(233),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on
trip=pos(tourmin',:); trip=[trip; trip(1,:)];
plot(trip(:,1),trip(:,2),'-')
title(['Genetic Search, cost = ' num2str(costmin)])
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lin and Kernighan hill climbing result
% Use local interchange to search for steepest descending local minimum solution
% 
disp('********** GREEDY SEARCH METHOD'); 

nn=input('local neighborhood size (default=2) = ');
if isempty(nn), nn=2; end
tour0=[1 randomize([2:N]')']; % random initial tour
locmin=N; 
converged = 0; iter=0;
while converged == 0, 
   iter=iter+1;
   % generate neighborhood trips
   tmp=nbgen(tour0(2:N),nn);
   nbtour=[tour0; [ones(size(tmp,1),1) tmp]]; % all the neighboring tours and initial tour
   ntour=size(nbtour,1);
   cost=zeros(1,ntour);
   for i=1:ntour, 
      cost(i)=sum(diag(citydistance(nbtour(i,:)',rshift(nbtour(i,:))')));
   end
   [locost,idx]=min(cost); % 1 step steepest descent search result
   tour=nbtour(idx,:); 
   if locost < locmin, % if still improvement, continue to next iteration
      locmin=locost; loctour=tour; tour0=tour;
   elseif locost == locmin, % a local minimum is reached
      converged=1; costmin=locost; tourmin=tour; 
   end
end

ncomplexity=iter*prod([N-1:-1:N-nn])/factorial(nn);
disp(['cost function evaluation: ' ...
      int2str(ncomplexity) ' times!'])
disp(['minmum trip length = ' num2str(costmin)])
disp('optimum tour = ')
disp(num2str(tourmin))

subplot(234),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on
trip=pos(tourmin',:); trip=[trip; trip(1,:)];
plot(trip(:,1),trip(:,2),'-')
title(['Greedy search, cost = ' num2str(costmin) ', nb size = ' int2str(nn)])
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% randomized search
disp('********** RANDOM SEARCH METHOD'); 
ntrial = input('max # of search steps = '); 

tour0=[1 randomize([2:N]')']; % initial tour
costmin=N;
for nt=1:ntrial,
   cost(nt)=sum(diag(citydistance(tour0',rshift(tour0)')));
   if cost(nt)<costmin,
      costmin=cost(nt);
      tourmin=tour0;
   end
   tournew=[1 randomize(tour0(2:N)')']; 
end

disp(['cost function evaluation: ' int2str(ntrial) ' times!'])
disp(['minmum trip length = ' num2str(costmin)])
disp('optimum tour = ')
disp(num2str(tourmin))

subplot(235),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on
trip=pos(tourmin',:); trip=[trip; trip(1,:)];
plot(trip(:,1),trip(:,2),'-')
title(['Random search, cost = ' num2str(costmin)])
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SA search
disp('********** Simulated Annealling SEARCH METHOD'); 

ntrial = input('max # of search steps (default 1000) = '); 
if isempty(ntrial), ntrial=1000; end
nn=input('neighborhood size (default 2) = ');
if isempty(nn), nn=2; end
% cooling schedule
T=1./(1+[1:ntrial]/1000); 

tour0=[1 randomize([2:N]')']; % initial tour
costmin=sum(diag(citydistance(tour0',rshift(tour0)'))); % initial cost
tourmin=tour0; 
converged = 0; 
nt=0; ccount=0; ccountmax=10;
while converged==0,
   nt=nt+1;
   % generate a new tour in the neighborhood of current tour0
   tmp=randomize(nbgen(tour0(2:N),nn)); 
   tour=[1 tmp(1,:)]; % pick only one
   cost=sum(diag(citydistance(tour',rshift(tour)'))); % cost of tour
   if cost <= costmin, % if reduction of cost,
      tour0=tour; costmin=cost; tourmin=tour0; 
      % count how long there is no further improvement
      % reset ccount when there is an improvement made
      if cost==costmin, ccount=ccount+1; else ccount=0; end
   elseif cost > costmin, % if no reduction
      % toss a dice
      dice=rand; % range between 0 and 1
      if dice <= T(nt), % accept higher cost
         tour0=tour; 
      end
   end
   % convergence test
   if nt == ntrial | ccount== ccountmax, 
      converged=1;
   end
end
disp(['cost function evaluation: ' int2str(nt) ' times!'])
disp(['minmum trip length = ' num2str(costmin)])
disp('optimum tour = ')
disp(num2str(tourmin))

subplot(236),plot(pos(:,1),pos(:,2),'o'),axis([0 1 0 1]), axis square, hold on
trip=pos(tourmin',:); trip=[trip; trip(1,:)];
plot(trip(:,1),trip(:,2),'-')
title(['SA search, cost = ' num2str(costmin)])
hold off