```% elin    Convex roof of the linear entropy (linear entropy of
%         entanglement)
%         See https://arxiv.org/abs/1409.3806
%         For an alternative realization see
%         http://de.mathworks.com/matlabcentral/fileexchange/47823-corona-convex-roof-numerical-analysis?requestedDomain=www.mathworks.com

function E_lin=elin(rho_full)

% TEST 1=ON,0=OFF
TEST=0;

% Treshold to take an eigenvalue zero
lambdatreshold=0.000001;

[sy,sx]=size(rho_full);

% Dimension of qudits
dq=floor(sqrt(sy)+0.5);

% Number of two-qudit units
% for the PPT symmetric extension
% N=2, first nontrivial level
N=2;

% Calculate the matrix for dimension reduction
% and the density matrix in the new basis
if 1 %rank(rho_full)<dq^2
[v,dd]=eig(rho_full);
dd=diag(dd);
% Tereshold to look for zero eigenvalue
index=find(abs(dd)>lambdatreshold);

Ureduce=[];
for n=1:length(index)
Ureduce=[Ureduce,v(:,index(n))];
end %for

% d=rank of the state
d=length(index);
rho=Ureduce'*rho_full*Ureduce;

% Just in case, normalize rho
rho=diag(diag(rho));
rho=real(rho);
rho=nm(rho);

%Ur=runitary(1,d);
%Ureduce=Ureduce*Ur;

%rho=diag(dd(index));
% Note: rho_full=Ureduce*rho*Ureduce'
else
d=dq^2;
rho=rho_full;
Ureduce=eye(d);
end %if

% Just to be sure, normalize it
rho=nm(rho);

% CHOOSE SOLVER
%yalmip('solver','sdpt3')
%yalmip('solver','sedumi')
sdpsettings('verbose',1);

% Observable in the expectation value
Obs=zeros(dq^4,dq^4);
oo=sud(dq);
for n=1:dq^2-1
Obs=Obs+mkron(oo(:,:,n),eye(dq),oo(:,:,n),eye(dq));
end %for

% ONLY FOR TEST
% Compute sum_k G_k^2
%sGk2=0;
%for n=1:dq^2-1
%    sGk2=sGk2+oo(:,:,n)^2;
%end %for
%trace_sGk2=trace(sGk2)

% Reduce to the range of rho
Obs0=Obs;
Obs=kron(Ureduce,Ureduce)'*Obs*kron(Ureduce,Ureduce);
Obs=(Obs+Obs')/2;

if d==2
rhols=sdpvar(N+1,N+1,'hermitian','complex');
%rhol=0*sdpvar(d^N,d^N,'hermitian','complex');
rhol=zeros(d^N,d^N);
for n=0:N
for m=0:N
rhol=rhol+rhols(n+1,m+1)*dstate(n,N)*dstate(m,N)';
end %for
end %for
elseif (d==3 || d==4 || d==9 || d==6) || N==2

Ps=proj_sym(N,d);
rPs=rank(Ps);
rhols=sdpvar(rPs,rPs,'hermitian','complex');
[v,dd]=eig(Ps);
index=find(diag(dd)>0.99);
v=v(:,index);

for n=1:rPs
for m=1:rPs
if TEST
disp(['General state-symmetric state ' num2str(n) ' ' num2str(m)]);
end %if
if n==1 && m==1
rhol=rhols(1,1)*v(:,1)*v(:,1)'; %Intialization
else
rhol=rhol+rhols(n,m)*v(:,n)*v(:,m)';
end %for

end %for
end %for

else

error(['Error with d and N: the procedure is not yet implemented for these paramaters: d=' num2str(d) '.'])

end %if

if TEST
disp(['1: general state-symmetric state conversion is done'])
end %if

% Computing the reduced states
[sx,sy]=size(rhol);
N=log2(sx)/log2(d);
N=floor(N+0.5);
listneg=1:2;
Nr=length(listneg);
if N==2
rho_twoqudit=rhol;
else
rho_twoqudit=sdpvar(d,d,'hermitian','complex');
% TEST
%[sy,sx]=size(rhol);
%R=randn(sx,sy);
%rho2=R'*rhol*R;
%disp(['11'])
for k=0:d^Nr-1
for l=0:d^Nr-1
if TEST
disp(['Two-body ' num2str(k) ' ' num2str(l)]);
end %if
rr=0*rhol(1,1);
i1=1+k*d^(N-Nr);
i2=1+l*d^(N-Nr);
for n=0:d^(N-Nr)-1
rr=rr+rhol(i1+n,i2+n);
end %for
rho_twoqudit(k+1,l+1)=rr;
end %for
end %for
end %if

if TEST
disp(['2: Computing two-body reduced state is done.'])
end %if

if N==2
rho_singlequdit=[];
for n=1:d
rho_singlequdit_row=[];
for m=1:d
rho_singlequdit_row=[rho_singlequdit_row,trace(rhol((n-1)*d+1:(n-1)*d+d,(m-1)*d+1:(m-1)*d+d))];
end %for
rho_singlequdit=[rho_singlequdit;rho_singlequdit_row];
end %for
else
listneg=1;
Nr=length(listneg);
rho_singlequdit=sdpvar(d,d,'hermitian','complex')
for k=0:d^Nr-1
for l=0:d^Nr-1
if TEST
disp(['Single-particle reduced state ' num2str(k) ' ' num2str(l)]);
end %if
rr=0*rhol(1,1);
i1=1+k*d^(N-Nr);
i2=1+l*d^(N-Nr);
for n=0:d^(N-Nr)-1
rr=rr+rhol(i1+n,i2+n);
end %for
rho_singlequdit(k+1,l+1)=rr;
end %for
end %for
end %if

if TEST
disp(['3: Computing single-particle reduced state is done.'])
end %if

if N==2
rhol_PT=[];
for n=1:d
rhol_PT_row=[];
for m=1:d
rhol_PT_row=[rhol_PT_row,rhol((n-1)*d+1:(n-1)*d+d,(m-1)*d+1:(m-1)*d+d).'];
end %for
rhol_PT=[rhol_PT;rhol_PT_row];
end %for
else
% Partial transpose
rhol_PT=rhol; % NOT VERY NICE, BUT WORKS.
if N==2
N1=1;
N2=1;
elseif N==3
N1=2;
N2=1;
elseif N==4
N1=1;
N2=3;
end %if
AM=kron(1:d^N,ones(d^N,1));
BM=kron(ones(1,d^N),(1:d^N).');
AT=pt_nonorm(AM,1:N1,d);
BT=pt_nonorm(BM,1:N1,d);
for m=1:d^N
for n=1:d^N
if TEST
disp(['Partial transpose ' num2str(m) ' ' num2str(n)]);
end %if
rhol_PT(m,n)=rhol(AT(m,n),BT(m,n));
end %for
end %for
end %if

if TEST
disp(['4: Comptuing the partial transpose is done'])
end %if

% Of course, we cannot juts write F=F+set(rho_singlequdit==rho);
F=set(rhols>=0)+set(trace(rhols)==1);
for n=1:d
for m=1:n
if TEST
disp(['Constraints ' num2str(n) ' ' num2str(m)]);
end %if
F=F+set(rho_singlequdit(m,n)==rho(m,n));
end %for
end %for

F=F+set(rhol_PT>=0);

if TEST
disp(['5: Contstaints are done'])
end %if

diagnostic=solvesdp(F,-real(trace(Obs*rho_twoqudit)));

% if numerics had problems, give back NaN
if diagnostic.problem~=0
E_Lin=NaN;
end %if

supremum=double(trace(Obs*rho_twoqudit));

E_lin=1-1/dq-(1/2)*supremum;

% Only for test
if TEST

diagnostic=diagnostic

supremum=supremum

supremum2_should_be_equal_to_supremum=double(trace(Obs0*kron(Ureduce,Ureduce)*rho_twoqudit*kron(Ureduce,Ureduce)'))

E_lin=E_lin

r2=double(rho_twoqudit);

r2=kron(Ureduce,Ureduce)*r2*kron(Ureduce,Ureduce)';

% A single copy of the bipartite state
r1=keep(r2,2:-1:1,dq);

difference_should_be_zero=max(max(abs(r1-rho_full)))

neg_rho=negativity(rho_full,1,dq)

end %if

if E_lin<-1e-5
error('Error: Negative enanglement!')
end %if

```