function MOO_CEM

%The workarea for the algorithm, including population of decision values,
%objective values and ranking
global WorkArea 
%Elite vector preserving good solutions and driving the CEM
global Elite;
%The means vector of the CEM; the standard deviation vector of the CEM
global mu; global sigma;
%Variable to collect each Elite mean and std dev after an iteration. 
%For info purposes, e.g. in thesis.
global SetMu; global SetSigma;

%The researcher can set this value to solve some standard MOO problems:
%1,2,3,4,6 in Coello Coello, and ZDT1, 2, 3.
MOP = 2;

%Problem parameters were coded in InitializeProblem to make the code
%more generic. "Limits" is used as a varying limit vector in the histogram
%assignments, while "L" is kept fixed. These initially have the same values
%(the limits of the decision variables).
[NumVars, NumObjectives, Limits, L, SheetName, ProblemN] = ...
    InitializeProblem(MOP);

%Clear all vectors:
WorkArea = []; mu=[]; Elite=[]; sigma=[]; SetMu=[]; SetSigma=[]; 

%The epsilon of the CEM. Set at 0.1-1.
eps = 0.1;

%Alpha of the CEM. Set at 0.7-0.8
alpha=0.7; %0.7

%Maximum number of outer loops. This loop supports the multi-objective
%search. In single-objective search, it would be one loop only.
NoOfLoops = 100;  

%Population size. Set at 30-100.
N = 100;  

%Probability to invert histograms to diversify search and to avoid trapping
Prob = 0.3;  
MaxEvaluations = 10000;  %Try to achieve good solutions with "few" tries.

%If the researcher wants to try different random numbers.
%rand('twister', 156089);
%rand('state', 34089);
%rand('seed', 14089);

NumEvaluations = 0;
Elite=[];
tic  %Start clock, not necessary.

