(*Some defs I use later*)
P0={{1,0},{0,1}};
PX={{0,1},{1,0}};
PY={{0,-I},{I,0}};
PZ={{1,0},{0,-1}};
BlockMatrix[bloecke_]:=ArrayFlatten[bloecke];
TP[m1_,m2_]:=BlockMatrix[Table[Table[m1[[i,j]]*m2,{j,1,Dimensions[m1][[2]]}],{i,1,Dimensions[m1][[1]]}]];

TPP[matrixliste_]:=Module[{ii,xx,imax},
xx=matrixliste[[1]];
imax=Dimensions[matrixliste][[1]]+1;
For[ii=2,ii<imax,
xx=TP[xx,matrixliste[[ii]]];
ii++;];
xx
];
makeDensity[st_]:=TP[{st},ConjugateTranspose[{st}]];
makeDensityReal[st_]:=TP[{st},Transpose[{st}]];

(*Sheet 6 Bell Inequalities*)
(Subscript[A, 1])=PX; (*setting Alice's first basis*)
(Subscript[A, 2])=PZ;(*setting Alice's Second basis*)
B[b_]:=Cos[b]PX+Sin[b]PZ;  (*Parametrising Bob's  bases*)
rho=makeDensityReal[{Sin[a],0,0,Cos[a]}] (*an arbitrary 2 qubit density matrix in a normal form*)

BI=Tr[(TP[(Subscript[A, 1]),B[(Subscript[b, 1])]]+TP[(Subscript[A, 1]),B[(Subscript[b, 2])]]+TP[(Subscript[A, 2]),B[(Subscript[b, 1])]]-TP[(Subscript[A, 2]),B[(Subscript[b, 2])]]).rho]//FullSimplify  (*Bell Inquality*)

(Cos[(Subscript[b, 1])]+Cos[(Subscript[b, 2])]) Sin[2 a]+Sin[(Subscript[b, 1])]-Sin[(Subscript[b, 2])]/.{(Subscript[b, 1])-> -(Subscript[b, 2])}//FullSimplify (*Just a guess, which appears to work*)

D[2 Cos[(Subscript[b, 2])] Sin[2 a]-2 Sin[(Subscript[b, 2])],(Subscript[b, 2])]//FullSimplify

sol=Solve[D[2 Cos[(Subscript[b, 2])] Sin[2 a]-2 Sin[(Subscript[b, 2])],(Subscript[b, 2])]==0,(Subscript[b, 2])]//FullSimplify

Assuming[(C[1])∈(Integers),FullSimplify[2 Cos[(Subscript[b, 2])] Sin[2 a]-2 Sin[(Subscript[b, 2])]/.sol[[2]]]]

Plot[{2 (Sqrt[1+((Sin[2 a])^2)]),2,2(Sqrt[2])},{a,0, Pi/4}] (*Plotting Violation*)

(*Sheet 9 Code*)

(*a state in the sheet Eq.1 *)
ρ[p_]:=p makeDensity[{1,0}]+(1-p)makeDensity[{0,1}];
(*quick check for the trace*)
Tr[ρ[p]]

(*Define vectors in Kraus Operators*)
αp[α_]:=(Sqrt[α]){1,0}+(Sqrt[1-α]){0,1};
αm[α_]:=(Sqrt[1-α]){1,0}-(Sqrt[α]){0,1};
(*check orthonormality*)
{αp[α].αp[α],αm[α].αm[α],αp[α].αm[α]}

(*Define Kraus Operators*)
Ap[α_,k_]:=(Sqrt[k])makeDensityReal[αp[α]]+(Sqrt[1-k])makeDensityReal[αm[α]];
Am[α_,k_]:=(Sqrt[1-k])makeDensityReal[αp[α]]+(Sqrt[k])makeDensityReal[αm[α]];
(*Do they sum up to Id?*)
Ap[a,k].Transpose[Ap[a,k]]+Am[a,k].Transpose[Am[a,k]]//FullSimplify//MatrixForm

(*here newρ is defined*)
tildeρ[p_,a_,k_]:=Ap[a,k].ρ[p].Transpose[Ap[a,k]]+Am[a,k].ρ[p].Transpose[Am[a,k]];
(*Check trace again*)
Tr[tildeρ[p,a,k]]//FullSimplify

(*Now part a) *)

Plot3D[{Tr[tildeρ[0,a,k].tildeρ[0,a,k]],1/2},{a,0,1},{k,0,1}, AxesLabel->Automatic]

Plot3D[{Tr[tildeρ[0.2,a,k].tildeρ[0.2,a,k]]},{a,0,1},{k,0,1}, AxesLabel->Automatic]

Plot3D[{Tr[tildeρ[0.4,a,k].tildeρ[0.4,a,k]]},{a,0,1},{k,0,1}, AxesLabel->Automatic]

(*b) Maximum average increase in Purity*)

(*Second derivative is always positive ==>  Critical point will correspond to Minimum*)
D[Tr[tildeρ[p,a,k].tildeρ[p,a,k]]-Tr[ρ[p].ρ[p]],{a,2}]//FullSimplify

(* Let's check: First derivative*)
Solve[D[Tr[tildeρ[p,a,k].tildeρ[p,a,k]]-Tr[ρ[p].ρ[p]],{a,1}]==0,a,Reals]//FullSimplify

Tr[tildeρ[p,a,k].tildeρ[p,a,k]]-Tr[ρ[p].ρ[p]]/.{a-> 1/2}//FullSimplify
(*This expression is negative for almost all p and k, so, we only get dicrease in purity*)

(*What do we do in this case? Check the boundaries! a=0 or a=1*)
Tr[tildeρ[p,a,k].tildeρ[p,a,k]]-Tr[ρ[p].ρ[p]]/.{a->0}//FullSimplify
Tr[tildeρ[p,a,k].tildeρ[p,a,k]]-Tr[ρ[p].ρ[p]]/.{a->1}//FullSimplify

Tr[tildeρ[p,a,k].tildeρ[p,a,k]/.{a-> 0}]//FullSimplify
Tr[ρ[p].ρ[p]]//FullSimplify

Tr[tildeρ[p,a,k].tildeρ[p,a,k]]-Tr[ρ[p].ρ[p]]/.{a-> 1}//FullSimplify