#include #include #include #include #include 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 DegreeSequence; typedef std::vector DiDegreeSequence; class Graph { public: Graph() {} Graph(int n) { adj.resize(n); } ~Graph() {} 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; } DegreeSequence getDegreeSequence() const { DegreeSequence d(adj.size()); std::transform(adj.begin(), adj.end(), d.begin(), [](const auto &u) { return u.size(); }); return d; } bool hasEdge(const Edge &e) const { for (int v : adj[e.u]) { if (v == e.v) return true; } return false; } bool addEdge(const Edge &e) { 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 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> adj; std::vector edges; }; class SwitchChain { public: SwitchChain() : mt(std::random_device{}()), 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.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) return false; g = gstart; edgeDistribution.param( std::uniform_int_distribution<>::param_type(0, g.edges.size() - 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() { 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}); SwitchChain chain; if (!chain.initialize(g)) { std::cerr << "Could not initialize Markov chain.\n"; return 1; } std::cout << "Starting switch Markov chain" << std::endl; int movesDone = 0; int movesTotal = 10000; for (int i = 0; i < movesTotal; ++i) if (chain.doMove()) ++movesDone; std::cout << movesDone << '/' << movesTotal << " moves succeeded." << std::endl; return 0; }