(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
Print["Qulib is now Melt! Visit the new page at https://melt1.notion.site"]
Print["This version of qulib will remain available, but will no longer be updated."]


(* ::Input::Initialization:: *)
(* Prints in matrix form but does not alter definition *)
Clear[mf];
mf[A_]:=(Print[MatrixForm@Chop@A];A)
mf::usage = "mf[A] will print the matrix in a pretty way using Mathematica's MatrixForm[] function. The printed version also chops any small numbers.";


(* ::Input::Initialization:: *)
(* FullSimplifty + ComplexExpand (to get rid of "Conjugate"s) *)
Clear[cf];
(*cf[expr_]:=FullSimplify@ComplexExpand@expr;*)
cf[expr_]:=FullSimplify[expr,Thread[DeleteDuplicates@Cases[Normal@expr,_Symbol,\[Infinity]]>0]]

cf::usage = "cf[expr] performs a FullSimplify on the expression, assuming that all symbols are real and positive.";


(* ::Input::Initialization:: *)
(* Outer product |v\[RightAngleBracket]\[LeftAngleBracket]u| between two vectors *)
(* Properly conjugates *)
Clear[out];
out[v_,u_]:=Outer[Times,v,Conjugate[u]];
out::usage = "Outer product: out[v,u] = |v\[RightAngleBracket]\[LeftAngleBracket]u|";


(* ::Input::Initialization:: *)
(* \[DoubleStruckCapitalC] writes as esc+dsC+esc *)
Clear[\[DoubleStruckCapitalC]];
\[DoubleStruckCapitalC][A_,B_]:=A . B-B . A;
\[DoubleStruckCapitalC][A_,B_,\[Zeta]_]:=A . B-\[Zeta] B . A

\[DoubleStruckCapitalC]::usage = "\[DoubleStruckCapitalC][A,B] = [A,B] = A.B-B.A. More generally \[DoubleStruckCapitalC][A,B,\[Zeta]] = A.B + \[Zeta] B.A. Use \[Zeta] = -1 for anti-commutator";


(* ::Input::Initialization:: *)
(* Quick way of checking what are the distinct entries of a matrix or tensor *)
Clear[ChopCheck];
ChopCheck[A_,tol_:10^-10]:=DeleteDuplicates@Flatten@Chop[A,tol]

ChopCheck::usage = "ChopCheck[A] gives 0 when the matrix is 0, to a tolerance of \!\(\*SuperscriptBox[\(10\), \(-10\)]\). I usually do ChopCheck[A-B] whenever I want to check numerically that A = B";


(* ::Input::Initialization:: *)
(* Use just like Range[] function *)
Clear[range];
range[xi_,xf_,df_:1]:=N@Rationalize@Range[xi,xf,df]; 
range::usage = "Same input/output as Range, but avoids a subtle floating point issue of the latter.";

(* Stolen from python/Matlab *)
Clear[linspace]
linspace[dd1_,dd2_,n_:100]:=With[{d1=N@dd1,d2=N@dd2},If[n==1,{d1},d1+(d2-d1)/(n-1) Range[0,n-1]]];
linspace::usage = "linspace[a,b,npts] yields npts linearly interpolated between a (inclusive) and b (inclusive).";

Clear[logspace];
logspace[dd1_,dd2_,n_:100]:=10^linspace[dd1,dd2,n];
logspace::usage = "Same as linspace, but in log scale. logscale[2,5,50] will give 50 pts between \!\(\*SuperscriptBox[\(10\), \(2\)]\) and \!\(\*SuperscriptBox[\(10\), \(5\)]\)";


(* ::Input::Initialization:: *)
Clear[PartitionIn];
PartitionIn[list_,x_,pad_:""]:=Partition[list,Ceiling[Length[list]/x],Ceiling[Length[list]/x],1,pad]

PartitionIn::usage = "PartitionIn[list,x] partitions list into x sublists, possibly with paddings at the end.";


(* ::Input::Initialization:: *)
Clear[PrettyTiming];
SetAttributes[PrettyTiming,HoldAll];
PrettyTiming[computation_]:=Module[{t1,t2,output,time,h,m,s},
t1 = UnixTime[];
output=ReleaseHold[computation];
t2 = UnixTime[];
time = t2-t1;
 h = Floor[time/3600];
m = Floor[Mod[time,3600]/60];
s = Round@Mod[time,60];
Print[ToString[h]<>"h : "<>ToString[m]<>"m : "<>ToString[s]<>"s"];
(*Print[PrettyTime[t2-t1]];*)
Return[output]];

PrettyTiming::usage = "PrettyTiming[computation] uses Mathematica's AbsoluteTiming[] to time a computation, but prints the result in a prettier form.";


(* ::Input::Initialization:: *)
Clear[do];
SetAttributes[do,HoldAll];
do[computation_,iterators__]:=Module[{iters = {iterators}[[All,1]],symb,list,i,d },
symb = ToString/@iters;
d=Length@symb;
list = Flatten[Table[iters,iterators],Length[iters]-1];

PrettyTiming@Monitor[Do[computation,iterators],
Check[
StringJoin@Table[{symb[[i]]," = ", ToString@StandardForm[iters[[i]]],"\n"},{i,d}]
ProgressIndicator[Position[list,Table[iters[[i]],{i,d}]][[1,1]],{1,Length@list}]
," "]
]
]


(* ::Input::Initialization:: *)
Clear[table];
SetAttributes[table,HoldAll];
table[computation_,iterators__]:=Module[{iters = {iterators}[[All,1]],symb,list,i,d },
symb = ToString/@iters;
d=Length@symb;
list = Flatten[Table[iters,iterators],Length[iters]-1];

PrettyTiming@Monitor[Table[computation,iterators],
Check[
StringJoin@Table[{symb[[i]]," = ", ToString@StandardForm[iters[[i]]],"\n"},{i,d}]
ProgressIndicator[Position[list,Table[iters[[i]],{i,d}]][[1,1]],{1,Length@list}]
," "]
]
]


(* ::Input::Initialization:: *)
Clear[nest];
nest[f_,init_,n_]:=Module[{it=0},PrintTemporary[ProgressIndicator[Dynamic[it/n]] Dynamic[ToString@Round[ 100 it/n]<>"%"]];
Nest[(it++;f[#])&,init,n]]


(* ::Input::Initialization:: *)
Clear[nestList];
nestList[f_,init_,n_]:=Module[{it=0},PrintTemporary[ProgressIndicator[Dynamic[it/n]] Dynamic[ToString@Round[ 100 it/n]<>"%"]];
NestList[(it++;f[#])&,init,n]]


(* ::Input::Initialization:: *)
Clear[makeGIF];
SetAttributes[makeGIF,HoldAll];
makeGIF[thing_,iterator_,filename_,duration_,fps_:20]:=Module[{ i0,if,\[CapitalDelta],\[Delta],frames,numberOfframes},

numberOfframes = duration fps;
i0 = iterator[[2]];
if = iterator[[3]];
\[Delta] = If[Length@iterator==4,iterator[[4]],1];
\[CapitalDelta] =  (if-i0)/numberOfframes ;   
If[IntegerQ[\[Delta]],\[CapitalDelta]=Ceiling[\[CapitalDelta]]];

If[Length@iterator==3,
frames = table[thing,Evaluate@Append[iterator, \[CapitalDelta]]],
frames = table[thing,Evaluate@ReplacePart[iterator,4-> \[CapitalDelta]]]
];
If[Length@frames < 0.8 numberOfframes,
Print["Could not match the desired framerate & duration. Reducing framerate to "<>ToString[(duration/Length@frames)^-1//Round]<>" fps"];
Print["Requires a longer list"]
];

Print["Exporting "<>filename<>".gif (can take a while...)"];
Export[filename <>".gif",frames,"DisplayDurations"->duration/Length@frames] ;
Animate[frames[[i]],{i,1,Length@frames,1}]
]


(* ::Input::Initialization:: *)
Clear[kron];
kron[mats__?MatrixQ]:=KroneckerProduct[mats];
kron[vecs__?VectorQ]:=Flatten@KroneckerProduct[vecs];

kron::usage ="kron[A,B,C...] = A \[CircleTimes] B \[CircleTimes] C \[CircleTimes] .... For matrices, it is exactly like the native KroneckerProduct[] from Mathematica. For vectors, kron[u,v] produces the proper |u\[RightAngleBracket] \[CircleTimes] |v\[RightAngleBracket]."


(* ::Input::Initialization:: *)
Clear[KronEye];
KronEye[A_,\[ScriptCapitalN]_,j_]:=kron[Eye[Length[A]^(j-1)],A,Eye[Length[A]^(\[ScriptCapitalN]-j)]]

KronEye::usage = "KronEye[A,\[ScriptCapitalN],j] = (\[DoubleStruckCapitalI] \[CircleTimes] ... \[DoubleStruckCapitalI] \[CircleTimes] A \[CircleTimes] \[DoubleStruckCapitalI] ... \[CircleTimes] \[DoubleStruckCapitalI], with A at site j and \[ScriptCapitalN] sites in total.";


(* ::Input::Initialization:: *)
Clear[Proj];
Proj[n_,i_]:=SparseArray[{{i,i}->1},{n,n}];
Proj[n_,i_,j_]:=SparseArray[{{i,j}->1},{n,n}];

Proj[{n_,m_},i_]:=SparseArray[{{i,i}->1},{n,m}];
Proj[{n_,m_},i_,j_]:=SparseArray[{{i,j}->1},{n,m}];


Proj::usage = "Proj[n,i,j] yields a n\[Times]n matrix with element (i,j) equal 1. Proj[n,i] = Proj[n,i,i].";


(* ::Input::Initialization:: *)
Clear[Basis];
Basis[n_,i_]:=SparseArray[{i->1},n]
Basis[n_,ilist_?ArrayQ]:=Total@Table[Basis[n,i],{i,ilist}]

Basis::usage = "Basis[n,i] = (0,...,0,1,0,...,0) with 1 at position i. Basis[n,{i1,i2,...}] yields the same, but with 1's at positions i1,i2,...";


(* ::Input::Initialization:: *)
Clear[Eye,NEye];
Eye[n_]:=SparseArray[{{i_,i_}->1},{n,n}];
Eye::usage = "Eye[n] is the n\[Times]n identity matrix as a sparse array. Same as NEye, but the 1s in the diagonal are integers";

NEye[n_]:=SparseArray[{{i_,i_}->1.},{n,n},0.];
Eye::usage = "NEye[n] is the n\[Times]n identity matrix as a sparse array. Same as Eye, but the 1.0s in the diagonal are in floating point.";


(* ::Input::Initialization:: *)
Clear[toeplitzMatrix];
toeplitzMatrix[a_,k_,L_]:=If[L==1,If[k==1,{{a}},{{0}}],SparseArray[{Band[{1,k}]->a,Band[{k,1}]->a},{L,L}]]

toeplitzMatrix::usage = "toeplitzMatrix[a,k,L] yields a L\[Times]L symmetric matrix with elements a along the kth diagonal (k=1 is the diagonal and k=2 is the first off-diagonal).";


(* ::Input::Initialization:: *)
Clear[StiffnessMatrix,CirculantMatrix];
StiffnessMatrix[n_]:=toeplitzMatrix[-1,2,n];
CirculantMatrix[n_]:=toeplitzMatrix[-1,2,n]+toeplitzMatrix[-1,-1,n];

StiffnessMatrix::usage = "StiffnessMatrix[n]: hopping matrix of the tight-binding model with open boundary conditions";
CirculantMatrix::usage = "CirculantMatrix[n]: hopping matrix of the tight-binding model with periodic boundary conditions";


(* ::Input::Initialization:: *)
Clear[FirstDerivative];
FirstDerivative[n_,\[CapitalDelta]x_:1,boundary_:None]:=
If[boundary===None,
1/(2\[CapitalDelta]x) SparseArray[{
Band[{1,2}]->1,
Band[{2,1}]->-1
},{n,n}],
1/(2\[CapitalDelta]x) SparseArray[{
{i_,j_}/;j==i+1&&i!=1->1,
{i_,j_}/;j==i-1&&i!=n->-1,
{1,2}->2,
{1,1}->-2,
{n,n}->2,
{n,n-1}->-2
},{n,n}]];

Clear[SecondDerivative]
SecondDerivative[n_,\[CapitalDelta]x_:1,boundary_:None]:=
If[boundary===None,
1/\[CapitalDelta]x^2 (2Eye[n]+StiffnessMatrix[n])]



(* ::Input::Initialization:: *)
Clear[Boxcar];
Boxcar[a_,b_][x_]:=UnitStep[x-a]-UnitStep[x-b]

Boxcar[a_,b_,smooth_][x_]:=If[Abs[smooth]<10^-14,Boxcar[a,b][x],1/2 (Tanh[(x-a)/smooth]-Tanh[(x-b)/smooth] )]

Boxcar::usage = "Boxcar[a,b][x] returns a boxcar function from a to b";


(* ::Input::Initialization:: *)
Clear[Eigvals,Eigvecs,Eigsys];

(*Eigvals[A_,neigs_:"all",\[CapitalDelta]_:10^5]:=Reverse@Eigenvalues[A + \[CapitalDelta] NEye@Length@A,If[neigs\[Equal]"all",Length@A,-neigs]]-\[CapitalDelta];*)
Eigvals[A_,neigs_:"all",\[CapitalDelta]_:10^5]:=If[neigs==="all",
Sort@Eigenvalues[N@A],
Reverse@Eigenvalues[A + \[CapitalDelta] NEye@Length@A,-neigs]-\[CapitalDelta]];

(*Eigvecs[A_,neigs_:"all",\[CapitalDelta]_:10^5]:=Transpose[(Reverse@Eigenvectors[A + \[CapitalDelta] NEye@Length@A,If[neigs\[Equal]"all",Length@A,-neigs]])];*)
Eigvecs[A_,neigs_:"all",\[CapitalDelta]_:10^5]:=
Module[{tmp,e,v,ord},
If[neigs=="all",(
{e,v}=Eigensystem[A];
ord=Ordering[e];
Transpose[v[[ord]]]
),
(
{e,v}=Eigensystem[A + \[CapitalDelta] NEye@Length@A,-neigs];
 Transpose@Reverse@v
)
]
]

(*Eigsys[A_, neigs_:"all",\[CapitalDelta]_:10^5]:=Module[{tmp},
tmp=Eigensystem[A + \[CapitalDelta] NEye@Length@A,If[neigs\[Equal]"all",Length@A,-neigs]];
 {Reverse@First@tmp-\[CapitalDelta],Transpose@Reverse@Last@tmp}
];*)

Clear[Eigsys];
Eigsys[A_, neigs_:"all",\[CapitalDelta]_:10^5]:=Module[{tmp,e,v,ord},

If[neigs=="all",(
{e,v}=Eigensystem[A];
ord=Ordering[e];
{e[[ord]],Transpose[v[[ord]]]}
),
(
{e,v}=Eigensystem[A + \[CapitalDelta] NEye@Length@A,-neigs];
 {Reverse@e-\[CapitalDelta],Transpose@Reverse@v}
)
]
]




(* ::Input::Initialization:: *)
Clear[GroundState];
GroundState[H_]:=Eigvecs[H,1][[All,1]]

GroundState[H_,n_]:=Eigvecs[H,n][[All,-1]];


(* ::Input::Initialization:: *)
(* y can be a ket or a density matrix *)
Clear[expect];
expect[A_,y_]:=If[Length@Dimensions@y==2,Tr[A . y],y\[Conjugate] . A . y];

expect::usage = "expect[A,y] computes \[LeftAngleBracket]A\[RightAngleBracket] in the state y, which can be either a ket or a density matrix";


(* ::Input::Initialization:: *)
(* Checks if H commutes with y, and y can be a pure state *)
Clear[checkCommute]
checkCommute[H_,y_]:=Module[{\[Rho]},
\[Rho] = If[Length@Dimensions@y==2,y,Outer[Times,y,Conjugate[y]]];
If[Norm[H . \[Rho]-\[Rho] . H]>10^-10,False,True]
]


(* ::Input::Initialization:: *)
Clear[DirectSum]
DirectSum[mats__]:=Module[{m=List@mats,tab, diag,a},
tab= Table[a[i],{i,Length@m}];
diag=DiagonalMatrix[Table[a[i],{i,Length@m}]];
diag/.Thread[tab->m]//ArrayFlatten
]

DirectSum::usage = "DirectSum[A,B,...] yields a diagonal block matrix diag(A,B,...)";


(* ::Input::Initialization:: *)
Clear[RescaleHamiltonian];
RescaleHamiltonian[H_]:=Module[{e1,ed},
e1 = Eigvals[H,1][[1]];
ed = Eigvals[H,-1][[1]];

{e1,ed,(e1+ed)/(e1-ed) Eye[Length@H] + (2 H)/(ed-e1) }
]


(* ::Input::Initialization:: *)
Clear[KPM,KPMfunc];
KPM[H_,Mcheb_,R_,\[ScriptCapitalO]_:"1"]:=Module[{h,\[DoubleStruckCapitalO],gs,rs,\[Nu]s,\[Mu]s,e1,ed,d=Length@H},
{e1,ed,h}=RescaleHamiltonian[H];

gs = Table[1/(Mcheb+1) ((Mcheb-m+1)Cos[(\[Pi] m)/(Mcheb+1.)] + Sin[(\[Pi] m)/(Mcheb+1)] Cot[\[Pi]/(Mcheb+1)] ),{m,0,Mcheb}]; 

\[DoubleStruckCapitalO] = SparseArray@If[\[ScriptCapitalO] ==="1",Eye[d],\[ScriptCapitalO]];

rs = Normalize/@(RandomVariate[NormalDistribution[0,1],{R,d}]+I RandomVariate[NormalDistribution[0,1],{R,d}]);
\[Nu]s = Table[With[{v01= {r,h . r}},NestList[{#[[2]],2 h . #[[2]]-#[[1]]}&,v01,Mcheb][[All,2]]],{r,rs}] ;
\[Mu]s = Prepend[Table[d/R Total@Table[rs[[i]]\[Conjugate] . \[DoubleStruckCapitalO] . \[Nu]s[[i,m]],{i,Length@rs}],{m,1,Mcheb}],1]; 

{e1,ed,gs,\[Mu]s}
]

KPMfunc[{e1_,ed_,gs_,\[Mu]s_},\[ScriptCapitalE]_]:=With[{\[Epsilon]=(e1+ed)/(e1-ed) + (2 \[ScriptCapitalE])/(ed-e1)}, 1/(\[Pi] Sqrt[1-\[Epsilon]^2]) (gs[[1]] \[Mu]s[[1]] + 2 Total@Table[gs[[m+1]] \[Mu]s[[m+1]] ChebyshevT[m,\[Epsilon]],{m,1,Length@\[Mu]s-1}])]


(* ::Input::Initialization:: *)
Clear[LoadPauliMatrices];
LoadPauliMatrices[]:=Module[{},
Clear[\[Sigma]0,\[Sigma]x,\[Sigma]y,\[Sigma]z,\[Sigma]p,\[Sigma]m,GMAT];
\[Sigma]0 = {{1,0},{0,1}};
\[Sigma]x = {{0,1},{1,0}};
\[Sigma]y = {{0,-I},{I,0}};
\[Sigma]z = {{1,0},{0,-1}};
\[Sigma]p = {{0,1},{0,0}}(*(\[Sigma]x + I \[Sigma]y)/2*);
\[Sigma]m = {{0,0},{1,0}}(*(\[Sigma]x - I \[Sigma]y)/2*);
(*nn = \[Sigma]p.\[Sigma]m; *)

GMAT[\[Theta]_,\[Phi]_] := MatrixExp[-I \[Phi] \[Sigma]z/2] . MatrixExp[-I \[Theta] \[Sigma]y/2];

"Matrices loaded: \[Sigma]0 (=1), \[Sigma]x, \[Sigma]y, \[Sigma]z, \[Sigma]p, \[Sigma]m, GMAT[\[Theta],\[Phi]]"
];


(* ::Input::Initialization:: *)
LoadPauliMatrices[]


(* ::Input::Initialization:: *)
LoadSpinChainMatrices[\[ScriptCapitalN]_,type_:"n"]:=Module[{one,\[Sigma]0,\[Sigma]x,\[Sigma]y,\[Sigma]z,\[Sigma]p,\[Sigma]m},

(* 2 x 2 Pauli matrices *)
one = If[type=="n",1.0,1];
\[Sigma]0 = SparseArray[one{{1,0},{0,1}}];
\[Sigma]x = SparseArray[one{{0,1},{1,0}}];
\[Sigma]y = SparseArray[one{{0,-I},{I,0}}];
\[Sigma]z =SparseArray[ one{{1,0},{0,-1}}];
\[Sigma]p =(\[Sigma]x + I \[Sigma]y)/2;
\[Sigma]m = (\[Sigma]x - I \[Sigma]y)/2;

Clear[eye,sp,sm,sx,sy,sz,c,cd];
eye = one Eye[2^\[ScriptCapitalN]];

If[\[ScriptCapitalN]>1,
Do[
eyeL = one Eye[2^(i-1)];
eyeR = one Eye[2^(\[ScriptCapitalN]-i)];
(* Spin matrices of site i *)
sx[i] = kron[eyeL,\[Sigma]x,eyeR];
sy[i] = kron[eyeL,\[Sigma]y,eyeR];
sz[i] = kron[eyeL,\[Sigma]z,eyeR];
sp[i] = kron[eyeL,\[Sigma]p,eyeR];
sm[i] = sp[i]\[Transpose];

(* Fermionic operators using the Jordan Wigner transformation *)
c[i] = kron@@Join[Table[(-\[Sigma]z),{j,1,i-1}],{\[Sigma]m},Table[\[Sigma]0,{j,i+1,\[ScriptCapitalN]}]];
cd[i] = c[i]\[ConjugateTranspose];

,{i,1,\[ScriptCapitalN]}],
(
sx[1] = \[Sigma]x;
sy[1] =\[Sigma]y;
sz[1] = \[Sigma]z;
sp[1] = \[Sigma]p;
sm[1] = \[Sigma]m;
c[1] = \[Sigma]m;
cd[1] = \[Sigma]p;)];

"Matrices loaded: \[Sigma]0 (=1), \[Sigma]x, \[Sigma]y, \[Sigma]z, \[Sigma]p, \[Sigma]m, eye, sm, sp, sx, sy, sz, c, cd"

];


(* ::Input::Initialization:: *)
Clear[LoadFermionicOperators];
LoadFermionicOperators[\[ScriptCapitalN]_,ord_:"local",type_:"n"]:=Module[{},

LoadSpinChainMatrices[\[ScriptCapitalN],type];
Clear[sp,sm,sx,sy,sz];
Clear[\[CapitalGamma]];
If[ord=="local",
Do[
\[CapitalGamma][2n-1] = c[n] + cd[n];
\[CapitalGamma][2n]= I (cd[n]-c[n]),
{n,1,\[ScriptCapitalN]}],
Do[
\[CapitalGamma][n] = c[n] + cd[n];
\[CapitalGamma][\[ScriptCapitalN]+n]= I (cd[n]-c[n]),
{n,1,\[ScriptCapitalN]}]
];

Clear[sp,sm,sx,sy,sz];
"Matrices loaded: eye, c, cd, \[CapitalGamma]"
]


(* ::Input::Initialization:: *)
LoadArbitrarySpinMatrices[S_,type_:"n"]:=Module[{it,srange,tmp},

it=If[type=="a",1,1.0];
srange = Range[S,-S,-it];

Clear[Sz,Sp,Sm,Sx,Sy,S0];

Sz= SparseArray@DiagonalMatrix[srange];
Sp = SparseArray@Table[Sqrt[(S-s)(S+s+1)]KroneckerDelta[r,s+1],{r,srange},{s,srange}];
Sm = Sp\[Transpose];
Sx = (Sp+Sm)/2;
Sy = (Sp-Sm)/(2I);

S0 =Eye[2S+1];


"Matrices loaded: S0 (=1), Sx, Sy, Sz, Sp, Sm"

];



(* ::Input::Initialization:: *)
LoadBosonicOperators[\[ScriptCapitalN]_,type_:"n"]:=Module[{one,tmp,n,m},
one = If[type=="a",1,1.0];
Clear[\[DoubleStruckCapitalI],a,X,P(*,CoherentState*)];
\[DoubleStruckCapitalI]= Eye[\[ScriptCapitalN]];
a=one DiagonalMatrix[SparseArray[Table[Sqrt[n],{n,1,\[ScriptCapitalN]-1}]],1];
\[DoubleStruckCapitalX] = 1/Sqrt[2] (a\[ConjugateTranspose]+a);
\[DoubleStruckCapitalP] = I/Sqrt[2] (a\[ConjugateTranspose]-a);
CoherentState[\[Alpha]_]:=Exp[-Abs[\[Alpha]]^2/2] Table[\[Alpha]^n/Sqrt[n!],{n,0,\[ScriptCapitalN]-1}];
ThermalState[n0_]:=DiagonalMatrix@Table[n0^n (1+n0)^(-1-n),{n,0,\[ScriptCapitalN]-1}];
"Matrices loaded: \[DoubleStruckCapitalI] (=1), a, \[DoubleStruckCapitalX], \[DoubleStruckCapitalP], CoherentState[\[Alpha]], ThermalState[n0]"
]


(* ::Input::Initialization:: *)
Clear[Trapz];
Trapz[fs_,\[CapitalDelta]x_]:=\[CapitalDelta]x/2 (fs[[1]] + 2 Total[fs[[2;;-2]]] + fs[[-1]]);

Trapz::usage = "Trapz[fs,\[CapitalDelta]x] integrates the list fs using the trapezoidal rule, with step \[CapitalDelta]x";


(* ::Input::Initialization:: *)
Clear[BracketRoot];
BracketRoot[f_,xmin_,xmax_,npts_:1000]:=Module[{xrange,tab,par,roots,pos,xtab,x},
xrange=linspace[xmin,xmax,npts];
tab = f/@xrange;
par = {Range[1,Length@tab-1],Most[tab] Rest[tab]}\[Transpose];
pos =Select[par,#[[2]]<= 0.0&][[All,1]];
Table[{xrange[[p]],xrange[[p+1]]},{p,pos}]
]

BracketRoot::usage = "BracketRoot[f,xmin,xmax] tries to bracket all roots of f(x) within the interval [xmin,xmax].";


(* ::Input::Initialization:: *)
Clear[FindAllRoots];
FindAllRoots[f_,xmin_,xmax_,npts_:1000]:=Module[{x,xtab},
xtab = BracketRoot[f,xmin,xmax,npts];
Table[x/.FindRoot[f[x],{x,xt[[1]],xt[[2]]}],{xt,xtab}]
]

FindAllRoots::usage = "FindAllRoots[f,xmin,xmax] finds all roots of f(x) within the interval [xmin,xmax].";


(* ::Input::Initialization:: *)
Clear[LegendreGauss];
LegendreGauss[n_,prec_]:=Block[{roots,p,x,w},
roots=Solve[LegendreP[n,x]==0,x];
p=N[\[Pi]/4(x+1)/.roots,prec];
w=(\[Pi]/4)N[2(1-x^2)/((n+1)LegendreP[n+1,x])^2/.roots,prec];
{Join[p,p+\[Pi]/2],Join[w Sin[p],w Sin[p+\[Pi]/2]]}
];

LegendreGauss[n\[Theta]_,n\[Phi]_,prec_]:=Block[{p\[Theta],w\[Theta],p\[Phi],w\[Phi],p,w},
{p\[Theta],w\[Theta]} = LegendreGauss[n\[Theta],prec];
p\[Phi] = linspace[0,2\[Pi],n\[Phi]];
w\[Phi] = ConstantArray[(2\[Pi])/( n\[Phi]-1),{n\[Phi]}];  
p = Tuples[{p\[Theta],p\[Phi]}];
w = Times@@#&/@Tuples[{w\[Theta],w\[Phi]}];
{p,w}
]

LegendreGauss::usage = "Generates grid of (points,weights) for numerical integration in a sphere.";


(* ::Input::Initialization:: *)
Clear[RK4,RK4t];
(* For time-indeendent RHS f(y) *)
RK4[f_,y0_,\[CapitalDelta]t_,tf_,onlyLastOutput_:False]:=Module[{F},
(* input function f[y] *)
F[y_]:=Block[{k1,k2,k3,k4},
k1 = f[y];
k2 = f[y + \[CapitalDelta]t k1/2];
k3 = f[y + \[CapitalDelta]t k2/2];
k4 = f[y+\[CapitalDelta]t k3];
y + \[CapitalDelta]t/6 (k1+2k2+2k3+k4)
];
(*If[onlyLastOutput===False,
Transpose[{Range[0,tf,\[CapitalDelta]t],NestList[F,y0,Round[tf/\[CapitalDelta]t]]}],
Nest[F,y0,Round[tf/\[CapitalDelta]t]]
]
*)
If[onlyLastOutput===False,
Transpose[{Range[0,tf,\[CapitalDelta]t],nestList[F,y0,Round[tf/\[CapitalDelta]t]]}],
nest[F,y0,Round[tf/\[CapitalDelta]t]]
]
]

(* for a time-dependent RHS f(t,y) -> y' *)
RK4t[f_,y0_,\[CapitalDelta]t_,tf_,onlyLastOutput_:False]:=Module[{F},
(* input function f[t, y] *)
F[{t_,y_}]:=Block[{k1,k2,k3,k4},
k1 = f[t,y];
k2 = f[t+\[CapitalDelta]t/2,y + \[CapitalDelta]t k1/2];
k3 = f[t+\[CapitalDelta]t/2,y + \[CapitalDelta]t k2/2];
k4 = f[t+\[CapitalDelta]t,y+\[CapitalDelta]t k3];
{t+\[CapitalDelta]t,y + \[CapitalDelta]t/6 (k1+2k2+2k3+k4)}
];
(*If[onlyLastOutput===False,
NestList[F,{0,y0},Round[tf/\[CapitalDelta]t]],
Nest[F,{0,y0},Round[tf/\[CapitalDelta]t]]\[LeftDoubleBracket]2\[RightDoubleBracket]
]*)

If[onlyLastOutput===False,
nestList[F,{0,y0},Round[tf/\[CapitalDelta]t]],
nest[F,{0,y0},Round[tf/\[CapitalDelta]t]][[2]]
]

]


(* ::Input::Initialization:: *)
GramChalierEdgeworth[cumulants_]:=Module[{d = Length@cumulants,\[Mu],\[Sigma],\[Kappa]list,Bn},
\[Mu] = cumulants[[1]];
\[Sigma] = Sqrt[cumulants[[2]]];
\[Kappa]list = Join[{0,0},cumulants[[3;;-1]]];
Bn = Table[Sum[BellY[n,k,\[Kappa]list],{k,1,n}],{n,3,d}];

Function[1/Sqrt[2\[Pi] \[Sigma]^2] Exp[-((#-\[Mu])^2/(2\[Sigma]^2))] (1 + Sum[1/(n! \[Sigma]^n) Bn[[n-2]] 2^(-n/2) HermiteH[n,(#-\[Mu])/(\[Sigma] Sqrt[2])],{n,3,d}])]

(*1/Sqrt[2\[Pi] \[Sigma]^2]Exp[-((x-\[Mu])^2/(2\[Sigma]^2))] (1 + Sum[1/(n!\[Sigma]^n)Bn\[LeftDoubleBracket]n-2\[RightDoubleBracket] 2^(-n/2)HermiteH[n,(x-\[Mu])/(\[Sigma] Sqrt[2])],{n,3,d}])
*)]


(* ::Input::Initialization:: *)
(* The Ps do not have to be normalised *)
Clear[SillySample];
SillySample[xs_?VectorQ,Ps_?VectorQ,nsamp_:1]:=If[nsamp===1,RandomChoice[Ps->xs],RandomChoice[Ps->xs,nsamp]]
SillySample[P_?MatrixQ,nsamp_:1]:=SillySample[P[[All,1]],P[[All,2]],nsamp]


(* ::Input::Initialization:: *)
Clear[TwoPointCorrelation];
TwoPointCorrelation[data_,hmax_,dt_:None,steps_:1]:=Module[{\[ScriptCapitalN]=Length@data,ave = Mean[data],vec1,M,corr},
vec1 = data[[1;;\[ScriptCapitalN]-hmax]];
M = Table[data[[1+i;;\[ScriptCapitalN]-hmax+i]],{i,1,hmax,steps}];
corr=1/\[ScriptCapitalN] (M . vec1-ave^2);
If[dt===None,corr,{dt Range[1, hmax,steps],corr}\[Transpose]]
]


(* ::Input::Initialization:: *)
Clear[FourierWithFreqs];
FourierWithFreqs[x_,dt_,transform_:Fourier]:=Module[{d,xf,\[Omega]s,\[Omega]s2,hamm},
d = Length@x;
If[OddQ[d],d=d-1];
hamm = Table[1/2 (1.0-Cos[(2\[Pi] j)/d]),{j,0,d-1}];
xf = transform[x[[1;;d]] hamm];
\[Omega]s = Range[0,d/2] (2\[Pi])/(d dt);
\[Omega]s2 = -Reverse@Drop[Drop[\[Omega]s,1],{-1}];
{Join[\[Omega]s2,\[Omega]s],Join[xf[[d/2+2;;d]],xf[[1;;d/2+1]]]}\[Transpose]
]


(* ::Input::Initialization:: *)
Clear[PowerSpectrum];
PowerSpectrum[x_,dt_,averaging_:1,overlap_:0.3]:=Module[{d,xf,\[Omega]s,hamm,S,k,xs,Ss},
d = Length@x;
If[OddQ[d],d=d-1];

k=If[averaging==1,d,(2d)/(2averaging-1)//Floor]; 
If[OddQ[k],k=k-1];

(*hamm = Table[1/2(1.0-Cos[(2\[Pi] j)/k]),{j,0,k-1}]; *)

\[Omega]s = Range[1,k/2] (2\[Pi])/(k dt); 
\[Omega]s = Join[-Reverse@Drop[\[Omega]s,{-1}],\[Omega]s];

xs  = Partition[x,k,Floor[overlap k]];

xf = Table[Fourier[(*hamm*)(xx[[1;;k]]-Mean[xx])],{xx,xs}]; 
S = Mean@Table[dt Abs[xxf]^2,{xxf,xf}];      

{\[Omega]s,Join[S[[k/2+2;;k]],S[[2;;k/2+1]]]}\[Transpose]  
]


(* ::Input::Initialization:: *)
Clear[butter];
butter[data_,\[Omega]c_,dt_:1,order_:2]:=Module[{f,b},
f=RecurrenceFilter[ToDiscreteTimeModel[ButterworthFilterModel[{"Lowpass",order,\[Omega]c }],dt],data];
b=RecurrenceFilter[ToDiscreteTimeModel[ButterworthFilterModel[{"Lowpass",order,\[Omega]c }],dt],Reverse[data]];
(f+Reverse[b])/2
]


(* ::Input::Initialization:: *)
Clear[DistributionFromData];
DistributionFromData[data_,\[CapitalDelta]\[CapitalDelta]x_:Automatic]:=Module[{\[CapitalDelta]x,max,min,d,x,p},
min = Min@data;
max = Max@data;
d = Length@data;
\[CapitalDelta]x=Which[
NumberQ[\[CapitalDelta]\[CapitalDelta]x],\[CapitalDelta]\[CapitalDelta]x,
matchString["discrete",\[CapitalDelta]\[CapitalDelta]x] &&Not[ Head[data[[1]]]===Real],First@DeleteCases[Sort@Differences@Sort@Y,0],
True,10 (max-min)/d];

x = N@Range[min,max,\[CapitalDelta]x];
p=BinCounts[data, {min,max+\[CapitalDelta]x,\[CapitalDelta]x}]/d  //N;    

{x,p}\[Transpose]
]


(* ::Input::Initialization:: *)
Clear[TightBindingHamiltonian];
TightBindingHamiltonian[L_,V_:1,J_:1]:=J StiffnessMatrix[L] + Which[
ArrayQ[V],DiagonalMatrix@SparseArray[Flatten[Table[V,{Ceiling[L/Length@V]}]][[1;;L]]],
True,toeplitzMatrix[V,1,L]
];


(* ::Input::Initialization:: *)
Clear[TightBindingHamiltonianPBC];
TightBindingHamiltonianPBC[L_,V_:1,J_:1]:=J CirculantMatrix[L] + Which[
ArrayQ[V],DiagonalMatrix@SparseArray[Flatten[Table[V,{Ceiling[L/Length@V]}]][[1;;L]]],
True,toeplitzMatrix[V,1,L]
];


(* ::Input::Initialization:: *)
Clear[MajoranaHamiltonian];
MajoranaHamiltonian[h_,G_:0,ord_:"global"]:=
If[ord==="local",If[MatrixQ[h],-kron[h,\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"0", 
RowBox[{"-", "I"}]},
{"I", "0"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\)],0] +If[MatrixQ[G],- I kron[(G-G\[ConjugateTranspose])/2,\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"0", "1"},
{"1", "0"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\)],0] +If[MatrixQ[G], kron[(G+G\[ConjugateTranspose])/2,\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"1", "0"},
{"0", 
RowBox[{"-", "1"}]}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\)],0] ,
If[MatrixQ[h],-kron[\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"0", 
RowBox[{"-", "I"}]},
{"I", "0"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\),h],0] +If[MatrixQ[G],- I kron[{{0,1},{1,0}},(G-G\[ConjugateTranspose])/2],0] +If[MatrixQ[G], kron[{{1,0},{0,-1}},(G+G\[ConjugateTranspose])/2],0]  ] 

Clear[MajoranaTightBinding];
MajoranaTightBinding[V_,L_,ord_:"global"]:=MajoranaHamiltonian[TightBindingHamiltonian[L,V],0,ord]

Clear[MajoranaXY];
MajoranaXY[L_,h0_,\[Kappa]_]:=Module[{\[CapitalDelta]},
\[CapitalDelta]=h0 Eye[L]+ SparseArray[{Band[{1,2}]->(1-\[Kappa])/2(*Jy*),Band[{2,1}]->(1+\[Kappa])/2(*Jx*)},{L,L}] ;
I ArrayFlatten[({
 {0, -\[CapitalDelta]},
 {\[CapitalDelta]\[ConjugateTranspose], 0}
})]]


(* ::Input::Initialization:: *)
(* Subscript[\[Lambda], c] = 2 *)
Clear[aubryAndreHarper];
aubryAndreHarper[L_,\[Phi]_:0,\[Alpha]_:0]:=Table[Cos[ 2.0\[Pi] GoldenRatio n + \[Phi]]/(1-\[Alpha] Cos[2\[Pi] GoldenRatio n + \[Phi]]),{n,1,L}];


(* ::Input::Initialization:: *)
Clear[fibonacciString];
fibonacciString[L_]:=Table[IntegerPart[(n+1) 1/GoldenRatio^2]-IntegerPart[n 1/GoldenRatio^2],{n,1,L}];

fibonacciString[L_,a_,b_]:=fibonacciString[L]/.{0->a,1->b}


(* ::Input::Initialization:: *)
sshHamiltonian[L_,\[Tau]0_,\[Delta]\[Tau]_,\[Epsilon]_:1]:= Module[{hop = Table[\[Tau]0-(-1)^j \[Delta]\[Tau],{j,1,L-1}]},
- \[Epsilon] DiagonalMatrix[ConstantArray[\[Epsilon],L]]  + DiagonalMatrix[hop,1] + DiagonalMatrix[hop,-1]//SparseArray
]


(* ::Input::Initialization:: *)
Clear[UnitaryDynamics]
UnitaryDynamics[H_?ListQ,dt_,\[Tau]_,y0_:"id",onlyLastOutput_:False]:=Module[{ntrotter,\[Psi]0,f,\[Psi],results},
ntrotter = Round[\[Tau]/dt];
\[Psi]0 = If[y0==="id",Eye[Length@H[[0]]],y0];
\[Psi] = \[Psi]0;
If[onlyLastOutput===False,
Transpose[{Range[0,(Length[H]-1) dt,dt],Table[\[Psi]=MatrixExp[-I h dt,\[Psi]],{h,H}]}],
Do[\[Psi]=MatrixExp[-I h dt,\[Psi]],{h,H}]; \[Psi]
]
]

UnitaryDynamics[H_,dt_,\[Tau]_,y0_:"id",onlyLastOutput_:False]:=Module[{ntrotter,\[Psi]0,f,\[Psi]},
ntrotter = Round[\[Tau]/dt];
\[Psi]0 = If[y0==="id",Eye[Length@H[0]],y0];
If[y0==="id",
f[{t_,\[Psi]_}]:={t+dt,MatrixExp[-I dt H[t]] . \[Psi]},
f[{t_,\[Psi]_}]:={t+dt,MatrixExp[-I dt H[t],\[Psi]]}];

If[onlyLastOutput===False,
NestList[f,{0,\[Psi]0},ntrotter],
Nest[f,{0,\[Psi]0},ntrotter][[2]]
]
]


(* ::Input::Initialization:: *)
Clear[CounterDiabaticHamiltonian];
CounterDiabaticHamiltonian[Hs_,tf_,dt_]:=Module[{Hbare,T,Vs,HCD,hcd,V},

Hbare = Table[Hs[t],{t,0,tf,dt}];
T = Length@Hbare;

Vs=Normal@Table[
V=Eigvecs[Hbare[[n]]];
(V\[Transpose] Exp[I Arg[V[[1]]]])\[Transpose]  (* fixes so that eigenvectors always have the same phases *)
,{n,T}];

HCD =  Table[
hcd=I ((Vs[[n+1]]-Vs[[n-1]])/(2dt)) . Vs[[n]]\[ConjugateTranspose]; 
(hcd+hcd\[ConjugateTranspose])/2,{n,2,T-1}];

HCD=Join[{0 Hbare[[1]]},HCD,{0 Hbare[[1]]}];
{Hbare,HCD}

]


(* ::Input::Initialization:: *)
Clear[PTr]
PTr[\[Rho]_,list_,locdimlist_]:=Module[{i,eye,v,L,n,tups,tab,\[ScriptCapitalN]},
n = Length@locdimlist;
Do[eye[l] = Eye[locdimlist[[l]]],{l,n}];

Do[v[l,i] = {eye[l][[i]]}\[Transpose],{l,n},{i,locdimlist[[l]]}];
tups = Tuples@Table[Range[1,l],{l,locdimlist[[list]]}];
tab = Table[
kron@@Table[If[MemberQ[list,l],v[l,j[[Position[list,l][[1,1]]]]],eye[l]],{l,n}],{j,tups}];
Total@Table[(tab[[i]])\[Transpose] . \[Rho] . tab[[i]],{i,Length@tab}]
];

PTr[\[Rho]_,list_,locdim_?NumberQ]:= PTr[\[Rho],list,ConstantArray[locdim,Log[locdim, Length@\[Rho]]]]

PTr[\[Rho]_,list_]:=PTr[\[Rho],list,2];

PTr::usage = "Computes the partial trace of \[Rho] over a list of subsystems";


(* ::Input::Initialization:: *)
Clear[vonNeumann];
vonNeumann[\[Rho]_,base_:E]:= Module[{eigs,func},
eigs = Eigenvalues[\[Rho]];
func[x_]:= If[NumberQ[x]&&Re@x<10^-12,0,-x Log[base,x]];
Chop@Sum[func[eigs[[i]]],{i,1,Length@eigs}]
]

vonNeumann::usage = "von Neumann entropy S(\[Rho])";


(* ::Input::Initialization:: *)
Renyi[\[Rho]_,\[Alpha]_,base_:E]:= Module[{eigs,func},
Which[
\[Alpha] ==1.0,vonNeumann[\[Rho]],
\[Alpha] ==\[Infinity],-Log[base,Max@Eigenvalues[\[Rho]]],
True,(eigs = Eigenvalues[\[Rho]];
1/(1-\[Alpha]) Log[base,Sum[(eigs[[i]])^\[Alpha],{i,1,Length@eigs}]])]//Chop
]

Renyi::usage = "Renyi entropy \!\(\*SubscriptBox[\(S\), \(\[Alpha]\)]\)(\[Rho])";


(* ::Input::Initialization:: *)
(* D(\[Rho]||\[Sigma]) *)
Clear[KullbackLeibler];
KullbackLeibler[\[Rho]_?MatrixQ,\[Sigma]_?MatrixQ]:=Module[{S1,S2,\[Lambda],Q,log},

S1 = vonNeumann[\[Rho]];
{\[Lambda],Q} = Eigensystem[\[Sigma]];
Q = Q\[Transpose]; 
log = DiagonalMatrix@Log[\[Lambda]];
S2 = - Tr[Q\[ConjugateTranspose] . \[Rho] . Q . log];
-S1+S2//Chop
]

KullbackLeibler[p_?VectorQ,q_?VectorQ]:=Total@Table[If[Abs[p[[i]]]<10^-12,0,p[[i]] Log[p[[i]]/q[[i]]]],{i,Length@p}]


(* ::Input::Initialization:: *)
(* Subscript[D, \[Alpha]](\[Rho]||\[Sigma]) *)
Clear[RenyiDivergence];
RenyiDivergence[\[Rho]_?MatrixQ,\[Sigma]_?MatrixQ,\[Alpha]_]:=Module[{S1,S2,\[Lambda],Q,log,m1},

If[Abs[\[Alpha]-1]<10^-12,KullbackLeibler[\[Rho],\[Sigma]],
m1 = MatrixPower[\[Sigma],(1-\[Alpha])/(2\[Alpha])];
-(1/(1-\[Alpha]))Log@Tr[MatrixPower[m1 . \[Rho] . m1,\[Alpha]]]]//Chop
]



(* ::Input::Initialization:: *)
Clear[MutualInformation];
MutualInformation[\[Rho]_,{A_,B_},locdimlist_]:=Module[{\[Rho]A,\[Rho]B,\[Rho]AB,n},
n= Length@locdimlist;(* Number of sub-systems in this case *)
\[Rho]A = PTr[\[Rho],Complement[Range[n],Flatten@{A}],locdimlist];
\[Rho]B = PTr[\[Rho],Complement[Range[n],Flatten@{B}],locdimlist];
\[Rho]AB = PTr[\[Rho],Complement[Range[n],Flatten@Join[Flatten@{A},Flatten@{B}]],locdimlist];
Chop[vonNeumann[\[Rho]A]+vonNeumann[\[Rho]B]-vonNeumann[\[Rho]AB],10^-12]
];

(* When all subsystems have different dimensions *)
MutualInformation[\[Rho]_,{A_,B_},locdim_?NumberQ]:= MutualInformation[\[Rho],{A,B},ConstantArray[locdim,Log[locdim, Length@\[Rho]]]]

(* When all subsystems are qubits *)
MutualInformation[\[Rho]_,{A_,B_}]:= MutualInformation[\[Rho],{A,B},ConstantArray[2,Log[2, Length@\[Rho]]]]


(* ::Input::Initialization:: *)
Concurrence[\[Rho]_]:=Module[{\[Sigma]y,\[Rho]t,\[Rho]sq,R,eigs},
\[Sigma]y = {{0,-I},{I,0}};
\[Rho]t = kron[\[Sigma]y,\[Sigma]y] . \[Rho]\[Conjugate] . kron[\[Sigma]y,\[Sigma]y];
\[Rho]sq = MatrixPower[\[Rho],1/2];
R = MatrixPower[\[Rho]sq . \[Rho]t . \[Rho]sq,1/2];
eigs = Chop@Reverse@Sort@Eigenvalues[R];
Max[0,Re[eigs[[1]]-eigs[[2]]-eigs[[3]]-eigs[[4]]]]
];


(* ::Input::Initialization:: *)
EntanglementOfFormation[\[Rho]_,b_:2]:=Module[{c,h},
h[x_]:=If[(x==0)||(x==1),0,-x Log[b,x]-(1-x)Log[b,1-x]];
c = Concurrence[\[Rho]];
h[(1+Sqrt[1-c^2])/2]
]



(* ::Input::Initialization:: *)
Clear[QuantumFidelity];
QuantumFidelity[\[Rho]_,\[Sigma]_]:=Tr[MatrixPower[MatrixPower[\[Rho],1/2] . \[Sigma] . MatrixPower[\[Rho],1/2],1/2]]^2


(* ::Input::Initialization:: *)
Clear[RelativeEntropyofCoherence];
RelativeEntropyofCoherence[\[Rho]_,X_]:=Module[{v,V,\[CapitalDelta]\[Rho]},
V = Normalize/@Eigenvectors[X];
\[CapitalDelta]\[Rho] = Total@Table[out[v,v] . \[Rho] . out[v,v],{v,V}];
Chop[vonNeumann[\[CapitalDelta]\[Rho]]-vonNeumann[\[Rho]],10^-12]
];


(* ::Input::Initialization:: *)
Clear[IntegratedCOVy];
IntegratedCOVy[A_,B_,\[Rho]_,tol_:10^-11]:=Module[{f,d=Length@\[Rho],ps,\[ScriptCapitalO],At,Bt},
{ps,\[ScriptCapitalO]} = Eigensystem[N@\[Rho]];
\[ScriptCapitalO] = \[ScriptCapitalO]\[Transpose];
(*ps = Abs[ps];*)

f=(DiagonalMatrix[ps]+Outer[Subtract,(1+tol)ps,ps] )/(Eye[d]+tol+Outer[Subtract,Log[ps],Log[ps]]); 

At = \[ScriptCapitalO]\[ConjugateTranspose] . A . \[ScriptCapitalO];
Bt = \[ScriptCapitalO]\[ConjugateTranspose] . B . \[ScriptCapitalO];
1/2 (Total@Flatten[f Bt At\[Transpose]]- Tr[A . \[Rho]]Tr[B . \[Rho]])//Chop
]


(* ::Input::Initialization:: *)
Clear[WignerYanaseDyson];
(*WignerYanaseDyson[A_,\[Rho]_,tol_:10^-11]:=1/2(Tr[A.A.\[Rho]]-Tr[A.\[Rho]]^2)-IntegratedCOVy[A,A,\[Rho],tol]//Chop*)

WignerYanaseDyson[A_,\[Rho]_,tol_:10^-11]:=Module[{f,d=Length@\[Rho],ps,\[ScriptCapitalO],At},
{ps,\[ScriptCapitalO]} = Eigensystem[N@\[Rho]];
\[ScriptCapitalO] = \[ScriptCapitalO]\[Transpose];

f=(DiagonalMatrix[ps]+Outer[Subtract,(1+tol)ps,ps] )/(Eye[d]+tol+Outer[Subtract,Log[ps],Log[ps]]); 
At = \[ScriptCapitalO]\[ConjugateTranspose] . A . \[ScriptCapitalO];
1/2 (Tr[A . A . \[Rho]]-Total@Flatten[(f At) At\[Transpose]])//Chop]


(* ::Input::Initialization:: *)
Clear[QuantumStateTomography];
QuantumStateTomography[\[Rho]_,nexperiments_,listOfObservables_]:=Module[{d=Length@\[Rho],\[Lambda],V,p,r,\[Rho]est},
\[Rho]est=Total@Table[

{\[Lambda],V} = Eigsys[A];
V = V\[Transpose];
p = Table[v\[Conjugate] . \[Rho] . v,{v,V}]//Chop;
r=SillySample[\[Lambda],p,nexperiments]//Mean//N;
r/Tr[A\[ConjugateTranspose] . A] A

,{A,listOfObservables}];

\[Rho]est = 1/d (1-Tr[\[Rho]est])Eye[d] + \[Rho]est
]

QuantumStateTomography[\[Rho]_,nexperiments_]:=Module[{d=Length@\[Rho],listOfObservables},

listOfObservables = Join[
Table[Proj[d,i,i],{i,d}],
Flatten[Table[Proj[d,i,j]+Proj[d,j,i],{i,d},{j,i+1,d}],1],
Flatten[Table[I(Proj[d,i,j]-Proj[d,j,i]),{i,d},{j,i+1,d}],1]
];
QuantumStateTomography[\[Rho],nexperiments,listOfObservables]]


(* ::Input::Initialization:: *)
Clear[QuantumProcessTomography];
QuantumProcessTomography[\[ScriptCapitalE]_,d_,nexperiments_,obs_]:=Module[{\[Lambda],\[Beta],\[Chi],\[ScriptCapitalE]reconstructed,\[Rho]s},
\[Rho]s =  Join[
Table[Proj[d,i,i],{i,d}],
Flatten[Table[out[(Basis[d,i]+Basis[d,j])/Sqrt[2],(Basis[d,i]+Basis[d,j])/Sqrt[2]],{i,d},{j,i+1,d}],1],
Flatten[Table[out[(Basis[d,i]+I Basis[d,j])/Sqrt[2],(Basis[d,i]+I Basis[d,j])/Sqrt[2]],{i,d},{j,i+1,d}],1]
];

\[Beta] = ArrayReshape[N@Table[Tr[\[Rho]k\[ConjugateTranspose] . Em . \[Rho]j . En\[ConjugateTranspose]]/Tr[\[Rho]k\[ConjugateTranspose] . \[Rho]k],{\[Rho]j,\[Rho]s},{\[Rho]k,\[Rho]s},{Em,obs},{En,obs}],{d^4,d^4}];

\[Lambda] = If[nexperiments==0,
Table[Tr[\[Rho]k\[ConjugateTranspose] . \[ScriptCapitalE][\[Rho]j]]/Tr[\[Rho]k\[ConjugateTranspose] . \[Rho]k],{\[Rho]j,\[Rho]s},{\[Rho]k,\[Rho]s}],
Table[Tr[\[Rho]k\[ConjugateTranspose] . QuantumStateTomography[\[ScriptCapitalE][\[Rho]j],nexperiments]]/Tr[\[Rho]k\[ConjugateTranspose] . \[Rho]k],{\[Rho]j,\[Rho]s},{\[Rho]k,\[Rho]s}]
];
\[Chi] = Unvec@LinearSolve[\[Beta],Vec[\[Lambda]]];

\[ScriptCapitalE]reconstructed = Function[\[Rho],Sum[\[Chi][[m,n]]obs[[m]] . \[Rho] . obs[[n]]\[ConjugateTranspose],{m,Length@obs},{n,Length@obs}]];
{obs,\[Chi],\[ScriptCapitalE]reconstructed}
];

QuantumProcessTomography[\[ScriptCapitalE]_,d_]:=QuantumProcessTomography[\[ScriptCapitalE],d,0]

QuantumProcessTomography[\[ScriptCapitalE]_,d_,nexperiments_]:=Module[{obs},
obs = Join[
Table[Proj[d,i,i],{i,d}],
Flatten[Table[Proj[d,i,j]+Proj[d,j,i],{i,d},{j,i+1,d}],1],
Flatten[Table[I(Proj[d,i,j]-Proj[d,j,i]),{i,d},{j,i+1,d}],1]
];
QuantumProcessTomography[\[ScriptCapitalE],d,nexperiments,obs]]



(* ::Input::Initialization:: *)
Clear[Negentropy];
Negentropy[cumulants_,nsigma_:7]:=Module[{pgce,w,pg,\[Mu],\[Sigma]}, 
pgce = GramChalierEdgeworth[cumulants];
\[Mu] = cumulants[[1]];
\[Sigma] = Sqrt[cumulants[[2]]];
pg = PDF[NormalDistribution[\[Mu],\[Sigma]],w];

NIntegrate[Re[pgce[w] Log[pgce[w]/pg]],{w,\[Mu]-nsigma \[Sigma],\[Mu]+nsigma \[Sigma]},MaxRecursion->\[Infinity]]
]


(* ::Input::Initialization:: *)
(* MatrixPower but with A^0=identity, even if A is singular *)
matrixpower[A_,n_]:=If[n==0,Eye[Length@A],MatrixPower[A,n]]
matrixpower[A_,n_,y_]:=If[n==0,y,MatrixPower[A,n,y]]


(* ::Input::Initialization:: *)
Clear[WorkMoments];
WorkMoments[nn_,y_,Hi_,Hf_]:=Module[{n,Y,result},
(*If[checkCommute[Hi,y]==False,Print["Initial state does not commute with Hi"]];*)
result=Which[
y[[1]]==="thermal",
Y = MatrixExp[-y[[2]] Hi];
Y = Y/Tr[Y];

Table[Sum[(-1)^(n-k) Binomial[n,k] Tr[matrixpower[Hf,k] . matrixpower[Hi,n-k] . Y],{k,0,n}],{n,Flatten[{nn}]}],

y[[1]]=="GS",
Y = GroundState[Hi];
Re@Table[Sum[(-1)^(n-k) Binomial[n,k] Conjugate[matrixpower[Hf,k,Y]] . matrixpower[Hi,n-k,Y],{k,0,n}],{n,Flatten[{nn}]}],

y[[1]]=="eigenstate",
Y = GroundState[Hi,y[[2]]];
Re@Table[Sum[(-1)^(n-k) Binomial[n,k] Conjugate[matrixpower[Hf,k,Y]] . matrixpower[Hi,n-k,Y],{k,0,n}],{n,Flatten[{nn}]}],

Length@Dimensions@y==2,
Table[
Sum[(-1)^(n-k) Binomial[n,k] Tr[matrixpower[Hf,k] . matrixpower[Hi,n-k] . y],{k,0,n}],{n,Flatten[{nn}]}],

True,Table[Sum[(-1)^(n-k) Binomial[n,k] Conjugate[matrixpower[Hf,k,y]] . matrixpower[Hi,n-k,y],{k,0,n}],{n,Flatten[{nn}]}]
]//Chop;
If[Head[nn]===List,result,First@result]
]
WorkMoments[n_,y_,Hi_,Hf_,U_]:=WorkMoments[n,y,Hi,U\[ConjugateTranspose] . Hf . U]



(* ::Input::Initialization:: *)
getCumulants[moms_]:=Module[{nmoms = Length@moms,rep,n},
rep = Thread[Table[Moment[j],{j,nmoms}]->moms];
Table[MomentConvert[Cumulant[n],Moment]/.rep,{n,1,nmoms}]]


(* ::Input::Initialization:: *)
ProbMerger[x_,P_,tol_:10^-11]:=List@@@Normal[GroupBy[Table[{N[Round[x[[i]],tol]],P[[i]]},{i,1,Length[x]}],First-> Last,Total]]


(* ::Input::Initialization:: *)
(* Full work distribution *)
Clear[WorkDistribution];
WorkDistribution[Hi_,Hf_,\[Rho]0type_,y_,U_:1,mergetol_:10^-11]:=Module[{ord,i,Ei,Vi,Ef,Vf,\[ScriptCapitalM],\[ScriptCapitalQ],d = Length@Hi,pi,works,probs},
{Ei,Vi} =Eigsys[Hi]; 
{Ef,Vf} =If[U===1,Eigsys[Hf],Eigsys[U\[ConjugateTranspose] . Hf . U]];
\[ScriptCapitalQ] = Abs[Vf\[ConjugateTranspose] . Vi]^2 (* The matrix |\[LeftAngleBracket]j|U|i\[RightAngleBracket]|^2 *);

(* Initial probabilities *)
Which[
(* ground-state *)
\[Rho]0type =="GS",
pi=IdentityMatrix[d][[1]],

(* Thermal: in this case y = \[Beta] = 1/T *)
\[Rho]0type=="thermal",
(pi = Table[Exp[-y(Ei[[i]]-Ei[[1]])],{i,1,d}];
pi = pi/Total@pi;
),

(* Eigenstate of Subscript[H, i]. In this case y labels the eigenstate, with y = 1 being the GS *)
\[Rho]0type =="eigenstate",
pi = IdentityMatrix[d][[y]],

(* Infinite temperature state *)
\[Rho]0type =="identity",
pi = ConstantArray[1/d,{d}],

(* Generic. In this case y is the initial density matrix *)
\[Rho]0type=="generic",
pi = Table[Vi[[All,i]]\[Conjugate] . y . Vi[[All,i]],{i,d}]
];

works = Flatten@Table[ef-ei,{ei,Ei},{ef,Ef}];
probs = Flatten@Table[\[ScriptCapitalQ][[j,i]] pi[[i]],{i,d},{j,d}];
ord = Ordering[works];

Select[ProbMerger[works[[ord]],probs[[ord]],mergetol],Abs[#[[2]]]>10^-10&]
]

(*WorkDistribution[Hi_,Hf_,\[Rho]0type_,y_,U_]:=WorkDistribution[Hi,U\[ConjugateTranspose].Hf.U,\[Rho]0type,y]*)


(* ::Input::Initialization:: *)
(* Given probabilities and energies, yields the points (x,y) of a thermo-majorization curve *)
ThermoMajorizationCurve[pps_,Es_,\[Beta]_]:=Module[{ord,x,y,ps},
ps = pps/Total[pps];
ord=Reverse@Ordering[ps Exp[\[Beta] Es]];
y=Prepend[Accumulate[ps[[ord]]],0];
x = Prepend[Accumulate[Exp[-\[Beta] Es][[ord]]],0]//N;
Transpose[{x,y}]
]


(* ::Input::Initialization:: *)
Ergotropy[H_,\[Rho]_]:=Module[{P,\[ScriptCapitalE],psiRho,psiH},
{P,psiRho}=Eigensystem[N@\[Rho]];
{\[ScriptCapitalE],psiH}=Eigensystem[1000 NEye[Length@H]-H];
\[ScriptCapitalE] = 1000 - \[ScriptCapitalE];
Sum[P[[i]] \[ScriptCapitalE][[j]](Abs[psiH[[j]]\[Conjugate] . psiRho[[i]]]^2-KroneckerDelta[i,j]),{i,Length@\[Rho]},{j,Length@H}]
]


(* ::Input::Initialization:: *)
QuantumFisherInformation[\[Rho]_,d\[Rho]_]:=2Vec[d\[Rho]\[Transpose]] . LinearSolve[kron[\[Rho]\[Transpose],Eye[Length@\[Rho]]]+kron[Eye[Length@\[Rho]],\[Rho]],Vec[d\[Rho]]]//Chop(*2Vec[d\[Rho]\[Transpose]].PseudoInverse[kron[\[Rho]\[Transpose],NEye[Length@\[Rho]]]+kron[NEye[Length@\[Rho]],\[Rho]]].Vec[d\[Rho]]*)


(* ::Input::Initialization:: *)
SymmetricLogarithmicDerivative[\[Rho]_,d\[Rho]_]:=Module[{n,A,b},
n = Length@\[Rho];
A = kron[\[Rho]\[Transpose],Eye[n]] + kron[Eye[n],\[Rho]];
b = 2 Vec[d\[Rho]];
LinearSolve[A,b]//Unvec
]


(* ::Input::Initialization:: *)
Clear[GammelmarkMolmerMFisher];
GammelmarkMolmerMFisher[\[Rho]_,\[ScriptCapitalL]_, H_,Hp_,cops_,copsprime_]:=Module[{op1,\[DoubleStruckCapitalI],op2,\[ScriptCapitalL]L,\[ScriptCapitalL]R,\[ScriptCapitalL]inv,outie},
op1 = -I Hp - 1/2 Sum[copsprime[[k]]\[ConjugateTranspose] . cops[[k]] + cops[[k]]\[ConjugateTranspose] . copsprime[[k]],{k,Length@cops}];
op2 = I Hp - 1/2 Sum[copsprime[[k]]\[ConjugateTranspose] . cops[[k]] + cops[[k]]\[ConjugateTranspose] . copsprime[[k]],{k,Length@cops}];
\[DoubleStruckCapitalI] = Eye[Length@H];
\[ScriptCapitalL]L=kron[\[DoubleStruckCapitalI], op1] + Sum[kron[cops[[k]]\[Conjugate],copsprime[[k]]],{k,Length@cops}];
\[ScriptCapitalL]R=kron[op2\[Transpose],\[DoubleStruckCapitalI]] + Sum[kron[copsprime[[k]]\[Conjugate],cops[[k]]],{k,Length@cops}];


(*outie = Eye[Length@\[ScriptCapitalL]]  - out[Vec[\[Rho]], Vec[Eye[Length@\[Rho]]]];
\[ScriptCapitalL]inv = outie.PseudoInverse[\[ScriptCapitalL]].outie;
4Re@(Sum[Tr[c.\[Rho].c\[ConjugateTranspose]],{c,copsprime}]- UnTr[\[ScriptCapitalL]L.\[ScriptCapitalL]inv.\[ScriptCapitalL]R.Vec[\[Rho]]]-UnTr[\[ScriptCapitalL]R.\[ScriptCapitalL]inv.\[ScriptCapitalL]L.Vec[\[Rho]]])*)
4Re@(Sum[Tr[c . \[Rho] . c\[ConjugateTranspose]],{c,copsprime}]- UnTr[\[ScriptCapitalL]L . DrazinApply2[\[ScriptCapitalL],Vec[\[Rho]],\[ScriptCapitalL]R . Vec[\[Rho]]]]-UnTr[\[ScriptCapitalL]R . DrazinApply2[\[ScriptCapitalL],Vec[\[Rho]],\[ScriptCapitalL]L . Vec[\[Rho]]]])
]


(* ::Input::Initialization:: *)
Clear[SequentialMaxLikelihood];
SequentialMaxLikelihood[loglike_,outcomes_,\[Theta]grid_,outputstep_:1]:=Module[{ll,posML,Lmax,\[Theta]max},
ll = Accumulate@Table[loglike[x,\[Theta]],{x,outcomes},{\[Theta],\[Theta]grid}];
posML = First@Ordering[#,-1]&/@ll;
(*Lmax=Max/@ll;*)
\[Theta]max = \[Theta]grid[[posML]];
\[Theta]max[[1;;-1;;outputstep]]
]


(* ::Input::Initialization:: *)
Clear[MeanSquaredError];
MeanSquaredError[\[Theta]est_?MatrixQ,\[Theta]real_]:=Mean/@Transpose[(\[Theta]est-\[Theta]real)^2]
MeanSquaredError[\[Theta]est_?VectorQ,\[Theta]real_]:=MeanSquaredError[{\[Theta]est},\[Theta]real]


(* ::Input::Initialization:: *)
Clear[AmplitudeDampingKrausOperators];
AmplitudeDampingKrausOperators[f_,\[Gamma]_]:=Module[{M},
M[0] = Sqrt[f]({
 {1, 0},
 {0, Sqrt[1-\[Gamma]]}
});
M[1] = Sqrt[f]({
 {0, Sqrt[\[Gamma]]},
 {0, 0}
});
M[2] = Sqrt[1-f]({
 {Sqrt[1-\[Gamma]], 0},
 {0, 1}
});
M[3] = Sqrt[1-f]({
 {0, 0},
 {Sqrt[\[Gamma]], 0}
});
{M[0],M[1],M[2],M[3]}

]


(* ::Input::Initialization:: *)
Clear[AmplitudeDamping];
AmplitudeDamping[f_,\[Gamma]_][\[Rho]_]:=Module[{M = AmplitudeDampingKrausOperators[f,\[Gamma]]},Total@Table[M[[i]] . \[Rho] . M[[i]]\[Transpose],{i,Length@M}]]

AmplitudeDamping[f_,\[Gamma]_,qubit_][\[Rho]_]:=Module[{M = AmplitudeDampingKrausOperators[f,\[Gamma]]},
Total@Table[KronEye[M[[i]],Log2[Length@\[Rho]],qubit] . \[Rho] . KronEye[M[[i]],Log2[Length@\[Rho]],qubit]\[Transpose],{i,Length@M}]

]


(* ::Input::Initialization:: *)
Clear[SWAP]
SWAP[a_,b_,L_]:=Module[{\[DoubleStruckCapitalI],\[Sigma]x,\[Sigma]y,\[Sigma]z},

\[DoubleStruckCapitalI] = Eye[2^L];
\[Sigma]x = {{0,1},{1,0}};
\[Sigma]y = {{0,-I},{I,0}};
\[Sigma]z = {{1,0},{0,-1}};

1/2 (\[DoubleStruckCapitalI] + KronEye[\[Sigma]x,L,a] . KronEye[\[Sigma]x,L,b]+ KronEye[\[Sigma]y,L,a] . KronEye[\[Sigma]y,L,b]+ KronEye[\[Sigma]z,L,a] . KronEye[\[Sigma]z,L,b])
]

SWAP[]:=SWAP[1,2,2];


(* ::Input::Initialization:: *)
Clear[PartialSWAP]
PartialSWAP[s_,a_,b_,L_]:=( -I Sqrt[1-s^2] Eye[2^L] +s SWAP[a,b,L])


(* ::Input::Initialization:: *)
Clear[RealPartialSWAP];
RealPartialSWAP[s_,a_,b_,L_]:=Module[{\[ScriptCapitalX],\[Sigma]p,\[Sigma]m},
\[Sigma]p = {{0,1},{0,0}};
\[Sigma]m = {{0,0},{1,0}};
\[ScriptCapitalX] = KronEye[\[Sigma]p,L,a] . KronEye[\[Sigma]m,L,b]-KronEye[\[Sigma]m,L,a] . KronEye[\[Sigma]p,L,b];
Eye[2^L] + s \[ScriptCapitalX] +(1-Sqrt[1-s^2])\[ScriptCapitalX] . \[ScriptCapitalX]
]


(* ::Input::Initialization:: *)
Clear[CNOT]
CNOT[a_,b_,L_]:=Module[{\[Sigma]0,\[Sigma]x,\[Sigma]p,\[Sigma]m},

\[Sigma]0 = {{1,0},{0,1}};
\[Sigma]x = {{0,1},{1,0}};
\[Sigma]p = {{0,1},{0,0}};
\[Sigma]m = {{0,0},{1,0}};

(* Partial SWAP *)
KronEye[\[Sigma]p . \[Sigma]m,L,a] . KronEye[\[Sigma]0,L,b]+ KronEye[\[Sigma]m . \[Sigma]p,L,a] . KronEye[\[Sigma]x,L,b]
]

(*CNOT[] := Normal@CNOT[1,2,2];*)
CNOT[]:=\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"1", "0", "0", "0"},
{"0", "1", "0", "0"},
{"0", "0", "0", "1"},
{"0", "0", "1", "0"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\);

CNOT::usage = "CNOT[] calls the standard 4x4 CNOT gate. 

CNOT[a,b,L] loads a CNOT gate in a space of L qubits, with a (=1,...,L) being the control qubit and b (=1,...,L) the target qubit. 

CNOT[1,2,2] is equivalent to CNOT[].";


(* ::Input::Initialization:: *)
Clear[Hadamard]
Hadamard[]:=1/Sqrt[2] ({
 {1, 1},
 {1, -1}
});
Hadamard[i_,\[ScriptCapitalN]_]:=kron[Eye[2^(i-1)],\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{
FractionBox["1", 
SqrtBox["2"]], 
FractionBox["1", 
SqrtBox["2"]]},
{
FractionBox["1", 
SqrtBox["2"]], 
RowBox[{"-", 
FractionBox["1", 
SqrtBox["2"]]}]}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\),Eye[2^(\[ScriptCapitalN]-i)]]
Hadamard::usage = "Hadamard[] calls the 2x2 Hadamard matrix. 

Hadamard[i,\[ScriptCapitalN]] is a Hadamard gate on qubit i, in a space of \[ScriptCapitalN] qubits";


(* ::Input::Initialization:: *)
Clear[Toffoli]
Toffoli[a_,b_,c_,L_]:=Module[{\[Sigma]0,\[Sigma]x,\[Sigma]p,\[Sigma]m},

\[Sigma]0 = {{1,0},{0,1}};
\[Sigma]x = {{0,1},{1,0}};
\[Sigma]p = {{0,1},{0,0}};
\[Sigma]m = {{0,0},{1,0}};

(* Partial SWAP *)
KronEye[\[Sigma]p . \[Sigma]m,L,a] . KronEye[\[Sigma]m . \[Sigma]p,L,b] . KronEye[\[Sigma]0,L,c]+KronEye[\[Sigma]m . \[Sigma]p,L,a] . KronEye[\[Sigma]p . \[Sigma]m,L,b] . KronEye[\[Sigma]0,L,c]+KronEye[\[Sigma]p . \[Sigma]m,L,a] . KronEye[\[Sigma]p . \[Sigma]m,L,b] . KronEye[\[Sigma]0,L,c]+ KronEye[\[Sigma]m . \[Sigma]p,L,a] . KronEye[\[Sigma]m . \[Sigma]p,L,b] . KronEye[\[Sigma]x,L,c]
]


(* ::Input::Initialization:: *)
(* Dephases A in the eigenbasis of X *)
Clear[FullDephasingMap]
FullDephasingMap[A_,X_,degeneracy_:False]:=Module[{\[Lambda],v,kd},

{\[Lambda],v}=Eigsys[X];
v=v\[Transpose];
kd[a_,b_,tol_:10^-10] := If[Abs[a-b]<tol,1,0];
If[degeneracy===False,
Total@Table[(v[[i]]\[Conjugate] . A . v[[i]]) out[v[[i]],v[[i]]],{i,Length@v}],
Total@Flatten[Table[(v[[i]]\[Conjugate] . A . v[[j]]) out[v[[i]],v[[j]]] kd[\[Lambda][[i]],\[Lambda][[j]]],{i,Length@\[Lambda]},{j,Length@\[Lambda]}],1]]
]


(* ::Input::Initialization:: *)
(* Not vectorized *)
Clear[LindbladDissipator]
LindbladDissipator[A_,\[Rho]_]:=A . \[Rho] . A\[ConjugateTranspose]-1/2 (A\[ConjugateTranspose] . A . \[Rho]+\[Rho] . A\[ConjugateTranspose] . A);
LindbladDissipator[A_,B_,\[Rho]_]:=A . \[Rho] . B-1/2 (B . A . \[Rho]+\[Rho] . B . A);


(* ::Input::Initialization:: *)
Vec[X_]:=Flatten[X\[Transpose]];
Unvec[x_]:=With[{d = Sqrt[Length@x]},Partition[x,d]\[Transpose]];
UnTr[X_]:=Tr@Unvec@X


(* ::Input::Initialization:: *)
LindbladToVec[L_]:=With[{\[ScriptCapitalI] = Eye[Length@L]},kron[Conjugate[L],L]-1/2 kron[(L\[ConjugateTranspose] . L)\[Transpose],\[ScriptCapitalI]]-1/2 kron[\[ScriptCapitalI],L\[ConjugateTranspose] . L]]
UnitaryToVec[H_]:=With[{\[ScriptCapitalI] = Eye[Length@H]},-I (kron[\[ScriptCapitalI],H]- kron[H\[Transpose],\[ScriptCapitalI]])]

Liouvillian[H_,cops_]:=Module[{c},
If[Length[H]!=Length[cops[[1]]],0,UnitaryToVec[H]]+ Sum[LindbladToVec[c],{c,cops}]]

Liouvillian[H_,cops_,rates_]:=Module[{c},
If[Length[H]!=Length[cops[[1]]],0,UnitaryToVec[H]] + Sum[rates[[i]] LindbladToVec[cops[[i]]],{i,Length@cops}]]


(* ::Input::Initialization:: *)
Clear[JumpOp];
(*JumpOp[c_,type_:"PD"]:=If[type==="PD",
kron[c\[Conjugate],c],
kron[Eye[Length@c],c]+kron[c\[Conjugate],Eye[Length@c]]
]*)

JumpOp[c_,type_:"PD"]:=
Which[
matchString["PD",type],kron[c\[Conjugate],c],
matchString["Homodyne",type],kron[Eye[Length@c],c]+kron[c\[Conjugate],Eye[Length@c]],
True,kron[Eye[Length@c],c]+kron[c\[Conjugate],Eye[Length@c]]
]


(* ::Input::Initialization:: *)
Clear[EffectiveNonHermitianHamiltonian];
EffectiveNonHermitianHamiltonian[H_,cops_]:=Module[{c},
H -I/2 Sum[c\[ConjugateTranspose] . c,{c,cops}]]

EffectiveNonHermitianHamiltonian[H_,cops_,rates_]:=Module[{c},
H -I/2 Sum[rates[[i]] cops[[i]]\[ConjugateTranspose] . cops[[i]],{i,Length@cops}]]


(* ::Input::Initialization:: *)
Clear[VecDiagsOrdering];
VecDiagsOrdering[d_]:=Module[{diags=Table[(d+1)i - d,{i,1,d}]},
Join[diags,Complement[Range[1,d^2],diags]]];

Clear[VecDiags];
VecDiags[X_]:=Vec[X][[VecDiagsOrdering[Length@X]]]
VecDiags[X_,S_]:=Chop@VecDiags[S\[ConjugateTranspose] . X . S]

Clear[UnvecDiags];
UnvecDiags[X_]:=X[[Ordering@VecDiagsOrdering[Sqrt[Length@X]]]]//Unvec;
UnvecDiags[X_,S_]:=S . UnvecDiags[X] . S\[ConjugateTranspose]//Chop;

Clear[LindbladToVecDiags]
LindbladToVecDiags[L_]:=Module[{ord = VecDiagsOrdering[Length@L]},LindbladToVec[L][[ord,ord]]];
LindbladToVecDiags[L_,S_]:=Module[{ord = VecDiagsOrdering[Length@L]},LindbladToVec[S\[ConjugateTranspose] . L . S][[ord,ord]]];

Clear[UnitaryToVecDiags];
UnitaryToVecDiags[H_]:=Module[{ord = VecDiagsOrdering[Length@H]},UnitaryToVec[H][[ord,ord]]];
UnitaryToVecDiags[H_,S_]:=Module[{ord = VecDiagsOrdering[Length@H]},UnitaryToVec[S\[ConjugateTranspose] . H . S][[ord,ord]]];

Clear[LiouvillianDiags];
LiouvillianDiags[H_,cops_,rates_]:=Module[{ord = VecDiagsOrdering[Length@H]},Liouvillian[H,cops,rates][[ord,ord]]]
LiouvillianDiags[H_,cops_,rates_,S_]:=Module[{ord = VecDiagsOrdering[Length@H]},Liouvillian[S\[ConjugateTranspose] . H . S,S\[ConjugateTranspose] . # . S&/@cops,rates][[ord,ord]]]


(* ::Input::Initialization:: *)
AnalyticalSteadyStateDiags[\[ScriptCapitalL]_]:=Module[{\[ScriptCapitalL]2,d,\[Rho]},d=Length@\[ScriptCapitalL];
(*\[ScriptCapitalL]2=Join[\[ScriptCapitalL],{ConstantArray[1,d]}];*)\[ScriptCapitalL]2=Normal@Join[\[ScriptCapitalL],{Vec@Eye[Sqrt[d]]}];
\[Rho]=UnvecDiags@LinearSolve[\[ScriptCapitalL]2,Eye[d+1][[-1]]];
\[Rho]]
(*Specify the vectorized identity matrix \[DoubleStruckCapitalI].Output in this case is not unvectorized.*)
AnalyticalSteadyStateDiags[\[ScriptCapitalL]_,\[DoubleStruckCapitalI]_]:=Module[{\[ScriptCapitalL]2,d},d=Length@\[ScriptCapitalL];
\[ScriptCapitalL]2=Normal@Join[\[ScriptCapitalL],Normal@{\[DoubleStruckCapitalI]}];
LinearSolve[\[ScriptCapitalL]2,Eye[d+1][[-1]]]]


(* ::Input::Initialization:: *)
(* Based on Norikazu Kamiya, Prog. Theor. Exp. Phys. 2015, 043A02 *)
(* Assumes \[ScriptCapitalL] was built using LiouvillianDiags[H_,cops_,rates_,S_] *)
KamiyaEffectiveLiouvillian[\[ScriptCapitalL]_]:=Module[{d,\[ScriptCapitalL]pp,\[ScriptCapitalL]cp,\[ScriptCapitalL]pc,\[ScriptCapitalL]cc},
d = Sqrt[Length@\[ScriptCapitalL]];
\[ScriptCapitalL]pp = \[ScriptCapitalL][[1;;d,1;;d]];
\[ScriptCapitalL]pc = \[ScriptCapitalL][[1;;d,d+1;;d^2]];
\[ScriptCapitalL]cp = \[ScriptCapitalL][[d+1;;d^2,1;;d]];
\[ScriptCapitalL]cc = \[ScriptCapitalL][[d+1;;d^2,d+1;;d^2]];
\[ScriptCapitalL]pp - \[ScriptCapitalL]pc . Inverse[\[ScriptCapitalL]cc] . \[ScriptCapitalL]cp
]

KamiyaTransformedObservable[\[ScriptCapitalL]_,A_]:=Module[{d,\[ScriptCapitalL]cp,\[ScriptCapitalL]cc},
d = Sqrt[Length@\[ScriptCapitalL]];
\[ScriptCapitalL]cp = \[ScriptCapitalL][[d+1;;d^2,1;;d]];
\[ScriptCapitalL]cc = \[ScriptCapitalL][[d+1;;d^2,d+1;;d^2]];
({
 {Eye[d], -\[ScriptCapitalL]cp\[ConjugateTranspose] . Inverse[\[ScriptCapitalL]cc]},
 {0, Eye[d]}
}) . A

]


(* ::Input::Initialization:: *)
(* Only meant for numerics *)
Clear[SteadyState];
SteadyState[A_]:=Module[{\[Rho],x},
x=SparseArray@First@Eigenvectors[A,-1];
\[Rho] = Unvec[x] ;
\[Rho] = (\[Rho]+\[Rho]\[HermitianConjugate]);
\[Rho] = \[Rho]/Tr[\[Rho]]//Chop
];


(* ::Input::Initialization:: *)
Clear[AnalyticalSteadyState];
AnalyticalSteadyState[\[ScriptCapitalL]_]:=Module[{\[ScriptCapitalL]2,d,\[Rho]},
d = Length@\[ScriptCapitalL];
(*\[ScriptCapitalL]2 = Join[\[ScriptCapitalL],{ConstantArray[1,d]}];*)
\[ScriptCapitalL]2 = Normal@Join[\[ScriptCapitalL],{Vec@Eye[Sqrt[d]]}];
\[Rho] = Unvec@LinearSolve[\[ScriptCapitalL]2,Eye[d+1][[-1]]];
\[Rho]
]

(* Specify the vectorized identity matrix \[DoubleStruckCapitalI]. Output in this case is not unvectorized. *)
AnalyticalSteadyState[\[ScriptCapitalL]_,\[DoubleStruckCapitalI]_]:=Module[{\[ScriptCapitalL]2,d},
d = Length@\[ScriptCapitalL];
\[ScriptCapitalL]2 = Normal@Join[\[ScriptCapitalL],Normal@{\[DoubleStruckCapitalI]}];
LinearSolve[\[ScriptCapitalL]2,Eye[d+1][[-1]]]
]


(* ::Input::Initialization:: *)
(* Kind of buggy: use with care *)
Clear[AnalyticalSteadyStateSequential]
AnalyticalSteadyStateSequential[\[ScriptCapitalL]_]:=Module[{B,nn,x,i,p,ph,phh,\[Rho],tr},
Clear[B];
B = Normal@\[ScriptCapitalL];
nn=Length@\[ScriptCapitalL];
Do[
Do[
B[[i]] = Together[B[[i]]-B[[j]]B[[i,j]]/B[[j,j]]]
,{i,j+1,nn}];
PrintTemporary[j];
(*Print[Normal@B\[LeftDoubleBracket]j+1\[RightDoubleBracket]];*)
,{j,1,nn}];
x=ConstantArray[0,nn];
x[[nn]]=1;
Do[
x[[i]]=-Together@Total@Table[B[[i,j]]x[[j]],{j,i+1,nn}]/B[[i,i]];
(*Print[x\[LeftDoubleBracket]i\[RightDoubleBracket]];*)
PrintTemporary[i]
,{i,nn-1,1,-1}];

Do[x[[i]] = Simplify[x[[i]]];PrintTemporary[i];PrintTemporary[x[[i]]];,{i,Length@x}];
p = Unvec[x];
ph = ConstantArray[0,{Length@p,Length@p}];
phh = p\[HermitianConjugate];
Do[
ph[[i,j]]=cf[phh[[i,j]]];
PrintTemporary[{i,j}];
PrintTemporary[ph[[i,j]]];
,{i,Length@p},{j,Length@p}];
\[Rho] = p + ph;
tr = Tr[\[Rho]]//Simplify;
\[Rho] = \[Rho]/tr;
Simplify@\[Rho]
]


(* ::Input::Initialization:: *)
Clear[WaitingTimeDistribution];
WaitingTimeDistribution[t_,\[ScriptCapitalL]0_,\[ScriptCapitalJ]i_,\[ScriptCapitalJ]f_,\[Rho]0_:None]:=Module[{\[Rho],d=Sqrt[Length@\[ScriptCapitalL]0]},
\[Rho]=If[\[Rho]0===None,Eye[d]/d,\[Rho]0];
UnTr[\[ScriptCapitalJ]f . MatrixExp[\[ScriptCapitalL]0 t,\[ScriptCapitalJ]i . Vec[\[Rho]]]]/UnTr[\[ScriptCapitalJ]i . Vec[\[Rho]]]]   


(* ::Input::Initialization:: *)
Clear[WaitingTimeDistributionClassical]
WaitingTimeDistributionClassical[t_,\[ScriptCapitalL]0_,\[ScriptCapitalJ]i_,\[ScriptCapitalJ]f_,\[Rho]0_:None]:=Module[{\[Rho],d=Length@\[ScriptCapitalL]0},
\[Rho]=If[\[Rho]0===None,ConstantArray[1,d]/d,\[Rho]0];
Total[\[ScriptCapitalJ]f . MatrixExp[\[ScriptCapitalL]0 t,\[ScriptCapitalJ]i . \[Rho]]]/Total[\[ScriptCapitalJ]i . \[Rho]]]      


(* ::Input::Initialization:: *)
Clear[WaitingTimeMoments];
WaitingTimeMoments[\[Rho]_?MatrixQ,\[ScriptCapitalL]_,\[ScriptCapitalL]1_,n_,reset_:True]:=Module[{f,\[Alpha],\[Rho]0},
f = LinearSolve[\[ScriptCapitalL]-\[ScriptCapitalL]1];
\[Rho]0 = If[reset==True,\[ScriptCapitalL]1 . Vec[\[Rho]]/UnTr[\[ScriptCapitalL]1 . Vec[\[Rho]]],Vec[\[Rho]]];
\[Alpha] = Nest[f,\[Rho]0,n];
(-1)^n n! UnTr[\[Alpha]]  
];

WaitingTimeMoments[\[Rho]_?VectorQ,\[ScriptCapitalL]_,\[ScriptCapitalL]1_,n_,reset_:True]:=Module[{f,\[Alpha],\[Rho]0},
f = LinearSolve[\[ScriptCapitalL]-\[ScriptCapitalL]1];
\[Rho]0 = If[reset==True,\[ScriptCapitalL]1 . \[Rho]/Total[\[ScriptCapitalL]1 . \[Rho]],\[Rho]];
\[Alpha] = Nest[f,\[Rho]0,n];
(-1)^n n! Total[\[Alpha]]  
]


(* ::Input::Initialization:: *)
(* The actual no-jump distribution is computed as Subscript[P, no]=tr(\[Rho]f), where f is the output of "NoJumpDistributionFactor" *)
(* This is what the function NoJumpDistribution does *)

Clear[NoJumpDistributionFactor];
NoJumpDistributionFactor[ts_,\[ScriptCapitalL]0_]:=Module[{d = (Length[\[ScriptCapitalL]0])^(1/2)},
Table[UnTr@MatrixExp[\[ScriptCapitalL]0 t, kronV[Basis[d,j],Basis[d,i]]],{j,d},{i,d},{t,ts}]]

NoJumpDistributionFactor[ts_,H_,cops_]:=Module[{d = H,He,\[ScriptCapitalE]},
He = H - I/2 Total@Table[c\[ConjugateTranspose] . c,{c,cops}];
Transpose[Table[
\[ScriptCapitalE] = MatrixExp[-I He t];
\[ScriptCapitalE]\[ConjugateTranspose] . \[ScriptCapitalE],{t,ts}],{3,1,2}]
]


(* ::Input::Initialization:: *)
Clear[NoJumpDistribution];
NoJumpDistribution[\[Rho]_,f_]:=Vec[\[Rho]] . Flatten[f,1]//Chop

NoJumpDistribution[\[Rho]_,ts_,\[ScriptCapitalL]0_]:=With[{f = NoJumpDistributionFactor[ts,\[ScriptCapitalL]0]},NoJumpDistribution[\[Rho],f]];

NoJumpDistribution[\[Rho]_,ts_,H_,cops_]:=With[{f = NoJumpDistributionFactor[ts,H,cops]},NoJumpDistribution[\[Rho],f]];


(* ::Input::Initialization:: *)
Clear[EmissionSpectrum];
EmissionSpectrum[\[Rho]_,\[Omega]s_,L_,c_]:=Module[{one,J,\[Alpha],\[Beta],tab,tmp,\[DoubleStruckCapitalI]=Eye[Length@L]},

tab = Table[
\[Alpha] = LinearSolve[L-I \[Omega] \[DoubleStruckCapitalI], Vec[c . \[Rho]]];
(*\[Beta] = LinearSolve[L+\[ImaginaryI] \[Omega] Eye[Length@L], Vec[\[Rho].c\[ConjugateTranspose]]];*)
tmp= -Tr[c\[ConjugateTranspose] . Unvec[\[Alpha]]](*-Tr[c.Unvec[\[Beta]]]*);
tmp+tmp\[Conjugate]
,{\[Omega],\[Omega]s}];

tab//Chop
]


(* ::Input::Initialization:: *)
Clear[AbsorptionSpectrum];
AbsorptionSpectrum[\[Rho]_,\[Omega]s_,L_,c_]:=Module[{one,J,\[Alpha],\[Beta],tab,tmp},

tab = Table[
\[Alpha] = LinearSolve[L-I \[Omega] Eye[Length@L], Vec[\[Rho] . c]];
(*\[Beta] = LinearSolve[L+\[ImaginaryI] \[Omega] Eye[Length@L], Vec[c\[ConjugateTranspose].\[Rho]]];*)
tmp= -Tr[c\[ConjugateTranspose] . Unvec[\[Alpha]]](*-Tr[c.Unvec[\[Beta]]]*);
tmp+tmp\[Conjugate]
,{\[Omega],\[Omega]s}];

tab//Chop
]


(* ::Input::Initialization:: *)
Clear[CollisionMap];
CollisionMap[\[Rho]s0_,U_,\[Rho]a_,n_,output_:"system"]:=Module[{ds=Length@\[Rho]s0,da=Length@\[Rho]a},

Which[
output=="system",
NestList[PTr[U . kron[#,\[Rho]a] . U\[ConjugateTranspose],{2},{ds,da}]&,\[Rho]s0,n],
output=="all",
NestList[
With[{\[Rho]sa=U . kron[#[[1]],\[Rho]a] . U\[ConjugateTranspose]} , 
{PTr[\[Rho]sa,{2},{ds,da}],PTr[\[Rho]sa,{1},{ds,da}]}]&,{\[Rho]s0,\[Rho]a},n]
]
]


(* ::Input::Initialization:: *)
Clear[CollisionModelSteadyState]
CollisionModelSteadyState[U_,\[Rho]a_]:=Module[{da=Length@\[Rho]a,ds,\[Rho],\[ScriptCapitalT],p},
ds = Length[U]/da;
\[Rho] = Table[p[i,j],{i,ds},{j,ds}];
\[ScriptCapitalT]=Last@CoefficientArrays[Vec[PTr[U . kron[\[Rho],\[Rho]a] . U\[ConjugateTranspose],{2},{ds,da}]],Vec[\[Rho]]];
(*SteadyState[\[ScriptCapitalT]-Eye[ds^2]]*)
AnalyticalSteadyState[\[ScriptCapitalT]-Eye[ds^2]] (* "AnalyticalSteadyState handles multiple steady-states, by selecting one that is guaranteed physical *)
]


(* ::Input::Initialization:: *)
Clear[CoherentQuantumAbsorber];
CoherentQuantumAbsorber[U_,\[Chi]_]:=Module[{da=Length@\[Chi],ds,\[Rho]a,\[Rho]ss,\[Lambda],f,g,V,\[Psi]t,\[Phi],support,i,f\[Chi],error},
ds = Length[U]/da;
\[Rho]a = out[\[Chi],\[Chi]];
\[Rho]ss = CollisionModelSteadyState[U,\[Rho]a];
{\[Lambda],f} = Eigsys[\[Rho]ss];
\[Psi]t = Sum[Sqrt[\[Lambda][[i]]]kron[f[[i]],f[[i]]],{i,ds}];

\[Phi] = kron[U,Eye[ds]] . Sum[Sqrt[\[Lambda][[i]]]kron[f[[i]],\[Chi],f[[i]]],{i,ds}];

support = Select[{Range[1,ds],\[Lambda]}\[Transpose],Abs[#[[2]]]>10^-12&][[All,1]];
g=Table[1/Sqrt[\[Lambda][[i]]] kron[{f[[i]]\[Conjugate]},Eye[ds da]] . \[Phi],{i,support}];
g = Join[g,NullSpace[g\[Conjugate]]];

f = Join[Table[f[[i]],{i,support}],Table[f[[i]],{i,Complement[Range[1,ds],support]}]];

f\[Chi] = Table[kron[\[Chi],f[[i]]],{i,support}];
f\[Chi] = Join[f\[Chi],NullSpace[f\[Chi]\[Conjugate]]];

V = Sum[kron[Eye[ds],out[f\[Chi][[i]],g[[i]]]],{i,ds^2}];

error=V . kron[U,Eye[ds]] . Sum[Sqrt[\[Lambda][[i]]]kron[f[[i]],\[Chi],f[[i]]],{i,ds}]-Sum[Sqrt[\[Lambda][[i]]]kron[f[[i]],\[Chi],f[[i]]],{i,ds}]//Chop;
If[Total@Abs[error]>10^-12,Print["Apparently it did not work. The error was " <>ToString[error]]];
{V,\[Psi]t}
]


(* ::Input::Initialization:: *)
DrazinApply[\[ScriptCapitalL]_,y_]:=LinearSolve[Join[\[ScriptCapitalL],{Vec@Eye[Sqrt[Length@\[ScriptCapitalL]]]}], Join[y ,{0}]]

DrazinApply2[\[ScriptCapitalL]_,\[Rho]v_,y_]:=Module[{z},
z = LinearSolve[\[ScriptCapitalL], y-\[Rho]v UnTr[y]];
z - \[Rho]v UnTr[z]
]


(* ::Input::Initialization:: *)
Clear[FCSAverage];
FCSAverage[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,res},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]];
\[Rho]v = Vec[\[Rho]];
res=Table[Total@Table[\[Mu][[\[Alpha],i]] UnTr[L[[i]] . \[Rho]v],{i,d}],{\[Alpha],Length@\[Mu]}]//Chop;
If[Length@res==1,res[[1]],res]
]


(* ::Input::Initialization:: *)
Clear[FCSDiffusionMatrix];
FCSDiffusionMatrix[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
Jk = Table[UnTr[lop . \[Rho]v],{lop,L}];

M = Which[
type=="PD",\[Mu] . DiagonalMatrix[Jk] . \[Mu]\[Transpose],
type=="Homodyne"||True,\[Mu] . \[Mu]\[Transpose]
]//Chop;

\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha]//Chop;
\[ScriptCapitalL]big =Join[\[ScriptCapitalL],{Vec@Eye[Sqrt[Length@\[ScriptCapitalL]]]}];
ws = Table[LinearSolve[\[ScriptCapitalL]big, Join[\[ScriptCapitalL]\[Alpha][[i]] . \[Rho]v-J\[Alpha][[i]] \[Rho]v ,{0}]],{i,R}];
\[ScriptCapitalD] = -Table[UnTr[\[ScriptCapitalL]\[Alpha][[i]] . ws[[j]]],{i,R},{j,R}];
\[ScriptCapitalD]=\[ScriptCapitalD] + \[ScriptCapitalD]\[Transpose] + M//Chop;
If[R===1,\[ScriptCapitalD][[1,1]],\[ScriptCapitalD]]
]


(* ::Input::Initialization:: *)
Clear[FCSPowerSpectrumEigen]
FCSPowerSpectrumEigen[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
Jk = Table[UnTr[lop . \[Rho]v],{lop,L}];

M = Which[
type=="PD",\[Mu] . DiagonalMatrix[Jk] . \[Mu]\[Transpose],
type=="Homodyne"||True,\[Mu] . \[Mu]\[Transpose]
]//Chop;

\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha]//Chop;
{\[Lambda],Q} = Eigsys[\[ScriptCapitalL]];
\[Lambda] = Drop[\[Lambda],-1];
x = (Q\[Transpose])[[1;;-2]];
y=Inverse[Q][[1;;-2]];
(*Defined so that \[ScriptCapitalL] = Sum[\[Lambda]\[LeftDoubleBracket]i\[RightDoubleBracket] out[x\[LeftDoubleBracket]i\[RightDoubleBracket],y\[LeftDoubleBracket]i\[RightDoubleBracket]\[Conjugate]],{i,Length[\[ScriptCapitalL]]-1}] *)

\[CapitalUpsilon] = Table[UnTr[\[ScriptCapitalL]\[Alpha][[\[Alpha]]] . out[x[[j]],y[[j]]\[Conjugate]] . \[ScriptCapitalL]\[Alpha][[\[Beta]]] . \[Rho]v],{\[Alpha],R},{\[Beta],R},{j,Length[\[Lambda]]}];

If[R===1,
Function[M[[1,1]] - Chop@Total@Table[1/(\[Lambda][[j]]-I #) \[CapitalUpsilon][[1,1,j]]+1/(\[Lambda][[j]]+I #) \[CapitalUpsilon][[1,1,j]],{j,Length[\[Lambda]]}] ],
Function[M -Table[Chop@Total@Table[1/(\[Lambda][[j]]-I #) \[CapitalUpsilon][[\[Beta],\[Alpha],j]]+1/(\[Lambda][[j]]+I #) \[CapitalUpsilon][[\[Alpha],\[Beta],j]],{j,Length[\[Lambda]]}] ,{\[Alpha],R},{\[Beta],R}]]
]

]


(* ::Input::Initialization:: *)
Clear[FCSPowerSpectrumLinearSys]
FCSPowerSpectrumLinearSys[\[Omega]s_,\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS],v1,v2,\[DoubleStruckCapitalI]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
Jk = Table[UnTr[lop . \[Rho]v],{lop,L}];

M = Which[
type=="PD",\[Mu] . DiagonalMatrix[Jk] . \[Mu]\[Transpose],
type=="Homodyne"||True,\[Mu] . \[Mu]\[Transpose]
]//Chop;

\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
(*J\[Alpha] = UnTr[#.\[Rho]v]&/@\[ScriptCapitalL]\[Alpha]//Chop;*)

\[DoubleStruckCapitalI] = Eye[Length@\[ScriptCapitalL]];
Table[
v1=Table[LinearSolve[\[ScriptCapitalL] - I \[Omega] \[DoubleStruckCapitalI], \[ScriptCapitalL]\[Alpha][[i]] . \[Rho]v],{i,R}];
v2=Table[LinearSolve[\[ScriptCapitalL] + I \[Omega] \[DoubleStruckCapitalI], \[ScriptCapitalL]\[Alpha][[i]] . \[Rho]v],{i,R}];
\[ScriptCapitalS]=M - Table[UnTr[\[ScriptCapitalL]\[Alpha][[\[Beta]]] . v1[[\[Alpha]]]] + UnTr[\[ScriptCapitalL]\[Alpha][[\[Alpha]]] . v2[[\[Beta]]]],{\[Alpha],R},{\[Beta],R}];
{\[Omega],If[R===1,\[ScriptCapitalS][[1,1]],\[ScriptCapitalS]]}
,{\[Omega],\[Omega]s}]//Chop
]


(* ::Input::Initialization:: *)
Clear[FCSCoherence];
FCSCoherence[\[ScriptCapitalS]_,\[Omega]_]:=((Abs[\[ScriptCapitalS][\[Omega]][[1,2]]]^2)/(\[ScriptCapitalS][\[Omega]][[1,1]] \[ScriptCapitalS][\[Omega]][[2,2]])) 


(* ::Input::Initialization:: *)
(* Without dividing by J^2 *)
Clear[FCSG2MEXP]
FCSG2MEXP[ts_,\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha]//Chop;

If[R===1,
Table[{t,(*1/J\[Alpha]\[LeftDoubleBracket]1\[RightDoubleBracket]^2*) UnTr[\[ScriptCapitalL]\[Alpha][[1]] . MatrixExp[\[ScriptCapitalL] t,\[ScriptCapitalL]\[Alpha][[1]] . \[Rho]v]]},{t,ts}],
Table[{t,Table[(*1/(J\[Alpha]\[LeftDoubleBracket]\[Alpha]\[RightDoubleBracket]J\[Alpha]\[LeftDoubleBracket]\[Beta]\[RightDoubleBracket])*) UnTr[\[ScriptCapitalL]\[Alpha][[\[Beta]]] . MatrixExp[\[ScriptCapitalL] t,\[ScriptCapitalL]\[Alpha][[\[Alpha]]] . \[Rho]v]],{\[Alpha],R},{\[Beta],R}]},{t,ts}]
]//Re
]


(* ::Input::Initialization:: *)
Clear[FCSg2Eigen]
FCSg2Eigen[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];

\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha];
{\[Lambda],Q} = Eigsys[\[ScriptCapitalL]];
\[Lambda] = Drop[\[Lambda],-1];
x = (Q\[Transpose])[[1;;-2]];
y=Inverse[Q][[1;;-2]];
(*Defined so that \[ScriptCapitalL] = Sum[\[Lambda]\[LeftDoubleBracket]i\[RightDoubleBracket] out[x\[LeftDoubleBracket]i\[RightDoubleBracket],y\[LeftDoubleBracket]i\[RightDoubleBracket]\[Conjugate]],{i,Length[\[ScriptCapitalL]]-1}] *)

\[CapitalUpsilon] = Table[UnTr[\[ScriptCapitalL]\[Alpha][[\[Beta]]] . out[x[[j]],y[[j]]\[Conjugate]] . \[ScriptCapitalL]\[Alpha][[\[Alpha]]] . \[Rho]v],{\[Alpha],R},{\[Beta],R},{j,Length[\[Lambda]]}];

If[R===1,
Function[1+1/J\[Alpha][[1]]^2 Total@Table[Exp[\[Lambda][[j]] #] \[CapitalUpsilon][[1,1,j]],{j,Length[\[Lambda]]}] ],
Function[Table[1+1/(J\[Alpha][[\[Alpha]]]J\[Alpha][[\[Beta]]]) Total@Table[Exp[\[Lambda][[j]] #]\[CapitalUpsilon][[\[Beta],\[Alpha],j]],{j,Length[\[Lambda]]}] ,{\[Alpha],R},{\[Beta],R}]]
]//Chop
]


(* ::Input::Initialization:: *)
Clear[FCSg2MEXP]
FCSg2MEXP[ts_,\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha];

If[R===1,
Table[{t,1/J\[Alpha][[1]]^2 Tr[\[ScriptCapitalL]\[Alpha][[1]] . MatrixExp[\[ScriptCapitalL] t,\[ScriptCapitalL]\[Alpha][[1]] . \[Rho]v]]},{t,ts}], 
Table[{t,Table[1/(J\[Alpha][[\[Alpha]]]J\[Alpha][[\[Beta]]]) UnTr[\[ScriptCapitalL]\[Alpha][[\[Beta]]] . MatrixExp[\[ScriptCapitalL] t,\[ScriptCapitalL]\[Alpha][[\[Alpha]]] . \[Rho]v]],{\[Alpha],R},{\[Beta],R}]},{t,ts}]
]//Re

]


(* ::Input::Initialization:: *)
Clear[FCSTwoPointFunctionEigen]
FCSTwoPointFunctionEigen[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
(*Jk = Table[UnTr[lop.\[Rho]v],{lop,L}];*)
(*M = Which[
type=="PD",\[Mu].DiagonalMatrix[Jk].\[Mu]\[Transpose],
type=="Homodyne"||True,\[Mu].\[Mu]\[Transpose]
]//Chop;

*)
\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha]//Chop;
{\[Lambda],Q} = Eigsys[\[ScriptCapitalL]];
\[Lambda] = Drop[\[Lambda],-1];
x = (Q\[Transpose])[[1;;-2]];
y=Inverse[Q][[1;;-2]];
(*Defined so that \[ScriptCapitalL] = Sum[\[Lambda]\[LeftDoubleBracket]i\[RightDoubleBracket] out[x\[LeftDoubleBracket]i\[RightDoubleBracket],y\[LeftDoubleBracket]i\[RightDoubleBracket]\[Conjugate]],{i,Length[\[ScriptCapitalL]]-1}] *)

\[CapitalUpsilon] = Table[UnTr[\[ScriptCapitalL]\[Alpha][[\[Beta]]] . out[x[[j]],y[[j]]\[Conjugate]] . \[ScriptCapitalL]\[Alpha][[\[Alpha]]] . \[Rho]v],{\[Alpha],R},{\[Beta],R},{j,Length[\[Lambda]]}];

If[R===1,
Function[ Chop@Total@Table[Exp[\[Lambda][[j]] #] \[CapitalUpsilon][[1,1,j]],{j,Length[\[Lambda]]}] ],
Function[Table[Chop@Total@Table[Exp[\[Lambda][[j]] #]\[CapitalUpsilon][[\[Beta],\[Alpha],j]],{j,Length[\[Lambda]]}] ,{\[Alpha],R},{\[Beta],R}]]
]
]


(* ::Input::Initialization:: *)
Clear[FCSTwoPointFunctionMEXP]
FCSTwoPointFunctionMEXP[ts_,\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],\[Lambda],Q,x,y,\[CapitalUpsilon],\[ScriptCapitalS]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
(*Jk = Table[UnTr[lop.\[Rho]v],{lop,L}];*)
(*M = Which[
type=="PD",\[Mu].DiagonalMatrix[Jk].\[Mu]\[Transpose],
type=="Homodyne"||True,\[Mu].\[Mu]\[Transpose]
]//Chop;

*)
\[ScriptCapitalL]\[Alpha] = \[Mu] . L;
J\[Alpha] = UnTr[# . \[Rho]v]&/@\[ScriptCapitalL]\[Alpha]//Chop;
(*{\[Lambda],Q} = Eigsys[\[ScriptCapitalL]];
\[Lambda] = Drop[\[Lambda],-1];
x = (Q\[Transpose])\[LeftDoubleBracket]1;;-2\[RightDoubleBracket];
y=Inverse[Q]\[LeftDoubleBracket]1;;-2\[RightDoubleBracket];
*)(*Defined so that \[ScriptCapitalL] = Sum[\[Lambda]\[LeftDoubleBracket]i\[RightDoubleBracket] out[x\[LeftDoubleBracket]i\[RightDoubleBracket],y\[LeftDoubleBracket]i\[RightDoubleBracket]\[Conjugate]],{i,Length[\[ScriptCapitalL]]-1}] *)

(*\[CapitalUpsilon] = Table[UnTr[\[ScriptCapitalL]\[Alpha]\[LeftDoubleBracket]\[Beta]\[RightDoubleBracket].out[x\[LeftDoubleBracket]j\[RightDoubleBracket],y\[LeftDoubleBracket]j\[RightDoubleBracket]\[Conjugate]].\[ScriptCapitalL]\[Alpha]\[LeftDoubleBracket]\[Alpha]\[RightDoubleBracket].\[Rho]v],{\[Alpha],R},{\[Beta],R},{j,Length[\[Lambda]]}];*)

If[R===1,
Chop@Table[{t,UnTr[\[ScriptCapitalL]\[Alpha][[1]] . MatrixExp[\[ScriptCapitalL] t,\[ScriptCapitalL]\[Alpha][[1]] . \[Rho]v]]-J\[Alpha][[1]]^2},{t,ts}],
Chop@Table[{t,Table[UnTr[\[ScriptCapitalL]\[Alpha][[\[Beta]]] . MatrixExp[\[ScriptCapitalL] t,\[ScriptCapitalL]\[Alpha][[\[Alpha]]] . \[Rho]v]]-J\[Alpha][[\[Alpha]]]J\[Alpha][[\[Beta]]],{\[Alpha],R},{\[Beta],R}]},{t,ts}]
]
]


(* ::Input::Initialization:: *)
Clear[FCSProbability]
FCSProbability[t_,\[Rho]_,\[ScriptCapitalL]_,cops_,\[Mu]\[Mu]_:1,npts_:1000,tol_:10^-8]:=Module[{L,\[Chi]s,\[Delta]\[Chi],\[ScriptCapitalL]s,\[Rho]v,Ms,n,\[Mu],\[Chi],nmax,prob,P},
\[Mu] = If[\[Mu]\[Mu]===1,
ConstantArray[1,Length@cops],
\[Mu]\[Mu]];

\[Chi]s = linspace[-\[Pi],\[Pi],npts];
\[Delta]\[Chi] = \[Chi]s[[2]]-\[Chi]s[[1]];
\[ScriptCapitalL]s = Table[\[ScriptCapitalL] +Total@Table[(E^(I \[Chi] \[Mu][[i]]) -1)JumpOp[cops[[i]]],{i,Length@cops}] ,{\[Chi],\[Chi]s}];

\[Rho]v = Vec[\[Rho]];
Ms = Chop@Table[Tr[MatrixExp[L t,\[Rho]v]],{L,\[ScriptCapitalL]s}];

prob[nns_]:=Table[{n,1/(2\[Pi]) Trapz[E^(-I n \[Chi]s) Ms,\[Delta]\[Chi]]},{n,nns}]//Chop;

nmax = 10;
P = prob[Range[-nmax,nmax]];
While[P[[1,2]]+P[[-1,2]] > tol,
P=Join[prob[Range[-nmax-10,-nmax-1]],P,prob[Range[nmax+1,nmax+10]]];
nmax = nmax+10;
];
P
]


(* ::Input::Initialization:: *)
Clear[HasegawaSteadyStateTUR];
HasegawaSteadyStateTUR[\[Rho]_,\[ScriptCapitalL]_,H_,cops_]:=GammelmarkMolmerMFisher[\[Rho],\[ScriptCapitalL],H,H,cops,#/2&/@cops]


(* ::Input::Initialization:: *)
Clear[HasegawaTransientTUR]
HasegawaTransientTUR[t_,\[Rho]_,H_,cops_]:=Module[{\[CapitalXi],Heff,V0},
Heff = EffectiveNonHermitianHamiltonian[H,cops];
V0 = MatrixExp[-I Heff t];
\[CapitalXi] = Tr[Inverse[(V0\[ConjugateTranspose] . V0)] . \[Rho]]- 1
]


(* ::Input::Initialization:: *)
Clear[ModifiedHasegawa];
ModifiedHasegawa[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M0,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],M,\[ScriptCapitalH],\[Psi]},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
Jk = Table[UnTr[lop . \[Rho]v],{lop,L}];

M0 = Which[
type=="PD",Total@Jk,
type=="Homodyne"||True,\[Mu] . \[Mu]\[Transpose]
]//Chop;

\[ScriptCapitalH] = \[ScriptCapitalL]-Sum[LindbladToVec[lop],{lop,LL}];
ws = Table[DrazinApply[\[ScriptCapitalL],\[ScriptCapitalH] . \[Rho]v],{r,R}];

\[Psi] = Table[1/\[Mu][[r]] . Jk Sum[\[Mu][[r,k]] UnTr[L[[k]] . ws[[r]]],{k,d}],{r,R}][[1]];
{M0,\[Psi],1/M0 (1+\[Psi])^2} 
]


(* ::Input::Initialization:: *)
(*Clear[ModifiedHasegawaOptimalDeformation];
ModifiedHasegawaOptimalDeformation[\[Rho]_,\[ScriptCapitalL]_,LL_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD]\[Mu],\[ScriptCapitalL]1,x,\[ScriptCapitalL]J},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
Jk = Table[UnTr[lop.\[Rho]v],{lop,L}];

M = Which[
type=="PD",\[Mu].DiagonalMatrix[Jk].\[Mu]\[Transpose],
type=="Homodyne"||True,\[Mu].\[Mu]\[Transpose]
]//Chop;

\[ScriptCapitalL]J = Sum[\[Mu]\[LeftDoubleBracket]1,k\[RightDoubleBracket] L\[LeftDoubleBracket]k\[RightDoubleBracket],{k,d}];
ws = Table[DrazinApply[\[ScriptCapitalL],LindbladToVec[LL\[LeftDoubleBracket]q\[RightDoubleBracket]].\[Rho]v],{i,R},{q,d}];
Total@Table[(\[Mu]\[LeftDoubleBracket]1,k\[RightDoubleBracket] Jk\[LeftDoubleBracket]k\[RightDoubleBracket] - UnTr[\[ScriptCapitalL]J.ws\[LeftDoubleBracket]1,k\[RightDoubleBracket]])^2/Jk\[LeftDoubleBracket]k\[RightDoubleBracket],{k,d}]
];*)


(* ::Input::Initialization:: *)
(* Apparently wrong *)
Clear[ModifiedHasegawaUnitary];
ModifiedHasegawaUnitary[\[Rho]_,\[ScriptCapitalL]_,LL_,H_,\[Mu]\[Mu]_:1,type_:"PD"]:=Module[{L,\[Mu],\[Rho]v,d,Jk,M0,R,\[ScriptCapitalL]big,\[ScriptCapitalL]\[Alpha],J\[Alpha],ws,\[ScriptCapitalD],M,\[ScriptCapitalH],\[Psi],wn,wd,\[ScriptCapitalL]J},

L = If[Length@Dimensions[LL]==2,{JumpOp[LL,type]},JumpOp[#,type]&/@LL]; (* List of operators *)
d = Length@L;
\[Mu] = Which[
Dimensions[\[Mu]\[Mu]]==={d},{\[Mu]\[Mu]},
Dimensions[\[Mu]\[Mu]]==={},{\[Mu]\[Mu] ConstantArray[1,d]},
True,\[Mu]\[Mu]]; (* R\[Times]d matrix *)
R = Length@\[Mu];
\[Rho]v = Vec[\[Rho]];
Jk = Table[UnTr[lop . \[Rho]v],{lop,L}];

M0 = Which[
type=="PD",Total@Jk,
type=="Homodyne"||True,\[Mu] . \[Mu]\[Transpose]
]//Chop;

wn = DrazinApply[\[ScriptCapitalL],Vec[-I(H . \[Rho]-\[Rho] . H)]];
wd = DrazinApply2[\[ScriptCapitalL], \[Rho]v,Vec[H . \[Rho]+\[Rho] . H]];
\[ScriptCapitalL]J = Sum[\[Mu][[1,k]] L[[k]],{k,d}];

-(UnTr[\[ScriptCapitalL]J . wn]^2/Tr[H . Unvec[wd]])
]


(* ::Input::Initialization:: *)
Gillespie[W_,x0_,njumps_]:=Module[{\[Chi],r1,\[ScriptCapitalW],x,\[Lambda],\[Tau],tab},
\[Chi] = Range[1,Length@W];
r1 = RandomReal[{0,1},njumps];
\[ScriptCapitalW] = W - DiagonalMatrix@Diagonal@W; (* Make sure the matrix has zero diagonals *)

x = x0;
\[Tau] = 0;
tab=Table[
\[Lambda] = Total[\[ScriptCapitalW][[All,x]]];
\[Tau] =\[Tau] -1/\[Lambda] Log[1-r1[[n]]];
x = RandomChoice[Drop[\[ScriptCapitalW][[All,x]],{x}]/\[Lambda]->Drop[\[Chi],{x}]];

{\[Tau],x},{n,njumps}];

Prepend[tab,{0,x0}]
];

GillespieMoment[result_,n_]:=Differences[result[[All,1]]] . (result[[1;;-2,2]])^n/result[[-1,1]]


(* ::Input::Initialization:: *)
Clear[PauliEquation];
PauliEquation[pup_?VectorQ,pdw_?VectorQ]:=Module[{W,d=Length@pup+1},
W = DiagonalMatrix[pup,-1] + DiagonalMatrix[pdw,1] - DiagonalMatrix[Append[pup,0]]-DiagonalMatrix[Prepend[pdw,0]]
]

PauliEquation[d_,pup_,pdw_]:=PauliEquation[ConstantArray[pup,{d-1}],ConstantArray[pdw,{d-1}]]



(* ::Input::Initialization:: *)
Clear[QuantumJumpUnravelling];
QuantumJumpUnravelling[H_,cops_,\[Psi]0_,ts_,seed_:0]:=Module[{r,\[CapitalDelta]t,npts,ncops = Length@cops,He,\[Psi],\[ScriptCapitalU],states,jump,jumpList,dyna,n,ss},
If[seed!=0,SeedRandom[seed]];
\[CapitalDelta]t = ts[[2]]-ts[[1]];
npts=Length[ts]-1;
r = RandomReal[{0,1},npts];
He = H - I/2 Sum[c\[ConjugateTranspose] . c,{c,cops}];
\[Psi] = N@\[Psi]0;
\[ScriptCapitalU] = Eye[Length@He]-I He \[CapitalDelta]t-(He . He \[CapitalDelta]t^2)/2+1/6 I He . He . He \[CapitalDelta]t^3+(He . He . He . He \[CapitalDelta]t^4)/24; 
jumpList = {};
states = {};
n=1;
While[n<npts,
ss=NestWhileList[\[ScriptCapitalU] . #&,\[Psi],Norm[#]^2>r[[n]]&,1,npts-n+1];
AppendTo[states,Normalize/@ss];
\[Psi] = ss[[-1]];
n = n + Length@ss;
jump=If[ncops==1,1,RandomChoice[Table[Re[\[Psi]\[Conjugate] . c\[ConjugateTranspose] . c . \[Psi]],{c,cops}]->Range[1,ncops]]];
If[n<=npts,AppendTo[jumpList,{n (*\[CapitalDelta]t*), jump}]];
\[Psi] = cops[[jump]] . ss[[-1]];
\[Psi] = \[Psi]/Norm[\[Psi]];

{states}];
{Flatten[states,1],jumpList}

]


(* ::Input::Initialization:: *)
Clear[QuantumJumpGillespie];
QuantumJumpGillespie[H_,cops_,\[Psi]0_,ts_,njumps_]:=Module[{He,\[ScriptCapitalJ],Qs,V,\[Psi],\[Tau],Ps,norm,Vs,tr,\[Delta]\[Tau]r,weights,jump,cr,tab},

\[ScriptCapitalJ] = Sum[c\[ConjugateTranspose] . c,{c,cops}];
He = H - I/2 \[ScriptCapitalJ];

tr = Range[1,Length@ts];
cr = Range[1,Length@cops];
Vs = Table[MatrixExp[-I He t],{t,ts}];
Qs = Table[ V\[ConjugateTranspose] . \[ScriptCapitalJ] . V,{V,Vs}];
norm = Norm[Qs[[-1]]];

\[Psi] = \[Psi]0;
\[Tau]= 0;
tab=Table[
Ps = Table[\[Psi]\[Conjugate] . Q . \[Psi],{Q,Qs}]//Chop;
\[Delta]\[Tau]r = SillySample[tr,Ps];
\[Tau] = \[Tau] + ts[[\[Delta]\[Tau]r]];
\[Psi] = Vs[[\[Delta]\[Tau]r]] . \[Psi];
weights = Table[\[Psi]\[Conjugate] . c\[ConjugateTranspose] . c . \[Psi],{c,cops}]//Chop;
jump=RandomChoice[weights->cr];
\[Psi] = cops[[jump]] . \[Psi];
\[Psi] = \[Psi]/Norm[\[Psi]];
{\[Tau],jump(*,\[Psi]*)},{n,njumps}];

(*{tab,norm,He}*)
tab

]


(* ::Input::Initialization:: *)
(* Assumes it is a renewal process *)
Clear[QuantumJumpGillespieRenewal];
QuantumJumpGillespieRenewal[H_,cops_,\[Psi]0_,ts_,njumps_]:=Module[{He,\[ScriptCapitalJ],Qs,V,\[Psi],\[Tau],Ps,norm,Vs,tr,\[Delta]\[Tau]r,weights,jump,cr,tab,\[Delta]\[Tau]rs},

\[ScriptCapitalJ] = Sum[c\[ConjugateTranspose] . c,{c,cops}];
He = H - I/2 \[ScriptCapitalJ];

tr = Range[1,Length@ts];
cr = Range[1,Length@cops];
Vs = Table[MatrixExp[-I He t],{t,ts}];
Qs = Table[ V\[ConjugateTranspose] . \[ScriptCapitalJ] . V,{V,Vs}];
norm = Norm[Qs[[-1]]];

Do[
\[Psi] = cops[[c]] . RandomReal[{0,1},Length@\[Psi]0];
\[Psi] = \[Psi]/Norm[\[Psi]];
Ps[c] = Table[\[Psi]\[Conjugate] . Q . \[Psi],{Q,Qs}]//Chop;
\[Delta]\[Tau]rs[c] = SillySample[tr,Ps[c],njumps];
,{c,cr}];

\[Tau]= 0;
jump = 1;
\[Psi] = cops[[jump]] . RandomReal[{0,1},Length@\[Psi]0];
\[Psi] = \[Psi]/Norm[\[Psi]];
tab=Table[
(*\[Delta]\[Tau]r = SillySample[Ps[jump],tr];*)
\[Delta]\[Tau]r = \[Delta]\[Tau]rs[jump][[n]];
\[Tau] = \[Tau] + ts[[\[Delta]\[Tau]r]];
\[Psi] = Vs[[\[Delta]\[Tau]r]] . \[Psi];
weights = Table[\[Psi]\[Conjugate] . c\[ConjugateTranspose] . c . \[Psi],{c,cops}]//Chop;

jump=RandomChoice[weights->cr];
\[Psi] = cops[[jump]] . \[Psi];
\[Psi] = \[Psi]/Norm[\[Psi]];
{\[Tau],jump},{n,njumps}];

(*{tab,norm,He}*)
tab

]


(* ::Input::Initialization:: *)
Clear[HomodyneSimulate];
HomodyneSimulate[y0_,H_,ccops_,dt_,tf_,ssops_:None,seed_:0]:=Module[{n,sops,nsteps,S,\[Rho],\[Psi],dW,M,xave,dY,tab,\[DoubleStruckCapitalI],M0,x,c2,xs,c2s,dc,tr,cops,\[CapitalUpsilon]},
(* cops are monitored. sops are not *)
If[seed!=0,SeedRandom[seed]];
cops = If[Length@Dimensions[ccops]===2,{ccops},ccops];
dc = Length@cops;

sops = If[ssops===None,
{ConstantArray[0,{Length@H,Length@H}]},
If[Length@Dimensions[sops]===2,
{ssops},
ssops]];

nsteps = Round[tf/dt];
S= Association[];

S["dt"]=dt;
S["tf"]=tf;
S["nsteps"]=nsteps;
\[DoubleStruckCapitalI] = NEye[Length@H];
M0 = \[DoubleStruckCapitalI] - (I H + 1/2 Sum[c\[ConjugateTranspose] . c,{c,cops}] + 1/2 Sum[c\[ConjugateTranspose] . c,{c,sops}])dt;
dW = RandomVariate[NormalDistribution[0,Sqrt[dt]],{nsteps,dc}];
\[Rho] = \[Rho]0;
xs = Table[c + c\[ConjugateTranspose],{c,cops}];
c2s = Table[c1 . c2,{c1,cops},{c2,cops}];
tr = Sum[c . c,{c,cops}];

If[MatrixQ[y0]||Not[sops===None],
(\[Rho] = If[MatrixQ[y0],y0,out[y0,y0]];
tab=Table[
xave = Table[Tr[x . \[Rho]],{x,xs}];
dY = xave dt + dW[[n]]; 

\[CapitalUpsilon] = dY . cops;
M = M0 + \[CapitalUpsilon] + 1/2 (\[CapitalUpsilon] . \[CapitalUpsilon]-dt tr);
\[Rho] = M . \[Rho] . M\[ConjugateTranspose] + dt Sum[c . \[Rho] . c\[ConjugateTranspose],{c,sops}];
\[Rho] = \[Rho]/Tr[\[Rho]];
{\[Rho],dY,xave},{n,nsteps}];
),
(\[Psi] = y0;
tab=Table[
xave = Table[\[Psi]\[Conjugate] . x . \[Psi],{x,xs}];
dY = xave dt + dW[[n]]; 

\[CapitalUpsilon] = dY . cops;
M = M0 + \[CapitalUpsilon] + 1/2 (\[CapitalUpsilon] . \[CapitalUpsilon]-dt tr);
\[Psi] = M . \[Psi]//Normalize;
{\[Psi],dY,xave},{n,nsteps}]
)];

S["state"]=tab[[All,1]]//Chop;
S["J"]=If[Length@Dimensions@ccops==3,tab[[All,2]],Flatten[tab[[All,2]]]]/dt//Chop;
S["xave"]=If[Length@Dimensions@ccops==3,tab[[All,3]],Flatten[tab[[All,3]]]]//Chop;
S
]


(* ::Input::Initialization:: *)
Clear[LyapunovEigen]
LyapunovEigen[W_,F_]:=Module[{\[ScriptCapitalI],\[CapitalLambda],S,Sinv},
{\[CapitalLambda],S}=Quiet@Eigensystem[W];
S = S\[Transpose];
Sinv = Inverse[S];

\[ScriptCapitalI] = 1/Outer[Plus,\[CapitalLambda],\[CapitalLambda]\[Conjugate]] Sinv . F . Sinv\[ConjugateTranspose];
S . \[ScriptCapitalI] . S\[ConjugateTranspose]//Chop
];


(* ::Input::Initialization:: *)
LyapunovDriven[W_,F_,\[ScriptCapitalC]0_,\[CapitalDelta]t_,tf_,sparse_:1]:=Module[{f,testW,testF},
(* Check if input is a function or a matrix *)
testW = Head[W] ===Symbol || Head[W] ===Function;
testF = Head[F] ===Symbol || Head[F] ===Function;

Which[
testW&&testF,
f[t_,\[ScriptCapitalC]_]:=-W[t] . \[ScriptCapitalC]-\[ScriptCapitalC] . W[t]\[ConjugateTranspose]+F[t],
testW && Not[testF],
f[t_,\[ScriptCapitalC]_]:=-W[t] . \[ScriptCapitalC]-\[ScriptCapitalC] . W[t]\[ConjugateTranspose]+F,
Not[testW] && testF,
f[t_,\[ScriptCapitalC]_]:=-W . \[ScriptCapitalC]-\[ScriptCapitalC] . W\[ConjugateTranspose]+F[t],
True,
f[t_,\[ScriptCapitalC]_]:=-W . \[ScriptCapitalC]-\[ScriptCapitalC] . W\[ConjugateTranspose]+F
];
RK4t[f,\[ScriptCapitalC]0,\[CapitalDelta]t,tf][[1;;-1;;sparse]]
]



(* ::Input::Initialization:: *)
Clear[BDMatrices]
BDMatrices[h_,\[Gamma]_,f_]:=Module[{L=Length@h,\[CapitalGamma],F0},
\[CapitalGamma] =If[Length@\[Gamma]===Length@h, SparseArray@DiagonalMatrix@\[Gamma],SparseArray[{{1,1}->\[Gamma],{L,L}->\[Gamma]},{L,L}]];
F0 =Which[
Length@f===Length@h, SparseArray@DiagonalMatrix@f ,
Length@f === 2,SparseArray[{{1,1}-> f[[1]],{L,L}-> f[[2]]},{L,L}],
True,f Eye[L]
];
{I h + \[CapitalGamma]/2, \[CapitalGamma] F0}
]

BDMatrices[V_,L_,\[Gamma]_,f_]:=Module[{\[CapitalGamma],A,W,F,g,h,i}, 
BDMatrices[TightBindingHamiltonian[L,V], \[Gamma],f]
];



(* ::Input::Initialization:: *)
Clear[BDCovMat];
BDCovMat[W_,F_]:=LyapunovEigen[W,F]
BDCovMat[h_,\[Gamma]_,{f1_,fL_}]:=LyapunovEigen@@BDMatrices[h,\[Gamma],{f1,fL}];
BDCovMat[V_,L_,\[Gamma]_,{f1_,fL_}]:=LyapunovEigen@@BDMatrices[V,L,\[Gamma],{f1,fL}]


(* ::Input::Initialization:: *)
(* Only holds for nearest-neighbor tight-binding models *)
Clear[BDCurrent];
BDCurrent[\[ScriptCapitalC]_]:=(*Im@\[ScriptCapitalC]\[LeftDoubleBracket]2,1\[RightDoubleBracket]*) -I (\[ScriptCapitalC][[2,1]]-\[ScriptCapitalC][[1,2]])//Chop;
BDCurrent[W_,F_]:=BDCurrent[LyapunovEigen[W,F]]
BDCurrent[h_,\[Gamma]_,{f1_,fL_}]:=BDCurrent[LyapunovEigen@@BDMatrices[h,\[Gamma],{f1,fL}]];
BDCurrent[V_,L_,\[Gamma]_,{f1_,fL_}]:=BDCurrent[LyapunovEigen@@BDMatrices[V,L,\[Gamma],{f1,fL}]];


(* ::Input::Initialization:: *)
Clear[BDPairingMatrices];
BDPairingMatrices[\[ScriptCapitalT]_,\[Gamma]_,f_]:=Module[{L,W,\[Gamma]m,\[Gamma]p,\[CapitalUpsilon],\[ScriptCapitalW],\[ScriptCapitalF]},

L=Length[\[ScriptCapitalT]]/2;
{W,\[Gamma]p} = BDMatrices[Eye[L],\[Gamma],f];
\[Gamma]m = W+W\[ConjugateTranspose]-\[Gamma]p;

\[CapitalUpsilon] = 1/4 ArrayFlatten[({
 {\[Gamma]p+\[Gamma]m, I(\[Gamma]p-\[Gamma]m)},
 {-I (\[Gamma]p-\[Gamma]m), \[Gamma]p+\[Gamma]m}
})];
\[ScriptCapitalW] = I \[ScriptCapitalT] + \[CapitalUpsilon]+\[CapitalUpsilon]\[Transpose];
\[ScriptCapitalF] = - 2(\[CapitalUpsilon]-\[CapitalUpsilon]\[Transpose]);
{\[ScriptCapitalW],\[ScriptCapitalF]}
]

BDPairingMatrices[h_,G_,\[Gamma]_,f_]:=BDPairingMatrices[MajoranaHamiltonian[h,G],\[Gamma],f]

(* XY chain, with z-field h0 and anisotropy \[Kappa] *)
BDPairingMatrices[L_,h0_,\[Kappa]_,\[Gamma]_,f_]:=Module[{\[CapitalDelta]},
\[CapitalDelta]=h0 Eye[L]+ SparseArray[{Band[{1,2}]->(1-\[Kappa])/2(*Jy*),Band[{2,1}]->(1+\[Kappa])/2(*Jx*)},{L,L}] ;
 BDPairingMatrices[I ArrayFlatten[({
 {0, -\[CapitalDelta]},
 {\[CapitalDelta]\[ConjugateTranspose], 0}
})],\[Gamma],f]
]


(* ::Input::Initialization:: *)
Clear[BDPairingCovMat];
BDPairingCovMat[\[ScriptCapitalW]_,\[ScriptCapitalF]_]:=LyapunovEigen[\[ScriptCapitalW],\[ScriptCapitalF]]
BDPairingCovMat[\[ScriptCapitalT]_,\[Gamma]_,f_]:=LyapunovEigen@@BDPairingMatrices[\[ScriptCapitalT],\[Gamma],f];
BDPairingCovMat[h_,G_,\[Gamma]_,f_]:=LyapunovEigen@@BDPairingMatrices[h,G,\[Gamma],f];
BDPairingCovMat[L_,h0_,\[Kappa]_,\[Gamma]_,f_]:=LyapunovEigen@@BDPairingMatrices[L,h0,\[Kappa],\[Gamma],f];


(* ::Input::Initialization:: *)
BDPairingCurrent[\[CapitalTheta]_,\[Gamma]_,f_]:=Module[{\[Gamma]\[Gamma] = First@Flatten@{\[Gamma]},ff=First@Flatten@{f}},
{First@Flatten@{\[Gamma]} (First@Flatten@{f}-FermionsReducedC[\[CapitalTheta]][[1,1]]),
Last@Flatten@{\[Gamma]} (Last@Flatten@{f}-FermionsReducedC[\[CapitalTheta]][[-1,-1]])}//Chop
]


(* ::Input::Initialization:: *)
Clear[DephasingSolve];
DephasingSolve[W_,F_,\[CapitalGamma]_]:=Module[{L = Length@W,vecW,vecRemoveDiag,vecF},
vecW = kron[W\[Conjugate],NEye[L]]+kron[NEye[L],W];
vecRemoveDiag = SparseArray[{i_,i_}/;Last@QuotientRemainder[i,L+1]!=1->1,{L^2,L^2}];

LinearSolve[vecW + \[CapitalGamma] vecRemoveDiag, Vec[F]]//Unvec//Chop
]


(* ::Input::Initialization:: *)
Clear[BDCovMatDephasing];
BDCovMatDephasing[W_,F_,\[CapitalGamma]_]:=DephasingSolve[W,F,\[CapitalGamma]];
BDCovMatDephasing[h_,\[Gamma]_,{f1_,fL_},\[CapitalGamma]_]:=Module[{W,F} ,
{W,F} = BDMatrices[h,\[Gamma],{f1,fL}];
DephasingSolve[W,F,\[CapitalGamma]]
];
BDCovMatDephasing[V_,L_,\[Gamma]_,{f1_,fL_},\[CapitalGamma]_]:=Module[{W,F} ,
{W,F} = BDMatrices[V,L,\[Gamma],{f1,fL}];
DephasingSolve[W,F,\[CapitalGamma]]
];



(* ::Input::Initialization:: *)
(* Only holds for nearest-neighbor tight-binding models *)
Clear[BDCurrentDephasing];
BDCurrentDephasing[\[ScriptCapitalC]_]:=(*Im@\[ScriptCapitalC]\[LeftDoubleBracket]2,1\[RightDoubleBracket]*) -I (\[ScriptCapitalC][[2,1]]-\[ScriptCapitalC][[1,2]])//Chop;
BDCurrentDephasing[W_,F_,\[CapitalGamma]_]:=BDCurrentDephasing[BDCovMatDephasing[W,F,\[CapitalGamma]]]
BDCurrentDephasing[h_,\[Gamma]_,{f1_,fL_},\[CapitalGamma]_]:=BDCurrentDephasing[BDCovMatDephasing[h,\[Gamma],{f1,fL},\[CapitalGamma]]];
BDCurrentDephasing[V_,L_,\[Gamma]_,{f1_,fL_},\[CapitalGamma]_]:=BDCurrentDephasing[BDCovMatDephasing[V,L,\[Gamma],{f1,fL},\[CapitalGamma]]];


(* ::Input::Initialization:: *)
Clear[CovMatRhoFromC];
CovMatRhoFromC[\[ScriptCapitalC]_,statistic_:"Fermions"]:=Module[{\[DoubleStruckCapitalI] = Eye[Length@\[ScriptCapitalC]],\[Lambda],V,Z,M},
{\[Lambda],V} = Eigsys[\[ScriptCapitalC]];
If[matchString["Fermions", statistic],
Z = 1/Det[\[DoubleStruckCapitalI]-\[ScriptCapitalC]];
M=V . DiagonalMatrix[Log[(1-\[Lambda])/\[Lambda]]] . V\[ConjugateTranspose];
,
Z = Det[\[DoubleStruckCapitalI]+\[ScriptCapitalC]];
M=V . DiagonalMatrix[Log[(1+\[Lambda])/\[Lambda]]] . V\[ConjugateTranspose];
];
{M,Z}//Chop
]


(* ::Input::Initialization:: *)
Clear[CovMatEntropy];
CovMatEntropy[\[ScriptCapitalC]_,\[Alpha]_:1,statistic_:"Fermions"]:=Module[{\[DoubleStruckCapitalI] = Eye[Length@\[ScriptCapitalC]],Z,M,\[Lambda],V},
{\[Lambda],V} = Eigsys[\[ScriptCapitalC]];
If[matchString["Fermions", statistic],
Z = 1/Det[\[DoubleStruckCapitalI]-\[ScriptCapitalC]];
M=V . DiagonalMatrix[Log[(1-\[Lambda])/\[Lambda]]] . V\[ConjugateTranspose];
,
Z = Det[\[DoubleStruckCapitalI]+\[ScriptCapitalC]];
M=V . DiagonalMatrix[Log[(1+\[Lambda])/\[Lambda]]] . V\[ConjugateTranspose];
];
Chop@If[Abs[\[Alpha]-1]<10^-12,
Log[Z] + Tr[M . \[ScriptCapitalC]],
-\[Alpha]/(1-\[Alpha]) Log[Z] + 1/(1-\[Alpha]) If[matchString["Fermions", statistic],
Total@Log[1+((1-\[Lambda])/\[Lambda])^-\[Alpha]],-Total@Log[1-((1+\[Lambda])/\[Lambda])^-\[Alpha]]]]
]


(* ::Input::Initialization:: *)
Clear[CovMatRelativeEntropy];
CovMatRelativeEntropy[\[ScriptCapitalC]1_,\[ScriptCapitalC]2_,\[Alpha]_:1,statistic_:"Fermions"]:=Module[{\[DoubleStruckCapitalI] = Eye[Length@\[ScriptCapitalC]1],Z1,M1,\[Lambda]1,V1,Z2,M2,\[Lambda]2,V2,s,\[Zeta],\[Xi]},
{\[Lambda]1,V1} = Eigsys[\[ScriptCapitalC]1];
{\[Lambda]2,V2} = Eigsys[\[ScriptCapitalC]2];
If[matchString["Fermions", statistic],
Z1 = 1/Det[\[DoubleStruckCapitalI]-\[ScriptCapitalC]1];
M1=V1 . DiagonalMatrix[Log[(1-\[Lambda]1)/\[Lambda]1]] . V1\[ConjugateTranspose];
Z2 = 1/Det[\[DoubleStruckCapitalI]-\[ScriptCapitalC]2];
M2=V2 . DiagonalMatrix[Log[(1-\[Lambda]2)/\[Lambda]2]] . V2\[ConjugateTranspose]; 
,
Z1 = Det[\[DoubleStruckCapitalI]+\[ScriptCapitalC]1];
M1=V1 . DiagonalMatrix[Log[(1+\[Lambda]1)/\[Lambda]1]] . V1\[ConjugateTranspose];
Z2 = Det[\[DoubleStruckCapitalI]+\[ScriptCapitalC]2];
M2=V2 . DiagonalMatrix[Log[(1+\[Lambda]2)/\[Lambda]2]] . V2\[ConjugateTranspose]; 
];
s = If[matchString["Fermions", statistic],1,-1];

Chop@If[Abs[\[Alpha]-1]<10^-12,
-(Log[Z1] + Tr[M1 . \[ScriptCapitalC]1]) + Log[Z2] + Tr[M2 . \[ScriptCapitalC]1]
,
\[Zeta] = DiagonalMatrix[(\[Lambda]2/(1-s \[Lambda]2))^((1-\[Alpha])/(2\[Alpha]))] . V2\[ConjugateTranspose] . V1 . DiagonalMatrix[(\[Lambda]1/(1-s \[Lambda]1))] . V1\[ConjugateTranspose] . V2 . DiagonalMatrix[(\[Lambda]2/(1-s \[Lambda]2))^((1-\[Alpha])/(2\[Alpha]))]; 
\[Xi] = Eigvals[\[Zeta]];
Log[Z2]+\[Alpha]/(1-\[Alpha]) Log[Z1] - 1/(1-\[Alpha]) s Total@Log[1+s \[Xi]^\[Alpha]]
]
]


(* ::Input::Initialization:: *)
Clear[ReactionCoordinateMappingNumerics];
ReactionCoordinateMappingNumerics[Js_,\[Omega]s_]:=Module[{\[CapitalLambda]1,\[CapitalOmega]1,P,\[Omega],\[Omega]min,\[Omega]max,d\[Omega],Jnew},

\[Omega]min = First@\[Omega]s;
\[Omega]max = Last@\[Omega]s;
d\[Omega] = \[Omega]s[[2]]-\[Omega]s[[1]];

\[CapitalLambda]1 = Sqrt[1/(2\[Pi]) Trapz[Js,d\[Omega]]];
\[CapitalOmega]1 = 1/(2\[Pi] \[CapitalLambda]1^2) Trapz[\[Omega]s Js,d\[Omega]];

P = 1/\[Pi] Table[Trapz[Drop[Js,{i}]/Drop[\[Omega]s-\[Omega]s[[i]],{i}],d\[Omega]],{i,Length@\[Omega]s}]; 
Jnew = (4 \[CapitalLambda]1^2 Js)/(P^2+Js^2);
{\[CapitalLambda]1,\[CapitalOmega]1,Jnew}
]


(* ::Input::Initialization:: *)
StarToChainMapping[\[ScriptCapitalJ]s_,\[Omega]s_,Lb_]:=Module[{\[CapitalLambda]1,\[CapitalOmega]1,P,\[Omega],\[Omega]min,\[Omega]max,d\[Omega],Jnew,Js},

\[Omega]min = First@\[Omega]s;
\[Omega]max = Last@\[Omega]s;
d\[Omega] = \[Omega]s[[2]]-\[Omega]s[[1]];
Js=\[ScriptCapitalJ]s;
Table[
\[CapitalLambda]1 = Sqrt[1/(2\[Pi]) Trapz[Js,d\[Omega]]];
\[CapitalOmega]1 = 1/(2\[Pi] \[CapitalLambda]1^2) Trapz[\[Omega]s Js,d\[Omega]];

P = 1/\[Pi] Table[Trapz[Drop[Js,{i}]/Drop[\[Omega]s-\[Omega]s[[i]],{i}],d\[Omega]],{i,Length@\[Omega]s}]; 
Js = (4 \[CapitalLambda]1^2 Js)/(P^2+Js^2);
{\[CapitalLambda]1^2,\[CapitalOmega]1}
,{Lb}]
]


(* ::Input::Initialization:: *)
LeadSamples[L_,W_,Wstar_,\[CapitalGamma]_]:=Module[{Llog,Llin,logsamp,samples,energies,spacings,gamma,couplings},
Llog = Round[0.2 L];
Llin=L -2 Round[0.2 L];
logsamp = logspace[Log10[Wstar],Log10[W],Llog+1];
samples=Join[
Most@Reverse[-logsamp],
linspace[-Wstar,Wstar,Llin+1],
Rest@logsamp];

(*samples = linspace[-W,W,L+1];*)
energies = (Most@samples+Rest@samples)/2;
spacings = Differences[samples];
gamma = spacings;
couplings = Sqrt[\[CapitalGamma] spacings/(2\[Pi])];
{energies,couplings,gamma}
];

LeadSamplesLinear[L_,W_,Wstar_,\[CapitalGamma]_]:=Module[{Llog,Llin,logsamp,samples,energies,spacings,gamma,couplings},
samples = linspace[-W,W,L+1];
energies = (Most@samples+Rest@samples)/2;
spacings = Differences[samples];
gamma = spacings;
couplings = Sqrt[\[CapitalGamma] spacings/(2\[Pi])];
{energies,couplings,gamma}
];


(* ::Input::Initialization:: *)
Clear[BDThermLeadMatrices];
BDThermLeadMatrices[h_,\[CapitalGamma]_,{TL_,TR_,\[Mu]L_,\[Mu]R_},Lb_,w_,disc_:"linlog"]:=Module[{energies,couplings,gamma,Ls=Length@h,HB1,HB2,VS1,VS2,H,W,F,Fermis,ws=w/2},
Which[
disc==="linlog",
{energies,couplings,gamma}=LeadSamples[Lb,w,ws,\[CapitalGamma]],
disc==="lin",
{energies,couplings,gamma}=LeadSamplesLinear[Lb,w,ws,\[CapitalGamma]]
];

HB1 = DiagonalMatrix[energies];
HB2 = HB1;

VS1=Sum[couplings[[n]] Proj[{Lb,Ls},n,1],{n,Length@couplings}];
VS2=Sum[couplings[[n]] Proj[{Ls,Lb},Ls,n],{n,Length@couplings}];

H=ArrayFlatten[({
 {HB1, VS1, 0},
 {VS1\[ConjugateTranspose], h, VS2},
 {0, VS2\[ConjugateTranspose], HB2}
})]; 
Fermis[T_,\[Mu]_,energies_]:=1/(Exp[(energies-\[Mu])/T]+1) ;

W = I H + 1/2 DiagonalMatrix@Join[gamma,ConstantArray[0,Ls],gamma];

F = DiagonalMatrix@Join[Fermis[TL,\[Mu]L,energies]gamma,ConstantArray[0,Ls],Fermis[TR,\[Mu]R,energies]gamma];

{W,F}
]


(* ::Input::Initialization:: *)
Clear[ThermleadsGetSysCM];
ThermleadsGetSysCM[\[ScriptCapitalC]_,Lb_]:=Module[{Ls = Length@\[ScriptCapitalC]-2Lb},\[ScriptCapitalC][[Lb+1;;Lb+Ls,Lb+1;;Lb+Ls]]]

Clear[ThermleadsDecomposeW];
ThermleadsDecomposeW[W_,Lb_,bath_]:=Module[{\[CapitalLambda],\[Gamma],\[Epsilon],Ls = Length@W-2Lb,rg},
rg = If[bath==="left",1;;Lb,Lb+Ls+1;;Ls+2Lb];
\[Gamma] = Diagonal[(W+W\[ConjugateTranspose])][[rg]];
\[Epsilon] = Im[Diagonal@W][[rg]];
\[CapitalLambda] = If[bath=="left",Im@W[[Lb+1,rg]],Im@W[[Lb+Ls,rg]]];
{\[Gamma],\[Epsilon],\[CapitalLambda]}//Chop
]

(*Clear[ThermleadsDecomposeWF];
ThermleadsDecomposeWF[W_,F_,Lb_,bath_]:=Module[{\[Gamma],\[Epsilon],\[CapitalLambda],f},

{\[Gamma],\[Epsilon],\[CapitalLambda]}=ThermleadsDecomposeW[W,Lb,bath];
f = Which[
bath=="left",Diagonal[F]\[LeftDoubleBracket]1;;Lb\[RightDoubleBracket]/\[Gamma],
bath=="right",Diagonal[F]\[LeftDoubleBracket]Lb+Ls+1;;-1\[RightDoubleBracket]/\[Gamma]];
{\[Gamma],\[Epsilon],\[CapitalLambda],f}
]*)


(* ::Input::Initialization:: *)
Clear[ThermleadsParticleCurrent];
ThermleadsParticleCurrent[\[ScriptCapitalC]_,W_,F_,Lb_,bath_:"left"]:=Module[{Ls = Length@\[ScriptCapitalC]-2Lb,\[Gamma],rg},
rg = If[bath==="left",
1;;Lb,
Lb+Ls+1;;Lb+Ls+Lb];

\[Gamma] = (W+W\[ConjugateTranspose])[[rg,rg]];
Tr[F[[rg,rg]]-\[Gamma] . \[ScriptCapitalC][[rg,rg]]]//Chop
]

Clear[ThermleadsEnergyCurrent];
ThermleadsEnergyCurrent[\[ScriptCapitalC]_,W_,F_,Lb_,bath_:"left"]:=Module[{Ls = Length@\[ScriptCapitalC]-2Lb,\[Gamma],\[Epsilon],\[CapitalLambda],JEprx,JEartur,\[CapitalOmega],Jnaive,Jfcs,term1,term2,term3,rg,JEarturRight,JEprxRight},
rg = If[bath==="left",
1;;Lb,
Lb+Ls+1;;Lb+Ls+Lb];

\[Gamma] = (W+W\[ConjugateTranspose]);
\[Epsilon]=DiagonalMatrix[Im[Diagonal@W]];

If[bath==="left",
(\[CapitalLambda] = SparseArray[{{Lb+1,j_}/;1<=j<=Lb:>Im@W[[Lb+1,j]] \[Gamma][[j,j]]},{Lb+Ls+Lb,Lb+Ls+Lb}];
\[CapitalOmega] = SparseArray[{{Lb+1,j_}/;1<=j<=Lb:>Im@W[[Lb+1,j]] \[Epsilon][[j,j]]},{Lb+Ls+Lb,Lb+Ls+Lb}];
),(
\[CapitalLambda] = SparseArray[{{Lb+Ls,j_}/;rg[[1]]<=j<=rg[[2]]:>Im@W[[Lb+Ls,j]] \[Gamma][[j,j]]},{Lb+Ls+Lb,Lb+Ls+Lb}];
\[CapitalOmega] = SparseArray[{{Lb+Ls,j_}/;rg[[1]]<=j<=rg[[2]]:>Im@W[[Lb+Ls,j]] \[Epsilon][[j,j]]},{Lb+Ls+Lb,Lb+Ls+Lb}];
)];

\[CapitalLambda]=\[CapitalLambda]+\[CapitalLambda]\[ConjugateTranspose];
\[CapitalOmega]=\[CapitalOmega]-\[CapitalOmega]\[ConjugateTranspose];
(* Relevant terms *)
(* tr(Subscript[H, L]\[ScriptCapitalL](\[Rho])) *)
term1 = Tr[\[Epsilon][[rg,rg]] . (F[[rg,rg]]-\[Gamma][[rg,rg]] . \[ScriptCapitalC][[rg,rg]])]//Chop;

(* tr(Subscript[V, SL]\[ScriptCapitalL](\[Rho])) *)
term2 = -(1/2)Tr[\[CapitalLambda] . \[ScriptCapitalC]]//Chop;

(* i \[LeftAngleBracket][Subscript[H, L],Subscript[V, SL]]\[RightAngleBracket] *)
term3 = -I Tr[\[CapitalOmega] . \[ScriptCapitalC]]//Chop;


(* Energy current in Marlon's PRX *)
(* tr(Subscript[H, L]\[ScriptCapitalL](\[Rho])) + tr(Subscript[V, SL]\[ScriptCapitalL](\[Rho])) *)
JEprx = term1+term2;

(* Change in energy of the bath (Artur's paper) *)
(* Meant to be the "correct" one; equal to JEprx in steady-state. *)
(* i \[LeftAngleBracket][Subscript[H, L],Subscript[V, SL]]\[RightAngleBracket] + tr(Subscript[V, SL]\[ScriptCapitalL](\[Rho])) *)
JEartur = term3+term2;

{{JEprx,JEartur},{term1,term2,term3}}
]


(* ::Input::Initialization:: *)
BDFCSAverage[W_,F_,{\[Mu]\[Mu]p_,\[Mu]\[Mu]m_}]:=Module[{\[Gamma]p,\[Gamma]m,\[DoubleStruckCapitalI] = Eye[Length@W],\[ScriptCapitalC],\[Mu]m,\[Mu]p},
\[Gamma]p = F;
\[Gamma]m = W+W\[ConjugateTranspose]-\[Gamma]p;
\[Mu]m = DiagonalMatrix[\[Mu]\[Mu]m];
\[Mu]p = DiagonalMatrix[\[Mu]\[Mu]p];
\[ScriptCapitalC] = BDCovMat[W,F];

Tr[\[Mu]m . \[Gamma]m . \[ScriptCapitalC] + \[Mu]p . \[Gamma]p . (\[DoubleStruckCapitalI]-\[ScriptCapitalC])]//Chop
]


(* ::Input::Initialization:: *)
BDFCSDiffusionMatrix[W_,F_,{\[Mu]\[Mu]p_,\[Mu]\[Mu]m_}]:=Module[{\[Gamma]p,\[Gamma]m,\[DoubleStruckCapitalI] = Eye[Length@W],\[ScriptCapitalC],\[CapitalOmega],\[ScriptCapitalC]t,\[CapitalUpsilon],\[Mu]m,\[Mu]p},
\[Gamma]p = F;
\[Gamma]m = W+W\[ConjugateTranspose]-\[Gamma]p;
\[Mu]m = DiagonalMatrix[\[Mu]\[Mu]m];
\[Mu]p = DiagonalMatrix[\[Mu]\[Mu]p];
\[ScriptCapitalC] = BDCovMat[W,F];
\[ScriptCapitalC]t = BDCovMat[W, -(*\[ImaginaryI]*) \[ScriptCapitalC] . \[Mu]m . \[Gamma]m . \[ScriptCapitalC]+(*\[ImaginaryI]*) (\[DoubleStruckCapitalI]-\[ScriptCapitalC]) . \[Mu]p . \[Gamma]p . (\[DoubleStruckCapitalI]-\[ScriptCapitalC])];
Tr[(\[Mu]m . \[Mu]m . \[Gamma]m - \[Mu]p . \[Mu]p . \[Gamma]p) . \[ScriptCapitalC]+\[Mu]p . \[Mu]p . \[Gamma]p]+- 2 (*\[ImaginaryI]*)(-1) Tr[(\[Mu]m . \[Gamma]m-\[Mu]p . \[Gamma]p) . \[ScriptCapitalC]t]//Chop
]


(* ::Input::Initialization:: *)
Clear[BDEmissionSpectra];
BDEmissionSpectra[W_,F_,\[Omega]s_]:=Module[{\[ScriptCapitalC],\[DoubleStruckCapitalI],inv1,inv2},
\[ScriptCapitalC] = BDCovMat[W,F];
\[DoubleStruckCapitalI] = Eye[Length@W];
Table[
inv1 = Inverse[W\[ConjugateTranspose]+I \[Omega] \[DoubleStruckCapitalI] ];
inv2 = Inverse[W-I \[Omega] \[DoubleStruckCapitalI]];
\[ScriptCapitalC] . inv1+inv2 . \[ScriptCapitalC],{\[Omega],\[Omega]s}]//Chop
]


(* ::Input::Initialization:: *)
(* Fermions *) 
Clear[BDAbsorptionSpectra];
BDAbsorptionSpectra[W_,F_,\[Omega]s_]:=Module[{\[ScriptCapitalC],\[DoubleStruckCapitalI],inv1,inv2,\[ScriptCapitalC]b},
\[DoubleStruckCapitalI] = Eye[Length@W];
\[ScriptCapitalC] = BDCovMat[W,F];
\[ScriptCapitalC]b = \[DoubleStruckCapitalI]-\[ScriptCapitalC];

Table[
inv1 = Inverse[W\[ConjugateTranspose]+I \[Omega] \[DoubleStruckCapitalI] ];
inv2 = Inverse[W-I \[Omega] \[DoubleStruckCapitalI]];
\[ScriptCapitalC]b . inv1+inv2 . \[ScriptCapitalC]b,{\[Omega],\[Omega]s}]//Chop
]


(* ::Input::Initialization:: *)
Clear[ThermleadsFCSGMatrix];
ThermleadsFCSGMatrix[W_,Lb_,type_:"particle",bath_:"left"]:=Module[{\[Gamma],\[Epsilon],\[CapitalLambda],G,Ls = Length@W-2Lb},
{\[Gamma],\[Epsilon],\[CapitalLambda]}= ThermleadsDecomposeW[W,Lb,bath];
G=Which[
type=="particle"&&bath=="left",
SparseArray[{{Lb+1,j_}/;1<=j<=Lb:>\[CapitalLambda][[j]]},{Lb+Ls+Lb,Lb+Ls+Lb}],
type=="energy"&&bath=="left",
SparseArray[{{Lb+1,j_}/;1<=j<=Lb:>(\[CapitalLambda](\[Epsilon]-I \[Gamma]/2))[[j]]},{Lb+Ls+Lb,Lb+Ls+Lb}],
type=="particle"&&bath=="right",
SparseArray[{{Lb+Ls,j_}/;Lb+Ls+1<=j<=2Lb+Ls:>\[CapitalLambda][[j-Lb-Ls]]},{Lb+Ls+Lb,Lb+Ls+Lb}],
type=="energy"&&bath=="right",
SparseArray[{{Lb+Ls,j_}/;Lb+Ls+1<=j<=2Lb+Ls:>(\[CapitalLambda](\[Epsilon]-I \[Gamma]/2))[[j-Lb-Ls]]},{Lb+Ls+Lb,Lb+Ls+Lb}],
type=="particle"&&bath=="difference",
SparseArray[{
{Lb+1,j_}/;1<=j<=Lb:>\[CapitalLambda][[j]],
{Lb+Ls,j_}/;Lb+Ls+1<=j<=2Lb+Ls:>-\[CapitalLambda][[j-Lb-Ls]]
},{Lb+Ls+Lb,Lb+Ls+Lb}]
];
G=G-G\[ConjugateTranspose]
]


(* ::Input::Initialization:: *)
Clear[ThermleadsFCSCurrent];
ThermleadsFCSCurrent[\[ScriptCapitalC]_,W_,F_,Lb_,bath_:"left"]:=Module[{G,Ls = Length@W-2Lb},
G=ThermleadsFCSGMatrix[W,Lb,"particle",bath];
-I Tr[G . \[ScriptCapitalC]]//Chop
]


(* ::Input::Initialization:: *)
Clear[ThermleadsFCSNoise]
ThermleadsFCSNoise[\[ScriptCapitalC]_,W_,F_,Lb_,bath_:"left"]:=Module[{\[DoubleStruckCapitalI] = Eye[Length@W],\[ScriptCapitalC]t,Ls = Length@W-2Lb,G,\[Psi]},
\[DoubleStruckCapitalI] = Eye[Lb+Ls+Lb];
G=ThermleadsFCSGMatrix[W,Lb,"particle",bath];
\[ScriptCapitalC]t = BDCovMat[W,- (1/2)((\[DoubleStruckCapitalI]-\[ScriptCapitalC]) . G . \[ScriptCapitalC]+ \[ScriptCapitalC] . G . (\[DoubleStruckCapitalI]-\[ScriptCapitalC]))];
2Tr[G . \[ScriptCapitalC]t]//Chop
]


(* ::Input::Initialization:: *)
Clear[ThermleadsFCSCurrentEnergyResolved];
ThermleadsFCSCurrentEnergyResolved[\[ScriptCapitalC]_,W_,F_,Lb_,bath_:"left"]:=Module[{\[CapitalLambda],Ls = Length@W-2Lb},
\[CapitalLambda]= Last@ThermleadsDecomposeW[W,Lb,bath];
If[bath=="left",
Table[2\[CapitalLambda][[j]] Im@\[ScriptCapitalC][[j,Lb+1]],{j,Lb}],
Table[2\[CapitalLambda][[j]] Im@\[ScriptCapitalC][[Lb+Ls+j,Lb+Ls]],{j,Lb}]
]//Chop
]


(* ::Input::Initialization:: *)
(* \[ScriptCapitalT]=tr(Subscript[\[CapitalGamma], 1]Subscript[G\[CapitalGamma], L]G\[ConjugateTranspose]) = \[LeftAngleBracket]1|S1/(\[CapitalLambda]-i \[Epsilon])S^-1|L\[RightAngleBracket]\[LeftAngleBracket]L| S^-1\[ConjugateTranspose]1/(\[CapitalLambda]\[ConjugateTranspose]+i \[Epsilon])S\[ConjugateTranspose]|1\[RightAngleBracket] = |\!\(
\*SubscriptBox[\(\[Sum]\), \(k\)]\ 
\*FractionBox[
SubscriptBox[\(
\*SubscriptBox[\(S\), \(1 k\)](
\*SuperscriptBox[\(S\), \(-1\)])\), \(kL\)], \(
\*SubscriptBox[\(\[CapitalLambda]\), \(k\)] - i\ \[Epsilon]\)]\)(|^2) *)
NEGFTransmissionFunction[h_,\[Gamma]_]:=Module[{L=Length@h,\[CapitalGamma]1,\[CapitalGamma]L,W,\[CapitalLambda],S,Sinv},
\[CapitalGamma]1 = \[Gamma] Proj[L,1];
\[CapitalGamma]L = \[Gamma] Proj[L,L];
(*W = \[ImaginaryI] h + (\[CapitalGamma]1 + \[CapitalGamma]L)/2;*)
W = -(I/2)(\[CapitalGamma]1+\[CapitalGamma]L)+ h;
{\[CapitalLambda],S}=Quiet@Eigensystem[W];
S = S\[Transpose]//SparseArray;
\[CapitalLambda] = I \[CapitalLambda]//SparseArray;
Sinv = Inverse[S]//SparseArray;

Function[With[{u = S[[1]],v=Sinv[[All,L]],\[Lambda] = \[CapitalLambda]},Abs[\[Gamma] u . (1/(\[Lambda]-I #) v)]^2]]

];

NEGFTransmissionFunction[V_,L_,\[Gamma]_]:=NEGFTransmissionFunction[TightBindingHamiltonian[V,L],\[Gamma]]
NEGFTransmissionFunction[V_,L_,\[Gamma]_,J_]:=NEGFTransmissionFunction[TightBindingHamiltonian[V,L,J],\[Gamma]]


(* ::Input::Initialization:: *)
Clear[NEGFTransmissionFunctionFiniteBand]
NEGFTransmissionFunctionFiniteBand[h_,\[Gamma]_,bandwidth_,npts_:10^2]:=Module[{L=Length@h,\[CapitalGamma]1,\[CapitalGamma]L,W0,\[CapitalLambda]0,\[Epsilon]min,\[Epsilon]max,G,\[Epsilon]s,\[DoubleStruckCapitalI] = NEye[Length@h],tab},
{\[Epsilon]min,\[Epsilon]max} = If[ArrayQ[bandwidth],bandwidth,{-bandwidth,bandwidth}];
\[Epsilon]s = Rest@Most@linspace[\[Epsilon]min,\[Epsilon]max,npts];
\[CapitalGamma]1 = \[Gamma] Proj[L,1];
\[CapitalGamma]L = \[Gamma] Proj[L,L];

W0 = h - I/2 (\[CapitalGamma]1+\[CapitalGamma]L);
\[CapitalLambda]0 = 1/(2\[Pi]) (\[CapitalGamma]1+\[CapitalGamma]L);
tab=Table[
G = I Inverse[\[Epsilon] \[DoubleStruckCapitalI] - W0 - \[CapitalLambda]0 Log[(\[Epsilon]-\[Epsilon]min)/(\[Epsilon]max-\[Epsilon])]]; 
{\[Epsilon],Tr[\[CapitalGamma]1 . G . \[CapitalGamma]L . G\[ConjugateTranspose]]}
,{\[Epsilon],\[Epsilon]s}];
AppendTo[tab,{\[Epsilon]max,0}];
PrependTo[tab,{\[Epsilon]min,0}];

tab//Interpolation
]


(* ::Input::Initialization:: *)
Clear[FermiEnergyRegion];
FermiEnergyRegion[T_?NumberQ,\[Mu]_?NumberQ,prec_:10^-8]:=Module[{\[Beta]=1/T,log,emin,emax,energ},
log = Log[prec];
emin = \[Mu] + log/\[Beta];
emax = \[Mu] - log/\[Beta];
{emin,emax}//N
]

FermiEnergyRegion[{TL_,TR_,\[Mu]L_,\[Mu]R_},prec_:10^-8]:=Module[{rL,rR},
rL = FermiEnergyRegion[TL,\[Mu]L,prec];
rR = FermiEnergyRegion[TR,\[Mu]R,prec];
{Min[rL[[1]],rR[[1]]],Max[rL[[2]],rR[[2]]]}
]


(* ::Input::Initialization:: *)
LandauerButtikerAllAvailableProperties = {
"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]","Shyp","IQL","IQR","P","\[Sigma]","\[CapitalDelta]P","\[CapitalDelta]QL","\[CapitalDelta]QR","\[CapitalDelta]\[Sigma]","turN","turE","turP","turQL","turQR","fcorr","SNRN","SNRE","SNRP","SNRQL","SNRQR","SNR\[Sigma]","a","b","Slhyp","turlhyp","turhyp"
};


(* ::Input::Initialization:: *)
Clear[LandauerButtiker];
LandauerButtiker[\[ScriptCapitalT]_,{TL_,TR_,\[Mu]L_,\[Mu]R_},pprops_:All,em_:\[Infinity],pres_:MachinePrecision,maxrec_:\[Infinity]]:=Module[{\[Beta]L = 1/TL,\[Beta]R=1/TR,\[Delta]\[Beta],\[Delta]\[Mu],\[Delta]\[Beta]\[Mu],fL,fR,IN,IE,IQL,IQR,P,\[Sigma],\[CapitalDelta]N,\[CapitalDelta]E,\[ScriptCapitalC],\[CapitalDelta]P,\[CapitalDelta]QL,\[CapitalDelta]QR,\[CapitalDelta]\[Sigma],turN,turE,turP,turQL,turQR,a,b,Jlhyp,\[CapitalDelta]lhyp,SNRN,SNRE,SNRP,SNRQL,SNRQR,Slhyp,Shyp,hyperparameters,fcorr,turlhyp,turhyp,SNR\[Sigma],g,G,integralsNeeded,intsNeeded,\[DoubleStruckCapitalI],result,LandauerButtikerIntegrals,emm,tab,props},

props = Flatten@{pprops};
(* DEFINE AUXILIARY VARIABLES & FUNCTIONS *)
\[Delta]\[Beta] = \[Beta]L-\[Beta]R;
\[Delta]\[Mu] = \[Mu]L-\[Mu]R;
\[Delta]\[Beta]\[Mu] = \[Beta]L \[Mu]L - \[Beta]R \[Mu]R;
fL[\[Epsilon]_]:=1/(Exp[Rationalize[\[Beta]L,10^-8](\[Epsilon]-Rationalize[\[Mu]L,10^-8])]+1); 
fR[\[Epsilon]_]:=1/(Exp[Rationalize[\[Beta]R,10^-8](\[Epsilon]-Rationalize[\[Mu]R,10^-8])]+1);  
g[\[Epsilon]_]:=fL[\[Epsilon]](1-fL[\[Epsilon]])+fR[\[Epsilon]](1-fR[\[Epsilon]]);
G[\[Epsilon]_]:=\[ScriptCapitalT][\[Epsilon]](g[\[Epsilon]] + (fL[\[Epsilon]]-fR[\[Epsilon]])^2 (1-\[ScriptCapitalT][\[Epsilon]]));

(* LIST ALL INTEGRALS NEEDED, BASED ON WHICH PROPS THE USER PASSED TO THE FUNCTION *)
integralsNeeded[prop_]:=Which[
prop==="IN",{"IN"},
prop==="IE",{"IE"},
prop==="\[CapitalDelta]N",{"\[CapitalDelta]N"},
prop==="\[CapitalDelta]E",{"\[CapitalDelta]E"},
prop==="\[ScriptCapitalC]",{"\[ScriptCapitalC]"},
prop==="Shyp",{"Shyp"},
prop==="IQL",{"IN","IE"},
prop==="IQR",{"IN","IE"},
prop==="P",{"IN"},
prop==="\[Sigma]",{"IN","IE"},
prop==="\[CapitalDelta]P",{"\[CapitalDelta]N"},
prop==="\[CapitalDelta]QL",{"\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="\[CapitalDelta]QR",{"\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="\[CapitalDelta]\[Sigma]",{"\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="turN",{"IN","IE","\[CapitalDelta]N"},
prop==="turE",{"IN","IE","\[CapitalDelta]E"},
prop==="turP",{"IN","IE","\[CapitalDelta]N"},
prop==="turQL",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="turQR",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="fcorr",{"\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="SNRN",{"IN","\[CapitalDelta]N"},
prop==="SNRE",{"IE","\[CapitalDelta]E"},
prop==="SNRP",{"IN","\[CapitalDelta]N"},
prop==="SNRQL",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="SNRQR",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="SNR\[Sigma]",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="a",{"IN","IE","\[CapitalDelta]N","\[ScriptCapitalC]"},
prop==="b",{"IN","IE","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="Slhyp",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="turlhyp",{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]"},
prop==="turhyp",{"IN","IE","Shyp"},
prop===All,{"IN","IE","\[CapitalDelta]N","\[CapitalDelta]E","\[ScriptCapitalC]","Shyp"}
];
intsNeeded=DeleteDuplicates@Flatten@Map[integralsNeeded,Flatten@{props}];

emm = If[ArrayQ[em],em,{-em,em}];
(* COMPUTE SAID INTEGRALS *)
LandauerButtikerIntegrals[int_,emm_]:=Module[{func,\[Epsilon]},
func = Which[
int==="IN",\[ScriptCapitalT][\[Epsilon]](fL[\[Epsilon]]-fR[\[Epsilon]]),
int==="IE",\[Epsilon] \[ScriptCapitalT][\[Epsilon]](fL[\[Epsilon]]-fR[\[Epsilon]]),
int==="\[CapitalDelta]N",G[\[Epsilon]],
int==="\[CapitalDelta]E",\[Epsilon]^2 G[\[Epsilon]],
int==="\[ScriptCapitalC]",\[Epsilon] G[\[Epsilon]],
int==="Shyp",(\[ScriptCapitalT][\[Epsilon]](fL[\[Epsilon]]-fR[\[Epsilon]])^2)/(fL[\[Epsilon]]+fR[\[Epsilon]]-2 fL[\[Epsilon]] fR[\[Epsilon]]-\[ScriptCapitalT][\[Epsilon]](fL[\[Epsilon]]-fR[\[Epsilon]])^2) 
];
1/(2\[Pi]) NIntegrate[func,{\[Epsilon],emm[[1]],0},WorkingPrecision->pres,MaxRecursion->maxrec]
+1/(2\[Pi]) NIntegrate[func,{\[Epsilon],0,emm[[2]]},WorkingPrecision->pres,MaxRecursion->maxrec]
];

Do[\[DoubleStruckCapitalI][int]=LandauerButtikerIntegrals[int,emm],{int,intsNeeded}];


(* USE INTEGRALS TO CONSTRUCT ALL DESIRED RESULTS *)
result[prop_]:=Which[
prop==="IN",\[DoubleStruckCapitalI]["IN"],
prop==="IE",\[DoubleStruckCapitalI]["IE"],
prop==="\[CapitalDelta]N",\[DoubleStruckCapitalI]["\[CapitalDelta]N"],
prop==="\[CapitalDelta]E",\[DoubleStruckCapitalI]["\[CapitalDelta]E"],
prop==="\[ScriptCapitalC]",\[DoubleStruckCapitalI]["\[ScriptCapitalC]"],
prop==="Shyp",\[DoubleStruckCapitalI]["Shyp"],
prop==="IQL",\[DoubleStruckCapitalI]["IE"] - \[Mu]L \[DoubleStruckCapitalI]["IN"],
prop==="IQR",\[DoubleStruckCapitalI]["IE"] - \[Mu]R \[DoubleStruckCapitalI]["IN"],
prop==="P",-\[Delta]\[Mu] \[DoubleStruckCapitalI]["IN"],
prop==="\[Sigma]",-\[Delta]\[Beta] \[DoubleStruckCapitalI]["IE"] + \[Delta]\[Beta]\[Mu] \[DoubleStruckCapitalI]["IN"],
prop==="\[CapitalDelta]P",\[Delta]\[Mu]^2 \[DoubleStruckCapitalI]["\[CapitalDelta]N"],
prop==="\[CapitalDelta]QL",\[DoubleStruckCapitalI]["\[CapitalDelta]E"] - 2 \[Mu]L \[DoubleStruckCapitalI]["\[ScriptCapitalC]"] + \[Mu]L^2 \[DoubleStruckCapitalI]["\[CapitalDelta]N"],
prop==="\[CapitalDelta]QR", \[DoubleStruckCapitalI]["\[CapitalDelta]E"] - 2 \[Mu]R \[DoubleStruckCapitalI]["\[ScriptCapitalC]"] + \[Mu]R^2 \[DoubleStruckCapitalI]["\[CapitalDelta]N"],
prop==="\[CapitalDelta]\[Sigma]",\[Delta]\[Beta]^2 \[DoubleStruckCapitalI]["\[CapitalDelta]E"] - 2 \[Delta]\[Beta] \[Delta]\[Beta]\[Mu] \[DoubleStruckCapitalI]["\[ScriptCapitalC]"] + \[Delta]\[Beta]\[Mu]^2 \[DoubleStruckCapitalI]["\[CapitalDelta]N"],
prop==="turN",\[DoubleStruckCapitalI]["\[CapitalDelta]N"]/\[DoubleStruckCapitalI]["IN"]^2 result["\[Sigma]"]//Quiet,
prop==="turE",\[DoubleStruckCapitalI]["\[CapitalDelta]E"]/\[DoubleStruckCapitalI]["IE"]^2 result["\[Sigma]"]//Quiet,
prop==="turP",result["\[CapitalDelta]P"]/result["P"]^2 result["\[Sigma]"]//Quiet, 
prop==="turQL",result["\[CapitalDelta]QL"]/result["IQL"]^2 result["\[Sigma]"]//Quiet, 
prop==="turQR",result["\[CapitalDelta]QR"]/result["IQR"]^2 result["\[Sigma]"]//Quiet, 
prop==="fcorr",\[DoubleStruckCapitalI]["\[ScriptCapitalC]"]/Sqrt[\[DoubleStruckCapitalI]["\[CapitalDelta]N"] \[DoubleStruckCapitalI]["\[CapitalDelta]E"]],
prop==="SNRN",\[DoubleStruckCapitalI]["IN"]^2/\[DoubleStruckCapitalI]["\[CapitalDelta]N"], 
prop==="SNRE",\[DoubleStruckCapitalI]["IE"]^2/\[DoubleStruckCapitalI]["\[CapitalDelta]E"],
prop==="SNRP",result["P"]^2/result["\[CapitalDelta]P"],
prop==="SNRQL",result["IQL"]^2/result["\[CapitalDelta]QL"],
prop==="SNRQR",result["IQR"]^2/result["\[CapitalDelta]QR"],
prop==="SNR\[Sigma]",result["\[Sigma]"]^2/result["\[CapitalDelta]\[Sigma]"],
(* Hyperaccurate results from arXiv:2108.05325 *)
prop==="a",\[DoubleStruckCapitalI]["IE"] \[DoubleStruckCapitalI]["\[CapitalDelta]N"] - \[DoubleStruckCapitalI]["IN"] \[DoubleStruckCapitalI]["\[ScriptCapitalC]"],
prop==="b",\[DoubleStruckCapitalI]["IN"] \[DoubleStruckCapitalI]["\[CapitalDelta]E"]-\[DoubleStruckCapitalI]["IE"] \[DoubleStruckCapitalI]["\[ScriptCapitalC]"], 
prop==="Slhyp",(result["a"] \[DoubleStruckCapitalI]["IE"] + result["b"] \[DoubleStruckCapitalI]["IN"])^2/(result["a"]^2 \[DoubleStruckCapitalI]["\[CapitalDelta]E"]+result["b"]^2 \[DoubleStruckCapitalI]["\[CapitalDelta]N"]+ 2 result["a"] result["b"] \[DoubleStruckCapitalI]["\[ScriptCapitalC]"]), 
prop==="turlhyp",result["\[Sigma]"]/result["Slhyp"],
prop==="turhyp",result["\[Sigma]"]/\[DoubleStruckCapitalI]["Shyp"]
];

tab=If[pprops ===All,
Chop@Table[p->result[p],{p,LandauerButtikerAllAvailableProperties}],
Chop@Table[p->result[p],{p,props}]];
If[Length@props==1 && Not[pprops ===All],
tab[[1,2]],
Association[tab]
]

]


(* ::Input::Initialization:: *)
(* Auxiliary function, pre-compiled for speed *)
Clear[NEGFAuxMat];
NEGFAuxMat= Compile[{{\[CapitalLambda],_Complex,1},{\[Beta],_Real},{\[Mu],_Real}}, 

Module[{x,y,\[Psi]p,\[Psi]m},

x = I \[CapitalLambda]\[Conjugate]-\[Mu];
y=-I \[CapitalLambda]-\[Mu];
\[Psi]p = Table[PolyGamma[(I \[Beta] y[[i]])/(2\[Pi])+1/2],{i,Length@\[CapitalLambda]}];
\[Psi]m = Table[PolyGamma[-((I \[Beta] x[[j]])/(2\[Pi]))+1/2],{j,Length@\[CapitalLambda]}];

Table[1/(x[[j]]-y[[i]] ) (I/2+1/(2\[Pi]) (\[Psi]m[[j]]-\[Psi]p[[i]])),{i,Length@\[CapitalLambda]},{j,Length@\[CapitalLambda]}]

](*,CompilationTarget\[Rule]"C"*)
];


(* ::Input::Initialization:: *)
NEGFCovMat[h_,\[Gamma]_,{TL_,TR_,\[Mu]L_,\[Mu]R_}]:=Module[{L = Length@h,\[CapitalGamma],W,\[CapitalLambda],S,Sinv,g,ppart,K,Y,\[ScriptCapitalI],amat,\[CapitalGamma]1,\[CapitalGamma]L,A1,AL},

\[CapitalGamma]1 = SparseArray[{{1,1}->\[Gamma]},{L,L}];
\[CapitalGamma]L = SparseArray[{{L,L}->\[Gamma]},{L,L}];
\[CapitalGamma] = \[CapitalGamma]1+\[CapitalGamma]L;
W = -I \[CapitalGamma]/2+ h;
{\[CapitalLambda],S}=Quiet@Eigensystem[W];
S = S\[Transpose]//SparseArray;
\[CapitalLambda] = I \[CapitalLambda]//SparseArray;
Sinv = Inverse[S]//SparseArray;
A1 =SparseArray@NEGFAuxMat[Normal@\[CapitalLambda],1/TL,\[Mu]L];
AL=SparseArray@NEGFAuxMat[Normal@\[CapitalLambda],1/TR,\[Mu]R];
Y = A1 (Sinv . \[CapitalGamma]1 . Sinv\[ConjugateTranspose])+AL (Sinv . \[CapitalGamma]L . Sinv\[ConjugateTranspose]);
S . Y . S\[ConjugateTranspose]//Chop
]


(* ::Input::Initialization:: *)
NEGFCurrents[V_,L_,\[Gamma]_,\[Beta]L_,\[Beta]R_,\[Mu]L_,\[Mu]R_]:=Module[{\[ScriptCapitalC],JN,JE,JQL,JQR,P,\[Sigma],v,i=2},

\[ScriptCapitalC] = NEGFCovMat[V,L,\[Gamma],\[Beta]L,\[Beta]R,\[Mu]L,\[Mu]R];
JN =  -I (\[ScriptCapitalC][[2,1]]-\[ScriptCapitalC][[1,2]])//Chop;
v = If[VectorQ[V],V[[1]],V]//Chop;
JE = I v ((\[ScriptCapitalC][[i-1,i]]-\[ScriptCapitalC][[i,i-1]])-(\[ScriptCapitalC][[i+1,i]]-\[ScriptCapitalC][[i,i+1]]))+ I (\[ScriptCapitalC][[i+1,i-1]]-\[ScriptCapitalC][[i-1,i+1]])//Chop;
JQL = JE - \[Mu]L JN//Chop;
JQR = JE - \[Mu]R JN//Chop;
P = -(\[Mu]L-\[Mu]R) JN//Chop;
\[Sigma] = -(\[Beta]L-\[Beta]R) JE + (\[Beta]L \[Mu]L - \[Beta]R \[Mu]R) JN//Chop;
<|"JN"->JN,"JE"->JE,"JQL"->JQL,"JQR"->JQR,"P"->P,"\[Sigma]"->\[Sigma]|>
]



(* ::Input::Initialization:: *)
Clear[FermionicChainWaitingTimeDistributionVacuum];
FermionicChainWaitingTimeDistributionVacuum[W_,F_,trange_,{i_,diri_},j_]:=Module[{L = Length@W,\[ScriptCapitalC],EM,EmM,Q,\[CapitalGamma],EX,EZ,P,\[Gamma]},

Q = W - F;
\[CapitalGamma] = Tr[F];
\[Gamma] = If[diri =="-",2 Re@W[[i,i]]-F[[i,i]],F[[i,i]]];

If[diri=="-",
Table[
EX = MatrixExp[-Q t];
EZ = EX\[ConjugateTranspose];
{t,\[Gamma] E^(-\[CapitalGamma] t) EX [[i,j]] EZ[[j,i]]},{t,trange}],

Table[
EX = MatrixExp[-Q t];
EZ = EX\[ConjugateTranspose];
{t,\[Gamma] E^(-\[CapitalGamma] t) ((EZ . EX)[[j,j]] - EX [[i,j]] EZ[[j,i]])},{t,trange}]]//Chop

]

Clear[FermionicChainWaitingTimeDistribution];
FermionicChainWaitingTimeDistribution[WW_,FF_,ttrange_,{i_,diri_},{j_,dirj_},CC_:"ss",pres_:MachinePrecision]:=Module[{L = Length@WW,W,F,\[ScriptCapitalC],EM,EmM,Q,X,\[CapitalGamma],c,\[DoubleStruckCapitalI],\[ScriptCapitalT],\[DoubleStruckCapitalD],\[Gamma],deno,trange,\[ScriptCapitalZ],EX,EmX,EZ,EmZ},

trange = SetPrecision[Rationalize@ttrange,pres];
W =SparseArray@SetPrecision[Rationalize@Normal@WW,pres];
F =SparseArray@SetPrecision[Rationalize@Normal@FF,pres];
\[DoubleStruckCapitalI] = Eye[L];

If[CC=="vac",Return[FermionicChainWaitingTimeDistributionVacuum[W,F,trange,{i,diri},j]]];

\[ScriptCapitalC]=If[CC==="ss",Chop@LyapunovSolve[W,F],CC];
EM = Inverse[\[ScriptCapitalC]]-\[DoubleStruckCapitalI];
EmM = Inverse[EM];

Q = W - F;
\[ScriptCapitalZ] = Det[\[DoubleStruckCapitalI] + EmM];
\[CapitalGamma] = Tr[F];

\[Gamma] = If[diri =="-",2 Re@W[[i,i]]-F[[i,i]],F[[i,i]]];
c = If[dirj=="-",\[ScriptCapitalC][[j,j]],1-\[ScriptCapitalC][[j,j]]];

Table[
X = -Q t;
(*Z = - Q\[ConjugateTranspose] t; (* X\[ConjugateTranspose] *)*)
EX = MatrixExp[X];
EmX = MatrixExp[-X];
EZ = EX\[ConjugateTranspose] (* \[ScriptCapitalE][Z]*);
EmZ=EmX\[ConjugateTranspose];

\[ScriptCapitalT] = Inverse[EmZ . EM . EmX+\[DoubleStruckCapitalI]];
\[DoubleStruckCapitalD] = Det[\[DoubleStruckCapitalI] + EX . EmM . EZ];


{t,\[Gamma]/c Exp[-\[CapitalGamma] t]/\[ScriptCapitalZ]  \[DoubleStruckCapitalD] Which[
(* \!\(
\*SubsuperscriptBox[\(c\), \(i\), \(\[ConjugateTranspose]\)]
\*SubscriptBox[\(c\), \(i\)]\ and\ 
\*SubsuperscriptBox[\(c\), \(j\), \(\[ConjugateTranspose]\)]
\*SubscriptBox[\(\[Rho]c\), \(j\)]\) *)diri=="-" && dirj=="+",  
((EM . EmX . \[ScriptCapitalT] . EX)[[j,j]] \[ScriptCapitalT][[i,i]] + (\[ScriptCapitalT] . EmZ . EM)[[i,j]] (EM . EmX . \[ScriptCapitalT])[[j,i]]),
(* Subscript[c, i]Subsuperscript[c, i, \[ConjugateTranspose]] and Subsuperscript[c, j, \[ConjugateTranspose]]Subscript[\[Rho]c, j] *)diri=="+" && dirj=="+", 
((EM . EmX . \[ScriptCapitalT] . EX)[[j,j]] (\[DoubleStruckCapitalI][[i,i]]-\[ScriptCapitalT][[i,i]]) - (\[ScriptCapitalT] . EmZ . EM)[[i,j]] (EM . EmX . \[ScriptCapitalT])[[j,i]]),
(* \!\(
\*SubsuperscriptBox[\(c\), \(i\), \(\[ConjugateTranspose]\)]
\*SubscriptBox[\(c\), \(i\)]\ and\ 
\*SubscriptBox[\(c\), \(j\)]
\*SubsuperscriptBox[\(\[Rho]c\), \(j\), \(\[ConjugateTranspose]\)]\) *)diri=="-" && dirj=="-", 
(\[ScriptCapitalT][[i,i]] (EmM[[j,j]]-(EmX . \[ScriptCapitalT] . EX . EmM)[[j,j]]) - (\[ScriptCapitalT] . EmZ)[[i,j]] (EmX . \[ScriptCapitalT])[[j,i]]),
(* Subscript[c, i]Subsuperscript[c, i, \[ConjugateTranspose]] and Subscript[c, j]Subsuperscript[\[Rho]c, j, \[ConjugateTranspose]] *)diri=="+" && dirj=="-", 
((\[DoubleStruckCapitalI][[i,i]]-\[ScriptCapitalT][[i,i]]) (EmM[[j,j]]-(EmX . \[ScriptCapitalT] . EX . EmM)[[j,j]]) + (\[ScriptCapitalT] . EmZ)[[i,j]] (EmX . \[ScriptCapitalT])[[j,i]])
]},
{t,trange}]//Chop

]

FermionicChainWaitingTimeDistribution[V_,L_,\[Gamma]_,{f1_,fL_},ttrange_,{i_,diri_},{j_,dirj_},CC_:"ss",pres_:MachinePrecision]:=Module[{W,F},
{W,F} = BDMatrices[V,L,\[Gamma],{f1,fL}];
FermionicChainWaitingTimeDistribution[W,F,ttrange,{i,diri},{j,dirj},CC,pres]]


(* ::Input::Initialization:: *)
Clear[FermionicChainNoJump];
FermionicChainNoJump[WW_,FF_,ttrange_,CC_:"ss",pres_:MachinePrecision]:=Module[{L = Length@WW,trange,W,F,\[DoubleStruckCapitalI],\[ScriptCapitalC],EM,EmM,Q,\[ScriptCapitalZ],\[CapitalGamma],EX},
trange = SetPrecision[Rationalize@ttrange,pres];
W =SparseArray@SetPrecision[Rationalize@Normal@WW,pres];
F =SparseArray@SetPrecision[Rationalize@Normal@FF,pres];
\[DoubleStruckCapitalI] = Eye[L];

\[ScriptCapitalC]=If[CC==="ss",Chop@LyapunovSolve[W,F],CC];
EM = Inverse[\[ScriptCapitalC]]-\[DoubleStruckCapitalI];
EmM = Inverse[EM];

Q = W - F;
\[ScriptCapitalZ] = Det[\[DoubleStruckCapitalI] + EmM];
\[CapitalGamma] = Tr[F];

Table[
EX = MatrixExp[-Q t];
{t,E^(-\[CapitalGamma] t)/\[ScriptCapitalZ] Det[\[DoubleStruckCapitalI] + EX . EmM . EX\[ConjugateTranspose]]},{t,trange}]//Chop
]

FermionicChainNoJump[V_,L_,\[Gamma]_,{f1_,fL_},ttrange_,CC_,pres_:MachinePrecision]:=Module[{W,F},
{W,F} = BDMatrices[V,L,\[Gamma],{f1,fL}];
FermionicChainNoJump[W,F,ttrange,CC,pres]]


(* ::Input::Initialization:: *)
Clear[RiccatiEigen];
RiccatiEigen[a_,q_,brb_]:=Module[{Z,U,d,U11,U21},

Z =ArrayFlatten[({
 {a, -brb},
 {-q, -a\[ConjugateTranspose]}
})];
(*U = Eigvecs[Z,Length[Z]/2];*)
U=Transpose[(Reverse@Eigenvectors[Z + 10^5 Eye@Length@Z,-Length[Z]/2])];
d=Length[Z]/2;
U11 = U[[1;;d,1;;d]];
U21 = U[[d+1;;2d,1;;d]];
U21 . Inverse[U11]
]



(* ::Input::Initialization:: *)
Clear[ConditionalFermionicDynamics];
ConditionalFermionicDynamics[H_,\[Gamma]_,f_,{\[Eta]m_,\[Eta]p_},tf_,CC_:"ss",solver_:"NDSolve",dt_:0.01]:=Module[{L,W,\[CapitalUpsilon],\[ScriptCapitalW],\[ScriptCapitalF],\[Gamma]p,\[Gamma]m,\[DoubleStruckCapitalC],\[ScriptCapitalC],\[ScriptCapitalC]0,\[ScriptCapitalA]m,\[ScriptCapitalA]p,\[CapitalTheta]0,\[CapitalTheta],\[DoubleStruckCapitalI],\[CapitalPhi]},

L=Length[\[Eta]m];
{W,\[Gamma]p} = BDMatrices[Eye[L],\[Gamma],f];
\[Gamma]m = W+W\[ConjugateTranspose]-\[Gamma]p;
\[ScriptCapitalA]m = Normal@DiagonalMatrix[\[Eta]m];
\[ScriptCapitalA]p = Normal@DiagonalMatrix[\[Eta]p];

If[L==Length@H,
W = I H + W;
(\[ScriptCapitalC]0=Normal@If[CC==="ss",Chop@LyapunovSolve[W,\[Gamma]p],CC];
If[solver ==="NDSolve",
NDSolveValue[{
\[ScriptCapitalC]'[t]==-(W . \[ScriptCapitalC][t] + \[ScriptCapitalC][t] . W\[ConjugateTranspose]) + \[Gamma]p+ \[ScriptCapitalC][t] . \[Gamma]m . \[ScriptCapitalA]m . \[ScriptCapitalC][t]-\[Gamma]p . \[ScriptCapitalA]p+ \[Gamma]p . \[ScriptCapitalA]p . \[ScriptCapitalC][t]+ \[ScriptCapitalC][t] . \[Gamma]p . \[ScriptCapitalA]p-\[ScriptCapitalC][t] . \[Gamma]p . \[ScriptCapitalA]p . \[ScriptCapitalC][t] ,
\[ScriptCapitalC][0]==\[ScriptCapitalC]0},\[ScriptCapitalC],{t,0,tf}],
Interpolation@RK4[-(W . # + # . W\[ConjugateTranspose]) + \[Gamma]p+ # . \[Gamma]m . \[ScriptCapitalA]m . #-\[Gamma]p . \[ScriptCapitalA]p+ \[Gamma]p . \[ScriptCapitalA]p . #+ # . \[Gamma]p . \[ScriptCapitalA]p-# . \[Gamma]p . \[ScriptCapitalA]p . #&,\[ScriptCapitalC]0,dt,tf]
]),
(
\[CapitalUpsilon] = 1/4 ArrayFlatten[({
 {\[Gamma]p+\[Gamma]m, I(\[Gamma]p-\[Gamma]m)},
 {-I (\[Gamma]p-\[Gamma]m), \[Gamma]p+\[Gamma]m}
})];
\[ScriptCapitalW] = I H + \[CapitalUpsilon]+\[CapitalUpsilon]\[Transpose];
\[ScriptCapitalF] = - 2(\[CapitalUpsilon]-\[CapitalUpsilon]\[Transpose]);

\[CapitalTheta]0 = Normal@If[CC==="ss",Chop@LyapunovSolve[\[ScriptCapitalW],\[ScriptCapitalF]],CC];
 \[CapitalPhi]=1/4 ArrayFlatten[({
 {\[Gamma]p . \[ScriptCapitalA]p+\[Gamma]m . \[ScriptCapitalA]m, I(\[Gamma]p . \[ScriptCapitalA]p-\[Gamma]m . \[ScriptCapitalA]m)},
 {-I (\[Gamma]p . \[ScriptCapitalA]p-\[Gamma]m . \[ScriptCapitalA]m), \[Gamma]p . \[ScriptCapitalA]p+\[Gamma]m . \[ScriptCapitalA]m}
})];
\[ScriptCapitalW] = \[ScriptCapitalW] -( \[CapitalPhi]+\[CapitalPhi]\[Transpose]);
\[ScriptCapitalF] = \[ScriptCapitalF] + \[CapitalPhi] - \[CapitalPhi]\[Transpose];
If[solver ==="NDSolve",
NDSolveValue[{
\[CapitalTheta]'[t]==-(\[ScriptCapitalW] . \[CapitalTheta][t] + \[CapitalTheta][t] . \[ScriptCapitalW]\[ConjugateTranspose]) + \[ScriptCapitalF]+\[CapitalTheta][t] . (\[CapitalPhi]-\[CapitalPhi]\[Transpose]) . \[CapitalTheta][t],
\[CapitalTheta][0]==\[CapitalTheta]0},\[CapitalTheta],{t,0,tf}],
Interpolation@RK4[-(\[ScriptCapitalW] . # + # . \[ScriptCapitalW]\[ConjugateTranspose]) + \[ScriptCapitalF]+# . (\[CapitalPhi]-\[CapitalPhi]\[Transpose]) . #&,\[CapitalTheta]0,dt,tf]
]
)]
];


(* ::Input::Initialization:: *)
Clear[ConditionalFermionicNoJumpDistribution];
ConditionalFermionicNoJumpDistribution[H_,\[Gamma]_,f_,{\[Eta]m_,\[Eta]p_},tf_,\[CapitalDelta]t_,\[ScriptCapitalC]_]:=Module[{\[Lambda],\[Gamma]p,\[Gamma]m,\[ScriptCapitalA]m,\[ScriptCapitalA]p,\[DoubleStruckCapitalI],Z0,\[Zeta],Zt,Pno,L,W},
L=Length[\[Eta]m];
{W,\[Gamma]p} = BDMatrices[Eye[L],\[Gamma],f];
\[Gamma]m = W+W\[ConjugateTranspose]-\[Gamma]p;
\[DoubleStruckCapitalI] = Eye[L];
\[ScriptCapitalA]m = Normal@DiagonalMatrix[\[Eta]m];
\[ScriptCapitalA]p = Normal@DiagonalMatrix[\[Eta]p];

\[Lambda]=If[L==Length@H,
Table[Tr[\[Gamma]p . \[ScriptCapitalA]p . (\[DoubleStruckCapitalI]-\[ScriptCapitalC][t]) + \[Gamma]m . \[ScriptCapitalA]m . \[ScriptCapitalC][t]],{t,0,tf,\[CapitalDelta]t}],Table[Tr[\[Gamma]p . \[ScriptCapitalA]p . (\[DoubleStruckCapitalI]-FermionsReducedC@\[ScriptCapitalC][t]) + \[Gamma]m . \[ScriptCapitalA]m . FermionsReducedC@\[ScriptCapitalC][t]],{t,0,tf,\[CapitalDelta]t}]
];
(*{Range[\[CapitalDelta]t,tf,\[CapitalDelta]t],Table[Exp[-Trapz[\[Lambda]\[LeftDoubleBracket]1;;j\[RightDoubleBracket],\[CapitalDelta]t]],{j,2,Length@\[Lambda]}]}\[Transpose]//Chop*)
{Range[\[CapitalDelta]t,tf,\[CapitalDelta]t],Table[-Trapz[\[Lambda][[1;;j]],\[CapitalDelta]t],{j,2,Length@\[Lambda]}]}\[Transpose]//Chop
]


(* ::Input::Initialization:: *)
Clear[ConditionalFermionicSteadyState];
ConditionalFermionicSteadyState[H_,\[Gamma]_,f_,{\[Eta]m_,\[Eta]p_},pres_:MachinePrecision]:=Module[{a,\[Gamma]p,\[Gamma]m,W,\[ScriptCapitalA]p,\[ScriptCapitalA]m,L,\[CapitalUpsilon],\[ScriptCapitalW],\[ScriptCapitalF],\[CapitalPhi]},
L=Length[\[Eta]m];
{W,\[Gamma]p} = BDMatrices[Eye[L],\[Gamma],f];
\[Gamma]m = W+W\[ConjugateTranspose]-\[Gamma]p;
\[ScriptCapitalA]m = Normal@DiagonalMatrix[\[Eta]m];
\[ScriptCapitalA]p = Normal@DiagonalMatrix[\[Eta]p];

If[L==Length@H,
W = I H + W;
(*RiccatiEigen[- W\[ConjugateTranspose]+\[Gamma]p.\[ScriptCapitalA]p, \[Gamma]p-\[Gamma]p.\[ScriptCapitalA]p, -(\[Gamma]m.\[ScriptCapitalA]m-\[Gamma]p.\[ScriptCapitalA]p)]*)
RiccatiSolve[{- W\[ConjugateTranspose]+\[Gamma]p . \[ScriptCapitalA]p,Chop@Sqrt[-(\[Gamma]m . \[ScriptCapitalA]m-\[Gamma]p . \[ScriptCapitalA]p)]},{ \[Gamma]p-\[Gamma]p . \[ScriptCapitalA]p,-NEye[L]} ]
,
\[CapitalUpsilon] = 1/4 ArrayFlatten[({
 {\[Gamma]p+\[Gamma]m, I(\[Gamma]p-\[Gamma]m)},
 {-I (\[Gamma]p-\[Gamma]m), \[Gamma]p+\[Gamma]m}
})];
\[ScriptCapitalW] = SetPrecision[I H + \[CapitalUpsilon]+\[CapitalUpsilon]\[Transpose],pres];
\[ScriptCapitalF] = SetPrecision[- 2(\[CapitalUpsilon]-\[CapitalUpsilon]\[Transpose]),pres];
 \[CapitalPhi]=SetPrecision[1/4 ArrayFlatten[({
 {\[Gamma]p . \[ScriptCapitalA]p+\[Gamma]m . \[ScriptCapitalA]m, I(\[Gamma]p . \[ScriptCapitalA]p-\[Gamma]m . \[ScriptCapitalA]m)},
 {-I (\[Gamma]p . \[ScriptCapitalA]p-\[Gamma]m . \[ScriptCapitalA]m), \[Gamma]p . \[ScriptCapitalA]p+\[Gamma]m . \[ScriptCapitalA]m}
})],pres];
(*RiccatiEigen[-\[ScriptCapitalW]\[ConjugateTranspose]+(\[CapitalPhi]+\[CapitalPhi]\[Transpose]),\[ScriptCapitalF]+\[CapitalPhi]-\[CapitalPhi]\[Transpose],-(\[CapitalPhi]-\[CapitalPhi]\[Transpose])]*)
RiccatiSolve[{-\[ScriptCapitalW]\[ConjugateTranspose]+(\[CapitalPhi]+\[CapitalPhi]\[Transpose]),Chop@Sqrt[-(\[CapitalPhi]-\[CapitalPhi]\[Transpose])]},{\[ScriptCapitalF]+\[CapitalPhi]-\[CapitalPhi]\[Transpose],kron[\[Sigma]x,NEye[L]]}]
]
]


(* ::Input::Initialization:: *)
Clear[GPTr];
GPTr[\[Sigma]_,list_]:=Module[{blist,clist,L = Length[\[Sigma]]/2},

blist = Riffle[2list-1,2list];
clist = Complement[Range[1,2L],blist];
If[2Length@list==Length@\[Sigma],\[Sigma],\[Sigma][[clist,clist]]]
]


(* ::Input::Initialization:: *)
Clear[SymplecticForm];
SymplecticForm[\[ScriptCapitalN]_]:=kron[Eye[\[ScriptCapitalN]],{{0,1},{-1,0}}];


(* ::Input::Initialization:: *)
Clear[SymplecticEigenvalues];
SymplecticEigenvalues[\[Sigma]_]:=Module[{\[ScriptCapitalN] = Length[\[Sigma]]/2,\[CapitalOmega]},
\[CapitalOmega] = SymplecticForm[\[ScriptCapitalN]];
Chop@Drop[Sort@Quiet@Eigenvalues[N[I \[CapitalOmega] . \[Sigma]]],\[ScriptCapitalN]]
]


(* ::Input::Initialization:: *)
(* Valid CM should have BonaFide \[GreaterEqual] 0 *)
Clear[BonaFide];
BonaFide[\[Sigma]_]:=Module[{\[CapitalOmega],\[ScriptCapitalN] = Length[\[Sigma]]/2},
\[CapitalOmega] = SymplecticForm[\[ScriptCapitalN]];
First@Sort@Eigenvalues[N[\[Sigma]+I \[CapitalOmega]/2]]//Chop
]


(* ::Input::Initialization:: *)
(* Transforms from (a,a\[ConjugateTranspose]) to quadratures *)
Clear[ModeQuadratureRotation];
ModeQuadratureRotation[\[ScriptCapitalN]_]:=kron[Eye[\[ScriptCapitalN]],1/Sqrt[2] ({
 {1, 1},
 {-I, I}
})];


(* ::Input::Initialization:: *)
Clear[GaussianPurity]
GaussianPurity[\[Sigma]_]:=1/Sqrt[Det[2\[Sigma]]];


(* ::Input::Initialization:: *)
Williamson[\[Sigma]_]:=Module[{\[ScriptCapitalN]=Length[\[Sigma]]/2,\[CapitalOmega],W,V,vals,vecs,list,perm,X,U,\[Gamma],\[CapitalGamma],S,P,error1,error2},
(*PhysRevA_79_052327*)
\[CapitalOmega]=SymplecticForm[\[ScriptCapitalN]];
\[CapitalGamma] = ModeQuadratureRotation[\[ScriptCapitalN]];
W =DiagonalMatrix@Abs@Eigenvalues[I \[CapitalOmega] . \[Sigma]];

V = MatrixPower[\[Sigma],-1/2];
X=V . \[CapitalOmega] . V;
{vals,vecs}=Chop@Eigensystem[X];

list=I Reverse[ Abs@vals Table[(-1)^i,{i,0,2\[ScriptCapitalN]-1}]];
perm=FindPermutation[vals,list];
U=Permute[vecs,perm]\[Transpose];
S=Inverse@(Sqrt[W] . (\[CapitalGamma] . U\[ConjugateTranspose]) . V);
error1 = Max@Flatten@Abs[S . W . S\[Transpose]-\[Sigma]]//Chop;
error2 = Max@Flatten@Abs[S . \[CapitalOmega] . S\[Transpose]-\[CapitalOmega]]//Chop;
{W,Chop@S,{error1,error2}}
]


(* ::Input::Initialization:: *)
GaussianVonNeumann[\[Sigma]_]:=Module[{h,g,\[ScriptCapitalN]=Length[\[Sigma]]/2,\[CapitalOmega],\[Nu]s},h[x_]:=If[x==0||x==1,0,x Log[x]];
g[x_]:=h[(x+1)/2]-h[(x-1)/2];
\[CapitalOmega]=kron[Eye[\[ScriptCapitalN]],{{0,1},{-1,0}}];
(*\[Nu]s=Chop@Drop[Sort@Quiet@Eigenvalues[I 2 \[CapitalOmega].\[Sigma]],\[ScriptCapitalN]];*)
\[Nu]s = 2 SymplecticEigenvalues[\[Sigma]];
Total@Table[g[\[Nu]],{\[Nu],\[Nu]s}]//Chop
];

GaussianRenyiEntropy[\[Sigma]_,\[Alpha]_]:=Module[{\[ScriptCapitalN]=Length[\[Sigma]]/2,\[CapitalOmega],\[Nu]s},
\[CapitalOmega]=kron[Eye[\[ScriptCapitalN]],{{0,1},{-1,0}}];
(*\[Nu]s=Chop@Drop[Sort@Quiet@Eigenvalues[I 2 \[CapitalOmega].\[Sigma]],\[ScriptCapitalN]];*)
\[Nu]s = 2 SymplecticEigenvalues[\[Sigma]];
If[\[Alpha]==1,GaussianVonNeumann[\[Sigma]],
1/(\[Alpha]-1) Total@Table[Log[((\[Nu]+1)/2)^\[Alpha]-((\[Nu]-1)/2)^\[Alpha]],{\[Nu],\[Nu]s}]]

]


(* ::Input::Initialization:: *)
GaussianMutualInformation[\[Sigma]_,Alist_,Blist_]:=GaussianVonNeumann[GPTr[\[Sigma],Blist]]+GaussianVonNeumann[GPTr[\[Sigma],Alist]]-GaussianVonNeumann[\[Sigma]]

GaussianMutualInformation[\[Sigma]_,Alist_,Blist_,\[Alpha]_]:=GaussianRenyiEntropy[GPTr[\[Sigma],Blist],\[Alpha]]+GaussianRenyiEntropy[GPTr[\[Sigma],Alist],\[Alpha]]-GaussianRenyiEntropy[\[Sigma],\[Alpha]]


(* ::Input::Initialization:: *)
nGmeasure[\[Rho]_]:=Module[{\[Sigma],\[ScriptCapitalN],a,X,P},
\[ScriptCapitalN]= Length@\[Rho];
a = SparseArray@Table[Sqrt[n]KroneckerDelta[m,n-1],{m,0,\[ScriptCapitalN]-1},{n,0,\[ScriptCapitalN]-1}];
X = 1/Sqrt[2] (a\[ConjugateTranspose]+a);
P = I/Sqrt[2] (a\[ConjugateTranspose]-a);

\[Sigma] = ({
 {Tr[\[Rho] . X . X]-Tr[\[Rho] . X]^2, 1/2 Tr[\[Rho] . (X . P+P . X)]-Tr[\[Rho] . X] Tr[\[Rho] . P]},
 {1/2 Tr[\[Rho] . (X . P+P . X)]-Tr[\[Rho] . X] Tr[\[Rho] . P], Tr[\[Rho] . P . P]-Tr[\[Rho] . P]^2}
});
GaussianVonNeumann[\[Sigma]]-vonNeumann[\[Rho]]
]


(* ::Input::Initialization:: *)
(* vacuum is 1/2 *)
TwoModeCanonicalForm[\[Sigma]_]:=Module[{\[Sigma]A,\[Sigma]B,\[Xi],DA,VA,ZA,DB,VB,ZB,M,\[CapitalSigma],UA,UB,S,\[Sigma]canon},
\[Sigma]A = \[Sigma][[1;;2,1;;2]];
\[Sigma]B = \[Sigma][[3;;4,3;;4]];
\[Xi] = \[Sigma][[1;;2,3;;4]];

{DA,VA} = Eigensystem[\[Sigma]A];
VA = VA\[Transpose];
ZA = DiagonalMatrix[{(DA[[2]]/DA[[1]])^(1/4),(DA[[2]]/DA[[1]])^(-1/4)}];

{DB,VB} = Eigensystem[\[Sigma]B];
VB = VB\[Transpose];
ZB = DiagonalMatrix[{(DB[[2]]/DB[[1]])^(1/4),(DB[[2]]/DB[[1]])^(-1/4)}];

M= ZA . VA\[Transpose] . \[Xi] . VB . ZB\[Transpose];

(*{UA,\[CapitalSigma]A,UB} = SingularValueDecomposition[ZA.VA\[Transpose].\[Xi].VB.ZB\[Transpose]];*)
(*UA = UA\[ConjugateTranspose]; *)
(*UB = UB\[ConjugateTranspose];*)

{\[CapitalSigma],UA} = Eigensystem[M . M\[Transpose]];
{\[CapitalSigma],UB} = Eigensystem[M\[Transpose] . M];
S = ArrayFlatten[({
 {UA . ZA . VA\[Transpose], 0},
 {0, UB . ZB . VB\[Transpose]}
})];
\[Sigma]canon = S . \[Sigma] . S\[Transpose];
{\[Sigma]canon,S}//Chop
]


(* ::Input::Initialization:: *)
Clear[GaussianEoF];
GaussianEoF[\[Sigma]bare_]:=Module[{\[Sigma],\[Theta],na,nb,kq,kp,L,Cp,Cq,X,sol,det1,det2,tab1,tab2,opt1,opt2,x0,x,\[Theta]min,detmin,Xmin,\[Gamma]p},
LoadPauliMatrices[];
\[Sigma] = First@TwoModeCanonicalForm[2\[Sigma]bare]; (* Factor of 2 because the calculations are done with vacuum being 1.*)

na = \[Sigma][[1,1]];
nb = \[Sigma][[3,3]];
kq = \[Sigma][[1,3]];
kp = \[Sigma][[2,4]];

L = \!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"1", "0", "0", "0"},
{"0", "0", "1", "0"},
{"0", "1", "0", "0"},
{"0", "0", "0", "1"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[SparseArray[Automatic, {4, 4}, 0, {1, {{0, 1, 2, 3, 4}, {{1}, {3}, {2}, {4}}}, {1, 1, 1, 1}}]]]]\);
Cq = ({
 {na, kq},
 {kq, nb}
});
Cp = ({
 {na, kp},
 {kp, nb}
});
X = x0 \[Sigma]0 + x Cos[\[Theta]] \[Sigma]z + x Sin[\[Theta]] \[Sigma]x;

sol=Quiet@Solve[{Det[Cq-X]==0,Det[X-Inverse[Cp]]==0},{x0,x}];
{det1,det2} = X[[1,1]] Inverse[X][[1,1]]/.sol//Simplify;

tab1 =Select[ Table[{\[Theta],det1},{\[Theta],linspace[10^-5,\[Pi]-10^-5,1000]}],Im[#[[2]]]==0.&];
opt1=First@MinimalBy[tab1,Last];

tab2 =Select[ Table[{\[Theta],det2},{\[Theta],linspace[10^-5,\[Pi]-10^-5,1000]}],Im[#[[2]]]==0&];
(*Print[tab2];*)
opt2=First@MinimalBy[tab2,Last];

(*Print[{opt1,opt2}];*)
If[opt1[[2]] < opt2[[2]],
(
{\[Theta]min,detmin} = opt1; 
Xmin = X/.sol[[1]]/.\[Theta]->\[Theta]min
),
(
{\[Theta]min,detmin} = opt2; 
Xmin = X/.sol[[2]]/.\[Theta]->\[Theta]min
)
];

\[Gamma]p = L\[Transpose] . ArrayFlatten[({
 {Xmin, 0},
 {0, Inverse[Xmin]}
})] . L//Chop;

{1/2 Log[Abs@detmin](*detmin*),\[Gamma]p/2}//Chop
(* 1st term has no 2 inside det because we were working with 2\[Sigma]*)
(* 2nd term is divided by 1/2 to get back to vacuum being 1/2 *)

]


(* ::Input::Initialization:: *)
Clear[SqueezedThermalState];
SqueezedThermalState[nmax_,\[ScriptCapitalN]_,r_,\[Theta]_:0]:=Module[{\[Beta],Z,b},
If[nmax>0,LoadBosonicOperators[nmax]];
\[Beta] =Log[(1+\[ScriptCapitalN])/\[ScriptCapitalN]];
Z=E^\[Beta]/(-1+E^\[Beta]);
b = a Cosh[r] + a\[ConjugateTranspose] E^(I \[Theta]) Sinh[r];
1/Z MatrixExp[-\[Beta] b\[ConjugateTranspose] . b]
]

(* Technical note: if input is nmax<0, the function will not load the bosonic operators. *)
(*This is useful when this function is being used within other calculations, which already used LoadBosonicOperators[nmax] *)


(* ::Input::Initialization:: *)
(* This is only a part of the SLD actually. It is the part called L^(2) in 1303.3682 *)
Clear[GaussianSLD]
GaussianSLD[\[Sigma]_,d\[Sigma]_]:=Module[{\[CapitalOmega]},
\[CapitalOmega] = SymplecticForm[Length[\[Sigma]]/2];
Unvec@LinearSolve[2 kron[\[Sigma]\[Transpose],\[Sigma]] + 1/2 kron[\[CapitalOmega]\[Transpose],\[CapitalOmega]],Vec[d\[Sigma]]]
]


(* ::Input::Initialization:: *)
GaussianQFI[\[Sigma]_,d\[Sigma]_]:=Module[{W,S,w,errors,\[CapitalOmega],M,vec,\[ScriptCapitalD]},
{W,S,errors}=Williamson[2*\[Sigma]];
\[CapitalOmega] = SymplecticForm[Length[W]/2];
M = \[CapitalOmega] . S\[Transpose] . \[CapitalOmega];
vec = Vec[M . d\[Sigma] . Transpose[M]];

w = SparseArray[W];
\[ScriptCapitalD] = 1/2 (kron[w,w]-kron[\[CapitalOmega],\[CapitalOmega]]);
vec\[Conjugate] . LinearSolve[\[ScriptCapitalD],vec]
]

GaussianQFI[\[Sigma]_,d\[Sigma]_,d\[Mu]_]:=d\[Mu] . LinearSolve[\[Sigma],d\[Mu]] + GaussianQFI[\[Sigma],d\[Sigma]]


(* ::Input::Initialization:: *)
FermionicFullCM[\[ScriptCapitalC]_,S_:0]:=Module[{L = Length[\[ScriptCapitalC]],vp={1,I},vm={1,-I}},
-kron[\[Sigma]y,Eye[L]] + kron[out[vp,vp],\[ScriptCapitalC]\[Transpose]]-kron[out[vm,vm],\[ScriptCapitalC]] +If[MatrixQ[S], kron[out[vm,vp],S]+kron[out[vp,vm],S\[ConjugateTranspose]],0]
]


(* ::Input::Initialization:: *)
(* Recovers the reduced (L\[Times]L) covariance matrix from the full (2L\[Times]2L) one *)
(* Still need to implement the local version *)
Clear[FermionsReducedC];
FermionsReducedC[\[CapitalTheta]_,ord_:"global"]:=Module[{L=Length[\[CapitalTheta]]/2},
If[ord==="global",
-(1/4)PTr[kron[\!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{"1", "I"},
{
RowBox[{"-", "I"}], "1"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\),Eye[L]] . (\[CapitalTheta]+kron[\[Sigma]y,Eye[L]]),{1},{2,L}]]
];

Clear[FermionsReducedS];
FermionsReducedS[\[CapitalTheta]_,ord_:"global"]:=Module[{L=Length[\[CapitalTheta]]/2},
If[ord==="global",
1/4 PTr[kron[\!\(\*
TagBox[
TagBox[
RowBox[{"(", "", GridBox[{
{"1", "I"},
{"I", 
RowBox[{"-", "1"}]}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "Rows" -> {{Baseline}}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.7]}, Offset[0.27999999999999997`]}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}}], "", ")"}],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]],
Function[BoxForm`e$, MatrixForm[BoxForm`e$]]]\),Eye[L]] . (\[CapitalTheta]+kron[\[Sigma]y,Eye[L]]),{1},{2,L}]]
];


(* ::Input::Initialization:: *)
Clear[grid];
grid[plots_,nrows_:1,lettpos_:{0.92,0.9},tight_:True]:=Module[{i,p2,LettList,tab,pads,L,R,T,B,ncols,imgsize},

ncols = Length[plots]/nrows;
imgsize=Max/@Transpose@Table[ImageSize/.AbsoluteOptions[p,ImageSize],{p,plots}];
tab=Table[Flatten[ImagePadding/.AbsoluteOptions[p,ImagePadding]],{p,plots}];
{L,R,B,T}=tab\[Transpose];
B = Max[B];
T = Max[T];
R = Max[R];

If[tight===False,
L = ConstantArray[Max[L],Length@plots],
(L = Transpose@Partition[L,ncols];
L = Flatten@Transpose@Table[ConstantArray[Max[L[[i]]],nrows],{i,ncols}])];

LettList  = "("<>#<>")"&/@CharacterRange["a","z"];
p2 = Table[Show[plots[[i]],
ImageSize->{Automatic,imgsize[[2]]},
ImagePadding->{{L[[i]],R},{B,T}},
Epilog->Text[LettList[[i]],Scaled[lettpos]]
],{i,Length@plots}];
Grid@PartitionIn[p2,nrows]
]




(* ::Input::Initialization:: *)
LoadPlotStyling[]:=Module[{},
Clear[LettList,blue,blue2,blue3,navy,yellow,purple,purple2,green,red,orange,Custom,opts];

(*LettList  = "("<>#<>")"&/@CharacterRange["a","z"];*)
blue = RGBColor[0.368417, 0.506779, 0.709798];
blue2 = RGBColor[68/256,151/256,151/256];
blue3 = RGBColor[205/256,205/256,222/256];
navy = RGBColor[{65,97,150}/255];

yellow=RGBColor[1, 0.75, 0];
purple = RGBColor[{158,137,183}/255];
purple2 = RGBColor[158/256,67/256,149/256];
green = Darker[RGBColor[{123,216,187}/255],.1];
red=RGBColor[{238,105,105}/255];
orange = RGBColor[224/256,148/256,57/256];
graphite = RGBColor[{95,95,95}/255];

(*nicePalette={Black,purple,navy,green,graphite,red};*)

opts = {
BaseStyle->{FontFamily->"Times",17},
PlotStyle->Automatic,
Axes->False,
ImageSize->330,
GridLinesStyle->Directive[Black,Dashed],
PlotRangePadding->None,
Frame->True,
FrameStyle->Black,
FrameTicksStyle->{FontFamily->"Times"(*,Black*)}
};
SetOptions[ListLinePlot,opts];
SetOptions[ListPlot,opts];
SetOptions[ListLogPlot,opts];
SetOptions[ListLogLinearPlot,opts];
SetOptions[ListLogLogPlot,opts];
SetOptions[Plot,opts];
SetOptions[LogPlot,opts];
SetOptions[LogLinearPlot,opts];
SetOptions[LogLogPlot,opts];
SetOptions[ParametricPlot,opts];
SetOptions[DiscretePlot,opts];
(*SetOptions[Graphics,opts];*)
];


(* ::Input::Initialization:: *)
Clear[legend];
legend[labels_,pos_,OptionsPattern[{
FontSize->16,
LegendMarkerSize->20,
LegendLayout->"Column",
framed->True,
BackgroundColor->Gray,
FrameColor->Gray
}]]:=Placed[LineLegend[Automatic,labels,
LabelStyle->{
FontFamily->"Times",
FontSize->OptionValue[FontSize]},
LegendLayout->OptionValue[LegendLayout],
LegendMarkerSize->{OptionValue[LegendMarkerSize],2},LegendFunction->If[OptionValue[framed],(Framed[#,RoundingRadius->4,FrameStyle->OptionValue[FrameColor],Background->Lighter[OptionValue[BackgroundColor],0.9]]&),None]],Scaled[pos]]


(* ::Input::Initialization:: *)
Clear[pointlegend];
pointlegend[labels_,pos_,OptionsPattern[{
FontSize->16,
LegendMarkerSize->10,
LegendLayout->"Column",
BackgroundColor->Gray,
FrameColor->Gray
}]]:=Placed[PointLegend[Automatic,labels,
LabelStyle->{
FontFamily->"Times",
FontSize->OptionValue[FontSize]},
Joined->False,
LegendLayout->OptionValue[LegendLayout],
LegendMarkerSize->OptionValue[LegendMarkerSize],LegendFunction->(Framed[#,RoundingRadius->4,FrameStyle->OptionValue[FrameColor],Background->Lighter[OptionValue[BackgroundColor],0.9]]&)],Scaled[pos]]


(* ::Input::Initialization:: *)
(* Ticks for logscale plots *)
(* Output appropriate for functions like FrameTicks or Ticks *)
(*logticks[i_,f_,\[Delta]_,log_:False]:=
Table[{k,Superscript["10",ToString[k]]},{k,i,f,\[Delta]}] *)

Clear[logticks];
logticks[i_,f_,\[Delta]_,ts_:0.8,label_:1]:=
Flatten[Table[{
{10^k,If[label==1,Superscript["10",ToString[k]],""],ts{0.026,0}},
{0.9 10^k,"",ts{0.012,0}},
{0.8 10^k,"",ts{0.011,0}},
{0.7 10^k,"",ts{0.010,0}},
{0.6 10^k,"",ts{0.009,0}},
{0.5 10^k,"",ts{0.008,0}},
{0.4 10^k,"",ts{0.007,0}},
{0.3 10^k,"",ts{0.006,0}},
{0.2 10^k,"",ts{0.005,0}}
},{k,i,f,\[Delta]}],1]


(* ::Input::Initialization:: *)
rasta[plot_,res_:100,sz_:400]:=Rasterize[plot,ImageResolution->res,ImageSize->sz];


(* ::Input::Initialization:: *)
BlochSphere[labels_:True]:=Module[{splineCircle,pointsAndConnection,surroundingCircles,texKet},
(*Needs["MaTeX`"];*)

splineCircle[m_List,r_,angles_List: {0,2 \[Pi]}]:=Module[{seg,\[Phi],start,end,pts,w,k},{start,end}=Mod[angles//N,2 \[Pi]];
If[end<=start,end+=2 \[Pi]];
seg=Quotient[end-start//N,\[Pi]/2];
\[Phi]=Mod[end-start//N,\[Pi]/2];
If[seg==4,seg=3;\[Phi]=\[Pi]/2];
pts=r RotationMatrix[start] . #&/@Join[Take[{{1,0},{1,1},{0,1},{-1,1},{-1,0},{-1,-1},{0,-1}},2 seg+1],RotationMatrix[seg \[Pi]/2] . #&/@{{1,Tan[\[Phi]/2]},{Cos[\[Phi]],Sin[\[Phi]]}}];
If[Length[m]==2,pts=m+#&/@pts,pts=m+#&/@Transpose[Append[Transpose[pts],ConstantArray[0,Length[pts]]]]];
w=Join[Take[{1,1/Sqrt[2],1,1/Sqrt[2],1,1/Sqrt[2],1},2 seg+1],{Cos[\[Phi]/2],1}];
k=Join[{0,0,0},Riffle[#,#]&@Range[seg+1],{seg+1}];
BSplineCurve[pts,SplineDegree->2,SplineKnots->k,SplineWeights->w]]/;Length[m]==2||Length[m]==3;

pointsAndConnection[points_]:=Sequence@@{Sequence@@Point/@#,Line@#}&@points;

surroundingCircles=GeometricTransformation[splineCircle[{0,0,0},1],{{RotationMatrix[0,{1,0,0}],{0,0,0}},{RotationMatrix[Pi/2,{1,0,0}],{0,0,0}},{RotationMatrix[Pi/2,{0,1,0}],{0,0,0}}}];

texKet[n_]:=Text@Style[StringTemplate["\!\(\*TemplateBox[{\"`1`\"},\n\"Ket\"]\)"][ToString@n] ,FontFamily->"Times",20];

Graphics3D[{White,Opacity@0.3,Sphere[{0,0,0},1],Opacity@1,Thickness@0.004,PointSize@0.02,Red,pointsAndConnection@{{0,0,1},{0,0,-1}},Blue,pointsAndConnection@{{1,0,0},{-1,0,0}},Darker@Green,pointsAndConnection@{{0,1,0},{0,-1,0}},Black,Point[{0,0,0}],
If[labels==True,
{Text[texKet[0],{0,0,1.2}],Text[texKet[1],{0,0,-1.2}],Text[texKet["+"],{1.2,0,0}],Text[texKet["-"],{-1.2,0,0}],Text[texKet["L"],{0,1.2,0}],Text[texKet["R"],{0,-1.2,0}]}],Gray,Thin,surroundingCircles},Boxed->False,PlotRange->ConstantArray[{-1.2,1.2},3],ImageSize->500,RotationAction->"Clip"]
]



(* ::Input::Initialization:: *)
Clear[COLORLIST];
COLORLIST=Join[ColorData[3,"ColorList"][[{1,2,4,6,8,9}]],ColorData[7,"ColorList"][[{-1}]],ColorData[4,"ColorList"][[{6}]],ColorData[3,"ColorList"][[{5}]],ColorData[17,"ColorList"][[{-1}]]]


(* ::Input::Initialization:: *)
(* Can also be adapted to show all defined functions *)
Clear[ListAllVariables]
ListAllVariables[]:=Module[{names = Names["Global`*"],variables,functions},
variables = Select[names,ToExpression[#,StandardForm,Function[sym,OwnValues[sym]=!={},HoldAll]]&];
functions = Complement[names,variables];
(*{variables,functions}*)
variables
]


(* ::Input::Initialization:: *)
matchString[str1_,str2_]:=StringMatchQ[str1,str2,SpellingCorrection->True]


(* ::Input::Initialization:: *)
matex[symbol_]:=MaTeX[symbol,FontSize->18,"Preamble"->{"\\usepackage{bm}"}];


(* ::Input::Initialization:: *)
LoadDataScienceFunctions[]:=Module[{},

(* Allows combining SortBy and ReverseSortBy in the same search entry *)
Clear[order,rev];
order[rev[x_],rev[y_]]:=-Order[x,y];
order[x_,y_]:=Order[x,y];

]



