diff --git a/cpp/Makefile b/cpp/Makefile
index 8ffa5bd74bcdec1dcdf1f8286bf75dada688fdfd..120aa577b78eb9f21a2335caa6508c465fd6bb03 100644
--- a/cpp/Makefile
+++ b/cpp/Makefile
@@ -10,20 +10,18 @@ CXXFLAGS += $(INCLUDES)
 CXXFLAGS += -DNDEBUG
 CXXFLAGS += -Wno-int-in-bool-context
 
-all: switchchain switchchain_exponent switchchain_initialtris switchchain_dsp
-
-
-switchchain:
-
-
-switchchain_exponent:
-
-
-switchchain_initialtris:
-
-
-switchchain_dsp:
-
+TARGETS += switchchain
+TARGETS += switchchain_dsp
+TARGETS += switchchain_exponent
+TARGETS += switchchain_initialtris
+TARGETS += switchchain_mixingtime
+TARGETS += switchchain_successrates
+TARGETS += switchchain_timeevol
+
+all: $(TARGETS)
+
+clean:
+	rm -f $(TARGETS)
 
 # target : dep1 dep2 dep3
 # 	$@ = target
diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp
index c28de798386c5f684874bda92e5ab545330f8104..30e4389eed3e925e9bbe5b4d4dd5b790503d056a 100644
--- a/cpp/switchchain.cpp
+++ b/cpp/switchchain.cpp
@@ -2,6 +2,7 @@
 #include "graph.hpp"
 #include "powerlaw.hpp"
 #include "graph_spectrum.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <array>
 #include <fstream>
@@ -10,61 +11,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
 void getTriangleDegrees(const Graph& g) {
     std::vector<std::array<std::size_t,3>> trids;
     const auto& adj = g.getAdj();
diff --git a/cpp/switchchain.hpp b/cpp/switchchain.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d574263da0d03e66ee0161ec70794773b19ed8e8
--- /dev/null
+++ b/cpp/switchchain.hpp
@@ -0,0 +1,59 @@
+#include "graph.hpp"
+#include <iostream>
+#include <random>
+
+// 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.5)
+    // 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() {
+        int e1index, e2index;
+        int timeout = 0;
+        // Keep regenerating while conflicting edges
+        do {
+            e1index = edgeDistribution(mt);
+            e2index = edgeDistribution(mt);
+            if (++timeout % 100 == 0) {
+                std::cerr << "Warning: sampled " << timeout
+                          << " random edges but they keep conflicting.\n";
+            }
+        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
+
+        // 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
+        bool switchType = permutationDistribution(mt);
+        return g.exchangeEdges(e1index, e2index, switchType);
+    }
+
+    Graph g;
+    std::mt19937 mt;
+    std::uniform_int_distribution<> edgeDistribution;
+    //std::uniform_int_distribution<> permutationDistribution;
+    std::bernoulli_distribution permutationDistribution;
+};
+
diff --git a/cpp/switchchain_dsp.cpp b/cpp/switchchain_dsp.cpp
index 268504b74542a89f3b293d016b0335c7511d9f3f..e7b052a0827d70c0d93e8499e01c935d770d18ae 100644
--- a/cpp/switchchain_dsp.cpp
+++ b/cpp/switchchain_dsp.cpp
@@ -1,6 +1,7 @@
 #include "exports.hpp"
 #include "graph.hpp"
 #include "powerlaw.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <fstream>
 #include <iostream>
