Files
@ 2c2c8fc81176
Branch filter:
Location: AENC/switchchain/cpp/switchchain.cpp
2c2c8fc81176
5.4 KiB
text/x-c++src
Add computation of degree-sequence-property and more
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 | #include "exports.hpp"
#include "graph.hpp"
#include "powerlaw.hpp"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <numeric>
#include <random>
#include <vector>
// 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);
}
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 Graph& gstart) {
if (gstart.edgeCount() == 0)
return false;
g = gstart;
edgeDistribution.param(
std::uniform_int_distribution<>::param_type(0, g.edgeCount() - 1));
return true;
}
bool doMove() {
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.getEdge(edgeDistribution(mt));
e2 = g.getEdge(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() {
// Generate a random degree sequence
std::mt19937 rng(std::random_device{}());
// Goal:
// Degrees follow a power-law distribution with some parameter tau
// Expect: #tri = const * n^{ something }
// The goal is to find the 'something' by finding the number of triangles
// for different values of n and tau
float tauValues[] = {2.2f, 2.35f, 2.5f, 2.65f, 2.8f};
Graph g;
std::ofstream outfile("graphdata.m");
outfile << '{';
bool outputComma = false;
for (int numVertices = 200; numVertices <= 1000; numVertices += 100) {
for (float tau : tauValues) {
DegreeSequence ds(numVertices);
powerlaw_distribution degDist(tau, 1, numVertices);
//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) {
// Generate a graph
// might require multiple tries
for (int i = 1; ; ++i) {
std::generate(ds.begin(), ds.end(),
[°Dist, &rng] { return degDist(rng); });
if (g.createFromDegreeSequence(ds))
break;
// When 10 tries have not worked, output a warning
if (i % 10 == 0) {
std::cerr << "Warning: could not create graph from "
"degree sequence. Trying again...\n";
}
}
SwitchChain chain;
if (!chain.initialize(g)) {
std::cerr << "Could not initialize Markov chain.\n";
return 1;
}
std::cout << "Starting switch Markov chain with n = "
<< numVertices << ", tau = " << tau << ". \t"
<< std::flush;
constexpr int mixingTime = 30000;
constexpr int measureTime = 20000;
constexpr int measureSkip =
200; // Take a sample every ... steps
constexpr int measurements =
(measureTime - 1) / measureSkip + 1;
int movesDone = 0;
int triangles[measurements];
for (int i = 0; i < mixingTime; ++i) {
if (chain.doMove())
++movesDone;
}
for (int i = 0; i < measureTime; ++i) {
if (chain.doMove())
++movesDone;
if (i % measureSkip == 0)
triangles[i / measureSkip] = chain.g.countTriangles();
}
std::cout << movesDone << '/' << mixingTime + measureTime
<< " moves succeeded." << std::endl;
if (outputComma)
outfile << ',';
outputComma = true;
std::sort(ds.begin(), ds.end());
outfile << '{' << '{' << numVertices << ',' << tau << '}';
outfile << ',' << triangles << ',' << ds << '}' << std::flush;
}
}
}
outfile << '}';
return 0;
}
|