cemcmc_satA(A , option)

function out = cemcmc_satA(A , option)
%
%
%  SAT solver by a Botev-Kroeze (BK) optimization
%  algorithm
%
%   S(X)=sum(j=1,...,m)C_j, where C_j=max(0 , (2X_i-1)*Aji
%
% 
%
% Usage
% ------
%
% out = cemcmc_satA(C , [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 0.5)
%    option.b             Number of gibbsampler burnin (default = 10*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)
%
%
% 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              = cemcmc_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        = 0.5;
    option.b          = 10*n;
    option.d          = 10;
    option.T_max      = 500;

end


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

St             = zeros(1 , option.T_max);

gammat         = zeros(1 , option.T_max);


ind_Nrho       = (1 : Nrho);

S_opt          = -inf;

S_old          = S_opt;

t              = 1;

cons           = 0;


%%%%%%%%%%%%% Initial run t = 0 %%%%%%%%%%%%%

X              = binornd(0.5 , n , option.N);

%%% Descending Sorting %%

[S , index]    = sort(cost_satA(X , A) , 2 , 'descend');

gammat_old     = S(Nrho);

%%% Index of Elite Samples %%

I              = index(ind_Nrho);
    
while((t < option.T_max) && (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 , option.N));
    
    Y               = Xtilde(: , ind_bootstrap);
    
    SY              = Stilde(ind_bootstrap);

     %%%%%%%%%%%%%%%% Gibbsampler Y = > X %%%%%%%%%%%%%%
     
    [X , SX]        = gibbsampler_satA(Y , SY , A , gammat_old , option.b); 
    
    [S , index]     = sort(SX , 2 , 'descend');
    

    %%% Index of Elite Samples %%

    I               = index(ind_Nrho);
       
    if(S(Nrho) > gammat_old)

        gammat(t)   = S(Nrho);

    else

        gammat(t)   = gammat_old;

    end
   
    %%% Best current Score %%

    S_new           = S(1);

    if (S_new == S_old );

        cons         = cons + 1;
             
    else 
        
        cons         = 0;

    end


    if(S_new > S_opt)
        
        S_opt      = S_new;

    end

    St(t)      = S_opt;

    S_old      = S_new;
    
    gammat_old = gammat(t);
    
    disp(sprintf('S_{opt} = %6.2f, S(%d) = %6.2f, gamma(%d) = %6.2f, d = %d' , S_opt, t, S(1) , t, gammat(t), cons))
    
    t          = t + 1;
    
    drawnow

end


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

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