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]

Cost1=CostfromObservable[{RandSelfadjMatrix[2],RandSelfadjMatrix[2],RandSelfadjMatrix[2]},True];

Cost2=CostfromObservable[{RandSelfadjMatrix[2],RandSelfadjMatrix[2],RandSelfadjMatrix[2]},True];

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

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

ρ01=RandState[2,2]; τ01=RandState[2,2];

ρ02=RandState[2,2]; τ02=RandState[2,2];

ρ03=RandState[2,2]; τ03=RandState[2,2];

ρ04=RandState[2,2]; τ04=RandState[2,2];

step=0.1;

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

ListPointPlot3D[coordinates,BoxRatios->{1,1,1}]

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

Off[SemidefiniteOptimization::parsuc]

ToExpression[ToString[Flatten[Table[{"time"<>ToString[i]<>ToString[j],"result"<>ToString[i]<>ToString[j]},{i,4},{j,4}],1]]<>"=Flatten[Table[AbsoluteTiming[TriIneq[ToExpression[\"ρ0\"<>ToString[i]],#,ToExpression[\"τ0\"<>ToString[i]],ToExpression[\"Cost\"<>ToString[j]],False,True]&/@ω2dim],{i,4},{j,4}],1];"]

Join[{{"time","result"}},ToExpression[ToString[Flatten[Table[{"time"<>ToString[i]<>ToString[j],"Min[result"<>ToString[i]<>ToString[j]<>"]"},{i,4},{j,4}],1]]]]//MatrixForm

{dim3time,dim3result}=AbsoluteTiming[Table[TriIneq[RandState[3,3],RandState[3,3],RandState[3,3],CostfromObservable[{RandSelfadjMatrix[3],RandSelfadjMatrix[3],RandSelfadjMatrix[3]},True],False,True],4000]];

dim3time/3600

Min[dim3result]

{dim4time,dim4result}=AbsoluteTiming[Table[TriIneq[RandState[4,4],RandState[4,4],RandState[4,4],CostfromObservable[{RandSelfadjMatrix[4],RandSelfadjMatrix[4],RandSelfadjMatrix[4]},True],False,True],4000]];

dim4time/3600

Min[dim4result]

{dim5time,dim5result}=AbsoluteTiming[Table[TriIneq[RandState[5,5],RandState[5,5],RandState[5,5],CostfromObservable[{RandSelfadjMatrix[5],RandSelfadjMatrix[5],RandSelfadjMatrix[5]},True],False,True],4000]];

dim5time/3600

Min[dim5result]