function best = GACE (dim, max_generations)

%% This is the algorithm presented to the Black Box Competition in GECCO 2015.
% This is an improved version of the method described in: 
% P. Lopez-Garcia, E. Onieva, E. Osaba, A. Masegosa, A. Perallos. Hybridizing Genetic Algorithm with Cross Entropy for Solving Continuous Functions. 
% Proceedings of the Companion Publication of the Genetic and Evolutionary Computation Conference, pp. 763-764, 2015.
%% ACM Digital Library URL: http://dl.acm.org/citation.cfm?id=2764881
%% Contact: p.lopez@deusto.es

%% Input parameters
%  dim = dimension used to optimize the function. 
%  max_generations = Number of maximum generations.

%% Output parameters
%  best = best individual found by the algorithm in max_generations.

%% Parameter initialization
    %% Default input options
    if nargin == 0
        dim = 10; 
        max_generations = 2500;
    end
    
	p_ga = 0.1; %GA size percentage
    p_up = 0.4; %Individuals taken to update CE part
    ind_size = dim; % Individual size
    
    %Get pop_size and calculate other population sizes
    
    if dim <= 10
        c = 5;
    else
        c = 2;
    end
    
    pop_size = dim*c;
    pop_ga = max(3, round(pop_size*p_ga));
    pop_ce = max(3, pop_size - pop_ga);
    n_up = round(pop_ce*p_up);
    
    %% GACE operators
    Selection       = @BinaryTournament;
    Crossover       = @BLX_05Crossover;
    Mutation        = @GaussianMutation_Variable;
    
    %% Common operators
    Initialization  = @InitializationDouble;
    Repair          = @RepairDouble;
    FitnessFunction = @generic; %Change to your fitness function
    
    %% GA parameters
     p_c = 0.9; %Crossover probability
     p_m = 1/ind_size; %Mutation probability
    
    %% CE parameters
    num_update      = max(3,n_up); %Number of individuals taken to update CE part
    learning_rate   = 0.7; 
    
    %% Common parameters (Problem-Dependent parameters) 
    ind_min  = -1.*ones(1,ind_size);
    ind_max  = 1.*ones(1,ind_size);
    
    
    %% Initialization
    Population = Initialization(pop_ga+pop_ce,ind_size,ind_min,ind_max);
    Means = Initialization(1,ind_size,ind_min,ind_max);
    Desvs = Initialization(1,ind_size,ind_min,ind_max);
    Desvs.x = (ind_max-ind_min)./2;
    Means = Desvs;
    best = Population(1);
	generation = 0;
	
	
    while generation < max_generations
    generation = generation+1;
        %% Repair and Evaluate Population
    Population = Repair(Population,ind_min,ind_max);
    
    for i = 1: pop_size 
            Population(i).f = FitnessFunction(Population(i).x);
    end
    
    %% Check the best individual     
    %  Elitism: Replace the last individual with the best, if the best does not exist in the population
    fitness = vertcat(Population.f);
    [a,b]=min(fitness);
    
    
    if a<best.f
        best = Population(b);
        gwc = 0;
    elseif a~=best.f    
        Population(end)=best;
    end
        
    if pop_ga > 0
        %% Selection
        Parents = Selection(pop_ga,Population);

        %% Crossover and Mutation
        Children_GA = Crossover(Parents,p_c);
        Children_GA = Mutation(Children_GA,p_m,ind_min,ind_max,0.005);
    end
    
    
    %% CE part
    if pop_ce > 0
        %% Select <num_update> best samples
        faux = vertcat(Population.f);
        [~,b] = sort(faux);    
        Selected = Population(b(1:num_update));
        
        %% Update <Means> and <Desvs>
        Means.x = learning_rate*Means.x + (1-learning_rate)*mean(vertcat(Selected.x));
        Desvs.x = learning_rate*Desvs.x + (1-learning_rate)*std(vertcat(Selected.x));
        

        %% Generate Samples
        Children_CE = struct([]);
        for i=1:pop_ce
            Children_CE(i).x = normrnd(Means.x,Desvs.x);                      
            Children_CE(i).f = 0; 
            % Resample out of range values (avoid round to min-max)
            Children_CE(i).x(Children_CE(i).x>ind_max) = normrnd(Means.x(Children_CE(i).x>ind_max),Desvs.x(Children_CE(i).x>ind_max)); 
            Children_CE(i).x(Children_CE(i).x<ind_min) = normrnd(Means.x(Children_CE(i).x<ind_min),Desvs.x(Children_CE(i).x<ind_min));    
        end
     
    end
    
    %% Replacement
    if pop_ce == 0
        Population = Children_GA;
    elseif pop_ga == 0
        Population = Children_CE;
    else
        Population = [Children_GA Children_CE];
    end
    
    
    end
	
	
	function parents = BinaryTournament(n_parents,Population)
    %%
		parents = [];
		pop_size=numel(Population);
		ffs = vertcat(Population.f);
		for padre=1:n_parents
			is = 1+round(rand(1,2)*(pop_size-1));
			[~,b]=min(ffs(is));
			parents = [parents Population(is(b))];
		end
	
    function off = BLX_05Crossover(par,p_c)
		off = par;
		dim = numel(par(1));
		ALPHA=0.5;
		for i=1:2:length(par)-1
			if rand()<p_c
				i1=i;i2=i+1;
				AMP = abs(off(i1).x-off(i2).x);
				MIN = min(off(i1).x,off(i2).x);
				off(i1).x = MIN-AMP*ALPHA + rand(1,dim)*(1+2*ALPHA)*AMP;
				off(i2).x = MIN-AMP*ALPHA + rand(1,dim)*(1+2*ALPHA)*AMP;
			end
		end
    function  off = GaussianMutation_Variable(par,p_m,ind_min,ind_max, value)
		%% Gaussian mutation per Variable
		off = par;    
		dim = numel(par(1).x);
		for i=1:length(off)
			for j=1:dim
				if rand()<p_m
					off(i).x(j) = normrnd(off(i).x(j), value*(ind_max(j)-ind_min(j)));  
				end
			end
		end
		
    function Population = InitializationDouble(pop_size,ind_size,ind_min,ind_max)
		Population = struct([]);    
		for i=1:pop_size
			Population(i).x = rand(1,ind_size);
			for j=1:ind_size
				Population(i).x(j)=ind_min(j) + Population(i).x(j)*(ind_max(j)-ind_min(j));
			end
			Population(i).f = inf;
		end
	
    function Population = RepairDouble(Population,ind_min,ind_max)
    
		for i=1:numel(Population)
			Population(i).x=min(Population(i).x,ind_max);
			Population(i).x=max(Population(i).x,ind_min);
		end
    
    function fitness = generic(x)
		fitness = sum(x.^2);