numPlots=20;

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

measureSkip=1;

minCount=Min[Map[Min[#[[2]]]&,selectedData]];

maxCount=Max[Map[Max[#[[2]]]&,selectedData]];

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

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

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

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

ListPlot[coarseData,Joined->True,PlotRange->{minCount,maxCount},DataRange->{0,measureSkip*maxTime},PlotLegends->labels]

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

(* ::Subsection:: *)

(*Compute 'mixing time'*)

(* Compute average of last part and check when the value drops below that for the first time *)

getMixingTime[values_]:=Module[{avg},

    (* average over the last 20 percent *)

    avg=Mean[values[[-Round[Length[values]/5];;-1]]];

    FirstPosition[values,_?(#<=avg&)][[1]]

]

(* gdata[[ tau index, n index, run index , {ntau, #tris, ds} ]] *)

measureSkip=1;

mixingTimes=Map[{#[[1,1]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]]}&,gdata,{3}];

mixingTimesBars=Map[

    {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]/Sqrt[Length[#]]]}&

,mixingTimes,{2}];

ErrorListPlot[mixingTimesBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","~~mixing time divided by n"},PlotLegends->taulabels]

(* For n fixed, look at function of tau *)

measureSkip=1;

mixingTimes=Map[{#[[1,2]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]]}&,gdata,{3}];

mixingTimesBars=Map[

    {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&

,mixingTimes[[All,-1]],{1}];

Show[

ErrorListPlot[mixingTimesBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"tau","~~mixing time divided by n, for n = 1000"},PlotRange->{{2,3},{0,30}}]

,Plot[(32-26(tau-2)),{tau,2,3}]]

(* ::Subsection:: *)

averages=Map[{#[[1]],Mean[#[[2,1;;-1]]]}&,gsraw];

(* averages=SortBy[averages,#[[1,1]]&]; (* Sort by n *) *)

averagesGrouped=GatherBy[averages,{#[[1,2]]&,#[[1,1]]&}]; (* Split by n,tau *)

(* averagesGrouped[[ tau index, n index, run index , {ntau, avgtri} ]] *)

nlabels=Map["n = "<>ToString[#]&,averagesGrouped[[1,All,1,1,1]]];

taulabels=Map["tau = "<>ToString[#]&,averagesGrouped[[All,1,1,1,2]]];

averagesErrorBars=Map[

{{#[[1,1,1]],Mean[#[[All,2]]]},

ErrorBar[StandardDeviation[#[[All,2]]]]

}&,averagesGrouped,{2}];

ErrorListPlot[averagesErrorBars,Joined->True,PlotMarkers->Automatic,PlotRange->All,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels]

ListLogLogPlot[averagesErrorBars[[All,All,1]],Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels]

(* ::Subsection:: *)

(*Fitting the log-log-plot*)

loglogdata=Log[averagesErrorBars[[All,All,1]]];

fits=Map[Fit[#,{1,logn},logn]&,loglogdata];

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

exponents=Transpose[{tauValues,fits[[All,2,1]]}];

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

(* ::Subsection:: *)

(*Plot #triangles distribution for specific (n,tau)*)

plotRangeByTau[tau_]:=Piecewise[{{{0,30000},tau<2.3},{{0,4000},2.3<tau<2.7},{{0,800},tau>2.7}},Automatic]

histograms=Map[Histogram[#[[All,2]],PlotRange->{plotRangeByTau[#[[1,1,2]]],Automatic}]&,averagesGrouped,{2}];

(* TableForm[histograms,TableHeadings->{taulabels,nlabels}] *)

TableForm[Transpose[histograms],TableHeadings->{nlabels,taulabels}]

(* ::Section:: *)

(*Greedy configuration model*)

(* ::Subsection:: *)

(*#triangles(GCM) distribution vs #triangles(SwitchChain)*)

timeWindow=Round[Length[gdata[[1,1,1,2]]]/10];

getStats[run_]:=Module[{avg,stddev},

    avg=N[Mean[run[[2,-timeWindow;;-1]]]];

    stddev=N[StandardDeviation[run[[2,timeWindow;;-1]]]];

    {run[[1]],stddev/avg,(run[[2,1]])/avg,Map[N[#/avg]&,run[[4]]]}

]

stats=Map[getStats,gdata,{3}];

histograms=Map[Histogram[{#[[1,4]]},PlotRange->{{0,2},Automatic},PlotLabel->"ErdosGallai="<>ToString[NumberForm[#[[1,3]],3]]<>"\[Cross]average. stddev="<>ToString[NumberForm[#[[1,2]],3]]<>"\[Cross]average"]&,stats,{2}];

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

(* ::Subsection:: *)

(*Greedy CM success rates*)

(* gdata[[ tau index, n index, run index , {ntau, #tris, ds, greedyTriangles} ]] *)

successrates=Map[{Length[#[[4]]],Length[#[[5]]]}&,gdata,{3}];

successrates=Map[Transpose,successrates,{2}];

successratesDelta=Map[Length[#[[5]]]-Length[#[[4]]]&,gdata,{3}];

rateHistograms=Map[Histogram[#,{10},PlotRange->{{0,100},Automatic}]&,successrates,{2}];

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

rateHistograms=Map[Histogram[#,{10},PlotRange->{{-100,100},Automatic}]&,successratesDelta,{2}];