ce_satA(A , option)

function out = ce_satA(A , option)

%  SAT solves by Cross-Entropy (CE) optimization
%  algorithm
%
%   S(X)=sum(j=1,...,m)C_j, where C_j=max(0 , (2X_i-1)*Aji
%
% Usage
% ------
%
% out = ce_satA(A , [option])
%
% 
%  Inputs
%  ------
%
%   A                     Sparse clauses matrix    (m x n)
%   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.theta_min     Minimum probability value (default 0.0)
%    option.theta_max     Mmaximul probability value (default 1.0)
%    option.T_max         Maximum number of iterations (default 50000)
%    option.verbose       Display iteration if = 1
%
%
% Outputs
% -------
%
% out
%    X                    Feasable solutions satisfying S_opt<=m clauses. a (n x FS(C)).
%    S_opt                Best clauses satisfied (S_opt <= m)
%    St                   Optimal cost versus time iteration of the algorithm
%    gammat               rho-Percentile versus time
%
%
% Example
% -------
%
% A                 = readsatA('uf20-01.cnf');
% out               = ce_satA(A);
% disp(sprintf('\n Pourcentage of Satisfied clauses = %5.2f, card(X) = %d' , 100*(out.S_opt/size(A , 1)) ,size(out.X , 2)))
%
% 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

n         = size(A , 2);

if ( (nargin < 2) || isempty(option))

    option.N          = 1000;
    option.rho        = 7*10e-3;
    option.alpha      = 0.6;
    option.d          = 10;
    option.theta_ini  = [];
    option.theta_min  = 0.0;
    option.theta_max  = 1.0;
    option.T_max      = 50000;    
    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) , 'theta_ini')))

    option.theta_ini              = [];

end

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

    option.theta_min              = 0.0;

end

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

    option.theta_max              = 1.0;

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;

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_satA(X , A) , 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    = min(option.theta_max , max(option.theta_min,(1 - option.alpha)*sum(X(: , ind_Nrho) , 2)*cteNrho + option.alpha*theta));
    
    if(S_new > S_opt)
        
        S_opt      = S_new;
        
    end
    
    if(option.verbose)

        disp(sprintf('S_{opt} = %6.2f, S(%d) = %6.2f, d = %d'  , S_opt, t  , S_new  , cons))

    end
    
    drawnow

    St(t)      = S_opt;

    S_old      = S_new;
    
    t          = t + 1;


end


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

out.S_opt               = S_opt;
out.X                   = unique(X(: , S == S_opt)' , 'rows')';
out.theta_opt           = theta;
out.St                  = St;