diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp index 6d02203eed4a5c9ae25660a37b4c29a3cbd7f315..b14a8ecfd6c5fdbcea2373c39d7f623ca350ea99 100644 --- a/cpp/switchchain.cpp +++ b/cpp/switchchain.cpp @@ -78,6 +78,9 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth degrees[i].second = i; } + std::vector available; + available.reserve(n); + // Clear the graph g.reset(n); @@ -99,6 +102,10 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth if (dmax > degrees.size() - 1) return false; + if (dmax == 0) { + std::cerr << "ERROR 1 in GCM.\n"; + } + unsigned int u = uIter->second; if (method2) { @@ -120,21 +127,28 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth vIter++; } } else { - // Pair with a random vertex that is not u itself - std::uniform_int_distribution<> distr(0, degrees.size() - 2); - auto vIter = degrees.begin() + distr(rng); - if (vIter == uIter) - vIter++; + // 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 in GCM1.\n"; - if (!g.addEdge({uIter->second, vIter->second})) + 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); @@ -144,6 +158,7 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth } } else { // Remove only v if it reaches zero + uIter->first--; vIter->first--; if (vIter->first == 0) degrees.erase(vIter); @@ -167,6 +182,7 @@ int main() { float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f}; Graph g; + Graph g1; Graph g2; std::ofstream outfile("graphdata.m"); @@ -181,8 +197,9 @@ int main() { //std::poisson_distribution<> degDist(12); // For a single n,tau take samples over several instances of - // the degree distribution - for (int degreeSample = 0; degreeSample < 500; ++degreeSample) { + // the degree distribution. + // 500 samples seems to give reasonable results + for (int degreeSample = 0; degreeSample < 200; ++degreeSample) { // Generate a graph // might require multiple tries for (int i = 1; ; ++i) { @@ -207,14 +224,23 @@ int main() { } // - // Test the greedy configuration model success rate + // Test the GCM1 and GCM2 success rate // - std::vector greedyTriangles; - int successrate = 0; + std::vector greedyTriangles1; + std::vector greedyTriangles2; + int successrate1 = 0; + int successrate2 = 0; for (int i = 0; i < 100; ++i) { - if (greedyConfigurationModel(ds, g2, rng, true)) { - ++successrate; - greedyTriangles.push_back(g2.countTriangles()); + Graph gtemp; + if (greedyConfigurationModel(ds, gtemp, rng, false)) { + ++successrate1; + greedyTriangles1.push_back(gtemp.countTriangles()); + g1 = gtemp; + } + if (greedyConfigurationModel(ds, gtemp, rng, true)) { + ++successrate2; + greedyTriangles2.push_back(gtemp.countTriangles()); + g2 = gtemp; } } @@ -224,6 +250,17 @@ int main() { 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; @@ -232,7 +269,7 @@ int main() { //constexpr int measureSkip = // 200; // Take a sample every ... steps int mixingTime = 0; - constexpr int measurements = 10000; + constexpr int measurements = 5000; constexpr int measureSkip = 1; @@ -255,7 +292,8 @@ int main() { << " moves succeeded (" << 100.0f * float(movesDone) / float(mixingTime + measurements * measureSkip) - << "%)." << std::endl; + << "%)."; + //std::cout << std::endl; if (outputComma) outfile << ',' << '\n'; @@ -263,8 +301,59 @@ int main() { std::sort(ds.begin(), ds.end()); outfile << '{' << '{' << numVertices << ',' << tau << '}'; - outfile << ',' << triangles << ',' << ds; - outfile << ',' << greedyTriangles << '}' << std::flush; + 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; } } }