function out = ce_tsp(L , option)
% Solve the classic Travel Sale Man problem with a Cross Entropy
% (CE) optimization algorithm.
%
% Usage :
%
% out = ce_tsp(L , [option])
%
% Inputs
% -------
%
% L : Distance/cost matrix (n x n)
% 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 solutions
% option.P_ini : Initial Probabilities Transition Matrix (n x n)
% option.digit : Truncation digit (default 5)
% option.verbose : 1 for real time plot, 2 for text display (default 2);
% option.V : Towns coordinates (n x 2). option.verbose = 1 is requiered;
%
% Outputs
%
% out
% X_opt Optimal solution. a (n + 1 x 1) index vector.
% S_opt Distance associated with the optimal path
% Pt_opt Transition Matrix (n x n) estimated at the end of the CE algorithm
%
% Example
%
% n = 20;
% V = 10*rand(n , 2);
% L = distmat(V);
% out = ce_tsp(L);
% plot(V(: , 1) , V(: , 2) , 'b+' , V(out.X_opt , 1) , V(out.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);
X1 = 1;
if((n ~= d))
error('L must be a square(n x n)');
end
if ( (nargin < 2) || 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.P_ini = [];
option.verbose = 2; % Display Best current solution of the CE algorithm
end
if(~any(strcmp(fieldnames(option) , 'N')))
option.N = round(5*n*n);
end
if(~any(strcmp(fieldnames(option) , 'rho')))
option.rho = 5*10e-3;
end
if(~any(strcmp(fieldnames(option) , 'alpha')))
option.alpha = 0.7;
end
if(~any(strcmp(fieldnames(option) , 'beta')))
option.beta = 0.7;
end
if(~any(strcmp(fieldnames(option) , 'd')))
option.d = 5;
end
if(~any(strcmp(fieldnames(option) , 'epsilon')))
option.epsilon = 100*eps;
end
if(~any(strcmp(fieldnames(option) , 'P_ini')))
option.P_ini = [];
end
if(~any(strcmp(fieldnames(option) , 'T_max')))
option.T_max = 10000;
end
if(~any(strcmp(fieldnames(option) , 'verbose')))
option.verbose = 2;
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
St = zeros(1 , option.T_max);
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 =====%%
if(~isempty(option.P_ini))
Pt = option.P_ini;
else
cten1 = 1/(n - 1);
Pt = cten1(ones(n));
Pt(1 : n + 1 : n*(n + 1)) = 0;
end
P0 = Pt;
%%===== 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 - option.gamma)*update_matrix(tp_X) + option.gamma*P0 + 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 == 1)
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)])
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 - option.gamma)*update_matrix(tp_X) + option.gamma*P0 + option.alpha*Pt ;
if (S(1) < S_opt)
X_opt = X(: , I(1));
S_opt = S(1);
end
Gt_old = Gt_new;
St(t) = S_opt;
X(: , ind_betaR1) = tp_X(: , ind_betaR1);
t = t + 1;
%%=========== Display ===========%%
if (option.verbose == 1)
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');
elseif(option.verbose == 2)
disp(sprintf('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));
end
drawnow
end
St(t:option.T_max) = [];
%gammat(t:option.T_max) = [];
out.S_opt = S_opt;
out.X_opt = X_opt;
out.Pt_opt = Pt;
out.St = St;
%out.gammat = gammat;
