%BUMDA=Boltzmann Univariate Marginal Distribution Algorithm
%S. Ivvan Valdez P.

%%% Please refer this code using the web page:
%%% http://personal.cimat.mx:8181/~ivvan/public/bumda.html
%%% Article:
%%% http://www.sciencedirect.com/science/article/pii/S0020025513001576
%%% A Boltzmann based estimation of distribution algorithm
%%% S. Ivvan Valdez, , Arturo Hernández , Salvador Botello 
%%% Information Sciences
%%% Volume 236, 1 July 2013, Pages 126–137
%%% And the web page:
%%% www.cimat.mx/~ivvan/bumda.m



%It is implemented for a maximization case.
%Arguments:
% npop    <- Population Size
% nvar    <- Number of decision variables
% xinf    <- Inferior limit for decision variables
% xsup    <- Superior limit for decision variables
% ngen    <- Maximum number of generations, stopping criterion
% problem <- matlab function, it is used via feval for computing the objective function values for the population
% minsd   <- most of the cases this is the best stopping criterion it is the standard deviation of the objective function values of the population


%Variables:
% fobj    <- Population, objective function 
% X       <- Population, decision variables
% ix      <- Selected set, indexes


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%BUMDA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xelite,felite,gen,F,X]=bumda(npop,nvar,xinf,xsup,ngen,problem,minsd)
%Settings
if (nargin<1)
npop=2000;                  %Population Size
end
if (nargin<2)
nvar=10;                   %Number of variables 
end
if (nargin<3)
xinf=-20.0*ones(nvar,1);  %Inferior search limits
end
if (nargin<4)
xsup=20.0*ones(nvar,1);   %Superior search limits
end
if (nargin<5)
ngen=200;                 %Number of generations
end
if (nargin<6)
problem=@xphere;          %Test problem
end
if (nargin<7)
minsd=1.0E-8;           %Stopping Criterion
end
gen=1;

theta=-inf;
X=initialize(xinf,xsup,npop,nvar);       %Generating the initial population
F=evaluate(problem,X);                   %Evaluate the objectie function
[ix,theta,ibest]=truncate(F,gen,theta);  %Truncating the population
xelite=X(ibest,:);                       %Copying the elite individual
felite=F(ibest,:);                       %Copying the elite individual
for gen=2:ngen
    G=translate_obj_func(F);                     %The objective function is traslated to get positive values
    [means,sds]=bumda_means_and_sds(G,X, ix,nvar); %Computing the means and standard deviations
    X=[xelite;regenerate(means,sds,npop)];       %Repopulating and preserving the elite individual
    F=evaluate(problem,X);                       %Evaluation
    [ix,theta,ibest]=truncate(F,gen,theta);      %Selection
    xelite=X(ibest,:);                           %Getting the elite
    felite=F(ibest,:);                           %Getting the elite
    if (std(F)<minsd) 
       break;
    end
end
end



%Test function: sphere objective function
%Note: the negative sign is for maximization
function [obj]=xphere(x)
nvar=length(x);
disp=0.1*(1:nvar); %An small displacement to check if everything is ok
obj=-sum((x-disp).^2);
end


%Generates the initial population
function [X]=initialize(xinf,xsup,npop,nvar)
   X=[];
   for i=1:nvar
      X=[X,unifrnd(xinf(i), xsup(i), npop, 1)];
   end
end

%Evaluates the objective function for the whole population,
%requires the function f and the population decision variables
function [F]=evaluate(f,X)
npop=size(X,1);
F=zeros(npop,1);
 for i=1:npop
    F(i)=feval(f,X(i,:));
 end
end

%Truncation method, returns the indexes of the selected set,
% the treshold "theta", and the elite individual index
function [ix,theta,ibest]=truncate(F,gen,theta)
npop=size(F,1);
[Ft,ix]=sort(F,'descend');
if (gen==1)
   theta=Ft(npop); 
else
   if(F(npop/2)>=theta)
      ix=ix(1:(npop/2));
      theta=F(ix(npop/2));
   else
      for i=((npop/2+1)-(1:npop/2))
         if (F(ix(i))>=theta)
            ix=ix(1:i);
            theta=F(ix(i));
            break;
         end
      end 
   end
end
ibest=ix(1);
end

%Computes the translated objective function, 
%in order to only have positive values.
function G=translate_obj_func(F)
mf=min(F);
G=F-mf+1.0e-6;
end

%Computes the mean and standad deviation per variable using the BUMDA formula
function [means,sds]=bumda_means_and_sds(G,X, ix,nvar)
means=1:nvar;
sds=means;
nsel=size(ix,1);
sumG=sum(G(ix));
   for i=1:nvar
       means(i)=sum(G(ix).*X(ix,i))/sumG;
       sds(i)=sqrt(sum(G(ix).*((X(ix,i)-means(i)).^2))/(1.0e-6+sumG));
   end
end

%Regenerate the population of size npop-1 considering the insertion of the elite 
%individual
function [X]=regenerate(means,sds,npop)
nvar=length(means);
X=[];
 for i=1:nvar
    X=[X,normrnd(means(i), sds(i), npop-1, 1)]; 
 end
end


