RandSelfadjMatrix[l_]:=Module[{A,i},
A=RandomVariate[NormalDistribution[],{l,l}]+I*RandomVariate[NormalDistribution[],{l,l}];
A=A+ConjugateTranspose[A];
A=Chop[A];
A]

RandPositive[l_,r_]:=Module[{A},
A=RandomVariate[NormalDistribution[],{l,r}]+I*RandomVariate[NormalDistribution[],{l,r}];
A=A.ConjugateTranspose[A];
A]

RandState[l_,r_]:=Module[{A},
A=RandPositive[l,r];
A=A/Tr[A];
A]

CostfromObservable[Observables_,transpose_]:=Sum[MatrixPower[KroneckerProduct[Observables[[k]],IdentityMatrix[Length[Observables[[1]]]]]-KroneckerProduct[IdentityMatrix[Length[Observables[[1]]]],If[transpose,Transpose[Observables[[k]]],Observables[[k]]]],2],{k,1,Length[Observables]}]

SaMCNB[dim_]:=SaMCNB[dim]=Flatten[{{Table[SparseArray[{{k,k}->1,{dim,dim}->0}],{k,1,dim}]},Table[Table[SparseArray[{{k,m}->Sqrt[2]/2,{m,k}->Sqrt[2]/2,{dim,dim}->0}],{m,k+1,dim}],{k,1,dim-1}],Table[Table[SparseArray[{{k,m}->-I*Sqrt[2]/2,{m,k}->I*Sqrt[2]/2,{dim,dim}->0}],{m,k+1,dim}],{k,1,dim-1}]},2];

SaMPB[2]={{{1,0},{0,1}},{{0,1},{1,0}},{{0,-I},{I,0}},{{1,0},{0,-1}}};
SaMPB[dim_]:=SaMPB[dim]=Flatten[Table[KroneckerProduct[SaMPB[dim/2][[j]],SaMPB[2][[k]]],{j,1,dim^2/4},{k,1,4}],1];

SaMPNB[dim_]:=SaMPNB[dim]=SaMPB[dim]/Sqrt[dim];

