bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool method2) {

    // Similar to Havel-Hakimi but instead of pairing up with the highest ones

    // that remain, simply pair up with random ones

    unsigned int n = ds.size();

    // degree, vertex index

    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);

    for (unsigned int i = 0; i < n; ++i) {

        degrees[i].first = ds[i];

        degrees[i].second = i;

    }

    // Clear the graph

    g.reset(n);

    while (!degrees.empty()) {

        std::shuffle(degrees.begin(), degrees.end(), rng);

        // Get the highest degree:

        // If there are multiple highest ones, we pick a random one,

        // ensured by the shuffle.

        // The shuffle is needed anyway for the remaining part

        unsigned int dmax = 0;

        auto uIter = degrees.begin();

        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {

            if (iter->first >= dmax) {

                dmax = iter->first;

                uIter = iter;

            }

        }

        if (dmax > degrees.size() - 1)

            return false;

        unsigned int u = uIter->second;

        if (method2) {

            // Take the highest degree out of the vector

            degrees.erase(uIter);

            // Now assign randomly to the remaining vertices

            // Since its shuffled, we can pick the first 'dmax' ones

            auto vIter = degrees.begin();

            while (dmax--) {

                if (vIter->first == 0)

                    std::cerr << "ERROR in GCM2.\n";

                if (!g.addEdge({u, vIter->second}))

                    std::cerr << "ERROR. Could not add edge in GCM2.\n";

                vIter->first--;

                if (vIter->first == 0)

                    vIter = degrees.erase(vIter);

                else

                    vIter++;

            }

        } else {

            // pair u to v

            if (vIter->first == 0)

                std::cerr << "ERROR. Could not add edge in GCM1.\n";

            // Purge anything with degree zero

            // Be careful with invalidating the other iterator!

            // Degree of u is always greater or equal to the degree of v

            if (dmax == 1) {

                // Remove both

                if (vIter > uIter) {

                    degrees.erase(vIter);

                    degrees.erase(uIter);

                } else {

                    degrees.erase(uIter);

                    degrees.erase(vIter);

                }

            } else {

                // Remove only v if it reaches zero

                vIter->first--;

                if (vIter->first == 0)

                    degrees.erase(vIter);

            }

            //degrees.erase(std::remove_if(degrees.begin(), degrees.end(),

            //                             [](auto x) { return x.first == 0; }));

        }

    }

    return true;

}

int main() {

    // Generate a random degree sequence

    std::mt19937 rng(std::random_device{}());

    // Goal:

    // Degrees follow a power-law distribution with some parameter tau

    // Expect:  #tri = const * n^{ something }

    // The goal is to find the 'something' by finding the number of triangles

    // for different values of n and tau

    float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f};

    Graph g;

    Graph g2;

    std::ofstream outfile("graphdata.m");

    outfile << '{';

    bool outputComma = false;

    for (int numVertices = 200; numVertices <= 1000; numVertices += 100) {

        for (float tau : tauValues) {

            DegreeSequence ds(numVertices);

            powerlaw_distribution degDist(tau, 1, numVertices);

            //std::poisson_distribution<> degDist(12);

            // For a single n,tau take samples over several instances of

                // Generate a graph

                // might require multiple tries

                for (int i = 1; ; ++i) {

                    std::generate(ds.begin(), ds.end(),

                                  [&degDist, &rng] { return degDist(rng); });

                    // First make the sum even

                    unsigned int sum = std::accumulate(ds.begin(), ds.end(), 0);

                    if (sum % 2) {

                        continue;

                        // Can we do this: ??

                        ds.back()++;

                    }

                    if (g.createFromDegreeSequence(ds))

                        break;

                    // When 10 tries have not worked, output a warning

                    if (i % 10 == 0) {

                        std::cerr << "Warning: could not create graph from "

                                     "degree sequence. Trying again...\n";

                    }

                }

                //

                //

                for (int i = 0; i < 100; ++i) {

                    }

                }

                SwitchChain chain;

                if (!chain.initialize(g)) {

                    std::cerr << "Could not initialize Markov chain.\n";

                    return 1;

                }

                std::cout << "Running n = " << numVertices << ", tau = " << tau

                          << ". \t" << std::flush;

                //int mixingTime = (32.0f - 26.0f*(tau - 2.0f)) * numVertices; //40000;

                //constexpr int measurements = 50;

                //constexpr int measureSkip =

                //    200; // Take a sample every ... steps

                int mixingTime = 0;

                constexpr int measureSkip = 1;

                int movesDone = 0;

                int triangles[measurements];

                for (int i = 0; i < mixingTime; ++i) {

                    if (chain.doMove())

                        ++movesDone;

                }

                for (int i = 0; i < measurements; ++i) {

                    for (int j = 0; j < measureSkip; ++j)

                        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)

                if (outputComma)

                    outfile << ',' << '\n';

                outputComma = true;

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

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

            }

        }

    }

    outfile << '}';

    return 0;

}

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

(**)

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

(**)

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

(*  *)

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

(*   *)

(*  *)

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

gsraw=Import[NotebookDirectory[]<>"data/graphdata.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"}]

histograms=Map[Histogram[{#[[1,4]]},PlotRange->{{0,2},Automatic},PlotLabel->"ErdosGallai="<>ToString[NumberForm[#[[1,3]],3]]<>"\[Cross]average. stddev="<>ToString[NumberForm[#[[1,2]],3]]<>"\[Cross]average"]&,stats,{2}];

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

(* ::Subsection:: *)

(*Greedy CM success rates*)

(* gdata[[ tau index, n index, run index , {ntau, #tris, ds, greedyTriangles} ]] *)

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

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

(*TableForm[Transpose[rateHistograms],TableHeadings->{nlabels,taulabels}]*)