% ~~~~~~~~~~  MAIN LOOP ~~~~~~~~~~  
for k=1:NoOfLoops 
    k
    %Init the CE vectors. sigmaeps keeps track of the consecutive sigma
    %differences for convergence checking.
    sigmaeps(1:NumVars)=Inf;
    %CEM requires an arbitrary large initial sigma
    sigma(1:NumVars) = 10*(L(1:NumVars,2) - L(1:NumVars,1));
    %Select a random mu vlaue for each DV in its definition range
    mu(1:NumVars) = (L(1:NumVars,1) + (L(1:NumVars,2) - ...
        L(1:NumVars,1)).*rand(NumVars,1))';

    t=0;
    WorkArea = []; %Must clear WorkArea in each main loop
 
    NotTerminate = true;

    while NotTerminate  %Will be governed by convergence
        t=t+1;
        bin_freq = []; 
        bin_edges = [];
        if size(Elite,1) > 0 && k>1 
           r = k + 2;
           for i=1:NumVars
               bin_edges(1:r+1) = 0;
               bin_edges(1) = L(i, 1);  
               bin_edges(r+1) = L(i, 2);  
               bin_edges(2) = min(Elite(:, i));
               bin_edges(r) = max(Elite(:, i));
               bin_edges(2:r)= bin_edges(2):(bin_edges(r) - ...
                   bin_edges(2))/(r-2):bin_edges(r); 
               bin_freq(1:r)   = histc(Elite(:,i), bin_edges(1:r));
               if rand<Prob
                  bin_freq(1:r) = max(bin_freq) - bin_freq(1:r);
               end
               bin_freq = floor(N*bin_freq./sum(bin_freq));
               s = sum(bin_freq(1:r));
               bin_freq(r) = N - sum(bin_freq(1:r)) + bin_freq(r);
               s=sum(bin_freq);
               UpTo=0;  %Indices into WorkArea
               %Now sample on each histogram bin range, with the number of
               %observations proportional to the histogram count
               for a=1:r 
                   Start = UpTo + 1; 
                   UpTo  = UpTo + bin_freq(a);
                   Limits(i, 1) = bin_edges(a); 
                   Limits(i, 2) = bin_edges(a+1);
                   if Start <= UpTo
                      %This is not the same mu as mu(i). 
                      %It moves with the histogram range
                      temp_mu = Limits(i,1) + rand*(Limits(i,2) - Limits(i,1));  
                      temp_s  = Limits(i,2) - Limits(i,1);
                      WorkArea(Start:UpTo, i) = ...
                          rand(UpTo-Start+1, 1)*(normcdf(Limits(i,2), ...
                         temp_mu, temp_s) - normcdf(Limits(i,1), ...
                         temp_mu, temp_s)) + normcdf(Limits(i,1), ...
                         temp_mu, temp_s);
                     %WorkArea(Start:UpTo,i) = log2(WorkArea(Start:UpTo,i));
                      WorkArea(Start:UpTo, i) = ...
                      norminv(WorkArea(Start:UpTo, i),...
                          temp_mu, temp_s);   
                   end
               end   %For each bin of the histogram of this i-th DV
           end       %For each DV
        else         % Build initial population
           for i=1:NumVars
               %Select random values from the definition ranges of the DVs
               WorkArea(1:N, i) = rand(N, 1)*(normcdf(L(i, 2), ...
                   mu(i), sigma(i)) - normcdf(L(i,1), mu(i), ...
                                 sigma(i))) + normcdf(L(i,1),...
                                 mu(i), sigma(i));
               WorkArea(1:N, i) = norminv(WorkArea(1:N, i), ...
                   mu(i), sigma(i));     
           end
        end    %Elite exists or initial population must be formed
        
        s = size(WorkArea, 1);
        % +2 is for assistance with ranking 
        WorkArea(1:s, NumVars+NumObjectives+2) = zeros(s, 1);  
        %Call function "f1" to evaluate the DV values:
        WorkArea(1:s, NumVars+1) = f1(WorkArea, NumVars, MOP);
        %Call function "f2" to evaluate the DV values:
        WorkArea(1:s, NumVars+2) = f2(WorkArea, NumVars, MOP);

        %Rank the solutions and remove any "good" solutions. 
        %Here, the threshold is th=2 (see text)    
        Temp = Rank(WorkArea, 2, NumVars, NumObjectives);
        %Add to Elite:
        Elite = vertcat(Elite, Temp); 

        %It's nice to see the progress, but it can slow down the algorithm:
        PlotDetailProgress(NumVars, WorkArea, ProblemN, NoOfLoops, t, k, N);

        %If good solutions were identified, Elite can be exploited.
        if size(Elite,1) > 1 
            AllEps = 0;
            for i=1:NumVars
                %Adjust the mu's according to the CEM:
                mu(i)       = (1-alpha)*mu(i)     + alpha*mean(Elite(:,i));
                sigmaeps(i) = sigma(i);
                %Adjust the sigmas according to the CEM:
                sigma(i)    = (1-alpha)*sigma(i) + alpha*std(Elite(:,i));
                %Calculate the consecutive changes, to determine
                %convergence:
                sigmaeps(i) = abs(sigma(i) - sigmaeps(i));
                %If this i-th DV has converged, mark it:
                AllEps = AllEps + (sigmaeps(i) > eps);  
            end
            %If they all converged, terminate inner loop:
            if AllEps == 0
                NotTerminate = false;
            end
            SetMu = vertcat(SetMu, mu);  %Record the current mu
            SetSigma = vertcat(SetSigma, sigma);  %and sigma for information
        end
        %Avoid trapping of the algorithm: allow no more than half the
        %maximum number of evaluations for the inner loop:
        NotTerminate = (NotTerminate && (N*t <= MaxEvaluations/2)); 
        if ~NotTerminate, break, end
        
        NumEvaluations = NumEvaluations + N;
        if  (NumEvaluations >= MaxEvaluations), break, end
    end  %while not Terminate

    NumEvaluations;
    %Rank again with th=1:
    Elite=Rank(Elite, 1, NumVars, NumObjectives);
    %Plot the intermediate Elite vector (approximation set), for
    %information purposes only:
    i = subplot(2,2,2);
    hold on
    scatter(Elite(:,NumVars+1), Elite(:, NumVars+2), 3, '*');
    xlabel('f1')
    ylabel('f2')
    title(['Current Pareto front, after iteration ' int2str(t)...
        '  ' int2str(NumEvaluations)])
    drawnow
    grid on
    hold off
    
    %Terminate outer loop if NumEvaluations was reached before NoOfLoops
    %was executed:
    if (NumEvaluations >= MaxEvaluations),  break, end  
