SetDirectory @ NotebookDirectory[];
Import["https://qtechtheory.org/QuESTlink.m"];
CreateDownloadedQuESTEnv[];

freqs={-0.7098301012032957`,-0.05177915010485945`,0.9065658350752286`,-0.9265872629174545`,0.09500423337423802`,-0.4959779888875864`,-0.16859651542448084`,0.8678962526289684`,-0.11732544489533625`,0.907000173126657`,-0.44233758592927064`,0.5741970506048069`,0.06929812708367722`,0.23796927148578906`,-0.9361278713972463`,0.4908133930451908`,0.9021460849017751`,0.7153291490064708`,-0.309302916408313`,-0.38233264968765734`};

spinChainHamil[J_,numQs_]:=Module[{hamil},
hamil=ExpandAll@Sum[J*((Subscript[X, k])(Subscript[X, k+1])+(Subscript[Y, k])(Subscript[Y, k+1])+(Subscript[Z, k])(Subscript[Z, k+1])),{k,0,numQs-2}]+ExpandAll[J*((Subscript[X, 0])(Subscript[X, numQs-1])+(Subscript[Y, 0])(Subscript[Y, numQs-1])+(Subscript[Z, 0])(Subscript[Z, numQs-1]))];
hamil+=Sum[freqs[[k]]*(Subscript[Z, k-1]),{k,1,numQs}];
Return@hamil;
];

calculate symbolic Hamiltonian

numQs=6;
hamil=Chop@spinChainHamil[.1,numQs]

code for generating all 3-local Paulis

pauli[type_, site_]:=(Subscript[type, site])
getListOf3LocalPaulis[numQs_]:=Join[
Flatten[Table[pauli[type, site],{site,0,numQs-1},{type,{X,Y,Z}}],1]
,
Flatten[
Table[pauli[type1, site1]pauli[type2, site2],{site1,0,numQs-1},{site2,site1+1,numQs-1},{type1,{X,Y,Z}},{type2,{X,Y,Z}}]
,3]
,
Flatten[
Table[pauli[type1, site1]pauli[type2, site2]pauli[type3, site3],{site1,0,numQs-1},{site2,site1+1,numQs-1},{site3,site2+1,numQs-1},{type1,{X,Y,Z}},{type2,{X,Y,Z}},{type3,{X,Y,Z}}]
,5]
]

getListOf2LocalPaulis[numQs_]:=Join[
Flatten[Table[{pauli[type, site]},{site,0,numQs-1},{type,{X,Y,Z}}],1]
,
Flatten[
Table[{pauli[type1, site1],pauli[type2, site2]},{site1,0,numQs-1},{site2,site1+1,numQs-1},{type1,{X,Y,Z}},{type2,{X,Y,Z}}]
,3]
]
Length@getListOf3LocalPaulis[numQs]

calculate Hamiltonian matrix and dt time-step propagator

Hmat=CalcPauliExpressionMatrix[hamil];
(* use all up to 3-local Paulis *)
listOfPaulis=getListOf3LocalPaulis[numQs];
eigs=Eigenvalues[Hmat,4];
vecs=Eigenvectors[Hmat,4];
dt=3;
Umat=MatrixExp[-I*dt*Hmat];
ll=200;(* number of total time steps *)
TT=ll*dt; (* total simulation time *)
Print["eigs:  ", eigs[[1;;3]]]

simulate dynamics and calculate expected values

psinull=(vecs[[1]]+0.1vecs[[2]]+0.1vecs[[3]]);psinull=psinull/Norm[psinull];
psi=psinull;
{ψ,ϕ}=CreateQuregs[numQs,2];
observables={};
observablesSHADOW={};
Do[
psi=Umat.psi;
SetQuregMatrix[ψ,psi];


AppendTo[observablesSHADOW,
(* these are built-in QuEST-link functions for simulating shadows *)
data = SampleClassicalShadow[ψ,1000];
CalcExpecPauliProdsFromClassicalShadow[data,listOfPaulis,3]
];

If[Mod[k,ll/50]==0,
Print["progress = ", N@k/ll];
(*Export["observables.m",Compress@observables];*)
Export["observablesSHADOW.m",Compress@observablesSHADOW];
];

,{k,1,ll}]

