From 0f3a4ccb62ea0c525cc49f885c8812e60964261e 2017-02-02 18:09:08
From: Tom Bannink <tombannink@gmail.com>
Date: 2017-02-02 18:09:08
Subject: [PATCH] Add generation of random degree distribution and split cpp file

---

diff --git a/cpp/exports.hpp b/cpp/exports.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8be142765b6e1d9e4d7aafa700f6d8b67d67655a
--- /dev/null
+++ b/cpp/exports.hpp
@@ -0,0 +1,23 @@
+#pragma once
+#include "graph.hpp"
+
+std::ostream &operator<<(std::ostream &s, const Edge &e) {
+    s << '{' << e.u << ',' << e.v << '}';
+    return s;
+}
+
+// Mathematica style export
+std::ostream &operator<<(std::ostream &s, const Graph &g) {
+    if (g.edgeCount() == 0) {
+        s << '{' << '}';
+        return s;
+    }
+    s << '{' << g.getEdge(0).u << '<' << '-' << '>' << g.getEdge(0).v;
+    for (unsigned int i = 1; i < g.edgeCount(); ++i) {
+        const Edge &e = g.getEdge(i);
+        s << ',' << e.u << '<' << '-' << '>' << e.v;
+    }
+    s << '}';
+    return s;
+}
+
diff --git a/cpp/graph.hpp b/cpp/graph.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6c5ea22dffe491607b800b87a7d540be013d5ea
--- /dev/null
+++ b/cpp/graph.hpp
@@ -0,0 +1,187 @@
+#pragma once
+#include <numeric>
+#include <vector>
+
+class Edge {
+  public:
+    unsigned int u, v;
+
+    bool operator==(const Edge &e) const { return u == e.u && v == e.v; }
+};
+
+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) { adj.resize(n); }
+
+    ~Graph() {}
+
+    void resize(unsigned int n) {
+        if (n < adj.size()) {
+            edges.clear();
+        }
+        adj.resize(n);
+    }
+
+    unsigned int edgeCount() const { return edges.size(); }
+
+    Edge &getEdge(unsigned int i) { return edges[i]; }
+    const Edge &getEdge(unsigned int i) const { return edges[i]; }
+
+    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;
+        }
+
+        edges.clear();
+        adj.resize(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()) {
+                edges.clear();
+                adj.clear();
+                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})) {
+                    edges.clear();
+                    adj.clear();
+                    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 {
+        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;
+        edges.push_back(e);
+        adj[e.u].push_back(e.v);
+        adj[e.v].push_back(e.u);
+        return true;
+    }
+
+    // There are two possible edge exchanges
+    // switchType indicates which one is desired
+    // Returns false if the switch is not possible
+    bool exchangeEdges(const Edge &e1, const Edge &e2, bool switchType) {
+        // 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 e1.v - e2.u
+        // First check if the move is possible
+        if (switchType) {
+            if (hasEdge({e1.u, e2.u}) || hasEdge({e1.v, e2.v}))
+                return false; // conflicting edges
+        } else {
+            if (hasEdge({e1.u, e2.v}) || hasEdge({e1.v, e2.u}))
+                return false; // conflicting edges
+        }
+
+        // Find the edges in the adjacency lists
+        unsigned int i1, j1, i2, j2;
+        for (i1 = 0; i1 < adj[e1.u].size(); ++i1) {
+            if (adj[e1.u][i1] == e1.v)
+                break;
+        }
+        for (j1 = 0; j1 < adj[e1.v].size(); ++j1) {
+            if (adj[e1.v][j1] == e1.u)
+                break;
+        }
+        for (i2 = 0; i2 < adj[e2.u].size(); ++i2) {
+            if (adj[e2.u][i2] == e2.v)
+                break;
+        }
+        for (j2 = 0; j2 < adj[e2.v].size(); ++j2) {
+            if (adj[e2.v][j2] == e2.u)
+                break;
+        }
+
+        // Remove the old edges
+        bool removedOne = false;
+        for (auto iter = edges.begin(); iter != edges.end();) {
+            if (*iter == e1) {
+                iter = edges.erase(iter);
+                if (removedOne)
+                    break;
+                removedOne = true;
+            } else if (*iter == e2) {
+                iter = edges.erase(iter);
+                if (removedOne)
+                    break;
+                removedOne = true;
+            } else {
+                ++iter;
+            }
+        }
+
+        // Add the new edges
+        if (switchType) {
+            adj[e1.u][i1] = e2.u;
+            adj[e1.v][j1] = e2.v;
+            adj[e2.u][i2] = e1.u;
+            adj[e2.v][j2] = e1.v;
+            edges.push_back({e1.u, e2.u});
+            edges.push_back({e1.v, e2.v});
+        } else {
+            adj[e1.u][i1] = e2.v;
+            adj[e1.v][j1] = e2.u;
+            adj[e2.u][i2] = e1.v;
+            adj[e2.v][j2] = e1.u;
+            edges.push_back({e1.u, e2.v});
+            edges.push_back({e1.v, e2.u});
+        }
+        return true;
+    }
+
+  private:
+    // Graph is saved in two formats for speed
+    // The two should be kept consistent at all times
+    std::vector<std::vector<unsigned int>> adj;
+    std::vector<Edge> edges;
+};
+
diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp
index 541119e116e1dfe1fe761fe234ebb333f4397f4d..436dc7a20821815e95ec973e42b11561b121fa49 100644
--- a/cpp/switchchain.cpp
+++ b/cpp/switchchain.cpp
@@ -4,13 +4,8 @@
 #include <numeric>
 #include <random>
 #include <vector>