@@ -8,61 +9,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
 double getProperty(const DegreeSequence& ds) {
     std::vector<std::vector<double>> vals(ds.size());
     for (auto& v : vals) {
@@ -104,8 +50,6 @@ int main() {
     float tauValues[] = {2.1f, 2.5f, 2.9f};
 
     Graph g;
-    Graph g1;
-    Graph g2;
 
     std::ofstream outfile("graphdata_dsp.m");
     outfile << '{';
diff --git a/cpp/switchchain_exponent.cpp b/cpp/switchchain_exponent.cpp
index a862a133c15ed69daaa554f23ea96850a24429a7..85d514b2ab9cbfc5f076c6c835e83996ddb9c59a 100644
--- a/cpp/switchchain_exponent.cpp
+++ b/cpp/switchchain_exponent.cpp
@@ -1,6 +1,7 @@
 #include "exports.hpp"
 #include "graph.hpp"
 #include "powerlaw.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <fstream>
 #include <iostream>
@@ -8,61 +9,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
 int main() {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());
@@ -75,8 +21,6 @@ int main() {
     float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f};
 
     Graph g;
-    Graph g1;
-    Graph g2;
 
     std::ofstream outfile("graphdata_exponent_hightau.m");
     outfile << '{';
diff --git a/cpp/switchchain_initialtris.cpp b/cpp/switchchain_initialtris.cpp
index 6daff04207b0bce9e7e21cd1abef76cfb3a38b37..d84d426b9f59c86cb73b6f2aba7341b905da2811 100644
--- a/cpp/switchchain_initialtris.cpp
+++ b/cpp/switchchain_initialtris.cpp
@@ -1,6 +1,7 @@
 #include "exports.hpp"
 #include "graph.hpp"
 #include "powerlaw.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <fstream>
 #include <iostream>
@@ -8,61 +9,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
 //
 // Assumes degree sequence does NOT contain any zeroes!
 //
diff --git a/cpp/switchchain_mixingtime.cpp b/cpp/switchchain_mixingtime.cpp
index 40ce0fc51672fa8726257558d51f4672966b5f9b..ba0ed098adafa94f400394f3e4e97db3953e2a23 100644
--- a/cpp/switchchain_mixingtime.cpp
+++ b/cpp/switchchain_mixingtime.cpp
@@ -1,6 +1,7 @@
 #include "exports.hpp"
 #include "graph.hpp"
 #include "powerlaw.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <array>
 #include <fstream>
@@ -9,192 +10,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
-void getTriangleDegrees(const Graph& g) {
-    std::vector<std::array<std::size_t,3>> trids;
-    const auto& adj = g.getAdj();
-    int triangles = 0;
-    for (auto& v : adj) {
-        for (unsigned int i = 0; i < v.size(); ++i) {
-            for (unsigned int j = i + 1; j < v.size(); ++j) {
-                if (g.hasEdge({v[i], v[j]})) {
-                    ++triangles;
-                    std::array<std::size_t, 3> ds = {{v.size(), adj[v[i]].size(),
-                                                     adj[v[j]].size()}};
-                    std::sort(ds.begin(), ds.end());
-                    trids.push_back(ds);
-                }
-            }
-        }
-    }
-    assert(triangles % 3 == 0);
-    return;
-}
-
-//
-// Assumes degree sequence does NOT contain any zeroes!
-//
-// method2 = true  -> take highest degree and finish its pairing completely
-// method2 = false -> take new highest degree after every pairing
-template <typename RNG>
-bool greedyConfigurationModel(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();
-
-    // degree, vertex index
-    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);
-    for (unsigned int i = 0; i < n; ++i) {
-        degrees[i].first = ds[i];
-        degrees[i].second = i;
-    }
-
-    std::vector<decltype(degrees.begin())> 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
-        unsigned int dmax = 0;
-        auto uIter = degrees.begin();
-        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
-            if (iter->first >= dmax) {
-                dmax = iter->first;
-                uIter = iter;
-            }
-        }
-
-        if (dmax > degrees.size() - 1)
-            return false;
-
-        if (dmax == 0) {
-            std::cerr << "ERROR 1 in GCM.\n";
-        }
-
-        unsigned int u = uIter->second;
-
-        if (method2) {
-            // Take the highest degree out of the vector
-            degrees.erase(uIter);
-
-            // Now assign randomly to the remaining vertices
-            // Since its shuffled, we can pick the first 'dmax' ones
-            auto vIter = degrees.begin();
-            while (dmax--) {
-                if (vIter->first == 0)
-                    std::cerr << "ERROR in GCM2.\n";
-                if (!g.addEdge({u, vIter->second}))
-                    std::cerr << "ERROR. Could not add edge in GCM2.\n";
-                vIter->first--;
-                if (vIter->first == 0)
-                    vIter = degrees.erase(vIter);
-                else
-                    vIter++;
-            }
-        } 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);
-            }
-            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 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);
-                } else {
-                    degrees.erase(uIter);
-                    degrees.erase(vIter);
-                }
-            } else {
-                // Remove only v if it reaches zero
-                uIter->first--;
-                vIter->first--;
-                if (vIter->first == 0)
-                    degrees.erase(vIter);
-            }
-            //degrees.erase(std::remove_if(degrees.begin(), degrees.end(),
-            //                             [](auto x) { return x.first == 0; }));
-        }
-    }
-    return true;
-}
-
 int main(int argc, char* argv[]) {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());
