(* ::Package:: *)

Quit[]


Needs["ErrorBarPlots`"]


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


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}],
Frame->True,FrameLabel->{MaTeX["\\tau"],"triangle powerlaw exponent"},
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]