end % for k

%Final ranking, non-dominated values, th=0:
Elite=Rank(Elite, 0, NumVars, NumObjectives);
size(Elite)

toc
%Show the known Pareto front and the approximate front:
Plot_WorkArea(Elite(:,NumVars+1), Elite(:,NumVars+2),...
    MOP, SheetName, NumEvaluations)

end  % Main function MOO_CEM
%------------------------------------------------------------------------


function Do_f1 = f1(X, NVars, MOP)
    if MOP == 1
       Do_f1=X(:,1).^2;
    elseif MOP == 2
       rt = 1/sqrt(NVars); 
       Do_f1=1-exp(-((X(:,1)-rt).^2+(X(:,2)-rt).^2+(X(:,3)-rt).^2));
    elseif MOP == 3
       A1 = 0.5*sin(1)-2*cos(1)+  sin(2) - 1.5*cos(2);
       A2 = 1.5*sin(1)-  cos(1)+2*sin(2) - 0.5*cos(2);
       B1 = 0.5*sin(X(:,1)) - 2* cos(X(:,1)) + sin(X(:,2)) - 1.5*cos(X(:,2));
       B2 = 1.5*sin(X(:,1)) -    cos(X(:,1)) + 2*sin(X(:,2)) - 0.5*cos(X(:,2));
       Do_f1 = -(1 + (A1 - B1).^2 + (A2 - B2).^2);
    elseif MOP == 4
       Do_f1 = -10*(exp(-0.2*sqrt(X(:,1).^2 + X(:,2).^2)) +...
           exp(-0.2*sqrt(X(:,2).^2 + X(:,3).^2))); 
    elseif MOP == 6
       Do_f1=X(:,1); 
    elseif MOP >= 8
       Do_f1 = X(:,1);  %ZDT1, 2 & 3 
    end   
end
function Do_f2 = f2(X, NVars, MOP)
    if MOP == 1
       Do_f2=(X(:,1)-2).^2; 
    elseif MOP == 2
        rt = 1/sqrt(NVars);
        Do_f2=1-exp(-((X(:,1)+rt).^2+(X(:,2)+rt).^2+(X(:,3)+rt).^2)); 
    elseif MOP == 3
       Do_f2 = -((X(:,1) + 3).^2 + (X(:,2) + 1).^2);
    elseif MOP == 4
       Do_f2   = abs(X(:,1)).^(0.8)+ 5*sin((X(:,1)).^3) + abs(X(:,2)).^(0.8)...
           + 5*sin((X(:,2)).^3) + abs(X(:,3)).^(0.8) + 5*sin((X(:,3)).^3);
    elseif MOP == 6
        x=X(:,1)./(1+10*X(:,2));
        y=x;
        x=x.^2;
        x=1-x;
        Do_f2=(1+10*X(:,2)).*(x - y.*sin(12*pi*X(:,1)));
    elseif MOP >= 8 && MOP <= 10 %ZDT1-3
        c = 9/(NVars-1);  
        x = transpose(sum(transpose(X(:,2:NVars))));
        gx = 1 + x.*c;
        gx_inv = 1./gx;
        if MOP == 8     %ZDT1
            Do_f2 = gx.*(1 - sqrt(gx_inv.*X(:,NVars+1)));
        elseif MOP == 9 %ZDT2
            Do_f2 = gx.*(1 - (gx_inv.*X(:,NVars+1)).^2);
        elseif MOP == 10 %ZDT3
            Ten_Pi = 10*pi;
            Do_f2 = gx.*(1 - sqrt(gx_inv.*X(:,1)) - ...
                gx_inv.*X(:,1).*sin(Ten_Pi*X(:,1)));
        end
    end
