close all; clear all;

J1 = @(x) 0.9 - (x.^2 - 0.1*cos(2*pi*8*x));
a = -1; b = 1; pts = 10000; stp = (b-a)/pts; x = a:stp:(b-stp);
% plot(x,J1(x)); hold on;

rand('state',0); randn('state',0);
NB = 10; mu = 20; R = 100; pm = 0.001; epsilon = 1e-6; b = mu;
G = 20; % Para a inicialização em 0, usar G = 150 e R = 200; para a inicialização aleatória, usar G = 20 e R = 100;
E = 2.^((NB-1):-1:0);
BKP = zeros(G,R);

for nr = 1:R,
for g = 1:G,
BKPTotal(g,nr).Dados = zeros(NB+1,mu);
end;
end;

for nr = 1:R,
% P = zeros(NB,mu); pm = 0.001; % plot((E*P-512)/512,J1((E*P-512)/512),'r.'); pause;
P = rand(NB,mu)>0.5; pm = 0.05; % plot((E*P-512)/512,J1((E*P-512)/512),'r.'); pause;
% pause;
Jmax = -1e12;
for g = 1:G,
% Custo
J = (J1((E*P-(2^(NB-1)))/2^(NB-1)))' + epsilon*randn(mu,1);
BKPTotal(g,nr).Dados = [P ; J'];
% [P ; J']
% pause;
% Controle do Melhor Resultado
if max(J)>Jmax, Jmax = max(J); Xmax = P(:,find(J==max(J))); end; BKP(g,nr) = Jmax;
% Selecao de Pais SUS
p = sortrows([(1:mu)' (J-min(J))/sum(J-min(J))],2); a = cumsum(p(:,2));
% a
MP2 = zeros(b,1);
r0 = rand(1,1)/mu; r = r0; k = 1; i = 1;
while k<b,
while r<=a(i),
MP2(k) = p(i,1);
r = r+1/mu;
k = k+1;
end;
i = i+1;
end;
MP2(find(MP2==0)) = p(mu,1);
% Recombinacao
PR = zeros(NB,b);
for k = 1:b/2,
a = 1 + floor(rand(1,1)*(NB-1));
PR(1:a,k) = P(1:a,MP2(k)); PR(a+1:end,k) = P(a+1:end,MP2(k+b/2));
PR(1:a,k+b/2) = P(1:a,MP2(k+b/2)); PR(a+1:end,k+b/2) = P(a+1:end,MP2(k));
end;
P = [P(:,MP2) PR];
% Mutacao
P = 0.5-2*(P-0.5).*((rand(size(P))<pm)-0.5);
% Selecao de Sobreviventes Ranking
J = (J1((E*P-(2^(NB-1)))/2^(NB-1)))';
p = sortrows([(1:length(J))' J],2);
P = P(:,p((end-mu+1):end,1));
% Jmax
% J
% P
% g
% plot(x,J1(x)); hold on; plot((E*P-512)/512,J1((E*P-512)/512),'g.'); hold off;
% pause;
end;
[nr Jmax]
end;

for nr=1:R,
for g=1:G,
P = BKPTotal(g,nr).Dados(1:NB,:);
J = BKPTotal(g,nr).Dados(NB+1,:);
if not(isempty(J(P(1,:) == 0))),
S0(g,nr) = mean(J(P(1,:) == 0));
end;
if not(isempty(J(P(1,:) == 1))),
S1(g,nr) = mean(J(P(1,:) == 1));
end;
end;
end;

for g=1:G,
MS0(g) = mean(S0(g,(S0(g,:)~=0)));
MS1(g) = mean(S1(g,(S1(g,:)~=0)));
end;

plot(MS0); hold on; plot(MS1,'r-');