Skip to content

Instantly share code, notes, and snippets.

@mitmul
Created April 15, 2013 10:30
Show Gist options
  • Save mitmul/5387198 to your computer and use it in GitHub Desktop.
Save mitmul/5387198 to your computer and use it in GitHub Desktop.
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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment