function out = cemcmc_knapsack(p , W , c , option)
%
%
% Multiple knapsack problem solves by Botev-Kroeze (BK) optimization
% algorithm
%
% max S(x)=(p^{t}x)
% st. Wx <= c
%
% equivalent to maximize
%
% max S(x)=C(x) + p^{t}x, where C(x)=alpha sum_{i=1,...,m}(min(ci -
% Wi^{t}xj , 0) and alpha = 1 + 1^{t}p/(max_{i,j}(ci-wij)
%
% Usage
% ------
%
% out = cemcmc_knapsack(p , W , c , [option])
%
% Inputs
% ------
%
% p Weights vector (n x 1)
% W Matrix Weights (m x n)
% c Costs vector (m x 1)
% option Structure for the option of the BK algorithm (optional)
%
% option.N Number of Bernouilli vectors to generate for each iteration (default 1000)
% option.rho Threshold for selectionning the best rho*N samples (default 0.1)
% option.b Number of gibbsampler burnin (default = 2*n)
% option.d Indicate how many BK iterations to wait until declaring a solution is found (default 10)
% option.T_max Maximum number of iterations (default 10000)
% option.epsilon Absolute relative difference between 2 consecutives solution
%
% Outputs
% -------
%
% out
% X_opt Optimal solution. a (n + 1 x 1) index vector.
% S_opt Cost associated with the optimal solution
% St Optimal cost versus time iteration of the algorithm
% gammat rho-Percentile versus time
%
%
%
% Example
% -------
%
% load weish01.mat
% out = cemcmc_knapsack(p , W , c);
% disp(sprintf('\nS_CEMCMC = %5.1f, S_opt = %5.1f' , out.S_opt , S_opt))
%
% 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
[m , n] = size(W);
if ( (nargin < 4) || isempty(option))
option.N = 1000;
option.rho = 0.5;
option.b = 10*n;
option.winjec = 5;
option.nbinjec = round(4*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 = 500;
option.verbose = 1;
end
if(~any(strcmp(fieldnames(option) , 'N')))
option.N = 10*round(log10(n))*n;
end
if(~any(strcmp(fieldnames(option) , 'b')))
option.b = round(5*log10(n))*n;
end
if(~any(strcmp(fieldnames(option) , 'rho')))
option.rho = log10(n)*0.2; %0.5;
end
if(~any(strcmp(fieldnames(option) , 'winjec')))
option.winjec = 5;
end
if(~any(strcmp(fieldnames(option) , 'nbinjec')))
option.nbinjec = round(4*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 = 1;
end
b = option.b;
rho = option.rho;
N = option.N;
Nrho = round(N*rho);
St = zeros(1 , option.T_max);
gammat = zeros(1 , option.T_max);
ind_Nrho = (1 : Nrho);
alpha = (1 + sum(p))./(max(max(c(: , ones(1 , n)) - W)));
%%%%%%%%%%%%% Initial run t = 0 %%%%%%%%%%%%%
if(~isempty(option.X_ini))
X = option.X_ini(: , ON);
else
X = binornd(0.5 , n , option.N);
end
%%% Descending Sorting %%
[S , index] = sort(cost_knapsack(X , p , W , c , alpha) , 2 , 'descend');
if(isempty(option.gamma_ini))
gammat_old = S(Nrho);
else
gammat_old = option.gamma_ini;
end
%%% Index of Elite Samples %%
I = index(ind_Nrho);
S_opt = -inf;
S_old = S_opt;
t = 1;
cons = 0;
injec = 0;
while((t < option.T_max) && (injec < option.nbinjec) )
%%%%%%%% 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));
Y = Xtilde(: , ind_bootstrap);
SY = Stilde(ind_bootstrap);
%%%%%%%%%%%%%%%% Gibbsampler Y = > X %%%%%%%%%%%%%%
[X , SX] = gibbsampler_knapsack(Y , SY , p , W , c , alpha , gammat_old , b);
[S , index] = sort(SX , 2 , 'descend');
% /// Matlab code of the gibbsampler_knapsack ///
%
% bino_mat = binornd(0.5 , option.b , option.N);
%
% indexgibbs = ceil(n*rand(1 , option.b)); % randomize update order of condition pdf's of the gibbsampler
%
% for i = 1:option.b
%
% ii = indexgibbs(i);
%
% Y(ii , :) = bino_mat(ii , :);
%
% Sd = cost_knapsack(Y , p , W , c , alpha);
%
% ind = (Sd < gammat_old);
%
% Y(ii , ind) = ~Y(ii , ind);
%
% SY(~ind) = Sd(~ind);
%
% end
%
% X = Y;
%
% [S , index] = sort(SY , 2 , 'descend');
%%% Index of Elite Samples %%
I = index(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);
index = [index , index(ind)];
I = index(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;
end
St(t) = S_opt;
S_old = S_new;
gammat_old = gammat(t);
if (option.verbose == 1)
disp(sprintf('N=%d, rho=%3.2f, S_{opt} = %6.2f, S(%d) = %6.2f, gamma(%d) = %6.2f, b=%d, d=%d, i=%d' , N , rho, S_opt, t, S(1) , t, gammat(t), b, cons, injec))
end
t = t + 1;
drawnow
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;