-
-class Edge {
-  public:
-    unsigned int u, v;
-
-    bool operator==(const Edge &e) const { return u == e.u && v == e.v; }
-};
+#include "graph.hpp"
+#include "exports.hpp"
 
 // Its assumed that u,v are distinct.
 // Check if all four vertices are distinct
@@ -18,202 +13,6 @@ 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);
 }
 
-std::ostream &operator<<(std::ostream &s, const Edge &e) {
-    s << '{' << e.u << ',' << e.v << '}';
-    return s;
-}
-
-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) { adj.resize(n); }
-
-    ~Graph() {}
-
-    void resize(unsigned int n) {
-        if (n < adj.size()) {
-            edges.clear();
-        }
-        adj.resize(n);
-    }
-
-    unsigned int edgeCount() const { return edges.size(); }
-
-    Edge &getEdge(unsigned int i) { return edges[i]; }
-    const Edge &getEdge(unsigned int i) const { return edges[i]; }
-
-    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;
-        }
-
-        edges.clear();
-        adj.resize(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()) {
-                edges.clear();
-                adj.clear();
-                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})) {
-                    edges.clear();
-                    adj.clear();
-                    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 {
-        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;
-        edges.push_back(e);
-        adj[e.u].push_back(e.v);
-        adj[e.v].push_back(e.u);
-        return true;
-    }
-
-    // There are two possible edge exchanges
-    // switchType indicates which one is desired
-    // Returns false if the switch is not possible
-    bool exchangeEdges(const Edge &e1, const Edge &e2, bool switchType) {
-        // 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 e1.v - e2.u
-        // First check if the move is possible
-        if (switchType) {
-            if (hasEdge({e1.u, e2.u}) || hasEdge({e1.v, e2.v}))
-                return false; // conflicting edges
-        } else {
-            if (hasEdge({e1.u, e2.v}) || hasEdge({e1.v, e2.u}))
-                return false; // conflicting edges
-        }
-
-        // Find the edges in the adjacency lists
-        unsigned int i1, j1, i2, j2;
-        for (i1 = 0; i1 < adj[e1.u].size(); ++i1) {
-            if (adj[e1.u][i1] == e1.v)
-                break;
-        }
-        for (j1 = 0; j1 < adj[e1.v].size(); ++j1) {
-            if (adj[e1.v][j1] == e1.u)
-                break;
-        }
-        for (i2 = 0; i2 < adj[e2.u].size(); ++i2) {
-            if (adj[e2.u][i2] == e2.v)
-                break;
-        }
-        for (j2 = 0; j2 < adj[e2.v].size(); ++j2) {
-            if (adj[e2.v][j2] == e2.u)
-                break;
-        }
-
-        // Remove the old edges
-        bool removedOne = false;
-        for (auto iter = edges.begin(); iter != edges.end();) {
-            if (*iter == e1) {
-                iter = edges.erase(iter);
-                if (removedOne)
-                    break;
-                removedOne = true;
-            } else if (*iter == e2) {
-                iter = edges.erase(iter);
-                if (removedOne)
-                    break;
-                removedOne = true;
-            } else {
-                ++iter;
-            }
-        }
-
-        // Add the new edges
-        if (switchType) {
-            adj[e1.u][i1] = e2.u;
-            adj[e1.v][j1] = e2.v;
-            adj[e2.u][i2] = e1.u;
-            adj[e2.v][j2] = e1.v;
-            edges.push_back({e1.u, e2.u});
-            edges.push_back({e1.v, e2.v});
-        } else {
-            adj[e1.u][i1] = e2.v;
-            adj[e1.v][j1] = e2.u;
-            adj[e2.u][i2] = e1.v;
-            adj[e2.v][j2] = e1.u;
-            edges.push_back({e1.u, e2.v});
-            edges.push_back({e1.v, e2.u});
-        }
-        return true;
-    }
-
-  private:
-    // Graph is saved in two formats for speed
-    // The two should be kept consistent at all times
-    std::vector<std::vector<unsigned int>> adj;
-    std::vector<Edge> edges;
-};
-
-// Mathematica style export
-std::ostream &operator<<(std::ostream &s, const Graph &g) {
-    if (g.edgeCount() == 0) {
-        s << '{' << '}';
-        return s;
-    }
-    s << '{' << g.getEdge(0).u << '<' << '-' << '>' << g.getEdge(0).v;
-    for (unsigned int i = 1; i < g.edgeCount(); ++i) {
-        const Edge &e = g.getEdge(i);
-        s << ',' << e.u << '<' << '-' << '>' << e.v;
-    }
-    s << '}';
-    return s;
-}
-
 class SwitchChain {
   public:
     SwitchChain() : mt(std::random_device{}()), permutationDistribution(0, 2) {
@@ -269,10 +68,34 @@ class SwitchChain {
 
 int main() {
     Graph g;
-    if (!g.createFromDegreeSequence({1, 2, 2, 2, 3, 3, 3})) {
+
+    // Generate a random degree sequence
+    std::mt19937 gen(std::random_device{}());
+
+    // 50 nodes with average degree 12
+    DegreeSequence ds(50);
+    std::poisson_distribution<> degDist(7);
+
+    // Try at most 10 times to generate a valid sequence
+    bool validGraph = false;
+    for (int i = 0; i < 10; ++i) {
+        std::generate(ds.begin(), ds.end(),
+                      [&degDist, &gen] { return degDist(gen); });
+        if (g.createFromDegreeSequence(ds)) {
+            validGraph = true;
+            break;
+        }
+    }
+    if (!validGraph) {
         std::cerr << "Could not create graph from degree sequence.\n";
         return 1;
     }
+    std::sort(ds.begin(), ds.end());
+
+    std::cout << "Degree sequence:";
+    for (auto i : ds)
+        std::cout << ' ' << i;
+    std::cout << std::endl;
 
     SwitchChain chain;
     if (!chain.initialize(g)) {
@@ -289,7 +112,7 @@ int main() {
     for (int i = 0; i < movesTotal; ++i) {
         if (chain.doMove())
             ++movesDone;
-        if (i % (movesTotal / 100) == (movesTotal / 100 - 1))
+        if (i % (movesTotal / 50) == (movesTotal / 50 - 1))
             outfile << ',' << chain.g;
     }
     outfile << '}';