skipPts=Max[1,Round[maxTime/5000]]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)

coarseData=Map[#[[2,1;;maxTime;;skipPts]]&,selectedData];

labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];

plotTimeEvol=ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},Frame->True,FrameLabel->{"timesteps","number of triangles"},ImageSize->300]

(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)

Export[NotebookDirectory[]<>"plots/timeevol.pdf",plotTimeEvol]

movingAvg=Map[MovingAverage[#[[2]],2000][[1;;-1;;skipPts]]-Mean[#[[2,-20000;;-1]]]&,selectedData[[1;;-1;;5]]];

plotMovingAvg=ListPlot[movingAvg,Joined->True,PlotRange->All,DataRange->{0,measureSkip*maxTime},Frame->True,FrameLabel->{"timesteps","number of triangles"}]

(* ::Subsection:: *)

(*Fit exponential to triangles time evolution*)

fitList=Map[NonlinearModelFit[#[[2]],Exp[-(t-t0)/tmix]+c,{{tmix,1000},{t0,10000},{c,2000}},t]&,selectedData];

(* tmix*Log[1/epsilon] gives the time it takes to get a factor epsilon close to the average *)

(* t0 gives the time it takes to be exactly 1 triangle (in absolute value) away from the average *)

(* Use fit["BestFitParameters"] to get parameters *)

(* Use fit[t] to get fit value *)

fitFuncsT=Map[#[t]&,fitList];

timeplot1=ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},PlotStyle->Opacity[0.5]];

Show[timeplot1,Plot[fitFuncsT,{t,1,maxTime},PlotRange->All]]

(* ::Subsection:: *)

(*Plot success rate over "time"*)

numPlots=10;

selectedData=gdata[[1,-1]][[-numPlots;;-1]];

measureSkip=100;

maxTime=Max[Map[Length[#[[4]]]&,selectedData]];

maxTime=10000;

coarseData=Map[#[[4,1;;maxTime/measureSkip]]&,selectedData];

labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];

ListPlot[coarseData,Joined->True,PlotRange->{0,100},DataRange->{0,maxTime},PlotLegends->labels]

(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)

(* ::Subsection:: *)

(*Correlation of avgsuccess rate vs other things*)

compare1=Map[{Mean[#[[4]]],Mean[#[[2]]]}&,gdata,{3}];

(* { GCM1 rate, GCM2 rate, mixing time from ErdosGallai } *)

scatterPlots=Map[ListPlot[#,PlotRange->{{0,100},All},PlotStyle->PointSize[Large]]&,compare1,{2}];

plot1a=Show[ListLogLogPlot[averagesErrorBars[[All,All,1]],Joined->False,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels],Plot[averagesFits,{logn,1,2000}]]

plot1b=Show[ListLogLogPlot[mediansErrorBars[[All,All,1]],Joined->False,PlotMarkers->Automatic,AxesLabel->{"n","median triangles"},PlotLegends->taulabels],Plot[mediansFits,{logn,1,2000}]]

Export[NotebookDirectory[]<>"plots/avgtris_n.pdf",plot1]

(* ::Subsection:: *)

(*Plot of T(\[Tau])*)

tauValues=averagesGrouped[[All,1,1,1,2]];

averagesExponents=Transpose[{tauValues,averagesFits[[All,2,1]]}];

mediansExponents=Transpose[{tauValues,mediansFits[[All,2,1]]}];

Show[ListPlot[{averagesExponents,mediansExponents},Joined->True,PlotMarkers->Automatic,AxesLabel->{"tau","exponent"},PlotRange->{{2,3},{0,1.6}}],Plot[3/2(3-tau),{tau,2,3}]]

(* ::Subsection:: *)

(*T(\[Tau]) including error bars*)

(* For visual, shift the tau values slightly left or right to distinguish the two datasets *)

tauValues=averagesGrouped[[All,1,1,1,2]];

averagesExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,

mediansExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,

Export[NotebookDirectory[]<>"plots/triangle_exponent.pdf",plot2]