function out = cemcmc_tsp(L , option)
% Solve the classic Travel Sale Man problem with Botev-Kroeze (BK) optimization
% algorithm
%
% Usage
% ------
%
% out = cemcmc_tsp(L , [option])
%
% Inputs
% -------
%
% L Distance/cost matrix (n x n)
% option Structure for the option of the BK algorithm (optional)
%
% option.N Number of paths to generate for each iteration (default 10*round(log10(n))*n)
% option.rho Threshold for selectionning the best rho*N paths (default log10(n)*0.25)
% option.b Number of gibbsampler burnin (default = round(10*log10(n)) )
% option.winjec Number of consecutive step to wait where gamma(t) not improved before injecting (default 5)
% option.nbinjec Number of injection (default 5*log10(n))
% option.Kgamma K Gain to add on gamma percentile (gamma(1+nb_inject*Kgamma) (default 1/100)
% option.Kb K Gain to add on b (b(1+nb_inject*Kb) (default 1/10)
% option.Krho K Gain to add on rho (rho(1+nb_inject*Krho) (default 1/50)
% option.KN K Gain to add on N (N(1+nb_inject*KN) (default 1/10)
% option.rhomax Maximum value of rho by injection (default 0.7)
% option.T_max Maximum number of iterations (default 10000)
% option.X_ini Initial vector (n + 1 x 1) (default = [])
% option.gamma_ini Initial percentile (default [])
% option.verbose 1 for real time plot, 2 for text display (default 2)
% option.V Cities 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
% St Optimal distance versus time iteration of the algorithm
% gammat rho-Percentile versus time
%
% Example
% ------
%
% n = 20;
% V = 10*rand(n , 2);
% L = distmat(V);
% out = cemcmc_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 (sebastien.paris@lsis.org)
% --------
%
% Ref : "An efficient Algorithm for Rare-Event Probability Estimation, Combinatorial, Optimization, and Counting",
% --- Z. I. Botev, D. P. Kroese, Methodology and Computing in Applied Probability, vol. 10
%
% 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(option))
option.N = 10*round(log10(n))*n; % number of tours to generate by MCWR
option.b = round(10*log10(n))*n;
option.rho = log10(n)*0.25; %0.5; % threshold for selectionning the best rho*N paths
option.winjec = 5;
option.nbinjec = round(5*log10(n));
option.Kgamma = 1/100;
option.Kb = 1/10;
option.Krho = 1/50;
option.rhomax = 0.7;
option.KN = 1/10;
option.X_ini = [];
option.gamma_ini = [];
option.T_max = 1000; % Maximum number of iteration
option.verbose = 2;
end
if(~any(strcmp(fieldnames(option) , 'N')))
option.N = 10*round(log10(n))*n;
end
if(~any(strcmp(fieldnames(option) , 'b')))
option.b = round(10*log10(n))*n;
end
if(~any(strcmp(fieldnames(option) , 'rho')))
option.rho = log10(n)*0.25; %0.5;
end
if(~any(strcmp(fieldnames(option) , 'winjec')))
option.winjec = 5;
end
if(~any(strcmp(fieldnames(option) , 'nbinjec')))
option.nbinjec = round(5*log10(n));
end
if(~any(strcmp(fieldnames(option) , 'Kgamma')))
option.Kgamma = 1/100;
end
if(~any(strcmp(fieldnames(option) , 'Kb')))
option.Kb = 1/10;
end
if(~any(strcmp(fieldnames(option) , 'Kb')))
option.Kb = 1/10;
end
if(~any(strcmp(fieldnames(option) , 'rhomax')))
option.rhomax = 0.7;
end
if(~any(strcmp(fieldnames(option) , 'KN')))
option.KN = 1/10;
end
if(~any(strcmp(fieldnames(option) , 'X_ini')))
option.X_ini = [];
end
if(~any(strcmp(fieldnames(option) , 'gamma_ini')))
option.gamma_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');
set(gca , 'nextplot' , 'replace');
str_n = num2str((1 : n)' , '%0.0f');
end
b = option.b;
rho = option.rho;
N = option.N;
Nrho = round(rho*N);
ind_Nrho = (1 : Nrho);
ON = ones(1 , option.N);
St = zeros(1 , option.T_max);
gammat = zeros(1 , option.T_max);
S_old = inf;
%%===== Initialization ========%%
if(~isempty(option.X_ini))
X = option.X_ini(: , ON);
else
[ignore,p] = sort(rand(n - 1 , N) , 1);
X = [ON ; p + 1 ; ON];
end
[S , indice] = sort(cost_TSP(L , X));
if(isempty(option.gamma_ini))
gammat_old = S(Nrho);
else
gammat_old = option.gamma_ini;
end
S_old = +inf;
I = indice(ind_Nrho);
X_opt = X(: , I(1));
S_opt = S(1);
t = 1;
cons = 0;
injec = 0;
if (option.verbose == 1)
maxi = max(option.V);
mini = min(option.V);
fig1 = figure(1);
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' , 12)
axis([mini(1) , maxi(1) , mini(2) , maxi(2)])
% title(sprintf('\\bf{n = %d, N = %d, \\rho=%3.2f, \\gamma_{%d}^{opt}=%6.2f , \\gamma_{%d}(%d) = %6.2f, d = %d}' , n , option.N , option.rho , t , S_opt , t , Nrho, gammat_old , cons) , 'fontsize' , 12 , 'fontname' , 'times');
title(sprintf('\\bf{n=%d, N=%d, \\rho=%3.2f, \\gamma_{%d}^{opt}=%6.2f , \\gamma_{%d}(%d)=%6.2f, b=%d, d=%d, i=%d}' , n , N , rho , t , S_opt , t , Nrho, gammat(t) , b, cons, injec) , 'fontsize' , 12 , 'fontname' , 'times');
elseif (option.verbose == 2)
disp(sprintf('n=%d, N=%d, \\rho=%3.2f, \\gamma_{%d}^{opt}=%6.2f , \\gamma_{%d}(%d)=%6.2f, b=%d, d=%d, i=%d' , n , N , rho , t , S_opt , t , Nrho, gammat(t) , b, cons, injec));
end
drawnow
%%========== Main Loop ==========%%
while( (t < option.T_max) && (injec <= option.nbinjec)) %(cons < option.d)
%%%%%%%% Best Elite Samples Xtilde ~ g*(x|gamma_{t-1}) %%%%%%%%%%
Xtilde = X(: , I);
Stilde = S(ind_Nrho);
%%%%%%%%%%%%%%%% Bootstrapping Xtilde => Y ~ g*(x|gamma_{t-1}) %%%%%%%%%%%%%%
ind_bootstrap = ceil(Nrho*rand(1 , N));
%ind_bootstrap = bootstrap(Stilde , option.N );
Y = Xtilde(: , ind_bootstrap);
SY = Stilde(ind_bootstrap);
%%%%%%%%%%%%%%%% Gibbsampler Y = > X %%%%%%%%%%%%%%
[X , SX] = gibbsampler_tsp(Y , SY , L , gammat_old , b);
%%% Index of Elite Samples %%
[S , indice] = sort(SX);
I = indice(ind_Nrho);
if(S(Nrho) < gammat_old)
gammat(t) = S(Nrho);
elseif(cons > option.winjec)
injec = injec + 1;
gammat(t) = round(gammat_old*(1 + injec*option.Kgamma));
b = round(option.b*(1 + injec*option.Kb));
rho = min(option.rhomax , option.rho*(1 + injec*option.Krho));
Nold = N;
N = round(option.N*(1 + injec*option.KN));
ind = ceil(Nold*rand(1 , (N - option.N))) ;
Nrho = round(rho*N);
ind_Nrho = (1 : Nrho);
indice = [indice , indice(ind)];
I = indice(ind_Nrho);
else
gammat(t) = gammat_old;
end
%%% Best current Score %%
S_new = S(1);
if ((gammat(t) == gammat_old) && (S_new == S_old));
cons = cons + 1;
else
cons = 0;
end
if(S_new < S_opt)
X_opt = X(: , I(1));
S_opt = S_new;
cons = 0;
end
St(t) = S_opt;
gammat_old = gammat(t);
S_old = S_new;
%%=========== Display ===========%%
if (option.verbose == 1)
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, \\rho=%3.2f, \\gamma_{%d}^{opt}=%6.2f , \\gamma_{%d}(%d)=%6.2f, b=%d, d=%d, i=%d}' , n , N , rho , t , S_opt , t , Nrho, gammat(t) , b, cons, injec) , 'fontsize' , 12 , 'fontname' , 'times');
elseif(option.verbose == 2)
disp(sprintf('n=%d, N=%d, \\rho=%3.2f, \\gamma_{%d}^{opt}=%6.2f , \\gamma_{%d}(%d)=%6.2f, b=%d, d=%d, i=%d' , n , N , rho , t , S_opt , t , Nrho, gammat(t) , b, cons, injec));
end
drawnow
t = t + 1;
end
St(t:option.T_max) = [];
gammat(t:option.T_max) = [];
out.S_opt = S_opt;
out.X_opt = X_opt;
out.St = St;
out.gammat = gammat;
