}

        return true;

    }

    DegreeSequence getDegreeSequence() const {

        DegreeSequence d(adj.size());

        std::transform(adj.begin(), adj.end(), d.begin(),

                       [](const auto &u) { return u.size(); });

        return d;

    }

    // Assumes valid vertex indices

        for (unsigned int v : adj[e.u]) {

            if (v == e.v)

                return true;

        }

        return false;

    }

    bool addEdge(const Edge &e) {

        if (e.u >= adj.size() || e.v >= adj.size())

            return false;

        if (hasEdge(e))

            return false;

#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:

        // random_device uses hardware entropy if available

        // std::random_device rd;

        // mt.seed(rd());

    }

    ~SwitchChain() {}

    bool initialize(const Graph& gstart) {

        if (gstart.edgeCount() == 0)

            return false;

        g = gstart;

        edgeDistribution.param(

            std::uniform_int_distribution<>::param_type(0, g.edgeCount() - 1));

        return true;

    }

    bool doMove() {

        int timeout = 0;

            ++timeout;

            if (timeout % 100 == 0) {

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

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

            }

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

        // Note that it might be that these new edges already exist

    }

    Graph g;

    std::mt19937 mt;

    std::uniform_int_distribution<> edgeDistribution;

};

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

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

                    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;

                }

                constexpr int measureTime = 20000;

                constexpr int measureSkip =

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

                constexpr int measurements =

                    (measureTime - 1) / measureSkip + 1;

                int movesDone = 0;

                int triangles[measurements];

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

                    if (chain.doMove())

                        ++movesDone;

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

(**)

(*- Improve runtime*)

(*   (a) Don't remove/add edges from the std::vector. Simply replace them*)

(*   (b) Better direct triangle counting? (I doubt it)*)

(*   (b') Better triangle counting by only keeping track of CHANGES in #triangles*)

(**)

(*- Experimental mixing time as function of n. At (n,tau)=(1000,2.5) it seems to be between 10.000 and 20.000 steps.*)

(**)

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

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

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

averagesErrorBars=Map[

{{#[[1,1,1]],Mean[#[[All,2]]]},

ErrorBar[StandardDeviation[#[[All,2]]]/Sqrt[Length[#]]]

}&,averagesGrouped,{2}];

ErrorListPlot[averagesErrorBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels]

ListLogLogPlot[averagesErrorBars[[All,All,1]],Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels]

(* ::Subsection:: *)

(*Fitting the log-log-plot*)

loglogdata=Log[averagesErrorBars[[All,All,1]]];

fits=Map[Fit[#,{1,logn},logn]&,loglogdata];

Show[ListLogLogPlot[averagesErrorBars[[All,All,1]],PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels],Plot[fits,{logn,1,2000}]]

tauValues=averagesGrouped[[All,1,1,1,2]];