close all; clear all; % clc; pack;

% Parte I - Ilustração dos Dados

% a = -5; b = 5; pts = 10000; stp = (b-a)/pts; x = a:stp:(b-stp); y = x.^2;
% a1 = -5; b1 = 5; pts1 = 1000; stp1 = (b1-a1)/pts1; x1 = a1:stp1:(b1-stp1); 
% a2 = -5; b2 = 5; pts2 = 1000; stp2 = (b2-a2)/pts2; x2 = a2:stp2:(b2-stp2); 
%    498   314   % Soluções: (500,314) e (500,688) => (0 ; -1.87) e (0 ; +1.87)
%    504   314
%    498   688
%    504   688
% [x1(500) x2(314) x2(688)]

a1 = -5; b1 = 5; pts1 = 50; stp1 = (b1-a1)/pts1; x1 = a1:stp1:(b1-stp1); 
a2 = -5; b2 = 5; pts2 = 50; stp2 = (b2-a2)/pts2; x2 = a2:stp2:(b2-stp2); 
y = zeros(pts1,pts2);
phi = y;

for g=1:pts1,
    for h=1:pts2,
        y(g,h) = x1(g)^2 + x2(h)^2;
        % if abs(7*x1(g)^2 + 3*x2(h)^2 - 10) <= 0.5,
        if (7*x1(g)^2 + 3*x2(h)^2) < 10,
            phi(g,h) = 1;
        end;
    end;
end;

% imagesc(y); colormap('gray'); axis square;
% figure; imagesc(phi); colormap('gray'); axis square;
% figure; imagesc(y.*phi); colormap('gray'); axis square;
figure(1); mesh(x1,x2,y.*phi); % [u,v] = find(y.*phi == max(max(y.*phi)));
% Solução exata: x1 = 0 e x2 = \pm sqrt(10/3) = 1.825;

% Parte II - Algoritmo Genético Simples com Peso de Restrição Adaptativo

J1 = @(u,v) u.^2+v.^2;
J2 = @(w,z) 7*w.^2+3*z.^2;

mu = 20;
W = -10;
Nruns = 100;
Ngeracoes = 1000;
rand('state',0); randn('state',0);
R = 2;
Jmax = 0;
BKP = [];

for run = 1:Nruns,
    P = [R*(rand(2,mu)-0.5) ; W*ones(1,20)];
    % P = [R*(rand(2,mu)-0.5)]; % Use esta linha para penalidade estática
    ngeracoes = 0;
    while ngeracoes < Ngeracoes,
        J = (J1(P(1,:),P(2,:)) + P(3,:).*(J2(P(1,:),P(2,:))>10))';
        % J = (J1(P(1,:),P(2,:)) -30*(J2(P(1,:),P(2,:))>10))'; % Use esta linha para penalidade estática
        % Critério de Parada e Controle do Melhor Resultado
        ngeracoes = ngeracoes + 1;
        if max(J)>Jmax,
            Jmax = max(J);
            Xmax = P(:,find(J==max(J)));
        end;
        % Seleção de Pais SUS
        p = sortrows([(1:mu)' (J-min(J))/sum(J-min(J))],2); a = cumsum(p(:,2)); MP2 = zeros(mu,1);
        r0 = rand(1,1)/mu; r = r0; k = 1; i = 1;
        while k<mu,
            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);
        % Recombinação
        P = [P(:,MP2) 0.5*P(:,MP2(1:2:end-1))+P(:,MP2(2:2:end))];
        % Mutação
        P = P+0.05*randn(size(P));
        % Seleção de Sobreviventes
        J = (J1(P(1,:),P(2,:)) + P(3,:).*(J2(P(1,:),P(2,:))>10))';
        % J = (J1(P(1,:),P(2,:)) -30*(J2(P(1,:),P(2,:))>10))'; % Use esta linha para penalidade estática
        p = sortrows([(1:length(J))' J],2);
        P = P(:,p((end-mu+1):end,1));      
    end;
    [run Jmax Xmax']
    BKP = [BKP ; run Jmax Xmax'];
    % Para obter gráficos mostrando a população final e a solução obtida, descomente o trecho abaixo:
    figure(2); plot(P(1,:),P(2,:),'r.'); axis([-2 2 -2 2]); grid on;
    a = -1.19; b = 1.19; pts = 100; stp = (b-a)/pts; x = a:stp:b-stp; y = sqrt((10-7*x.^2)/3);
    hold on; plot(x,y,'b-'); plot(x,-y,'b-'); plot(Xmax(1),Xmax(2),'g.');
    pause;
    hold off;
end;

