function [PP,Y,SS,SSw,eval]=MOEAD_in(w,A,L,max_eval,N,refs)
% w: weight
% A: adjacency matrix
% L: Lambda
% max_eval: evaluation budget
% N: neighborhood size
% ref: reference points (targets)
% PP: final population
% Y: penalized weights of PP
% SS: final output
% SSw: penalized weights of SS
% eval: number of evaluation after termination
l=size(L,1);
[k,m]=size(w);
n=int32(size(A,1));
% 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
w_max=max(max(w));size_max=ceil(max(sum(w,2))/l);
if nargin<5
    N=1;% Default neighborhood size
else
    N=min(max(N,1),l);
end
scalar_w=L*w;
eval=0;
neighbor=cell(1,l);
P=cell(1,l);P_size=int32(ones(1,l));
fit=zeros(l,1);
weights=cell(1,l);
S=false(size_max*l,m);Sw=zeros(size_max*l,k);
for i=1:l% Initialize population + set up neighbors
    temp=sqrt(sum((L-repmat(L(i,:),l,1)).^2,2));
    [~,temp]=sort(temp);
    neighbor{i}=temp(1:N);
    x=randi([0,1],1,m)>0;
    Px=false(size_max,m);Px(1,:)=x;
    P{i}=Px;
    S(i,:)=x;
    penalty=m*w_max*(n-rank_graphic_in(x)-1);
    fit(i)=penalty+scalar_w(i,:)*(x');
    t_w=(w*(x'))';
    weights{i}=zeros(size_max,k);weights{i}(1,:)=t_w;
    Sw(i,:)=double(penalty)+t_w;
    eval=eval+1;
end
temp=true(l,1);
for i=1:l% Check dominance
    if temp(i)
        t=repmat(Sw(i,:),l,1)-Sw(1:l,:);
        s=all(t<=0,2);
        temp(s)=false;
        temp(i)=true;
    end
end
S_size=sum(temp);S(1:S_size,:)=S(temp,:);Sw(1:S_size,:)=Sw(temp,:);
if nargin>=6% If reference points are supplied, set up termination flag
    lref=size(refs,1);
    satflag=false(lref,1);
    for i=1:S_size;check(Sw(i,:));end
    if all(satflag);max_eval=eval;end
end
mutate_strength=1/m;
while eval+l<=max_eval
    for j=1:l
        x=P{j}(randi([1,P_size(j)]),:);% Selection
        y=xor(x,rand(1,m)<=mutate_strength);% Mutation
        z=m*w_max*(n-rank_graphic_in(y)-1);zw=w*(y');zy=(zw+double(z))';
        eval=eval+1;
        neigh=neighbor{j};
        for h=1:numel(neigh)% Loop through neighborhood
            neigh_ind=neigh(h);
            zz=z+L(neigh_ind,:)*zw;% Penalized scalar fitness
            if zz<fit(neigh_ind)
                fit(neigh_ind)=zz;
                P{neigh_ind}(1,:)=y;
                weights{neigh_ind}(1,:)=zy;
                P_size(neigh_ind)=1;
            elseif zz==fit(neigh_ind)
                P_size(neigh_ind)=P_size(neigh_ind)+1;
                P{neigh_ind}(P_size(neigh_ind),:)=y;
                weights{neigh_ind}(P_size(neigh_ind),:)=zy;
            end
        end
        t=repmat(zy,S_size,1)-Sw(1:S_size,:);s=all(t<=0,2);ss=sum(s);
        if ss>0% If dominating some solutions in S
            S_size=S_size-ss+1;S(1:S_size,:)=[S(~s,:);y];Sw(1:S_size,:)=[Sw(~s,:);zy];
            if nargin>=6;check(zy);if all(satflag);max_eval=eval;break;end;end
        else if ~any(all(t>=0,2))% If not dominated by any solution in S
                S_size=S_size+1;S(S_size,:)=y;Sw(S_size,:)=zy;
                if nargin>=6;check(zy);if all(satflag);max_eval=eval;break;end;end
            end
        end
    end
    for j=1:l% Remove solutions with duplicate images
        [temp,ia]=unique(weights{j}(1:P_size(j),:),'rows');
        P_size(j)=size(temp,1);
        P{j}(1:P_size(j),:)=P{j}(ia,:);
        weights{j}(1:P_size(j),:)=temp;
    end
end
P_final_size=[0,cumsum(P_size)];
PP=false(P_final_size(end),m);Y=zeros(P_final_size(end),k);
for i=1:l
    PP((P_final_size(i)+1):P_final_size(i+1),:)=P{i}(1:P_size(i),:);
    Y((P_final_size(i)+1):P_final_size(i+1),:)=weights{i}(1:P_size(i),:);
end
[~,I]=sortrows(Y);
PP=PP(I,:);Y=Y(I,:);
SS=S(1:S_size,:);SSw=Sw(1:S_size,:);
[~,I]=sortrows(SSw);
SS=SS(I,:);SSw=SSw(I,:);

function check(f_val)% Set flag to check if each ref point is dominated
satflag=satflag|all(repmat(f_val,lref,1)<=refs,2);
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