QOT[ρ_,ω_,C_,dual_,transpose_]:=Module[{dim, sol,Π,x,y,X,Y},
dim=Length[ρ];
If[
dual,{x=Table[Symbol["x"<>ToString[n]],{n,dim^2}];
y=Table[Symbol["y"<>ToString[n]],{n,dim^2}];
sol=SemidefiniteOptimization[-Tr[(y.SaMCNB[dim]).ω+(x.SaMCNB[dim]).ρ],{VectorGreaterEqual[{C-KroneckerProduct[(y.SaMCNB[dim]),IdentityMatrix[dim]]-KroneckerProduct[IdentityMatrix[dim],If[transpose,Transpose[(x.SaMCNB[dim])],x.SaMCNB[dim]]],0},{"SemidefiniteCone", dim^2}]} ,
Flatten[{x,y}]];
X=Chop[sol[[1;;dim^2,2]].SaMCNB[dim],10^-3]; Y=Chop[sol[[dim^2+1;;2*dim^2,2]].SaMCNB[dim],10^-3];
Chop[Sqrt[Tr[X.ω+Y.ρ]],10^-3],X//MatrixForm,Y//MatrixForm},

{x=Table[Symbol["x"<>ToString[n]],{n,dim^4}];
sol=SemidefiniteOptimization[Chop[Simplify[Tr[(x.SaMCNB[dim^2]).C]],10^-3],
{ResourceFunction["MatrixPartialTrace"][x.SaMCNB[dim^2],2,{dim,dim}]==ω,ResourceFunction["MatrixPartialTrace"][x.SaMCNB[dim^2],1,{dim,dim}]==If[transpose,Transpose[ρ],ρ],
 VectorGreaterEqual[{x.SaMCNB[dim^2],0}, {"SemidefiniteCone", dim^2}] },
x];
Π=Chop[Sum[sol[[n,2]]SaMCNB[dim^2][[n]],{n,1,dim^4}],10^-3];
Chop[Sqrt[Tr[C.Π]],10^-3],Π//MatrixForm}]]

ModQOT [ρ_,ω_,cost_,dual_,transpose_]:=Sqrt[QOT[ρ,ω,cost,dual,transpose][[1]]^2-(QOT[ρ,ρ,cost,dual,transpose][[1]]^2+QOT[ω,ω,cost,dual,transpose][[1]]^2)/2]

TriIneq[ρ_,ω_,τ_,C_,dual_,transpose_]:=ModQOT[ρ,ω,C,dual,transpose]+ModQOT[ω,τ,C,dual,transpose]-ModQOT[ρ,τ,C,dual,transpose]

ρ01=1/2*{1,1/Sqrt[2],1/Sqrt[3],0}.SaMPB[2]

τ01=1/2*{1,0,1/3,1/4}.SaMPB[2]

step1=0.01;

coordinates1=Select[Flatten[Table[{x,y},{x,-1,1,step1},{y,-1,1,step1}],1],PositiveSemidefiniteMatrixQ[ 1/2*(Flatten[{1,#,1/Sqrt[2]}]).SaMPB[2]]&];

ListPlot[coordinates1,AspectRatio->1]

Cost1=CostfromObservable[SaMPB[2][[{2,4}]],True]

ω1=1/2*(Flatten[{1,#,1/Sqrt[2]}]).SaMPB[2]&/@coordinates1;

Off[SemidefiniteOptimization::parsuc]

{time1,result1}=AbsoluteTiming[TriIneq[ρ0,#,τ0,Cost1,False,True]&/@ω1];

time1/3600

ListPlot3D[Join[coordinates1,Transpose@{result1},2],AspectRatio->1]

ListDensityPlot[Join[coordinates1,Transpose@{result1},2],AspectRatio->1]

Max[Abs[Im[result1]]]

Min[result1]

ρ02=1/4*{1,0,0,0,0,0.1,0,0,0.2,0,0,0,0.3,0,0,0}.SaMPB[4]

τ02=1/4*{1,0,0,0.3,0,0,0,0.2,0,0,0,0,0.1,0,0,0}.SaMPB[4]

step2=0.01;

coordinates2=Select[Flatten[Table[{x,y},{x,-1,1,step2},{y,-1,1,step2}],1],PositiveSemidefiniteMatrixQ[ 1/4*(Flatten[{1,#,0,0.1,0.1,0.1,0,0.3,0,0.2,0,0,0,0,0}]).SaMPB[4]]&];

ListPlot[coordinates2,AspectRatio->1]

Cost2=CostfromObservable[SaMPB[4][[2;;16]],True]

Cost2//MatrixForm

ω2=1/4*(Flatten[{1,#,0,0.1,0.1,0.1,0,0.3,0,0.2,0,0,0,0,0}]).SaMPB[4]&/@coordinates2;

Off[SemidefiniteOptimization::parsuc]

{time2,result2}=AbsoluteTiming[TriIneq[ρ02,#,τ02,Cost2,False,True]&/@ω2];

time2/3600

ListPlot3D[Join[coordinates2,Transpose@{result2},2],AspectRatio->1]

ListDensityPlot[Join[coordinates2,Transpose@{result2},2],AspectRatio->1]

Max[Abs[Im[result2]]]

Min[result2]

ρ03=RandState[2,2]

%//MatrixForm

τ03=RandState[2,2]

%//MatrixForm

step3=0.01;

coordinates3=Select[Flatten[Table[{x,y},{x,-1,1,step3},{y,-1,1,step3}],1],PositiveSemidefiniteMatrixQ[ 1/2*(Flatten[{1,#,0.2}]).SaMPB[2]]&];

ListPlot[coordinates3,AspectRatio->1]

Cost3=CostfromObservable[{RandSelfadjMatrix[2],RandSelfadjMatrix[2],RandSelfadjMatrix[2]},True]

%//MatrixForm

ω3=1/2*(Flatten[{1,#,0.2}]).SaMPB[2]&/@coordinates3;

Off[SemidefiniteOptimization::parsuc]

{time3,result3}=AbsoluteTiming[TriIneq[ρ03,#,τ03,Cost3,False,True]&/@ω3];

time3/3600

ListPlot3D[Join[coordinates3,Transpose@{result3},2]]

ListDensityPlot[Join[coordinates3,Transpose@{result3},2]]

Max[Abs[Im[result3]]]

Min[result3]

ρ04=RandState[4,4]

%//MatrixForm

τ04=RandState[4,4]

%//MatrixForm

step4=0.01;

coordinates4=Select[Flatten[Table[{x,y},{x,-1,1,step4},{y,-1,1,step4}],1],PositiveSemidefiniteMatrixQ[ 1/4*(Flatten[{1,#,0,0.1,0.1,0.1,0,0.3,0,0.2,0,0,0,0,0}]).SaMPB[4]]&];

ListPlot[coordinates4,AspectRatio->1]

Cost4=CostfromObservable[{RandSelfadjMatrix[4],RandSelfadjMatrix[4],RandSelfadjMatrix[4],RandSelfadjMatrix[4],RandSelfadjMatrix[4]},True];

Round[Cost4,1]//MatrixForm

ω4=1/4*(Flatten[{1,#,0,0.1,0.1,0.1,0,0.3,0,0.2,0,0,0,0,0}]).SaMPB[4]&/@coordinates4;

Off[SemidefiniteOptimization::parsuc]

{time4,result4}=AbsoluteTiming[TriIneq[ρ04,#,τ04,Cost4,False,True]&/@ω4];

time4/3600

ListPlot3D[Join[coordinates4,Transpose@{result4},2],AspectRatio->1]

ListDensityPlot[Join[coordinates4,Transpose@{result4},2],AspectRatio->1]

Max[Abs[Im[result4]]]

Min[result4]