(* ::Package:: *)

Needs["ErrorBarPlots`"]


(* ::Section:: *)
(*Triangle exponent*)


(* ::Text:: *)
(*Expected number of triangles is the following powerlaw*)


(* ::DisplayFormula:: *)
(*#triangles = const\[Cross]n^(T(\[Tau]))      where    T(\[Tau])=3/2 (3-\[Tau])*)


(* When importing from exponent-only-data file or property-data file *)
(* graphdata_exponent_mix32.m *)
(* graphdata_exponent_highN.m *)
(* graphdata_properties2.m *)
(* graphdata_canonical_properties.m *)
(* graphdata_exponent_accurate1.m *)
gsraw=Import[NotebookDirectory[]<>"data/graphdata_exponent_hightau.m"];
gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *)
averagesGrouped=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];


gsraw2=Import[NotebookDirectory[]<>"data/graphdata_canonical_properties2.m"];
gsraw2=SortBy[gsraw2,#[[1,1]]&]; (* Sort by n *)
averagesGrouped2=GatherBy[gsraw2,{#[[1,2]]&,#[[1,1]]&}];
canonicalDatapoints=Map[{#[[1,1,1]],Mean[#[[All,2]]]}&,averagesGrouped2,{2}];


(* 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]]];
mediansErrorBars=Map[
{{#[[1,1,1]],Median[#[[All,2]]]},
ErrorBar[MedianDeviation[#[[All,2]]]]
}&,averagesGrouped,{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]
ErrorListPlot[mediansErrorBars,Joined->True,PlotMarkers->Automatic,PlotRange->All,AxesLabel->{"n","median triangles"},PlotLegends->taulabels]


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


(* ::Subsection:: *)
(*Fitting the log-log-plot*)


nRange=5;;-1;
nRange=All;

averagesLoglogdata=Log[averagesErrorBars[[All,nRange,1]]];
mediansLoglogdata=Log[mediansErrorBars[[All,nRange,1]]];
averagesFits=Map[Fit[#,{1,logn},logn]&,averagesLoglogdata];
averagesFitsExtra=Map[LinearModelFit[#,logn,logn]&,averagesLoglogdata];
mediansFits=Map[Fit[#,{1,logn},logn]&,mediansLoglogdata];
mediansFitsExtra=Map[LinearModelFit[#,logn,logn]&,mediansLoglogdata];


canonicalloglogdata=Log[canonicalDatapoints[[All,nRange]]];
canonicalFitsExtra=Map[LinearModelFit[#,logn,logn]&,canonicalloglogdata];
canonicalNLMs=Map[NonlinearModelFit[#,factor*n^exponent,{factor,exponent},n]&,canonicalDatapoints];


averagesFitsExtra[[1]]["ParameterConfidenceIntervalTable"]
averagesFitsExtra[[1]]["BestFitParameters"]
averagesFitsExtra[[1]]["ParameterErrors"]
averagesFitsExtra[[1]]["ParameterConfidenceIntervals"]


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}]]


canonicalNLMs[[1]]["ParameterConfidenceIntervals"]


(* ::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]];
canonicalTauValues=averagesGrouped2[[All,1,1,1,2]];
averagesExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,
Transpose[{tauValues-0.003,averagesFitsExtra}]];
mediansExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,
Transpose[{tauValues+0.003,mediansFitsExtra}]];
canonicalExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,
Transpose[{canonicalTauValues+0.000,canonicalFitsExtra}]];
canonicalExponentsErrorBars2=Map[{{#[[1]],exponent/.#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-exponent/.#[[2]]["BestFitParameters"]]}&,
Transpose[{canonicalTauValues+0.000,canonicalNLMs}]];
plot2=Show[
ErrorListPlot[{averagesExponentsErrorBars,mediansExponentsErrorBars,canonicalExponentsErrorBars},
Joined->True,PlotMarkers->Automatic,
PlotLegends->Placed[{"average","median","canonical"},{Left,Bottom}],
Frame->True,FrameLabel->{"tau","triangle exponent"},
PlotRange->{{2,3},{0,1.6}},
ImageSize->300],
Plot[3/2(3-tau),{tau,2,3},PlotStyle->{Black,Dashed}]]


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