cma_es.m

strfitnessfct = 'frastrigin10';
N = 10;
xmean = rand(N,1);
init_xmean = xmean;
sigma = 0.5;
stopfitness = 1e-10;
stopeval = 1e3*N^2;

lambda = 4+floor(300*log(N));
mu = lambda/2;
weights = log(mu+1/2)-log(1:mu)';
mu = floor(mu);
weights = weights/sum(weights);
mueff=sum(weights)^2/sum(weights.^2);

cc = (4 + mueff/N) / (N+4 + 2*mueff/N);
cs = (mueff+2) / (N+mueff+5);
c1 = 2 / ((N+1.3)^2+mueff);
cmu = min(1-c1, 2 * (mueff-2+1/mueff) / ((N+2)^2+mueff));
damps = 1 + 2*max(0, sqrt((mueff-1)/(N+1))-1) + cs;

pc = zeros(N,1); ps = zeros(N,1);
B = eye(N,N);
D = ones(N,1);
C = B * diag(D.^2) * B';
invsqrtC = B * diag(D.^-1) * B';
eigeneval = 0;
chiN=N^0.5*(1-1/(4*N)+1/(21*N^2));

out.dat = []; out.datx = [];

function f=frosenbrock(x)
  if size(x,1) < 2 error('dimension must be greater one'); end
  f = 100*sum((x(1:end-1).^2 - x(2:end)).^2) + sum((x(1:end-1)-1).^2);
endfunction

function f=frastrigin10(x)
  N = size(x,1); if N < 2 error('dimension must be greater one'); end
  scale=10.^((0:N-1)'/(N-1));
  f = 10*size(x,1) + sum((scale.*x).^2 - 10*cos(2*pi*(scale.*x)));
endfunction

function f=fschwefel(x)
  f = 0;
  for i = 1:size(x,1),
    f = f+sum(x(1:i))^2;
  end
endfunction

function f=fgriewank(x)
  f = 1 + sum(x.^2)/4000 - prod(cos(x./sqrt(1:size(x,1))'));
endfunction

counteval = 0;
while counteval < stopeval

  for k=1:lambda,
    arx(:,k) = xmean + sigma * B * (D .* randn(N,1));
    arfitness(k) = feval(strfitnessfct, arx(:,k));
    counteval = counteval+1;
  end

  [arfitness, arindex] = sort(arfitness);
  xold = xmean;
  xmean = arx(:,arindex(1:mu)) * weights;

  ps = (1-cs) * ps ...
        + sqrt(cs*(2-cs)*mueff) * invsqrtC * (xmean-xold) / sigma;
  hsig = sum(ps.^2)/(1-(1-cs)^(2*counteval/lambda))/N < 2 + 4/(N+1);
  pc = (1-cc) * pc ...
        + hsig * sqrt(cc*(2-cc)*mueff) * (xmean-xold) / sigma;

  artmp = (1/sigma) * (arx(:,arindex(1:mu)) - repmat(xold,1,mu));

  C = (1-c1-cmu) * C ...
       + c1 * (pc * pc' ...
               + (1-hsig) * cc*(2-cc) * C) ...
       + cmu * artmp * diag(weights) * artmp';

  sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1));

  if counteval - eigeneval > lambda/(c1+cmu)/N/10
    eigeneval = counteval;
    % C = triu(C) + triu(C,1)';
    [B,D] = eig(C);
    D = sqrt(diag(D));
    invsqrtC = B * diag(D.^-1) * B';
  end

  if arfitness(1) <= stopfitness || max(D) > 1e7 * min(D)
    break;
  end

  % more off;  % turn pagination off in Octave
  % disp([num2str(counteval) ': ' num2str(arfitness(1)) ' ' ...
  %       num2str(sigvalma*sqrt(max(diag(C)))) ' ' ...
  %       num2str(max(D) / min(D))]);
  % init_xmean'
  % xmean'

  out.dat = [out.dat; arfitness(1) sigma 1e5*D' ];
  out.datx = [out.datx; xmean'];
end

disp([num2str(counteval) ': ' num2str(arfitness(1))]);
xmin = arx(:, arindex(1));  % Return best point of last iteration.
                            % Notice that xmean is expected to be even
                            % better.
figure(1); hold off; semilogy(abs(out.dat)); hold on;  % abs for negative fitness
semilogy(out.dat(:,1) - min(out.dat(:,1)), 'k-');  % difference to best ever fitness, zero is not displayed
title('fitness, sigma, sqrt(eigenvalues)'); grid on; xlabel('iteration');
figure(2); hold off; plot(out.datx);
title('Distribution Mean'); grid on; xlabel('iteration')