Graph() {}

    Graph(unsigned int n) { reset(n); }

    ~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

        // Based on Erdos-Gallai theorem

        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;

#include "exports.hpp"

#include "graph.hpp"

#include "powerlaw.hpp"

#include <algorithm>

#include <fstream>

#include <iostream>

#include <numeric>

#include <random>

#include <vector>

// Its assumed that u,v are distinct.

// Check if all four vertices are distinct

bool edgeConflicts(const Edge& e1, const Edge& e2) {

    return (e1.u == e2.u || e1.u == e2.v || e1.v == e2.u || e1.v == e2.v);

}

class SwitchChain {

  public:

    SwitchChain()

        : mt(std::random_device{}()), permutationDistribution(0.5)

    // permutationDistribution(0, 2)

    {

        // random_device uses hardware entropy if available

        // std::random_device rd;

        // mt.seed(rd());

    }

    ~SwitchChain() {}

        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;

};

//

// Assumes degree sequence does NOT contain any zeroes!

//

// method2 = true  -> take highest degree and finish its pairing completely

// method2 = false -> take new highest degree after every pairing

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;

    }

    std::vector<decltype(degrees.begin())> available;

    available.reserve(n);

    // Clear the graph

    g.reset(n);

    while (!degrees.empty()) {

    // 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, 2.9f};

    Graph g;

    Graph g1;

    Graph g2;

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

    outfile << '{';

    bool outputComma = false;

    for (int numVertices = 200; numVertices <= 2000; numVertices += 400) {

        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.

            // 500 samples seems to give reasonable results

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

                    }

                }

                //

                SwitchChain chain1, chain2;

                bool do1 = true, do2 = true;

                if (!chain1.initialize(g1)) {

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

                    do1 = false;

                }

                if (!chain2.initialize(g2)) {

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

                    do2 = false;

                }

                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 measurements = 50000;

                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)

                          << "%).";

                //std::cout << std::endl;

                if (outputComma)

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