From 446bcd991614eefca806092578c40cd31909062a 2017-01-26 17:16:41
From: Tom Bannink <tom.bannink@cwi.nl>
Date: 2017-01-26 17:16:41
Subject: [PATCH] Add initial cpp source with skeleton

---

diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..25b3fe70a86d760de72c089966eb1d7c90724457
--- /dev/null
+++ b/cpp/switchchain.cpp
@@ -0,0 +1,206 @@
+#include <iostream>
+#include <numeric>
+#include <random>
+#include <vector>
+
+class Edge {
+  public:
+    int u, v;
+
+    bool operator==(const Edge &e) const { return u == e.u && v == e.v; }
+};
+
+// 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);
+}
+
+std::ostream &operator<<(std::ostream &s, const Edge &e) {
+    s << '{' << e.u << ',' << e.v << '}';
+    return s;
+}
+
+class DiDegree {
+  public:
+    int in;
+    int out;
+};
+
+typedef std::vector<int> DegreeSequence;
+typedef std::vector<DiDegree> DiDegreeSequence;
+
+class Graph {
+  public:
+    Graph() : edgeCount(0) {}
+
+    Graph(int n) : edgeCount(0) { adj.resize(n); }
+
+    ~Graph() {}
+
+    bool createFromSequence(const DegreeSequence &d) {
+        edgeCount = std::accumulate(d.begin(), d.end(), 0);
+
+        //
+        // TODO
+        //
+
+        return false;
+    }
+
+    bool hasEdge(const Edge &e) const {
+        for (int v : adj[e.u]) {
+            if (v == e.v)
+                return true;
+        }
+        return false;
+    }
+
+    // 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
+        int i1, j1, i2, j2;
+        for (i1 = 0; i1 < (int)adj[e1.u].size(); ++i1) {
+            if (adj[e1.u][i1] == e1.v)
+                break;
+        }
+        for (j1 = 0; j1 < (int)adj[e1.v].size(); ++j1) {
+            if (adj[e1.v][j1] == e1.u)
+                break;
+        }
+        for (i2 = 0; i2 < (int)adj[e2.u].size(); ++i2) {
+            if (adj[e2.u][i2] == e2.v)
+                break;
+        }
+        for (j2 = 0; j2 < (int)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;
+    }
+
+    // Graph is saved in two formats for speed
+    // The two should be kept consistent at all times
+    std::vector<std::vector<int>> adj;
+    std::vector<Edge> edges;
+    int edgeCount;
+};
+
+class SwitchChain {
+  public:
+    SwitchChain() : permutationDistribution(0, 2) {
+        // random_device uses hardware entropy if available
+        std::random_device rd;
+        mt.seed(rd());
+    }
+    ~SwitchChain() {}
+
+    bool initialize(const DegreeSequence &d) {
+        if (!g.createFromSequence(d))
+            return false;
+        edgeDistribution.param(
+            std::uniform_int_distribution<>::param_type(0, g.edgeCount - 1));
+        return true;
+    }
+
+    bool doMove() {
+        Edge e1 = g.edges[edgeDistribution(mt)];
+        Edge e2 = g.edges[edgeDistribution(mt)];
+        // Keep regenerating while conflicting edges
+        int timeout = 0;
+        while (edgeConflicts(e1, e2)) {
+            e1 = g.edges[edgeDistribution(mt)];
+            e2 = g.edges[edgeDistribution(mt)];
+            ++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
+        // in which case we also reject the move
+        // This is checked in exchangeEdges
+
+        int perm = permutationDistribution(mt);
+        if (perm == 0) // Original permutation
+            return false;
+        return g.exchangeEdges(e1, e2, perm == 1);
+    }
+
+    Graph g;
+    std::mt19937 mt;
+    std::uniform_int_distribution<> edgeDistribution;
+    std::uniform_int_distribution<> permutationDistribution;
+};
+
+int main() {
+    SwitchChain chain;
+    if (!chain.initialize({3, 2, 4, 2, 2, 1})) {
+        std::cerr << "Could not initialize Markov chain.\n";
+        return 1;
+    }
+
+    std::cout << "Starting switch Markov chain" << std::endl;
+    int movesDone = 0;
+    for (int i = 0; i < 100; ++i)
+        if (chain.doMove())
+            ++movesDone;
+
+    std::cout << movesDone << "/100 moves succeeded." << std::endl;
+
+    return 0;
+}