function [PP,Pw,eval]=GSEMO(w,A,max_eval,refs)
% w: weight
% A: adjacency matrix
% max_eval: evaluation budget
% ref: reference points (targets)
% PP: final population
% Pw: penalized weights of PP
% eval: number of evaluation after termination
[k,m]=size(w);n=int32(size(A,1));
w_max=max(max(w));
max_size=ceil(max(sum(w,2)));
% Prepare graph for rank calculation
Ind=int32(zeros(m,2));ind=1;
for i=1:(n-1)
    for j=(i+1):n
        if A(i,j)
            Ind(ind,:)=[i,j];
            ind=ind+1;
        end
    end
end
P=false(max_size,m);
x=randi([0,1],1,m)>0;
P(1,:)=x;
eval=0;
Y=zeros(max_size,k);
y=double(m*w_max*(n-rank_graphic_in(x)-1))+(w*(x'))';
Y(1,:)=y;
eval=eval+1;
% If reference points are supplied, set up termination flag
if nargin>=4;l=size(refs,1);satflag=false(l,1);if check(y);max_eval=eval;end;end
pop_size=1;mutate_strength=1/m;
while eval<max_eval
    x=xor(P(randi([1,pop_size]),:),rand(1,m)<=mutate_strength);% Selection
    y=double(m*w_max*(n-rank_graphic_in(x)-1))+(w*(x'))';% Mutation
    eval=eval+1;
    t=repmat(y,pop_size,1)-Y(1:pop_size,:);
    s=all(t<=0,2);ss=sum(s);
    if ss>0% If dominating some solutions in P
        pop_size=pop_size-ss+1;P(1:pop_size,:)=[P(~s,:);x];Y(1:pop_size,:)=[Y(~s,:);y];
        if nargin>=4;if check(y);break;end;end
    else if ~any(all(t>=0,2))% If not dominated by any solution in P
            pop_size=pop_size+1;P(pop_size,:)=x;Y(pop_size,:)=y;
            if nargin>=4;if check(y);break;end;end
        end
    end
end
PP=P(1:pop_size,:);Pw=Y(1:pop_size,:);
[Pw,I]=sortrows(Pw);
PP=PP(I,:);
function b=check(f_val)% Set flag to check if each ref point is dominated
satflag=satflag|all(repmat(f_val,l,1)<=refs,2);
b=all(satflag);
end
function [r,is_independent]=rank_graphic_in(X)
ll=size(X,1);
r=zeros(ll,1);is_independent=false(ll,1);
for ii=1:ll
    YY=false(n,n);
    for jj=1:m
        YY(Ind(jj,1),Ind(jj,2))=X(ii,jj);
        YY(Ind(jj,2),Ind(jj,1))=X(ii,jj);
    end
    r(ii)=n-count_cc(YY);
    is_independent(ii)=(r(ii)==sum(X(ii,:)));
end

function cc=count_cc(edge_ind)
ver=false(1,n);cc=0;
for J=1:n
    if ver(J);continue;end
    ver(J)=true;cc=cc+1;
    new_ver=false(1,n);new_ver(J)=true;
    while sum(new_ver)>0
        new_ver=any(edge_ind(new_ver,:),1);
        new_ver(ver)=false;
        ver(new_ver)=true;
    end
end
end
end
end