From 1bdb42fca62f7146758a54969040a88df9b7c69e 2017-08-14 15:48:05 From: Tom Bannink Date: 2017-08-14 15:48:05 Subject: [PATCH] Fix Constrained Configuration Model sampling --- diff --git a/cpp/bruteforce_triangles.cpp b/cpp/bruteforce_triangles.cpp index 9b0eb2f4fcf045fffc28d85f07fc1685f8cf14d3..5d840773326130a124fe2b16d51971967662e187 100644 --- a/cpp/bruteforce_triangles.cpp +++ b/cpp/bruteforce_triangles.cpp @@ -39,20 +39,24 @@ int main() { std::cout << "Comparing with Erdos-Gallai triangles." << std::endl; int total = 0; int unequal = 0; + int differenceSum = 0; for (auto &keyvalue : optimalTriangles) { auto ds = keyvalue.first; auto optimalTris = keyvalue.second; g.createFromDegreeSequence(ds); auto EGtris = g.countTriangles(); if (EGtris < optimalTris) { + int difference = optimalTris - EGtris; std::cout << "Erdos-Gallai triangles: " << EGtris << "; Optimal triangles: " << optimalTris - << "; Difference: " << (optimalTris - EGtris) + << "; Difference: " << difference << "; ds = " << ds << std::endl; unequal++; + differenceSum += difference; } total++; } std::cout << "Done." << std::endl; std::cout << unequal << " out of " << total << " degree sequences gave non-optimal triangles.\n"; + std::cout << "Average difference of those: " << float(differenceSum) / float(unequal) << std::endl; } diff --git a/cpp/graph_ccm.hpp b/cpp/graph_ccm.hpp index c3033c2136b8e5af284533d20fa6298a5aa42848..e9844e398e31ed5f4f0b8c75f384c5375239951c 100644 --- a/cpp/graph_ccm.hpp +++ b/cpp/graph_ccm.hpp @@ -10,7 +10,8 @@ // method2 = true -> take highest degree and finish its pairing completely // method2 = false -> take new highest degree after every pairing template -bool constrainedConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool method2) { +bool constrainedConfigurationModel(DegreeSequence &ds, Graph &g, RNG &rng, + bool method2) { // Similar to Havel-Hakimi but instead of pairing up with the highest ones // that remain, simply pair up with random ones unsigned int n = ds.size(); @@ -22,18 +23,16 @@ bool constrainedConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool degrees[i].second = i; } - std::vector available; + // remaining half-edges , iterator into `degrees` + std::vector> available; available.reserve(n); // Clear the graph g.reset(n); while (!degrees.empty()) { - std::shuffle(degrees.begin(), degrees.end(), rng); // Get the highest degree: - // If there are multiple highest ones, we pick a random one, - // ensured by the shuffle. - // The shuffle is needed anyway for the remaining part + // If there are multiple highest ones, we pick the first one unsigned int dmax = 0; auto uIter = degrees.begin(); for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) { @@ -48,45 +47,79 @@ bool constrainedConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool if (dmax == 0) { std::cerr << "ERROR 1 in CCM.\n"; + return false; } unsigned int u = uIter->second; - if (method2) { - // Take the highest degree out of the vector - degrees.erase(uIter); + // Pair with a random vertex that is not u itself and to which + // u has not paired yet! + // Uniform over all available half-edges, not available vertices + unsigned int availableEdges = 0; + available.clear(); + for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) { + if (iter->second != u && !g.hasEdge({u, iter->second})) { + availableEdges += iter->first; + available.push_back(std::make_pair(iter->first, iter)); + } + } + if (availableEdges == 0) + return false; + if (method2) { // Now assign randomly to the remaining vertices - // Since its shuffled, we can pick the first 'dmax' ones - auto vIter = degrees.begin(); + // weighted by the number of half-edges they have left + auto vIter = uIter; while (dmax--) { + std::uniform_int_distribution<> distr(1, availableEdges); + // Go from edge index to number. We should really have a proper + // data structure like some cumulative tree + unsigned int halfEdge = distr(rng); + unsigned int cumulative = 0; + for (auto iter = available.begin(); iter != available.end(); + ++iter) { + cumulative += iter->first; + if (halfEdge <= cumulative) { + vIter = iter->second; + availableEdges -= iter->first; + available.erase(iter); + break; + } + } + if (vIter == uIter) + std::cerr << "ERROR 1 in CCM2.\n"; if (vIter->first == 0) - std::cerr << "ERROR in CCM2.\n"; + std::cerr << "ERROR 2 in CCM2.\n"; if (!g.addEdge({u, vIter->second})) - std::cerr << "ERROR. Could not add edge in CCM2.\n"; + std::cerr << "ERROR 3 in CCM2.\n"; + uIter->first--; vIter->first--; - if (vIter->first == 0) - vIter = degrees.erase(vIter); - else - vIter++; } + for (auto iter = degrees.begin(); iter != degrees.end();) + if (iter->first == 0) + iter = degrees.erase(iter); + else + ++iter; } else { - // 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); + std::uniform_int_distribution<> distr(1, availableEdges); + // Go from edge index to number. We should really have a proper + // data structure like some cumulative tree + auto vIter = uIter; + unsigned int halfEdge = distr(rng); + unsigned int cumulative = 0; + for (auto iter = available.begin(); iter != available.end(); + ++iter) { + cumulative += iter->first; + if (halfEdge <= cumulative) { + vIter = iter->second; + break; + } } - 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 2 in CCM1.\n"; + std::cerr << "ERROR 4 in CCM1.\n"; if (!g.addEdge({u, vIter->second})) - std::cerr << "ERROR. Could not add edge in CCM1.\n"; + std::cerr << "ERROR 5. Could not add edge in CCM1.\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 @@ -107,11 +140,11 @@ bool constrainedConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool if (vIter->first == 0) degrees.erase(vIter); } - //degrees.erase(std::remove_if(degrees.begin(), degrees.end(), - // [](auto x) { return x.first == 0; })); + // degrees.erase(std::remove_if(degrees.begin(), degrees.end(), + // [](auto x) { return x.first == 0; + // })); } } return true; } -