(* ::Package:: *)

Needs["ErrorBarPlots`"]


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


(* ::Text:: *)
(*- 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.*)
(*   (c) Compare two greedy ways: (c1) first take highest and finish all its pairings  (c2) take new highest after every single pairing*)
(*   (d) Time evolution for GCM on top of Erdos-Gallai time evolution.*)
(**)
(*- Count #moves that result in +-k triangles (one move could create many triangles at once!)*)
(**)
(*- For a graph snapshot: for all V shapes, compute the number of ways to make it into a triangle:*)
(*  Let u1,u2 be the endpoints of the V. For all neighbors v1 of u1 and v2 of u2, see of v1,v2 has an edge. Meaning, if we were to select randomly an u1 edge and an u2 edge, then whats the probability that it can be used to switch the V into a triangle.*)
(**)
(*- Why does GCM-2 start with very low #triangles*)
(**)
(*- Improve runtime*)
(*   (a) Better direct triangle counting? (I doubt it)*)
(*   (b) Better triangle counting by only keeping track of CHANGES in #triangles*)


(* ::Subsection:: *)
(*Done*)


(* ::Text:: *)
(*- Do a single very long run: nothing weird seems to happen with the triangle counts. Tried 10 million steps.*)
(**)
(*- 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)]) .*)
(*  Computing the f(i,j) = (1-Exp[..]) terms is fine, but then computing Sum[ f(i,j) f(j,k) f(i,k) ) ] over n^3 entries is very slow.*)
(*  *)
(*  - Improve runtime*)
(*   (a) Don't remove/add edges from the std::vector. Simply replace them. Done, is way faster for large n.*)
(*   (b) Do not choose the three permutations with 1/3 probability: choose the "staying" one with zero probability. Should still be a valid switch chain?*)
(*   *)
(*   - Experimental mixing time as function of n. At (n,tau)=(1000,2.5) it seems to be between 10.000 and 20.000 steps.*)
(*     Done. Seems to be something like  (1/2)(32-26(tau-2))n  so we run it for that time without the factor (1/2).*)
(*  *)
(*  *)


(* ::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[]<>"data/graphdata_partial.m"];
gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *)


gdata=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];
(* gdata[[ tau index, n index, run index , {ntau, #tris, ds} ]] *)
nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
taulabels=Map["tau = "<>ToString[#]&,gdata[[All,1,1,1,2]]];


newData=Import[NotebookDirectory[]<>"data/graphdata_3.m"];
mergedData=Import[NotebookDirectory[]<>"data/graphdata_merged.m"];
Export[NotebookDirectory[]<>"data/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 #trianges vs some degree-sequence-property*)


getProperty[ds1_]:=Module[{ds,n=Length[ds1],tmp=ConstantArray[0,{Length[ds1],Length[ds1]}]},
ds=N[ds1/Sqrt[N[Total[ds1]]]]; (* scale degrees by 1/Sqrt[total] *)
(* The next table contains unneeded entries, but its faster to have a square table for the sum *)
tmp=Table[1.-Exp[-ds[[i]]ds[[j]]],{i,1,n},{j,1,n}];
Sum[tmp[[i,j]]*tmp[[j,k]]*tmp[[i,k]],{i,3,n},{j,2,i-1},{k,1,j-1}] (* somehow i>j>k is about 60x faster than doing i<j<k !!! *)
(* This sparser table is slower
tmp=Table[1.-Exp[-ds[[i]]ds[[j]]],{i,1,n-1},{j,i+1,n}];
(* tmp[[a,b]] is now with  ds[[a]]*ds[[a+b]] *)
Sum[tmp[[i,j-i]]*tmp[[j,k-j]]*tmp[[i,k-i]],{i,1,n-2},{j,i+1,n-1},{k,j+1,n}]
*)
];


(* gdata[[ tau index, n index, run index , {ntau, #tris, ds} ]] *)
avgAndProp=ParallelMap[{getProperty[#[[3]]],Mean[#[[2,1;;-1]]]}&,gdata[[2,2,1;;100]]];


Show[ListPlot[avgAndProp,AxesOrigin->{0,0},AxesLabel->{"degree-sequence-property","<#triangles>"},AspectRatio->1],Plot[x,{x,1,1000}]]


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


numPlots=20;
selectedData=gdata[[5,-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:: *)
(*Plot average #triangles vs n*)


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


Show[ListLogLogPlot[averagesErrorBars[[All,All,1]],PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels],Plot[fits,{logn,1,2000}]]


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]],(run[[2,1]]-avg)/stddev,Map[N[(#-avg)/stddev]&,run[[4]]]}
]
stats=Map[getStats,gdata,{3}];


histograms=Map[Histogram[#[[1,3]],PlotRange->{{-8,8},Automatic},PlotLabel->"ErdosGallai deviation: "<>ToString[#[[1,2]]]]&,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]]]&,gdata,{3}];

rateHistograms=Map[Histogram[#,{10},PlotRange->{{0,100},Automatic}]&,successrates,{2}];
TableForm[rateHistograms,TableHeadings->{taulabels,nlabels}]
(*TableForm[Transpose[rateHistograms],TableHeadings->{nlabels,taulabels}]*)