plot data -- row index is time variable, column index is observable

SetDirectory[NotebookDirectory[]];
observablesSHADOW=Uncompress@Import["observablesSHADOW.m"];
ll=Dimensions[observablesSHADOW][[1]];
dt=3;TT=ll*dt; 
observablesSHADOW=Transpose@observablesSHADOW;
(*  above is the non-square matrix that holds expected values as a function of time stepss  *)

observablesSHADOW=Standardize/@observablesSHADOW;
plotD=Rasterize[MatrixPlot[Transpose@observablesSHADOW,AspectRatio->1/2],RasterSize->{Automatic,1000}, ImageSize->{Automatic,200}]
(* With enough classical shadow data the column dimension should be much larger than the row dimension *)

perform a Ljung-box pre-screening of observables

dd=observablesSHADOW;
pValues=(AutocorrelationTest[#,ll-1,"LjungBox"]&)/@dd;

(* take only the best NMAX signals *)
NMAX=100;
ListLinePlot[Sort[pValues][[1;;NMAX]],PlotRange->All]
pos=Ordering[pValues][[1;;NMAX]];
data=observablesSHADOW[[pos]];

calculate the quadratic “covariance matrix”

(*data=observablesSHADOW;*)
covariance=(Transpose[data].data)/(Length@observablesSHADOW);
Dimensions@covariance
plotC=Rasterize[MatrixPlot[covariance,AspectRatio->1],RasterSize->{Automatic,1000}, ImageSize->{Automatic,200}]

plot decaying eigenvalues of “covariance matrix” and calculate eigenvectors

ee=Eigenvalues[covariance];
ListPlot[Chop@ee,PlotMarkers->Automatic,Joined->True,PlotRange->{{0,50},All}]
eigvs=Eigenvectors[covariance,4];

Calculate spectral cross-correlation between the different eigenvectors

(* define window function *)
win=HannWindow[Range[-0.5,0.5-1/ll,1/ll]]

xcorr[x_,y_]:=Module[{ll},
(*see def at https://www.mathworks.com/help/matlab/ref/xcorr.html *)
(*sx=Standardize[x];
sy=Standardize[y];*)
ll=Length@x;
xc=Table[Mean[Drop[x,-m]*Drop[y,m]],{m,1,ll-1}];
Return@xc;
];
llH=ll-1; dt=3;TTH=llH*dt;
win=HannWindow[Range[-0.5,0.5-1/llH,1/llH]];
data=Table[
Fourier[win*xcorr[eigvs[[k]],eigvs[[l]]]]
,{k,1,Length@eigvs},{l,1,Length@eigvs}];

solution=Table[
aa=Chop@data[[;;,;;,k]];
Max@SingularValueList[aa]
,{k,1,llH}];

spectrum=ListLinePlot[solution,PlotRange->{{0,2π*llH/(TTH)/2},All}
,PlotStyle->Black
,DataRange->{0,2π*llH/(TTH)}
,GridLines->{{Abs[eigs[[1]]-eigs[[2]]],
Abs[eigs[[1]]-eigs[[3]]],
Abs[eigs[[2]]-eigs[[3]]]},None}
,GridLinesStyle->Directive[Gray, Dashed,Thick]
, ImageSize->{Automatic,200}]

pos=Reverse[Ordering[Abs@solution]][[1]];
N@(llH-pos)/llH(2π*llH/(TTH))
N@(pos)/llH(2π*llH/(TTH))
Abs[eigs[[1]]-eigs[[2]]]
Abs[eigs[[1]]-eigs[[3]]]