if (chain.doMove())

                            ++movesDone;

                    triangles[i] = chain.g.countTriangles();

                }

                std::cout << movesDone << '/' << mixingTime + measurements * measureSkip

                          << " moves succeeded ("

                          << 100.0f * float(movesDone) /

                                 float(mixingTime + measurements * measureSkip)

                          << "%)." << std::endl;

                if (outputComma)

                outputComma = true;

                std::sort(ds.begin(), ds.end());

                outfile << '{' << '{' << numVertices << ',' << tau << '}';

                outfile << ',' << triangles << ',' << ds;

                outfile << ',' << greedyTriangles << '}' << std::flush;

            }

        }

    }

    outfile << '}';

    return 0;

}

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

(**)

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

(**)

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

(* ::Section:: *)

(*Plot triangle counts*)

(* ::Subsection:: *)

(*Data import and data merge*)

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

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

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

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

]

stats=Map[getStats,gdata,{3}];

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