~Graph() {}

    // Clears any previous edges and create

    // an empty graph on n vertices

    void reset(unsigned int n) {

        edges.clear();

        adj.resize(n);

        for (auto &v : adj)

            v.clear();

        badj.resize(n);

        for (auto &v : badj) {

            v.resize(n);

            v.assign(n, false);

        }

    }

    unsigned int edgeCount() const { return edges.size(); }

    const Edge &getEdge(unsigned int i) const { return edges[i].e; }

    // When the degree sequence is not graphics, the Graph can be

    // in any state, it is not neccesarily empty

    bool createFromDegreeSequence(const DegreeSequence &d) {

        // Havel-Hakimi algorithm

        unsigned int n = d.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 = d[i];

            degrees[i].second = i;

        }

        // Clear the graph

        reset(n);

        while (!degrees.empty()) {

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

            // Highest degree is at back of the vector

            // Take it out

            unsigned int degree = degrees.back().first;

            unsigned int u = degrees.back().second;

            degrees.pop_back();

            if (degree > degrees.size()) {

                return false;

            }

            // Now loop over the last 'degree' entries of degrees

        do {

            e1index = edgeDistribution(mt);

            e2index = edgeDistribution(mt);

            if (++timeout % 100 == 0) {

                std::cerr << "Warning: sampled " << timeout

                          << " random edges but they keep conflicting.\n";

            }

        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));

        // Consider one of the three possible permutations

        // 1) e1.u - e1.v and e2.u - e2.v (original)

        // 2) e1.u - e2.u and e1.v - e2.v

        // 3) e1.u - e2.v and e1.v - e2.u

        bool switchType = permutationDistribution(mt);

        return g.exchangeEdges(e1index, e2index, switchType);

    }

    Graph g;

    std::mt19937 mt;

    std::uniform_int_distribution<> edgeDistribution;

    //std::uniform_int_distribution<> permutationDistribution;

    std::bernoulli_distribution permutationDistribution;

};

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;

    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

            // the degree distribution

            for (int degreeSample = 0; degreeSample < 500; ++degreeSample) {

                // Generate a graph

                // might require multiple tries

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

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

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

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

                    }

                }

                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 movesDone = 0;

                int triangles[measurements];

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

(**)

(*- Count #moves that result in +-k triangles (one move could create many triangles at once!)*)

(**)

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

(*   *)

(*  *)

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

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

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

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