(* ::Package:: *)

(* ::Section:: *)
(*Data import*)


gsraw=Import[NotebookDirectory[]<>"data/graphdata_timeevol.m"];
(* gsraw=SortBy[gsraw,{#[[1,1]]&,#[[1,2]]&}]; (* Sort by n and then by tau. The {} forces a *stable* sort because otherwise Mathematica sorts also on triangle count and other things. *) *)


gdata=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];
(* Data format: *)
(* gdata[[ tau index, n index, run index , datatype index ]] *)
(* datatype index:
1: {n,tau}
2: #triangles time sequence
3: degree sequence
4: GCM1 starting triangle counts
5: GCM2 starting triangle counts
6: GCM1 time sequence
7: GCM2 time sequence
*)
nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
taulabels=Map["tau = "<>ToString[#]&,gdata[[All,1,1,1,2]]];


(* ::Section:: *)
(*Triangle counts autocorrelation*)


doAutoCorrelate[dataset_,maxlag_]:=Module[{avg,data2,result},
avg=Mean[dataset];
data2=dataset-avg;
result=ParallelTable[Mean[data2[[1+lag;;]]*data2[[;;-1-lag]]],{lag,0,maxlag}];
result=result/result[[1]]
]


doPeriodicAutocorrelate[dataset_,maxlag_]:=Module[{f,result},
f=Fourier[dataset-Mean[dataset]];
result=InverseFourier[f*Conjugate[f]];
result=result[[1;;Min[Floor[Length[result]/2],maxlag]]];
result=Chop[result];
result=result/result[[1]]
]


dataX=gdata[[1,1]][[3,2]];
ListPlot[dataX,Joined->True,PlotRange->All]
dataX=dataX[[15000;;50000]];


cor1=doAutoCorrelate[dataX,2000];


cor2=doPeriodicAutocorrelate[dataX,20000];


ListPlot[{cor1,cor2},Joined->True,PlotRange->{-0.5,1}]