function [xopt, fopt, eval_count] = dxnes(f, dim, pop_size, initMu, initSigma, maxGens, stopfx)
%[xopt, fopt] = dxnes(@f, dim, pop, initMu, initSigma)
%INPUTS:
% @f: objective function to be minimized
% dim: dimension of the objective function
% pop_size: number of the individuals generated in one iteration
% initMu: initial mean vector
% initSigma: initial stepsize
% **maxGens: maximum generation
% **stopfx: the minimum objective function value
%OUTPUTS:
% xopt: best individuals
% fopt: best objective function value
% eval_count: number of the evaluations
%
% @N. Fukushima.

    if (nargin < 7)
        stopfx = 10^(-10);
        if (nargin < 6) 
            maxGens = 5 .* 10.^5 .* dim ./ pop_size;
        end
    end
    % -- Algorithm parameters ---------------------------------
    % dimension and population size
    n = dim; 
    L = 2 * floor(pop_size / 2);
    
    % learning rates
    etaMu = 1.0;               
    etaSigmaset = [1.0 0.5*(1.0 + L ./ (L + 2.*n)) 1.0 + L ./ (L + 2.*n)];
    etaBset = [(L + 2.*n) ./ (L + 2.*n.^2 + 100) .* min(1, sqrt(L/n)) L ./ (L + 2.*n.^2 + 100)];

    % utility function
    wi = max(0, log(L/2 + 1) - log(1:1:L));
    uRank = wi ./ sum(wi) - 1 ./ L;
    alpha = 0.9 + 0.15 * log(n);
    
    % evolution path
    base = sqrt(n) * (1 - 1/(4*n) + 1/(21*n^2));
    mueff = 1./sum((uRank + 1./L).^2);
    csigma = (mueff + 2.0) ./ (n + mueff + 5.0) / sqrt(n);

    % -- Algorithm starts -------------------------------------
    % initialize
    mu = initMu; sigma = initSigma; I = eye(n, n); B = I;
    psigma = repmat(0, n, 1);
    xopt = mu; fopt = f(mu);
    g = 0;
    
    condition = true;
    while condition
        g = g + 1;
        A = sigma * B;
        
        % sampling the individuals with the antithetic variates method
        z = randn(n, L/2); z = cat(2, z, -z);
        x = A * z + repmat(mu, 1, L);
        
        % evaluate and sort the individuals
        fx = f(x); [fx idxfx] = sort(fx, 2, 'ascend');
        z = z(:,idxfx);
        
        if fx(1) <= fopt
            xopt = x(:, idxfx(1));
            fopt = fx(1);
        end
        
        % update the evolution path
        sumUZ = sum(repmat(uRank, n, 1) .* z, 2);
        psigma = (1 - csigma) .* psigma + sqrt( csigma .* (2 - csigma) .* mueff ) * sumUZ;
        rate = norm(psigma) / base;
        
        % calculate the gradients with the detection of the center's movement
        % and adapt the learning rates by the identification of the search phases
        GM = repmat(0, n, n);
        if rate >= 1.0
            expZ = exp(alpha * sqrt(sum(z.^2, 1)));
            uDist = (wi .* expZ) ./ sum(wi .* expZ, 2) - 1 ./ L;
            Gdelta = sum(repmat(uDist, n, 1) .* z, 2);
            GM = (repmat(uDist, n, 1) .* z) * z' - sum(uDist) * I;
            etaSigma = etaSigmaset(1); etaB = etaBset(1);     
        elseif rate < 1.0        
            Gdelta = sumUZ;
            GM = (repmat(uRank, n, 1).*z) * z' - sum(uRank) * I;
            etaSigma = etaSigmaset(2); etaB = etaBset(2);
            if rate < 0.1
                etaSigma = etaSigmaset(3);
            end 
        end
 
        % update the parameters
        mu = mu + etaMu * A * Gdelta;
        Gsigma = trace(GM) ./ n; GB = GM - Gsigma * I;        
        sigma = sigma * exp(etaSigma./2 .* Gsigma);
        B = B * expm(etaB./2 .* GB);
        
        if (g >= maxGens) || (fopt <= stopfx) || (sigma < 10e-15)
            condition = false;
        end
    end
    eval_count = g * L;
end