diff --git a/cpp/graph.hpp b/cpp/graph.hpp
index f171ce2ce002b72a8c98f59079433d0496429a7a..f233d4c98ea6c3df4926a36a6c22d2b1a36917e6 100644
--- a/cpp/graph.hpp
+++ b/cpp/graph.hpp
@@ -60,6 +60,7 @@ class Graph {
     // in any state, it is not neccesarily empty
     bool createFromDegreeSequence(const DegreeSequence &d) {
         // Havel-Hakimi algorithm
+        // Based on Erdos-Gallai theorem
 
         unsigned int n = d.size();
 
diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp
index 23d96eda05bceffb77cd0a7ed52c0f6a575ae204..ee360658dcd8328924f7605304c2ddf2eda8a255 100644
--- a/cpp/switchchain.cpp
+++ b/cpp/switchchain.cpp
@@ -63,6 +63,63 @@ class SwitchChain {
     std::bernoulli_distribution permutationDistribution;
 };
 
+bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng) {
+    // Same as 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);
+
+    std::uniform_int_distribution<> distr(1, n - 1);
+
+    while (!degrees.empty()) {
+        // Get the highest degree:
+        unsigned int dmax = 0;
+        unsigned int u = 0;
+        auto maxIter = degrees.begin();
+        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
+            if (iter->first >= dmax) {
+                dmax = iter->first;
+                u = iter->second;
+                maxIter = iter;
+            }
+        }
+        // Take the highest degree out of the vector
+        degrees.erase(maxIter);
+
+        available.clear();
+        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
+            if (iter->first)
+                available.push_back(iter);
+        }
+
+        if (dmax > available.size())
+            return false;
+
+        std::shuffle(available.begin(), available.end(), rng);
+
+        // Now assign randomly to the remaining vertices
+        // Pick 'cmax' distinct integers between 1 and degrees.size()
+        for (unsigned int i = 0; i < dmax; ++i) {
+            if (!g.addEdge({u,available[i]->second})) {
+                std::cerr << "ERROR. Could not add edge in greedy configuration model.\n";
+            }
+            available[i]->first--;
+        }
+    }
+    return true;
+}
+
 int main() {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());
@@ -95,8 +152,28 @@ int main() {
                 for (int i = 1; ; ++i) {
                     std::generate(ds.begin(), ds.end(),
                                   [&degDist, &rng] { return degDist(rng); });
-                    if (g.createFromDegreeSequence(ds))
+                    // First make the sum even
+                    unsigned int sum = std::accumulate(ds.begin(), ds.end(), 0);
+                    if (sum % 2) {
+                        continue;
+                        // Can we do this: ??
+                        ds.back()++;
+                    }
+
+                    // Option 1:
+                    if (g.createFromDegreeSequence(ds)) {
+                        // Test option 2:
+                        //int good = 0;
+                        //for (int i = 0; i < 100; ++i) {
+                        //    if (greedyConfigurationModel(ds, g, rng))
+                        //        ++good;
+                        //}
+                        //std::cerr << "Greedy configuration model success rate: " << good << "%\n";
+
+                        //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 "
diff --git a/showgraphs.m b/showgraphs.m
index 2d71c7db010b24428e83e09aab7709ae7920c6c3..42beccd4041b9dccdceb090f5778542b37bf61fb 100644
--- a/showgraphs.m
+++ b/showgraphs.m
@@ -8,8 +8,6 @@ Needs["ErrorBarPlots`"]
 
 
 (* ::Text:: *)
-(*- Experimental mixing time as function of n. At (n,tau)=(1000,2.5) it seems to be between 10.000 and 20.000 steps.*)
-(**)
 (*- Use different starting point for switch chain that is closer to uniform:*)
 (*   Do configuration model, starting with the vertex with highest degree and keeping track of a "forbidden list" meaning dont pair something that is not allowed*)
 (*   (a) How close is this to uniform ? At least w.r.t. the measure of #triangles*)
@@ -35,6 +33,9 @@ Needs["ErrorBarPlots`"]
 (*  - Improve runtime*)
 (*   (a) Don't remove/add edges from the std::vector. Simply replace them. Done, is way faster for large n.*)
 (*   (b) Do not choose the three permutations with 1/3 probability: choose the "staying" one with zero probability. Should still be a valid switch chain?*)
+(*   *)
+(*   - Experimental mixing time as function of n. At (n,tau)=(1000,2.5) it seems to be between 10.000 and 20.000 steps.*)
+(*     Done. Seems to be something like  (1/2)(32-26(tau-2))n  so we run it for that time without the factor (1/2).*)
 (*  *)
 (*  *)
 
@@ -64,7 +65,7 @@ Grid[Partition[gs,10],Frame->All]
 (*Data import and data merge*)
 
 
-gsraw=Import[NotebookDirectory[]<>"data/graphdata2.m"];
+gsraw=Import[NotebookDirectory[]<>"data/graphdata.m"];
 gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *)