ce_knapsack(p , W , c , option)

function out = ce_knapsack(p , W , c , option)

% Multiple knapsack problem solves by cross-entropy optimization
%
%  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)
%  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 7*10e-3)
%    option.alpha         Smoothing coefficient (0<alpha<1) for updating Bernouilli parameter
%    option.d             Indicate how many CE iterations to wait until declaring a solution is found (default 10)
%    option.theta_ini     Initial parameter values (default [])
%    option.T_max         Maximum number of iterations (default 50000)
%    option.verbose       Display iteration if = 1
%
%
% Outputs
% -------
%
% out
%    X_opt                Optimal solution. a (n + 1 x 1) index vector.
%    S_opt                Cost associated with the optimal solution
%    theta_opt            Optimal Bernouilli parameter
%    St                   Optimal cost versus time iteration of the algorithm
%
%
% Example
% -------
%
% load weish01.mat
% out = ce_knapsack(p , W , c);
% disp(sprintf('\nS_CE = %5.1f, S_opt = %5.1f' , out.S_opt , S_opt))
%
% Author : Sbastien PARIS ([e-mail])
% --------
%
%
% 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          = 500;
    option.rho        = 7*10e-3;
    option.alpha      = 0.6;
    option.d          = 10;
    option.T_max      = 50000;
    option.theta_ini  = [];
    option.verbose    = 1;
    

end

if(~any(strcmp(fieldnames(option) , 'N')))

    option.N                     = 1000;

end

if(~any(strcmp(fieldnames(option) , 'rho')))

    option.rho                   = 7*10e-3;

end

if(~any(strcmp(fieldnames(option) , 'alpha')))

    option.alpha                  = 0.6;

end

if(~any(strcmp(fieldnames(option) , 'd')))

    option.d                      = 10;

end

if(~any(strcmp(fieldnames(option) , 'tjeta_ini')))

    option.theta_ini              = [];

end

if(~any(strcmp(fieldnames(option) , 'T_max')))

    option.T_max                = 50000;

end

if(~any(strcmp(fieldnames(option) , 'verbose')))

    option.verbose               = 1;

end


if(~isempty(option.theta_ini))

    theta     = option.theta_ini;

else

    theta     = 0.5*ones(n , 1);

end

Nrho      = round(option.N*option.rho);

cteNrho   = 1/Nrho;

St        = zeros(1 , option.T_max);

ind_Nrho  = (1 : Nrho);

ON        = ones(1 , option.N);

S_opt     = -inf;

S_old     = -inf;

t         = 1;

cons      = 0;

alpha     = (1 + sum(p))./(max(max(c(: , ones(1 , n)) - W)));

while((t < option.T_max) && (cons < option.d) )

    %%%%%%%% t = 1, T ~ Ber(theta) %%%%%%

    T               = theta(: , ON);

    X               = binornd(T);

    %%% Descending Sorting %%

    [S , index]     = sort(cost_knapsack(X , p , W , c , alpha) , 2 , 'descend');

    X               = X(: , index);

    %%% Best current Score %%

    S_new           = S(1);

    if (S_new == S_old);

        cons         = cons + 1;
        
    else
        
        cons         = 0;
        
    end

    %%% Update pdf parameter theta with the best samples %%

    theta    = (1 - option.alpha)*sum(X(: , ind_Nrho) , 2)*cteNrho + option.alpha*theta;

    if(S_new > S_opt)
        
        X_opt      = X(: , 1);

        S_opt      = S_new;

    end

    St(t)      = S_opt;

    S_old      = S_new;

    t          = t + 1;

    if(option.verbose == 1)
        
        disp(sprintf('S_{opt} = %6.2f, S(%d) = %6.2f, d = %d', S_opt , t , S(1)  , cons))
        
    end

    drawnow

end


St(t:option.T_max)      = [];

out.S_opt               = S_opt;
out.X_opt               = X_opt;
out.theta_opt           = theta;
out.St                  = St;