#pragma once

#include <algorithm>

#include <cassert>

#include <numeric>

#include <vector>

#include <iostream>

class Edge {

  public:

    unsigned int u, v;

    bool operator==(const Edge &e) const { return u == e.u && v == e.v; }

};

class StoredEdge {

  public:

    Edge e;

    // indices into adjacency lists

    // adj[u][u2vindex] = v;

    // adj[v][v2uindex] = u;

    unsigned int u2vindex, v2uindex;

};

class DiDegree {

  public:

    unsigned int in;

    unsigned int out;

};

typedef std::vector<unsigned int> DegreeSequence;

typedef std::vector<DiDegree> DiDegreeSequence;

class Graph {

  public:

    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;

            degrees.pop_back();

            if (degree > degrees.size()) {

                return false;

            }

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

            auto rit = degrees.rbegin();

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

                if (rit->first == 0 || !addEdge({u, rit->second})) {

                    return false;

                }

                rit->first--;

                ++rit;

            }

        }

        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

    bool hasEdge(const Edge& e_) const {

        return badj[e_.u][e_.v];

        //Edge e;

        //if (adj[e_.u].size() <= adj[e_.v].size()) {

        //    e = e_;

        //} else {

        //    e.u = e_.v;

        //    e.v = e_.u;

        //}

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

        StoredEdge se;

        se.e = e;

        se.u2vindex = adj[e.u].size();

        se.v2uindex = adj[e.v].size();

        adj[e.u].push_back(e.v);

        adj[e.v].push_back(e.u);

        edges.push_back(se);

        badj[e.u][e.v] = 1;

        badj[e.v][e.u] = 1;

        return true;

    }

    // There are two possible edge exchanges

    // switchType indicates which one is desired

    // Returns false if the switch is not possible

    bool exchangeEdges(unsigned int e1index, unsigned int e2index, bool switchType) {

        StoredEdge &se1 = edges[e1index];

        StoredEdge &se2 = edges[e2index];

        const Edge &e1 = se1.e;

        const Edge &e2 = se2.e;

        // The new edges configuration is one of these two

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

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

        // Note that to do (B) instead of (A), simply swap e2.u <-> e2.v

        // Now we can just consider switch type (A)

        if (switchType) {

            std::swap(se2.e.u, se2.e.v);

            std::swap(se2.u2vindex, se2.v2uindex);

        }

        // First check if the move is possible

        if (hasEdge({e1.u, e2.u}) || hasEdge({e1.v, e2.v}))

            return false; // conflicting edges

        // Clear old edges

        badj[e1.u][e1.v] = false;

        badj[e1.v][e1.u] = false;

        badj[e2.u][e2.v] = false;

        badj[e2.v][e2.u] = false;

        adj[e1.u][se1.u2vindex] = e2.u;

        adj[e1.v][se1.v2uindex] = e2.v;

        adj[e2.u][se2.u2vindex] = e1.u;

        adj[e2.v][se2.v2uindex] = e1.v;

        // Carefull: when updating se1,se2 also e1 and 2e change

        std::swap(se1.e.v, se2.e.u);

        std::swap(se1.v2uindex, se2.u2vindex);

        // e1 and e2 now contain the NEW edges!!

        badj[e1.u][e1.v] = true;

        badj[e1.v][e1.u] = true;

        badj[e2.u][e2.v] = true;

        badj[e2.v][e2.u] = true;

        return true;

    }

    int countTriangles() const {

        int triangles = 0;

        for (auto& v : adj) {

            for (unsigned int i = 0; i < v.size(); ++i) {

                for (unsigned int j = i + 1; j < v.size(); ++j) {

                    if (hasEdge({v[i], v[j]})) {

                        ++triangles;

                    }

                }

            }

        }

        assert(triangles % 3 == 0);

        return triangles / 3;

    }

    // Should return zero

    int consistencyCheck() const {

        // Check if info in 'edges' is present

        // in adj and badj

        for (auto &se : edges) {

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

                return 1;

            if (!badj[se.e.u][se.e.v])

                return 2;

            if (!badj[se.e.v][se.e.u])

                return 3;

            if (se.u2vindex >= adj[se.e.u].size())

                return 4;

            if (se.v2uindex >= adj[se.e.v].size())

                return 5;

            if (adj[se.e.u][se.u2vindex] != se.e.v)

                return 6;

            if (adj[se.e.v][se.v2uindex] != se.e.u)

                return 7;

        }

        // Check if info in 'adj' is present

        // in badj and edges

        for (unsigned int u = 0; u < adj.size(); ++u) {

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

                if (!badj[u][v])

                    return 8;

                if (!badj[v][u])

                    return 9;

                // Check if it appears in edges

                bool found = false;

                for (auto &se : edges) {

                    if ((se.e.u == u && se.e.v == v) ||

                        (se.e.u == v && se.e.v == u)) {

                        found = true;

                        break;

                    }

                }

                if (!found)

                    return 10;

            }

        }

        // Check if info in 'badj' is present

        // in adj and edges

        // TODO

        return 0;

    }

  private:

    // Graph is saved in three formats for speed

    // They should be kept consistent at all times

    std::vector<std::vector<unsigned int>> adj;

