gsraw=Import[NotebookDirectory[]<>"data/graphdata_ecm_initialtris.m"];

gsraw2=Import[NotebookDirectory[]<>"data/graphdata_ecm_initialtris_extra.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]]&}];

(* Data format: *)

(* gdata[[ tau index, n index, datatype index ]] *)

(* datatype index:

1: {n,tau}

2: {uniform triangle samples}

3: {ECM triangle samples}

*)

(* ::Section:: *)

(*Erased configuration model*)

(* ::Subsection:: *)

(*Distribution of initial #triangles for ECM compared to uniform triangle distribution*)

getHistogram[run_]:=Histogram[{run[[2]],run[[3]]},Automatic,"Probability",

ChartLegends->Placed[{"Uniform","ECM"},Bottom],

ImageSize->250,

Frame->True,

FrameLabel->{"Triangles","Probability"},

PlotLabel->("n = "<>ToString[run[[1,1]]]<>", \[Tau] = "<>ToString[run[[1,2]]])

];

histograms=Map[getHistogram,gdata,{2}];

TableForm[histograms]

(* ::Subsubsection:: *)

(*Exporting plots*)

Needs["MaTeX`"]

getHistogram[run_,bins_,plotrange_,paddings_,tickDelta_,bottomTicks_,textpos_,framelabel_]:=Histogram[{run[[3]],run[[2]]},bins,"Probability",

ImageSize->150+paddings[[1,1]]+paddings[[1,2]],

ImagePadding->paddings,

AspectRatio->4/6,

PlotRange->plotrange,

Frame->True,

FrameLabel->framelabel,

FrameTicks->{{{#,NumberForm[#,{2,2}]}&/@Range[0,1,tickDelta],Automatic},{bottomTicks,Automatic}},

Epilog->Text[MaTeX["n = "<>ToString[run[[1,1]]]<>",\\; \\tau = "<>ToString[run[[1,2]]]],textpos],

ChartStyle->{ColorData[97][1],ColorData[97][2]}

];

textpos=Scaled[{0.65,0.90}];

ticks1={#,ToString[#]}&/@Range[800,5000,800];

ticks2={#*10^4,ToString[#\[CenterDot]Superscript[10,4],TraditionalForm]}&/@Range[2,10,2];

{h1,h2,h3}={

getHistogram[gdata[[2,1]],{50},{0,0.25},{{40,0},{12,0}},0.10,ticks1,textpos,{None,"Probability"}],

getHistogram[gdata[[6,1]], {3},{0,0.10},{{40,0},{12,0}},0.05,Automatic,textpos,{None,"Probability"}],

getHistogram[gdata[[10,1]],{1},{0,0.18},{{40,0},{32,0}},0.05,Automatic,textpos,{"Triangles","Probability"}]

};

(*plotgrid1=Column[{h1,h2,h3}]*)

{h4,h5,h6}={

getHistogram[gdata[[2,-1]],{800},{0,0.50},{{20,5},{12,0}},0.10,ticks2,textpos,None],

getHistogram[gdata[[6,-1]], {10},{0,0.10},{{20,5},{12,0}},0.05,Automatic,textpos,None],

getHistogram[gdata[[10,-1]], {1},{0,0.10},{{20,5},{32,0}},0.05,Automatic,textpos,{"Triangles"}]

};

(*SwatchLegend[{ColorData[97][1],ColorData[97][2]},{"ECM","Uniform"},LegendLayout\[Rule]"Row"]*)

plotgrid2=Grid[Transpose[{{SwatchLegend[{ColorData[97][1]},{"ECM"}],h1,h2,h3},{SwatchLegend[{ColorData[97][2]},{"Uniform"}],h4,h5,h6}}]]

Export[NotebookDirectory[]<>"plots/ecm_initialtris2.pdf",plotgrid2]

(* Dataset with extra samples *)

getHistogram[gsraw2[[1]], {10},{0,0.10},{{20,5},{12,0}},0.05,Automatic,textpos,None]

(* ::Section:: *)

(*As function of n*)

dataPointsUniform=Map[{#[[1,1]],Mean[#[[2]]]}&,gdata,{2}];

dataPointsECM=Map[{#[[1,1]],Mean[#[[3]]]}&,gdata,{2}];

taulabels=Map["\[Tau] = "<>ToString[#[[1,1,2]]]&,gdata];

(* Standard Deviation with division by N instead of N-1 *)

mySD[xs_]:=Sqrt[Total[(xs-Mean[xs])^2]/Length[xs]];

(* { {x,y}, ErrorBar[err] } *)

getErrorBars[run_]:=Module[{n,tau,avgUni,avgECM,sdUni,sdECM},

{n,tau}=run[[1]];

avgUni=Mean[run[[2]]];

sdUni=mySD[run[[2]]];

avgECM=Mean[run[[3]]];

sdECM=mySD[run[[3]]];

{

{{n,avgUni},ErrorBar[sdUni]},

{{n,avgECM},ErrorBar[sdECM]}

}

]

allErrorBars=Map[getErrorBars,gdata,{2}];

Map[ErrorListPlot[Transpose[#],

ImageSize->500,

Frame->True,

PlotMarkers->Automatic

]&,allErrorBars[[{1,6,11}]]]

(* ::Subsection:: *)

(*Fitting the log-log-plot*)

nRange=All;

(* Weight: 1/err^2 *)

(* Weight: N/Total[(xs-Mean[xs])^2] *)

getWeight[xs_]:=1/(Log[Mean[xs]]-Log[Mean[xs]-Sqrt[Total[(xs-Mean[xs])^2]/Length[xs]]])^2;

(* Several runs for fixed tau but different n *)

getFitData[runs_,index_]:=Map[{Log[#[[1,1]]],Log[Mean[#[[index]]]],getWeight[#[[index]]]}&,runs];

uniformFitData=Map[getFitData[#,2]&,gdata[[All,nRange]]];

ECMFitData=Map[getFitData[#,3]&,gdata[[All,nRange]]];

uniformFits=Map[LinearModelFit[#[[All,{1,2}]],logn,logn(*,Weights\[Rule]#[[All,3]]*)]&,uniformFitData];

ECMFits=    Map[LinearModelFit[#[[All,{1,2}]],logn,logn(*,Weights\[Rule]#[[All,3]]*)]&,ECMFitData];

(*

uniformloglog=Log[dataPointsUniform[[All,nRange]]];

ECMloglog=Log[dataPointsECM[[All,nRange]]];

uniformFits=Map[LinearModelFit[#,logn,logn]&,uniformloglog];

ECMFits=Map[LinearModelFit[#,logn,logn]&,ECMloglog];

*)

uniformFuncs=Map[#[logn]&,uniformFits];

ECMFuncs=Map[#[logn]&,ECMFits];

uniformFits[[1]]["ParameterTable"] (* Get `Standard Error' by "ParameterErrors" *)

uniformFits[[1]]["ParameterConfidenceIntervalTable"] (* Get confidence by "ParameterConfidenceIntervals *)

(* estimate +- standard error *)

uniformFits[[1]]["BestFitParameters"]-uniformFits[[1]]["ParameterErrors"]

uniformFits[[1]]["BestFitParameters"]+uniformFits[[1]]["ParameterErrors"]

tauChoices={1,4,6,8,11};

taulabels=Map["\[Tau] = "<>ToString[NumberForm[#[[1,1,2]],{3,2}]]&,gdata[[tauChoices]]];

repeatColors[n_,k_]:=Table[ColorData[97][Mod[i,k,1]],{i,1,n}]

plot1=Show[ListLogLogPlot[Evaluate[dataPointsUniform[[tauChoices]]~Join~dataPointsECM[[tauChoices]]],

Frame->True,

FrameTicks->{{Table[{10^k,Superscript[10,k]},{k,0,6}],Automatic},{{1000,2000,5000,10000},Automatic}},

FrameLabel->{"n","triangles"},

ImageSize->250,

PlotMarkers->Automatic,

PlotLegends->taulabels,

PlotStyle->repeatColors[10,5]

],

Plot[Evaluate[uniformFuncs[[tauChoices]]],{logn,1,20000}],

Plot[Evaluate[ECMFuncs[[tauChoices]]],{logn,1,20000}]

]

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

(* ::Subsection:: *)

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

gsraw2=Import[NotebookDirectory[]<>"data/graphdata_exponent_hightau.m"];

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

averagesGrouped=GatherBy[gsraw2,{#[[1,2]]&,#[[1,1]]&}];

averagesLoglogdata=Map[{Log[#[[1,1,1]]],Log[Mean[#[[All,2]]]]}&,averagesGrouped[[All,nRange]],{2}];

averagesFitsExtra=Map[LinearModelFit[#,logn,logn]&,averagesLoglogdata];

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

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

Transpose[{avgTauValues-0.000,averagesFitsExtra}]];

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

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

uniformExponents=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, Transpose[{tauValues+0.000,uniformFits}]];

ECMExponents    =Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, Transpose[{tauValues+0.000,ECMFits}]];

Needs["MaTeX`"]

plot2=Show[

ErrorListPlot[{ECMExponents,uniformExponents,averagesExponentsErrorBars},

Joined->True,PlotMarkers->Automatic,

PlotLegends->Placed[{"ECM","canonical","average"},{Left,Bottom}],

PlotRange->{{2,3},{0,1.6}},

ImageSize->300],

Plot[3/2(3-tau),{tau,2,3},PlotStyle->{Black,Dashed},PlotLegends->Placed[LineLegend[{MaTeX["\\frac{3}{2}(3-\\tau)"]},LegendMarkerSize->20],{Left,Bottom}]]]

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