From a79267af1717d3a82f640c21e571936023913be9 2017-01-27 16:09:48 From: Tom Bannink Date: 2017-01-27 16:09:48 Subject: [PATCH] Add Havel-Hakimi algorithm and mathematica output --- diff --git a/cpp/showgraphs.m b/cpp/showgraphs.m new file mode 100644 index 0000000000000000000000000000000000000000..6aa966bf83cd84e5e5e239bea18f967b6b50f0ab --- /dev/null +++ b/cpp/showgraphs.m @@ -0,0 +1,9 @@ +(* ::Package:: *) + +gsraw=Import[NotebookDirectory[]<>"graphdata.m"]; + + +gs=Graph/@gsraw; + + +Grid[Partition[gs,10],Frame->All] diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp index 2b2a5efafb28ed999295cd523ccbd8ab3b981c0e..541119e116e1dfe1fe761fe234ebb333f4397f4d 100644 --- a/cpp/switchchain.cpp +++ b/cpp/switchchain.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -6,7 +7,7 @@ class Edge { public: - int u, v; + unsigned int u, v; bool operator==(const Edge &e) const { return u == e.u && v == e.v; } }; @@ -24,30 +25,72 @@ std::ostream &operator<<(std::ostream &s, const Edge &e) { class DiDegree { public: - int in; - int out; + unsigned int in; + unsigned int out; }; -typedef std::vector DegreeSequence; +typedef std::vector DegreeSequence; typedef std::vector DiDegreeSequence; class Graph { public: Graph() {} - Graph(int n) { adj.resize(n); } + 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) { - // - // TODO - // - // See - // http://stackoverflow.com/questions/13102738/creating-graphs-using-a-degree-sequence - // See https://en.wikipedia.org/wiki/Havel%E2%80%93Hakimi_algorithm - (void)d; - return false; + // Havel-Hakimi algorithm + + unsigned int n = d.size(); + + // degree, vertex index + std::vector> 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 { @@ -57,8 +100,9 @@ class Graph { return d; } + // Assumes valid vertex indices bool hasEdge(const Edge &e) const { - for (int v : adj[e.u]) { + for (unsigned int v : adj[e.u]) { if (v == e.v) return true; } @@ -66,6 +110,8 @@ class Graph { } 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); @@ -91,20 +137,20 @@ class Graph { } // Find the edges in the adjacency lists - int i1, j1, i2, j2; - for (i1 = 0; i1 < (int)adj[e1.u].size(); ++i1) { + 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 < (int)adj[e1.v].size(); ++j1) { + for (j1 = 0; j1 < adj[e1.v].size(); ++j1) { if (adj[e1.v][j1] == e1.u) break; } - for (i2 = 0; i2 < (int)adj[e2.u].size(); ++i2) { + for (i2 = 0; i2 < adj[e2.u].size(); ++i2) { if (adj[e2.u][i2] == e2.v) break; } - for (j2 = 0; j2 < (int)adj[e2.v].size(); ++j2) { + for (j2 = 0; j2 < adj[e2.v].size(); ++j2) { if (adj[e2.v][j2] == e2.u) break; } @@ -146,12 +192,28 @@ class Graph { return true; } + private: // Graph is saved in two formats for speed // The two should be kept consistent at all times - std::vector> adj; + std::vector> adj; std::vector 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) { @@ -161,31 +223,23 @@ class SwitchChain { } ~SwitchChain() {} - bool initialize(const DegreeSequence &d) { - if (!g.createFromDegreeSequence(d)) - return false; - edgeDistribution.param( - std::uniform_int_distribution<>::param_type(0, g.edges.size() - 1)); - return true; - } - bool initialize(const Graph &gstart) { - if (gstart.edges.size() == 0) + if (gstart.edgeCount() == 0) return false; g = gstart; edgeDistribution.param( - std::uniform_int_distribution<>::param_type(0, g.edges.size() - 1)); + 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)]; + Edge e1 = g.getEdge(edgeDistribution(mt)); + Edge e2 = g.getEdge(edgeDistribution(mt)); // Keep regenerating while conflicting edges int timeout = 0; while (edgeConflicts(e1, e2)) { - e1 = g.edges[edgeDistribution(mt)]; - e2 = g.edges[edgeDistribution(mt)]; + e1 = g.getEdge(edgeDistribution(mt)); + e2 = g.getEdge(edgeDistribution(mt)); ++timeout; if (timeout % 100 == 0) { std::cerr << "Warning: sampled " << timeout @@ -214,13 +268,11 @@ class SwitchChain { }; int main() { - Graph g(6); - g.addEdge({0, 1}); - g.addEdge({0, 2}); - g.addEdge({0, 3}); - g.addEdge({2, 3}); - g.addEdge({3, 4}); - g.addEdge({3, 5}); + Graph g; + if (!g.createFromDegreeSequence({1, 2, 2, 2, 3, 3, 3})) { + std::cerr << "Could not create graph from degree sequence.\n"; + return 1; + } SwitchChain chain; if (!chain.initialize(g)) { @@ -228,14 +280,22 @@ int main() { return 1; } + std::ofstream outfile("graphdata.m"); + outfile << '{' << g; + std::cout << "Starting switch Markov chain" << std::endl; int movesDone = 0; - int movesTotal = 10000; - for (int i = 0; i < movesTotal; ++i) + int movesTotal = 100000; + for (int i = 0; i < movesTotal; ++i) { if (chain.doMove()) ++movesDone; + if (i % (movesTotal / 100) == (movesTotal / 100 - 1)) + outfile << ',' << chain.g; + } + outfile << '}'; - std::cout << movesDone << '/' << movesTotal << " moves succeeded." << std::endl; + std::cout << movesDone << '/' << movesTotal << " moves succeeded." + << std::endl; return 0; }