#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() {}

    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 e1index, e2index;

        int timeout = 0;

        // Keep regenerating while conflicting edges

        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()) {

        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;

        if (dmax == 0) {

            std::cerr << "ERROR 1 in GCM.\n";

        }

        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 with a random vertex that is not u itself and to which

            // u has not paired yet!

            available.clear();

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

                if (iter->second != u && !g.hasEdge({u, iter->second}))

                    available.push_back(iter);

            }

            if (available.empty())

                return false;

            std::uniform_int_distribution<> distr(0, available.size() - 1);

            auto vIter = available[distr(rng)];

            // pair u to v

            if (vIter->first == 0)

                std::cerr << "ERROR 2 in GCM1.\n";

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

                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

                // Erasure invalidates all iterators AFTER the erased one

                if (vIter > uIter) {

                    degrees.erase(vIter);

                    degrees.erase(uIter);

                } else {

                    degrees.erase(uIter);

                    degrees.erase(vIter);

                }

            } else {

                // Remove only v if it reaches zero

                uIter->first--;

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

                    }

                }

                //

                // Test the GCM1 and GCM2 success rate

                //

                std::vector<int> greedyTriangles1;

                std::vector<int> greedyTriangles2;

                int successrate1 = 0;

                int successrate2 = 0;

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

                    Graph gtemp;

                    // Take new highest degree every time

                    if (greedyConfigurationModel(ds, gtemp, rng, false)) {

                        ++successrate1;

                        greedyTriangles1.push_back(gtemp.countTriangles());

                        g1 = gtemp;

                    }

                    // Finish all pairings of highest degree first

                    if (greedyConfigurationModel(ds, gtemp, rng, true)) {

                        ++successrate2;

                        greedyTriangles2.push_back(gtemp.countTriangles());

                        g2 = gtemp;

                    }

                }

                SwitchChain chain;

                if (!chain.initialize(g)) {

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

                    return 1;

                }

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

                outputComma = true;

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

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

                outfile << ',' << triangles;

                outfile << ',' << ds;

                outfile << ',' << greedyTriangles1;

                outfile << ',' << greedyTriangles2;

                if (do1) {

                    movesDone = 0;

                    SwitchChain& c = chain1;

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

                        if (c.doMove())

                            ++movesDone;

                    }

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

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

                            if (c.doMove())

                                ++movesDone;

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

                    }

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

                        << " moves succeeded ("

                        << 100.0f * float(movesDone) /

                        float(mixingTime + measurements * measureSkip)

                        << "%).";

                    outfile << ',' << triangles;

                }

                if (do2) {

                    movesDone = 0;

                    SwitchChain& c = chain2;

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

                        if (c.doMove())

                            ++movesDone;

                    }

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

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

                            if (c.doMove())

                                ++movesDone;

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

                    }

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

                        << " moves succeeded ("

                        << 100.0f * float(movesDone) /

                        float(mixingTime + measurements * measureSkip)

                        << "%).";

                    outfile << ',' << triangles;

                }

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

                std::cout << std::endl;

            }

        }

    }

    outfile << '}';

    return 0;

}