(* ::Package:: *)

Needs["ErrorBarPlots`"]


gsraw=Import[NotebookDirectory[]<>"data/graphdata_etmt.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: etmt
*)
tauvalues=gdata[[All,1,1,1,2]];
nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
taulabels=Map["\[Tau] = "<>ToString[#]&,gdata[[All,1,1,1,2]]];


etmtMean=Map[Mean[N[#[[All,2]]]]&,gdata,{2}];
etmtSD=Map[StandardDeviation[N[#[[All,2]]]]&,gdata,{2}];
etmtQuantile=Map[Quantile[#[[All,2]],99/100]&,gdata,{2}];


histograms=Map[Histogram[#[[All,2]]]&,gdata,{2}];


histogramsWithLines=MapIndexed[
Histogram[#[[All,2]],PlotRange->All,
Epilog->Line[{
{{etmtMean[[#2/.List->Sequence]],0},
{etmtMean[[#2/.List->Sequence]],500}},
{{etmtMean[[#2/.List->Sequence]]+etmtSD[[#2/.List->Sequence]],0},
{etmtMean[[#2/.List->Sequence]]+etmtSD[[#2/.List->Sequence]],500}},
{{etmtQuantile[[#2/.List->Sequence]],0},
{etmtQuantile[[#2/.List->Sequence]],500}}
}]
]
&,gdata,{2}];


TableForm[histogramsWithLines,TableHeadings->{taulabels,nlabels}]


gdataSwitched=Transpose[gdata];
tauChoices={1,3,5,7,9};
combiHistograms=Map[SmoothHistogram[#[[tauChoices,All,2]],PlotRange->All,PlotLegends->taulabels[[tauChoices]],PlotLabel->"n = "<>ToString[#[[1,1,1,1]]]]&,gdataSwitched,{1}];


combiHistograms[[9]]


Export[NotebookDirectory[]<>"plots/ETMTdistribution.pdf",combiHistograms[[9]]]


mixingTimesBars=Map[{{#[[1,1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&,gdata,{2}];
mixingTimesQuantiles=Map[{#[[1,1,1]],Quantile[#[[All,2]],99/100]}&,gdata,{2}];


plot1=ErrorListPlot[mixingTimesBars[[{1,2,3,5,8}]],Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"n","ETMT"},PlotLegends->taulabels];
plot2=ListPlot[mixingTimesQuantiles[[{1,2,3,5,8}]],Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"n","ETMT"},PlotLegends->taulabels];
Show[plot2,plot1]


mixingTimesDivN=Map[{#[[1,1]],#[[1,2]]/(#[[1,1]])}&,mixingTimesBars,{2}];
mixingTimesDivNlogN=Map[{#[[1,1]],#[[1,2]]/(#[[1,1]]*Log[#[[1,1]]])}&,mixingTimesBars,{2}];
etmtQuantileDivN=Map[{#[[1]],#[[2]]/(#[[1]])}&,mixingTimesQuantiles,{2}];
etmtQuantileDivNlogN=Map[{#[[1]],#[[2]]/(#[[1]]*Log[#[[1]]])}&,mixingTimesQuantiles,{2}];


plotN=ListPlot[mixingTimesDivN,Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"n","\[LeftAngleBracket]ETMT\[RightAngleBracket]/n"},PlotLegends->taulabels,ImageSize->300]
plotNlogN=ListPlot[mixingTimesDivNlogN,Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"n","\[LeftAngleBracket]ETMT\[RightAngleBracket]/(n log n)"},PlotLegends->taulabels,ImageSize->300]
plotQuantileN=ListPlot[etmtQuantileDivN,Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"n","q(ETMT,99%)/n"},PlotLegends->taulabels,ImageSize->300]
plotQuantileNlogN=ListPlot[etmtQuantileDivNlogN,Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"n","q(ETMT,99%)/(n log n)"},PlotLegends->taulabels,ImageSize->300]


Export[NotebookDirectory[]<>"plots/ETMTdivN.pdf",plotN]


etmtQuantileDivNmax=Map[Max[#[[All,2]]]&,etmtQuantileDivN];
etmtQuantileDivNmax=Transpose[{tauvalues,etmtQuantileDivNmax}];

mixingTimesDivNmax=Map[Max[#[[All,2]]]&,mixingTimesDivN];
mixingTimesDivNmax=Transpose[{tauvalues,mixingTimesDivNmax}];


Show[
Plot[{(50-30(tau-2)),32-26(tau-2)},{tau,2,3},AxesOrigin->{2,0}],
ListPlot[{etmtQuantileDivNmax,mixingTimesDivNmax},Joined->True,PlotMarkers->Automatic,PlotRange->{{2,3},All}]
]