function [X_opt , S_opt , Pt , T] = TSP_CE(L , X1 , option)
% Solve the classic Travel Sale Man problem with a Cross Entropy
% (CE) optimization algorithm.
%
% Usage :
%
% [X_opt , S_opt , Pt , T] = TSP_CE(L , [X1] , [option])
%
% Inputs
% -------
%
% L : Distance/cost matrix (n x n)
% X1 : (optional) start city (default X1 = 1);
% option : (optional) structure for the option of the CE algorithm
%
% option.N : Number of paths to generate for each iteration (default 3*n*n)
% option.rho : Threshold for selectionning the best rho*N paths (default 5*10e-3)
% option.alpha : Weight for updating Transition Matrix (default 0.8)
% option.beta : Determine the fraction of the Nrho*N sample who will be kept for each iteration (default 0.7)
% option.d : Indicate how many CE iterations to wait until declaring a solution is found (default 5)
% option.T_max : Maximum number of iterations (default 10000)
% option.epsilon : Absolute relative difference between 2 consecutives solution
% option.digit : Truncation rank (default 5)
% option.verbose : 1 for real time plot;
% option.V : Towns coordinates (n x 2). option.verbose = 1 is requiered;
%
% Outputs
%
% X_opt : optimal solution. a (n + 1 x 1) index vector.
% S_opt : Distance associated with the optimal path
% Pt : Transition Matrix (n x n) estimated at the end of the CE algorithm
% T : Number of CE iteration needed( < option.T_max)
%
% Example
%
% n = 20;
% V = 10*rand(n , 2);
% L = distmat(V);
% X_opt = TSP_CE(L);
% plot(V(: , 1) , V(: , 2) , 'b+' , V(X_opt , 1) , V(X_opt , 2) , 'r' , 'markersize' , 8 , 'linewidth' , 2)
% text(V(: , 1) , V(: , 2) , num2str((1 : n)' , '%0.0f') , 'fontname' , 'times' , 'fontsize' , 13);
%
% Author : Sbastien PARIS
%
% Ref : "The Cross Entropy method", P-T de Boer, D. P. Kroese, S.
% Mannor, R. Y. Rubinstein.
%
% see http://iew3.technion.ac.il/CE/about.php
[n , d] = size(L);
if((n ~= d))
error('L must be a square(n x n)');
end
if ((nargin < 2) || isempty(X1) )
X1 = 1;
end
if ( (X1 < 1) || (X1 > n))
error('Start city index must be within [1, ... , n]');
end
if ( (nargin < 3) || isempty(option))
option.N = round(5*n*n);
option.rho = 5*10e-3; % threshold for selectionning the best rho*N paths
option.alpha = 0.7; % weight for updating Transition Matrix with new paths candidates
option.beta = 0.7; % Determine the fraction of the Nrho*N sample who will be keep in each iteration
option.d = 5; % indicate how many iterations to wait until declaring the solution is reached
option.T_max = 10000; % Maximum number of CE iterations
option.epsilon = 100*eps;
option.digit = 5;
option.verbose = 0; % Display Best current solution of the CE algorithm
end
if (option.verbose == 1)
if (~any(strcmp(fieldnames(option) , 'V')))
disp('option.V must be provided tp verbose, switching verbose = 0');
option.verbose = 0;
end
fig1 = figure(1);
set(fig1 , 'renderer' , 'zbuffer');
set(fig1 , 'doublebuffer' , 'on');
subplot(121)
set(gca , 'nextplot' , 'replace');
subplot(122)
set(gca , 'nextplot' , 'replace');
str_n = num2str((1 : n)' , '%0.0f');
end
Nrho = round(option.rho*option.N);
betaR1 = round(option.beta*Nrho);
ind_Nrho = (1 : Nrho);
ind_betaR1 = (1 : betaR1);
ind_betaR2 = (betaR1 + 1 : option.N);
ind_n = (n + 1: -1 : 1);
ctelog = 1/log(10);
NN = option.N - betaR1;
%%==== Initial Transition Matrix =====%%
cten1 = 1/(n - 1);
Pt = cten1(ones(n));
Pt(1 : n + 1 : n*(n + 1)) = 0;
%%===== Initialization iteration ========%%
X = generation_TSP(Pt , option.N , X1);
[S , indice] = sort(cost_TSP(L , X));
Gt = S(1);
I = indice(ind_Nrho);
tp_X = X(: , I);
%a) troncate to the N digits the Nrho first values of S %
x = abs(S(ind_Nrho));
m = ceil(log(x)*ctelog);
k = (x <= 10.^(m - 1));
m(k) = m(k) - 1;
f = 10.^(option.digit - m);
x = fix(x.*f)./f;
%b) Find sequential identical values
ind2 = find((x(2 : Nrho) - x(1 : Nrho - 1)) == 0) + 1;
if ( ~isempty(ind2) )
ind1 = [ind2(1) - 1 , ind2((ind2(2 : end) - ind2(1 : end - 1) ~= 1)) + 1];
ind3 = ind2((tp_X(2 , ind2) < tp_X(n , ind2)));
tp_X(: , ind3) = tp_X(ind_n , ind3);
ind4 = ind1(((tp_X(2 , ind1) ~= tp_X(2 , ind1 + 1))));
tp_X(: , ind4) = tp_X(ind_n , ind4);
end
Pt = (1 - option.alpha)*update_matrix(tp_X) + option.alpha*Pt ;
X_opt = X(: , I(1));
S_opt = S(1);
t = 1;
cons = 0;
Gt_old = Gt;
X(: , ind_betaR1) = tp_X(: , ind_betaR1);
if (option.verbose)
maxi = max(option.V);
mini = min(option.V);
fig1 = figure(1);
subplot(121);
p1 = plot(option.V(: , 1) , option.V( : , 2) , '+' , 'markersize' , 7);
text(option.V(: , 1) , option.V(: , 2) , str_n , 'fontname' , 'times' , 'fontsize' , 12)
axis([mini(1) , maxi(1) , mini(2) , maxi(2)])
drawnow
subplot(122);
bar3(Pt);
axis([0 , n+1 , 0 , n+1 , 0 , 1])
view([-50 50]);
title(sprintf('\\bf{P_{%d}}' ,t ) , 'fontsize' , 12 , 'fontname' , 'times');
drawnow
end
%%========== Main Loop ==========%%
while( (t < option.T_max) && (cons < option.d))
X(: , ind_betaR2) = generation_TSP(Pt , NN , X1);
[S , indice] = sort(cost_TSP(L , X));
Gt_new = S(1);
if ((abs(Gt_new - Gt_old) < option.epsilon) && (Gt_new == Gt_old) );
cons = cons + 1;
else
cons = 0;
end
I = indice(ind_Nrho);
tp_X = X(: , I);
%a) troncate to the N digits the Nrho first values of S %
x = abs(S(ind_Nrho));
m = ceil(log(x)*ctelog);
k = (x <= 10.^(m - 1));
m(k) = m(k) - 1;
f = 10.^(option.digit - m);
x = fix(x.*f)./f;
%b) Find sequential identical values
ind2 = find((x(2 : Nrho) - x(1 : Nrho - 1)) == 0) + 1;
if ( ~isempty(ind2) )
ind3 = ind2((tp_X(2 , ind2) < tp_X(n , ind2)));
tp_X(: , ind3) = tp_X(ind_n , ind3);
ind1 = [ind2(1) - 1 , ind2((ind2(2 : end) - ind2(1 : end - 1) ~= 1)) + 1];
ind4 = ind1(((tp_X(2 , ind1) ~= tp_X(2 , ind1 + 1))));
tp_X(: , ind4) = tp_X(ind_n , ind4);
end
Pt = (1 - option.alpha)*update_matrix(tp_X) + option.alpha*Pt ;
if (S(1) < S_opt)
X_opt = X(: , I(1));
S_opt = S(1);
end
Gt_old = Gt_new;
X(: , ind_betaR1) = tp_X(: , ind_betaR1);
t = t + 1;
%%=========== Display ===========%%
if (option.verbose)
subplot(121);
p1 = plot(option.V(X_opt , 1) , option.V(X_opt , 2) , 'r' , option.V(: , 1) , option.V( : , 2) , '+' , 'markersize' , 8 , 'linewidth' , 2);
text(option.V(: , 1) , option.V(: , 2) , str_n , 'fontname' , 'times' , 'fontsize' , 13);
axis([mini(1) , maxi(1) , mini(2) , maxi(2)])
title(sprintf('\\bf{n = %d, N = %d, \\gamma_{%d}^{opt}=%6.2f , \\gamma_{%d} = %6.2f, cons = %d}' , n , option.N , t , S_opt , t , Gt_new , cons) , 'fontsize' , 12 , 'fontname' , 'times');
subplot(122);
bar3(Pt);
axis([0 , n+1 , 0 , n+1 , 0 , 1])
view([-50 50]);
title(sprintf('\\bf{P_{%d}}' , t ) , 'fontsize' , 12 , 'fontname' , 'times');
drawnow
end
end
T = t - 1;