end

function PPF =Plot_WorkArea(x,y, MOP, xlSheetName, NEval)
subplot(2,2,4);
hold on;
if MOP <= 4 || (MOP >= 6 && MOP < 11)
    z=[]; %The true Pareto fronts provided by Coello Coello are in Excel
    z = xlsread('True_PFs_Coello.xls', xlSheetName);
    scatter(z(:,1), z(:,2), 5, 'v', 'filled');
end 
scatter(x, y, 3, 'o');
hold off
grid on
xlabel('f1');
ylabel('f2');
title(['Final Pareto Front of MOP', int2str(MOP), ...
    ' after ', int2str(NEval), ' evaluations']);
end
function RankIt = Rank(Pop, Threshold, NVars, NObj)  

K = NVars+NObj;
Pop(:,K+1)=0;

F=[];
Sp=[];
signe = -1; %+1 for MOP3

for z=NVars+1:K - 1
    Pop=sortrows(Pop, signe*(z));
    for p=1:size(Pop,1)-1
        for q=p+1:size(Pop,1)
            if Pop(p,z+1) >  Pop(q, z+1)  %Turn around for MOP3, maxim.
               %Rank is in last col.
               Pop(p, K+1) = Pop(p, K+1) + 1;   
                %No need to look further, this candidate is not making it
               if Pop(p, K+1) > (NObj-1)*Threshold, break, end  
            end
        end
        if Pop(p, K+1) <= Threshold
           F = vertcat(F, Pop(p, :));
        end
    end    
    F = vertcat(F, Pop(size(Pop,1), :));   %Add the last vector to the set
end
RankIt  = F;

end    

function [NumVars, NumObjectives, Limits, L, SheetName, ProblemN] = ...
    InitializeProblem(MOP)

MOP_Config = [1 1 2 -1E5 1E5, %MOP1
              2 3 2 -4 4,     %MOP2
              3 2 2 -pi pi,   %MOP3
              4 3 2 -4 4,     %MOP4
              5 0 0 0 0,
              6 2 2 0 1,      %MOP6   
              7 0 0 0 0,
              8 30 2 0 1,     %ZDT1
              9 30 2 0 1,     %ZDT2
              10 30 2 0 1     %ZDT3
              ]; 
          
NumVars       = MOP_Config(MOP, 2);
NumObjectives = MOP_Config(MOP, 3);

for i=1:NumVars    %Set problem boundaries
    L(i,1) = MOP_Config(MOP, 4);  
    L(i,2) = MOP_Config(MOP, 5); 
end

ProblemName = ['MOP1 ', 'MOP2 ', 'MOP3 ', 'MOP4 ', 'MOP5 ', 'MOP6 ',...
    'MOP7 ', 'ZDT1 ', 'ZDT2 ', 'ZDT3 '];

for i=1:NumVars
    Limits(i,1) = L(i,1);
    Limits(i,2) = L(i,2);
end;    

ProblemN = ProblemName((MOP-1)*5+1:5*MOP);
SheetName = ProblemN(1:4);
end %function InitializeProblem

%PlotDetailProgress
function PlotDetailProgress(NumVars, WorkArea, ProblemN, NoOfLoops, t, ...
    k, N)
  h = subplot(2,2,[1 3]);
  hold on    
  scatter(WorkArea(1:N, NumVars+1), WorkArea(1:N, NumVars+2), 3, '*');
  xlabel('f1')
  ylabel('f2')
  title(['Search space for ' ProblemN ': k=' int2str(k)...
      ' of ' int2str(NoOfLoops) ', t=' int2str(t)]);
  drawnow
  grid on
  hold off
end



