(* ::Package:: *)

Quit[]


Needs["ErrorBarPlots`"]
Needs["MaTeX`"]


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


(* ::Text:: *)
(*- Triangle law exponent: gather more data*)
(**)
(*- Why does GCM-2 start with very low #triangles*)
(*  Do not only consider number of standard deviations but also relative number of triangles.*)
(*  Look at the following: for all triangles (v1, v2, v3) consider the degrees d1<d2<d3 and make a scatter plot of di vs dj. Make such a scatter plot for the initial GCM-2 graph and for a mixed graph and see how it changes.*)
(**)
(*- GCM success rates: for the degree sequences where it "always fails", look at the degree sequence. Does it have a low/high number of degree 1 nodes? Is the maximum degree very low/high?*)
(**)
(*- Does GCM start closer to uniform?*)
(*   (a) How close 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.*)
(*   (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.*)
(**)
(*- 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).*)


(* ::Subsection:: *)
(*Observations on `Greedy Configuration Model'*)


(* ::Text:: *)
(*- GCM1: Greedy Configuration Model 1: take highest degree and do a single pairing, then take new highest degree*)
(*- GCM2: Greedy Configuration Model 2: take highest degree and completely do all its pairings (at random)*)
(*The difference does not matter if one node is by far the highest.*)
(**)
(*The success rates, conditioned on the degree sequence being graphical, is almost always higher using GCM2. For certain degree sequences the success rate of GCM2 can be 0.9 higher than that of GCM1. (i.e. amost always works vs almost always fails).*)
(*For tau > ~2.3 the success rate of GCM2 seems to be higher than 80% for most sequences.*)
(*For tau < ~2.3 the success rate of GCM2 can drop to less than 10% for some sequences but for many sequences it is still larger than 80%.*)
(**)
(*Success rate of GCM seems to be correlated with mixing time from Erdos-Gallai: higher success rate implies lower mixing time.  *)
(**)
(*Initial #triangles in both GCM1 and GCM2 is always below the average #triangles whereas Erdos-Gallai is usually many times higher than average.*)


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


gsraw=Import[NotebookDirectory[]<>"data/graphdata_timeevol.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]]&,#[[1,1]]&}];
(* Data format: *)
(* gdata[[ tau index, n index, run index , datatype index ]] *)
(* datatype index:
1: {n,tau}
2: #triangles time sequence
3: degree sequence
4: GCM1 starting triangle counts
5: GCM2 starting triangle counts
6: GCM1 time sequence
7: GCM2 time sequence
*)
nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
taulabels=Map["tau = "<>ToString[#]&,gdata[[All,1,1,1,2]]];


(* ::Subsection:: *)
(*Merge data*)


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


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


(* ::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[[4,-1,All,2]],PlotRange->{{0,2000},{0,100}},AxesLabel->{"d_max","frequency"}]
Histogram[maxDegrees[[-1,-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[[1,1]][[-numPlots;;-1]];
measureSkip=1;
minCount=Min[Map[Min[#[[2]]]&,selectedData]];
maxCount=Max[Map[Max[#[[2]]]&,selectedData]];
maxTime=Max[Map[Length[#[[2]]]&,selectedData]];
(* maxTime=30000; *)
skipPts=Max[1,Round[maxTime/500]]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)
coarseData=Map[#[[2,1;;maxTime;;skipPts]]&,selectedData];
labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];
ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},PlotLegends->labels]
(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)


selectedData=gdata[[1,1]];
measureSkip=1;
minCount=Min[Map[Min[#[[2]]]&,selectedData]];
maxCount=Max[Map[Max[#[[2]]]&,selectedData]];
maxTime=Max[Map[Length[#[[2]]]&,selectedData]];
maxTime=30000;
skipPts=Max[1,Round[maxTime/100]]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)
coarseData=Map[#[[2,1;;maxTime;;skipPts]]&,selectedData];
labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];
plot1=ListPlot[coarseData[[1;;5]],Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime}]
plot2=ListPlot[coarseData[[6;;10]],Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime}]
plot3=ListPlot[coarseData[[11;;15]],Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime}]
plot4=ListPlot[coarseData[[16;;20]],Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime}]


(* For export *)
numPlots=25;
selectedPlots={6,7,8,11,12,13,16,17,18,21,22,23};
selectedData=gdata[[2,1]][[selectedPlots]];
measureSkip=1;
minCount=Min[Map[Min[#[[2]]]&,selectedData]];
maxCount=Max[Map[Max[#[[2]]]&,selectedData]];
maxTime=Max[Map[Length[#[[2]]]&,selectedData]];
(* maxTime=30000; *)
skipPts=Max[1,Round[maxTime/5000]]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)
coarseData=Map[#[[2,1;;maxTime;;skipPts]]&,selectedData];
labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];
plotTimeEvol=ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},
Frame->True,FrameLabel->{MaTeX["\\text{timesteps}"],MaTeX["\\text{number of triangles}"]},
PlotLabel->MaTeX["n=1000,\\; \\tau = 2.2"],
ImageSize->250]
(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)


Export[NotebookDirectory[]<>"plots/timeevol.pdf",plotTimeEvol]


movingAvg=Map[MovingAverage[#[[2]],2000][[1;;-1;;skipPts]]-Mean[#[[2,-20000;;-1]]]&,selectedData[[1;;-1;;5]]];
plotMovingAvg=ListPlot[movingAvg,Joined->True,PlotRange->All,DataRange->{0,measureSkip*maxTime},Frame->True,FrameLabel->{"timesteps","number of triangles"}]


(* ::Subsection:: *)
(*Fit exponential to triangles time evolution*)


fitList=Map[NonlinearModelFit[#[[2]],Exp[-(t-t0)/tmix]+c,{{tmix,1000},{t0,10000},{c,2000}},t]&,selectedData];
(* tmix*Log[1/epsilon] gives the time it takes to get a factor epsilon close to the average *)
(* t0 gives the time it takes to be exactly 1 triangle (in absolute value) away from the average *)
(* Use fit["BestFitParameters"] to get parameters *)
(* Use fit[t] to get fit value *)
fitFuncsT=Map[#[t]&,fitList];
tmixList=Map[tmix/.#["BestFitParameters"]&,fitList];


timeplot1=ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},PlotStyle->Opacity[0.5]];
Show[timeplot1,Plot[fitFuncsT,{t,1,maxTime},PlotRange->All]]


(* Log version of exponential fits *)
fitAverages=Map[c/.#["BestFitParameters"]&,fitList];
shiftedFitFuncsT=MapIndexed[#1[t]-fitAverages[[#2[[1]]]]&,fitList];
shiftedCoarseData=MapIndexed[MovingAverage[#1[[2]],1000][[1;;-1;;skipPts]]-fitAverages[[#2[[1]]]]&,selectedData];


(* Plot log version *)
timeplot2=ListLogPlot[shiftedCoarseData[[1;;5]],Joined->True,PlotRange->{0*minCount+0.1,maxCount},DataRange->{0,measureSkip*maxTime},PlotStyle->Opacity[0.5]];
Show[timeplot2,LogPlot[Evaluate[shiftedFitFuncsT[[1;;5]]],{t,1,maxTime},PlotRange->All,PlotStyle->Dotted]]


(* ::Subsection:: *)
(*Plot success rate over "time"*)


numPlots=10;
selectedData=gdata[[1,-1]][[-numPlots;;-1]];
measureSkip=100;
maxTime=Max[Map[Length[#[[4]]]&,selectedData]];
maxTime=10000;
coarseData=Map[#[[4,1;;maxTime/measureSkip]]&,selectedData];
labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];
ListPlot[coarseData,Joined->True,PlotRange->{0,100},DataRange->{0,maxTime},PlotLegends->labels]
(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)


(* ::Subsection:: *)
(*Correlation of avgsuccess rate vs other things*)


compare1=Map[{Mean[#[[4]]],Mean[#[[2]]]}&,gdata,{3}];
(* { GCM1 rate, GCM2 rate, mixing time from ErdosGallai } *)


scatterPlots=Map[ListPlot[#,PlotRange->{{0,100},All},PlotStyle->PointSize[Large]]&,compare1,{2}];
TableForm[scatterPlots,TableHeadings->{taulabels,nlabels}]


(* ::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]]
]
(* Get fit of Exp[-t/tmix] *)
getMixingTime2[values_]:=Module[{avg,etmt,fitVals},
    (* average over the last 20 percent *)
    avg=Mean[values[[-Round[Length[values]/5];;-1]]];
    etmt=FirstPosition[values,_?(#<=avg&)][[1]];
    fitVals=FindFit[values,Exp[-(t-t0)/tmix]+tavg,{{tmix,etmt/4},{t0,2*etmt},{tavg,avg}},t];
    tmix/.fitVals
(* tmix*Log[1/epsilon] gives the time it takes to get a factor epsilon close to the average *)
(* t0 gives the time it takes to be exactly 1 triangle (in absolute value) away from the average *)
]
(* 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[(PrintTemporary[#[[1]]];
{#[[1,2]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]],(1/#[[1,1]])measureSkip * getMixingTime2[#[[2]]]}
)&,gdata,{3}];


mixingTimesBars=Map[
    {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&
,mixingTimes[[All,-1,All,{1,2}]],{1}];
mixingTimesBars2=Map[
    {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&
,mixingTimes[[All,-1,All,{1,3}]],{1}];


Show[
ErrorListPlot[{mixingTimesBars,mixingTimesBars2},Joined->True,PlotMarkers->Automatic,
AxesLabel->{"tau","~~mixing time divided by n, for n = 1000"},
PlotRange->{{2,3},{0,30}}]
,Plot[(32-20(tau-2)),{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:: *)
(*Distribution of initial #triangles for GCM1,GCM2,EG compared to average #triangles.*)


(* Data format: *)
(* gdata[[ tau index, n index, run index , datatype index ]] *)
(* datatype index:
1: {n,tau}
2: #triangles time sequence
3: degree sequence
4: GCM1 starting triangle counts
5: GCM2 starting triangle counts
6: GCM1 time sequence
7: GCM2 time sequence
*)


 (* Stats for a single run at every (n,tau) *)
timeWindow=Round[Length[gdata[[1,1,1,2]]]/10];
skipPts=Max[1,Round[timeWindow/100]];
getSingleStats[runs_]:=Module[{run,avg,stddev},
    run=runs[[1]]; (* Select some run *)
    avg=N[Mean[run[[2,-timeWindow;;-1]]]];
    stddev=N[StandardDeviation[run[[2,timeWindow;;-1]]]];
    {run[[1]], (* {n,tau} *)
    stddev/avg,
    (run[[2,1]])/avg, (* EG starting point *)
    Map[N[#/avg]&,run[[4]]],  (* GCM1 counts *)
    Map[N[#/avg]&,run[[5]]],  (* GCM2 counts *)
    (run[[2,-timeWindow;;-1;;skipPts]])/avg  (* counts in uniform distribution *)
    }
]
singleStats=Map[getSingleStats,gdata,{2}];


(* Yellow: GCM1 (take new highest everytime *)
(* Blue: GCM2 (finish highest first, more similar to EG) *)
(* Green: Actual uniform distribution *)
histogramsSingle=Map[Histogram[{#[[4]],#[[5]],#[[6]]},PlotRange->{{0,3},Automatic},ImageSize->300,PlotLabel->"ErdosGallai="<>ToString[NumberForm[#[[3]],3]]<>"\[Cross]average"]&,singleStats,{2}];


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


 (* Consider all runs *)
timeWindow=Round[Length[gdata[[1,1,1,2]]]/10];
skipPts=Max[1,Round[timeWindow/100]];
getAverage[run_]:=Module[{avg,stddev},
    avg=N[Mean[run[[2,-timeWindow;;-1]]]];
    {
    Mean[run[[4]]]/avg,(* GCM1 counts *)
    Mean[run[[5]]]/avg, (* GCM2 counts *)
    (run[[2,1]])/avg, (* EG starting point *)
    }
]
getTotalStats[runs_]:=Transpose[Map[getAverage,runs]];
totalStats=Map[getTotalStats,gdata,{2}];


(* Yellow: GCM1 (take new highest everytime *)
(* Blue: GCM2 (finish highest first, more similar to EG) *)
histogramsTotal=Map[Histogram[#,{0.1},PlotRange->{{0,5},Automatic},ImageSize->300]&,totalStats,{2}];


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


(* ::Subsection:: *)
(*GCM1 vs GCM2 success rates*)


(* gdata[[ tau index, n index, run index , {ntau, #tris, ds, greedyTriangles} ]] *)
successrates=Map[{Length[#[[4]]],Length[#[[5]]]}&,gdata,{3}];
successrates=Map[Transpose,successrates,{2}];
successratesDelta=Map[Length[#[[5]]]-Length[#[[4]]]&,gdata,{3}];

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

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


(* ::Subsection:: *)
(*Compare success rate with mixing time*)


successrates2=Map[{Length[#[[4]]],Length[#[[5]]],getMixingTime[#[[2]]]}&,gdata,{3}];
(* { GCM1 rate, GCM2 rate, mixing time from ErdosGallai } *)


scatterPlots=Map[ListPlot[#[[All,{1,3}]],PlotRange->{All,All},PlotStyle->PointSize[Large]]&,successrates2,{2}];
TableForm[scatterPlots,TableHeadings->{taulabels,nlabels}]