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=20;

selectedData=gdata[[2,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/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->{"timesteps","number of triangles"},ImageSize->300]

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

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

(* ::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.*)

(* ::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 *)

gsraw=Import[NotebookDirectory[]<>"data/graphdata_exponent_mix32.m"];

gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *)

averagesGrouped=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];

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

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

(* ::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]];

averagesExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,

mediansExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,

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