diff --git a/cpp/switchchain_successrates.cpp b/cpp/switchchain_successrates.cpp
index 4f8df72a039202ab01001d3dd58bd88d50dee767..38df1c97466f53f30cfcfd77d6e32d51b769d325 100644
--- a/cpp/switchchain_successrates.cpp
+++ b/cpp/switchchain_successrates.cpp
@@ -1,6 +1,7 @@
 #include "exports.hpp"
 #include "graph.hpp"
 #include "powerlaw.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <array>
 #include <fstream>
@@ -9,192 +10,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
-void getTriangleDegrees(const Graph& g) {
-    std::vector<std::array<std::size_t,3>> trids;
-    const auto& adj = g.getAdj();
-    int triangles = 0;
-    for (auto& v : adj) {
-        for (unsigned int i = 0; i < v.size(); ++i) {
-            for (unsigned int j = i + 1; j < v.size(); ++j) {
-                if (g.hasEdge({v[i], v[j]})) {
-                    ++triangles;
-                    std::array<std::size_t, 3> ds = {{v.size(), adj[v[i]].size(),
-                                                     adj[v[j]].size()}};
-                    std::sort(ds.begin(), ds.end());
-                    trids.push_back(ds);
-                }
-            }
-        }
-    }
-    assert(triangles % 3 == 0);
-    return;
-}
-
-//
-// Assumes degree sequence does NOT contain any zeroes!
-//
-// method2 = true  -> take highest degree and finish its pairing completely
-// method2 = false -> take new highest degree after every pairing
-template <typename RNG>
-bool greedyConfigurationModel(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();
-
-    // degree, vertex index
-    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);
-    for (unsigned int i = 0; i < n; ++i) {
-        degrees[i].first = ds[i];
-        degrees[i].second = i;
-    }
-
-    std::vector<decltype(degrees.begin())> 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
-        unsigned int dmax = 0;
-        auto uIter = degrees.begin();
-        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
-            if (iter->first >= dmax) {
-                dmax = iter->first;
-                uIter = iter;
-            }
-        }
-
-        if (dmax > degrees.size() - 1)
-            return false;
-
-        if (dmax == 0) {
-            std::cerr << "ERROR 1 in GCM.\n";
-        }
-
-        unsigned int u = uIter->second;
-
-        if (method2) {
-            // Take the highest degree out of the vector
-            degrees.erase(uIter);
-
-            // Now assign randomly to the remaining vertices
-            // Since its shuffled, we can pick the first 'dmax' ones
-            auto vIter = degrees.begin();
-            while (dmax--) {
-                if (vIter->first == 0)
-                    std::cerr << "ERROR in GCM2.\n";
-                if (!g.addEdge({u, vIter->second}))
-                    std::cerr << "ERROR. Could not add edge in GCM2.\n";
-                vIter->first--;
-                if (vIter->first == 0)
-                    vIter = degrees.erase(vIter);
-                else
-                    vIter++;
-            }
-        } 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);
-            }
-            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 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);
-                } else {
-                    degrees.erase(uIter);
-                    degrees.erase(vIter);
-                }
-            } else {
-                // Remove only v if it reaches zero
-                uIter->first--;
-                vIter->first--;
-                if (vIter->first == 0)
-                    degrees.erase(vIter);
-            }
-            //degrees.erase(std::remove_if(degrees.begin(), degrees.end(),
-            //                             [](auto x) { return x.first == 0; }));
-        }
-    }
-    return true;
-}
-
 int main(int argc, char* argv[]) {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());
@@ -208,8 +23,6 @@ int main(int argc, char* argv[]) {
     float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f};
 
     Graph g;
-    Graph g1;
-    Graph g2;
 
     std::ofstream outfile;
 
diff --git a/cpp/switchchain_timeevol.cpp b/cpp/switchchain_timeevol.cpp
index a69b23d9008771c33259ce7be0ba2cef391cd807..374a466dbe551eda889568394fd95173969b6334 100644
--- a/cpp/switchchain_timeevol.cpp
+++ b/cpp/switchchain_timeevol.cpp
@@ -1,6 +1,7 @@
 #include "exports.hpp"
 #include "graph.hpp"
 #include "powerlaw.hpp"
