(* ::Package:: *)

Needs["ErrorBarPlots`"]


(* ::Section:: *)
(*TODO*)


(* ::Text:: *)
(*- Compute  Sum over i<j<k  of  (1-Exp[- d_i d_j / (2E)]) * (1 - Exp[-d_j d_k / (2E)]) * (1 - Exp[-d_k d_i / (2E)]) .*)
(*   Only depends on degree sequence. Can be done with current data.*)
(*- Do a single very long run to see if there are any weird things after many steps*)
(*- Use different starting point for switch chain that is closer to uniform:*)
(*   Do configuration model, starting with the vertex with highest degree and keeping track of a "forbidden list" meaning dont pair something that is not allowed*)
(*   (a) How close is this to uniform ? At least w.r.t. the measure of #triangles*)
(*   (b) How often does this procedure work/fail. Might still be faster to do switchings from Erdos-Gallai.*)
(*- Improve runtime*)
(*   (a) Don't remove/add edges from the std::vector. Simply replace them*)
(*   (b) Better triangle counting? (I doubt it)*)
(*   (c) Do not choose the three permutations with 1/3 probability: choose the "staying" one with a very low probability. Should still be a valid switch chain?*)
(*- ?*)


(* ::Section:: *)
(*Visualize graphs*)


gsraw=Import[NotebookDirectory[]<>"graphdata.m"];


ListPlot[gsraw[[2]],Joined->True,PlotRange->All,AxesLabel->{"Step","Triangles"}]


gs=Map[Graph[#,GraphLayout->"CircularEmbedding"]&,gsraw[[1]]];
gs2=Map[Graph[#,GraphLayout->Automatic]&,gsraw[[1]]];


Grid[Partition[gs,10],Frame->All]


(* ::Section:: *)
(*Plot triangle counts*)


(* ::Subsection:: *)
(*Data import and data merge*)


gsraw=Import[NotebookDirectory[]<>"graphdata_merged.m"];


newData=Import[NotebookDirectory[]<>"graphdata_tau_multi4.m"];
mergedData=Import[NotebookDirectory[]<>"graphdata_merged.m"];
Export[NotebookDirectory[]<>"graphdata_merged_new.m",Join[mergedData,newData]]


(* ::Subsection:: *)
(*Plot empirical distribution of maximum degree*)


maxDegrees=Map[{#[[1]],Max[#[[3]]]}&,gsraw];
maxDegrees=GatherBy[maxDegrees,{#[[1,2]]&,#[[1,1]]&}];
(* maxDegrees[[ tau index, n index, run index,  ntau or dmax ]] *)


Histogram[maxDegrees[[1,-1,All,2]],PlotRange->{{0,2000},{0,100}},AxesLabel->{"d_max","frequency"}]
Histogram[maxDegrees[[2,-1,All,2]],PlotRange->{{0,2000},{0,100}},AxesLabel->{"d_max","frequency"}]
Histogram[maxDegrees[[3,-1,All,2]],PlotRange->{{0,2000},{0,100}},AxesLabel->{"d_max","frequency"}]


(* ::Subsection:: *)
(*Plot triangle count over "time" in Markov chain instances*)


numPlots=15;
selectedData=gsraw[[-numPlots-1;;-1]];
minCount=Min[Map[Min[#[[2]]]&,selectedData]];
maxCount=Max[Map[Max[#[[2]]]&,selectedData]];
maxTime=Max[Map[Length[#[[2]]]&,selectedData]];
skipPts=Round[maxTime/100]; (* 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,maxTime},PlotLegends->labels]
(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)


(* ::Subsection:: *)
(*Plot average #triangles vs n*)


averages=Map[{#[[1]],Mean[#[[2,1;;-1]]]}&,gsraw];
(* Sort by n *)
averages=SortBy[averages,#[[1,1]]&];
(* Split by n,tau *)
averagesGrouped=GatherBy[averages,{#[[1,2]]&,#[[1,1]]&}];
(* averagesGrouped[[ tau index, n index, run index , {ntau, tri, ds} ]] *)
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]]]/Sqrt[Length[#]]]
}&,averagesGrouped,{2}];


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


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


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


histograms=Map[Histogram[#[[All,2]]]&,averagesGrouped,{2}];
nlabels=Map["n = "<>ToString[#]&,averagesGrouped[[1,All,1,1,1]]];
taulabels=Map["tau = "<>ToString[#]&,averagesGrouped[[All,1,1,1,2]]];


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