+#include "switchchain.hpp"
 #include <algorithm>
 #include <array>
 #include <fstream>
@@ -9,192 +10,6 @@
 #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.5)
-    // 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() {
-        int e1index, e2index;
-        int timeout = 0;
-        // Keep regenerating while conflicting edges
-        do {
-            e1index = edgeDistribution(mt);
-            e2index = edgeDistribution(mt);
-            if (++timeout % 100 == 0) {
-                std::cerr << "Warning: sampled " << timeout
-                          << " random edges but they keep conflicting.\n";
-            }
-        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));
-
-        // 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
-        bool switchType = permutationDistribution(mt);
-        return g.exchangeEdges(e1index, e2index, switchType);
-    }
-
-    Graph g;
-    std::mt19937 mt;
-    std::uniform_int_distribution<> edgeDistribution;
-    //std::uniform_int_distribution<> permutationDistribution;
-    std::bernoulli_distribution permutationDistribution;
-};
-
-void getTriangleDegrees(const Graph& g) {
-    std::vector<std::array<std::size_t,3>> trids;
-    const auto& adj = g.getAdj();
-    int triangles = 0;
-    for (auto& v : adj) {
-        for (unsigned int i = 0; i < v.size(); ++i) {
-            for (unsigned int j = i + 1; j < v.size(); ++j) {
-                if (g.hasEdge({v[i], v[j]})) {
-                    ++triangles;
-                    std::array<std::size_t, 3> ds = {{v.size(), adj[v[i]].size(),
-                                                     adj[v[j]].size()}};
-                    std::sort(ds.begin(), ds.end());
-                    trids.push_back(ds);
-                }
-            }
-        }
-    }
-    assert(triangles % 3 == 0);
-    return;
-}
-
-//
-// Assumes degree sequence does NOT contain any zeroes!
-//
-// method2 = true  -> take highest degree and finish its pairing completely
-// method2 = false -> take new highest degree after every pairing
-template <typename RNG>
-bool greedyConfigurationModel(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();
-
-    // degree, vertex index
-    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);
-    for (unsigned int i = 0; i < n; ++i) {
-        degrees[i].first = ds[i];
-        degrees[i].second = i;
-    }
-
-    std::vector<decltype(degrees.begin())> 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
-        unsigned int dmax = 0;
-        auto uIter = degrees.begin();
-        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
-            if (iter->first >= dmax) {
-                dmax = iter->first;
-                uIter = iter;
-            }
-        }
-
-        if (dmax > degrees.size() - 1)
-            return false;
-
-        if (dmax == 0) {
-            std::cerr << "ERROR 1 in GCM.\n";
-        }
-
-        unsigned int u = uIter->second;
-
-        if (method2) {
-            // Take the highest degree out of the vector
-            degrees.erase(uIter);
-
-            // Now assign randomly to the remaining vertices
-            // Since its shuffled, we can pick the first 'dmax' ones
-            auto vIter = degrees.begin();
-            while (dmax--) {
-                if (vIter->first == 0)
-                    std::cerr << "ERROR in GCM2.\n";
-                if (!g.addEdge({u, vIter->second}))
-                    std::cerr << "ERROR. Could not add edge in GCM2.\n";
-                vIter->first--;
-                if (vIter->first == 0)
-                    vIter = degrees.erase(vIter);
-                else
-                    vIter++;
-            }
-        } 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);
-            }
-            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 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);
-                } else {
-                    degrees.erase(uIter);
-                    degrees.erase(vIter);
-                }
-            } else {
-                // Remove only v if it reaches zero
-                uIter->first--;
-                vIter->first--;
-                if (vIter->first == 0)
-                    degrees.erase(vIter);
-            }
-            //degrees.erase(std::remove_if(degrees.begin(), degrees.end(),
-            //                             [](auto x) { return x.first == 0; }));
-        }
-    }
-    return true;
-}
-
 int main(int argc, char* argv[]) {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());
@@ -208,8 +23,6 @@ int main(int argc, char* argv[]) {
     float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f};
 
     Graph g;
-    Graph g1;
-    Graph g2;
 
     std::ofstream outfile;