diff --git a/appveyor.yml b/appveyor.yml
index 67009131..10e0db39 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -49,6 +49,10 @@ before_build:
   - git submodule update --init --recursive
 
 install:
+  - bash -lc "curl -O http://repo.msys2.org/msys/x86_64/msys2-keyring-r21.b39fb11-1-any.pkg.tar.xz"
+  - bash -lc "curl -O http://repo.msys2.org/msys/x86_64/msys2-keyring-r21.b39fb11-1-any.pkg.tar.xz.sig"
+  - bash -lc "pacman-key --verify msys2-keyring-r21.b39fb11-1-any.pkg.tar.xz.sig"
+  - bash -lc "pacman -U --noconfirm msys2-keyring-r21.b39fb11-1-any.pkg.tar.xz"
   - bash -lc "pacman --needed --noconfirm -Sy pacman-mirrors"
   - bash -lc "pacman --noconfirm -Sy"
   - bash -lc "pacman --noconfirm -S zstd"
diff --git a/include/GraphHelper.h b/include/GraphHelper.h
index fd44b623..7b67c356 100644
--- a/include/GraphHelper.h
+++ b/include/GraphHelper.h
@@ -100,8 +100,6 @@ class Graph
     vector<size_t> const& get_neighbours(size_t v, igraph_neimode_t mode);
     size_t get_random_neighbour(size_t v, igraph_neimode_t mode, igraph_rng_t* rng);
 
-    pair<size_t, size_t> get_endpoints(size_t e);
-
     inline size_t get_random_node(igraph_rng_t* rng)
     {
       return get_random_int(0, this->vcount() - 1, rng);
@@ -113,7 +111,7 @@ class Graph
     inline size_t ecount() { return igraph_ecount(this->_graph); };
     inline double total_weight() { return this->_total_weight; };
     inline size_t total_size() { return this->_total_size; };
-    inline int is_directed() { return igraph_is_directed(this->_graph); };
+    inline int is_directed() { return this->_is_directed; };
     inline double density() { return this->_density; };
     inline int correct_self_loops() { return this->_correct_self_loops; };
     inline int is_weighted() { return this->_is_weighted; };
@@ -128,12 +126,15 @@ class Graph
       return this->_edge_weights[e];
     };
 
+    inline void edge(size_t eid, size_t &from, size_t &to) {
+      from = IGRAPH_FROM(this->get_igraph(), eid);
+      to = IGRAPH_TO(this->get_igraph(), eid);
+    }
+
     inline vector<size_t> edge(size_t e)
     {
-      igraph_integer_t v1, v2;
-      igraph_edge(this->_graph, e, &v1, &v2);
       vector<size_t> edge(2);
-      edge[0] = v1; edge[1] = v2;
+      this->edge(e, edge[0], edge[1]);
       return edge;
     }
 
@@ -147,7 +148,7 @@ class Graph
 
     inline size_t degree(size_t v, igraph_neimode_t mode)
     {
-      if (mode == IGRAPH_IN)
+      if (mode == IGRAPH_IN || !this->is_directed())
         return this->_degree_in[v];
       else if (mode == IGRAPH_OUT)
         return this->_degree_out[v];
@@ -159,7 +160,7 @@ class Graph
 
     inline double strength(size_t v, igraph_neimode_t mode)
     {
-      if (mode == IGRAPH_IN)
+      if (mode == IGRAPH_IN || !this->is_directed())
         return this->_strength_in[v];
       else if (mode == IGRAPH_OUT)
         return this->_strength_out[v];
@@ -173,6 +174,7 @@ class Graph
 
   private:
     igraph_t* _graph;
+    igraph_vector_t _temp_igraph_vector;
 
     // Utility variables to easily access the strength of each node
     vector<double> _strength_in;
@@ -199,6 +201,7 @@ class Graph
     double _total_weight;
     size_t _total_size;
     int _is_weighted;
+    bool _is_directed;
 
     int _correct_self_loops;
     double _density;
diff --git a/include/MutableVertexPartition.h b/include/MutableVertexPartition.h
index 9aec67fd..caa57f2f 100644
--- a/include/MutableVertexPartition.h
+++ b/include/MutableVertexPartition.h
@@ -77,10 +77,11 @@ class MutableVertexPartition
     inline Graph* get_graph() { return this->graph; };
 
     void renumber_communities();
-    vector<size_t> renumber_communities(map<size_t, size_t> const& original_fixed_membership);
+    void renumber_communities(map<size_t, size_t> const& original_fixed_membership);
     void renumber_communities(vector<size_t> const& new_membership);
     void set_membership(vector<size_t> const& new_membership);
-    vector<size_t> static renumber_communities(vector<MutableVertexPartition*> partitions);
+    void relabel_communities(vector<size_t> const& new_comm_id);
+    vector<size_t> static rank_order_communities(vector<MutableVertexPartition*> partitions);
     size_t get_empty_community();
     size_t add_empty_community();
     void from_coarse_partition(vector<size_t> const& coarse_partition_membership);
@@ -97,8 +98,36 @@ class MutableVertexPartition
     inline double total_weight_in_all_comms()         { return this->_total_weight_in_all_comms; };
     inline size_t total_possible_edges_in_all_comms() { return this->_total_possible_edges_in_all_comms; };
 
-    double weight_to_comm(size_t v, size_t comm);
-    double weight_from_comm(size_t v, size_t comm);
+    inline double weight_to_comm(size_t v, size_t comm)
+    {
+      if (this->_current_node_cache_community_to != v)
+      {
+        this->cache_neigh_communities(v, IGRAPH_OUT);
+        this->_current_node_cache_community_to = v;
+      }
+
+      if (comm < this->_cached_weight_to_community.size())
+        return this->_cached_weight_to_community[comm];
+      else
+        return 0.0;
+    }
+
+    inline double weight_from_comm(size_t v, size_t comm)
+    {
+      if (!this->graph->is_directed())
+        return weight_to_comm(v, comm);
+
+      if (this->_current_node_cache_community_from != v)
+      {
+        this->cache_neigh_communities(v, IGRAPH_IN);
+        this->_current_node_cache_community_from = v;
+      }
+
+      if (comm < this->_cached_weight_from_community.size())
+        return this->_cached_weight_from_community[comm];
+      else
+        return 0.0;
+    }
 
     vector<size_t> const& get_neigh_comms(size_t v, igraph_neimode_t);
     set<size_t> get_neigh_comms(size_t v, igraph_neimode_t mode, vector<size_t> const& constrained_membership);
diff --git a/src/GraphHelper.cpp b/src/GraphHelper.cpp
index 4c5b72c4..35996c48 100644
--- a/src/GraphHelper.cpp
+++ b/src/GraphHelper.cpp
@@ -97,6 +97,7 @@ Graph::Graph(igraph_t* graph,
   this->_node_self_weights = node_self_weights;
 
   this->_correct_self_loops = correct_self_loops;
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
 }
 
@@ -120,6 +121,7 @@ Graph::Graph(igraph_t* graph,
   this->_correct_self_loops = this->has_self_loops();
 
   this->_node_self_weights = node_self_weights;
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
 }
 
@@ -140,6 +142,7 @@ Graph::Graph(igraph_t* graph,
   this->_node_sizes = node_sizes;
 
   this->_correct_self_loops = correct_self_loops;
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -161,6 +164,7 @@ Graph::Graph(igraph_t* graph,
 
   this->_correct_self_loops = this->has_self_loops();
 
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -175,6 +179,7 @@ Graph::Graph(igraph_t* graph, vector<double> const& edge_weights, int correct_se
   this->_edge_weights = edge_weights;
   this->_is_weighted = true;
   this->set_default_node_size();
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -191,6 +196,7 @@ Graph::Graph(igraph_t* graph, vector<double> const& edge_weights)
   this->_correct_self_loops = this->has_self_loops();
 
   this->set_default_node_size();
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -207,6 +213,7 @@ Graph::Graph(igraph_t* graph, vector<size_t> const& node_sizes, int correct_self
 
   this->set_default_edge_weight();
   this->_is_weighted = false;
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -225,6 +232,7 @@ Graph::Graph(igraph_t* graph, vector<size_t> const& node_sizes)
 
   this->_correct_self_loops = this->has_self_loops();
 
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -236,6 +244,7 @@ Graph::Graph(igraph_t* graph, int correct_self_loops)
   this->_correct_self_loops = correct_self_loops;
   this->set_defaults();
   this->_is_weighted = false;
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -249,6 +258,7 @@ Graph::Graph(igraph_t* graph)
 
   this->_correct_self_loops = this->has_self_loops();
 
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -260,6 +270,7 @@ Graph::Graph()
   this->set_defaults();
   this->_is_weighted = false;
   this->_correct_self_loops = false;
+  igraph_vector_init(&this->_temp_igraph_vector, this->vcount());
   this->init_admin();
   this->set_self_weights();
 }
@@ -271,6 +282,7 @@ Graph::~Graph()
     igraph_destroy(this->_graph);
     delete this->_graph;
   }
+  igraph_vector_destroy(&this->_temp_igraph_vector);
 }
 
 int Graph::has_self_loops()
@@ -364,81 +376,64 @@ void Graph::init_admin()
 {
 
   size_t m = this->ecount();
+  size_t n = this->vcount();
+  this->_is_directed = igraph_is_directed(this->_graph);
+
+  this->_strength_in.clear();
+  this->_strength_in.resize(n, 0.0);
+
+  this->_degree_in.clear();
+  this->_degree_in.resize(n, 0.0);
+
+  if (this->_is_directed) {
+    this->_strength_out.clear();
+    this->_strength_out.resize(n, 0.0);
+
+    this->_degree_out.clear();
+    this->_degree_out.resize(n, 0);
+
+    this->_degree_all.clear();
+    this->_degree_all.resize(n, 0);
+  }
 
   // Determine total weight in the graph.
   this->_total_weight = 0.0;
-  for (size_t e = 0; e < m; e++)
-    this->_total_weight += this->edge_weight(e);
+  for (size_t e = 0; e < m; e++) {
+    double w = this->edge_weight(e);
+    this->_total_weight += w;
+
+    size_t from, to;
+    this->edge(e, from, to);
+
+    if (this->is_directed()) {
+      this->_strength_in[to] += w;
+      this->_strength_out[from] += w;
+
+      this->_degree_in[to]++;
+      this->_degree_out[from]++;
+      this->_degree_all[to]++;
+      this->_degree_all[from]++;
+    } else {
+      // we only compute strength_in and degree_in for undirected graphs
+      this->_strength_in[to] += w;
+      this->_strength_in[from] += w;
+
+      // recall that igraph ignores the mode for undirected graphs
+      this->_degree_in[to]++;
+      this->_degree_in[from]++;
+    }
+  }
 
   // Make sure to multiply by 2 for undirected graphs
   //if (!this->is_directed())
   //  this->_total_weight *= 2.0;
 
-  size_t n = this->vcount();
-
   this->_total_size = 0;
   for (size_t v = 0; v < n; v++)
     this->_total_size += this->node_size(v);
 
-  igraph_vector_t weights;
-  igraph_vector_t res;
-
-  // Strength IN
-  igraph_vector_init(&res, n);
-  // Copy weights to an igraph_vector_t
-  igraph_vector_init_copy(&weights, &this->_edge_weights[0], this->ecount());
-  // Calculcate strength
-  igraph_strength(this->_graph, &res, igraph_vss_all(), IGRAPH_IN, true, &weights);
-  igraph_vector_destroy(&weights);
-
-  // Assign to strength vector
-  this->_strength_in.clear();
-  this->_strength_in.resize(n);
-  for (size_t v = 0; v < n; v++)
-    this->_strength_in[v] = VECTOR(res)[v];
-  igraph_vector_destroy(&res);
-
-  // Strength OUT
-  igraph_vector_init(&res, n);
-  // Copy weights to an igraph_vector_t
-  igraph_vector_init_copy(&weights, &this->_edge_weights[0], this->ecount());
-  // Calculcate strength
-  igraph_strength(this->_graph, &res, igraph_vss_all(), IGRAPH_OUT, true, &weights);
-  igraph_vector_destroy(&weights);
-
-  // Assign to strength vector
-  this->_strength_out.clear();
-  this->_strength_out.resize(n);
-  for (size_t v = 0; v < n; v++)
-    this->_strength_out[v] = VECTOR(res)[v];
-  igraph_vector_destroy(&res);
-
-  // Degree IN
-  igraph_vector_init(&res, n);
-  igraph_degree(this->_graph, &res, igraph_vss_all(), IGRAPH_IN, true);
-  this->_degree_in.clear();
-  this->_degree_in.resize(n);
-  for (size_t v = 0; v < n; v++)
-    this->_degree_in[v] = VECTOR(res)[v];
-  igraph_vector_destroy(&res);
-
-  // Degree OUT
-  igraph_vector_init(&res, n);
-  igraph_degree(this->_graph, &res, igraph_vss_all(), IGRAPH_OUT, true);
-  this->_degree_out.clear();
-  this->_degree_out.resize(n);
-  for (size_t v = 0; v < n; v++)
-    this->_degree_out[v] = VECTOR(res)[v];
-  igraph_vector_destroy(&res);
-
-  // Degree ALL
-  igraph_vector_init(&res, n);
-  igraph_degree(this->_graph, &res, igraph_vss_all(), IGRAPH_ALL, true);
-  this->_degree_all.clear();
-  this->_degree_all.resize(n);
-  for (size_t v = 0; v < n; v++)
-    this->_degree_all[v] = VECTOR(res)[v];
-  igraph_vector_destroy(&res);
+  // this is initialized in the constructors
+  igraph_vector_t *res = &this->_temp_igraph_vector;
 
   // Calculate density;
   double w = this->total_weight();
@@ -477,9 +472,8 @@ void Graph::cache_neighbour_edges(size_t v, igraph_neimode_t mode)
     cerr << "Degree: " << degree << endl;
   #endif
 
-  igraph_vector_t incident_edges;
-  igraph_vector_init(&incident_edges, degree);
-  igraph_incident(this->_graph, &incident_edges, v, mode);
+  igraph_vector_t *incident_edges = &this->_temp_igraph_vector;
+  igraph_incident(this->_graph, incident_edges, v, mode);
 
   vector<size_t>* _cached_neigh_edges = NULL;
   switch (mode)
@@ -497,14 +491,12 @@ void Graph::cache_neighbour_edges(size_t v, igraph_neimode_t mode)
       _cached_neigh_edges = &(this->_cached_neigh_edges_all);
       break;
   }
-  _cached_neigh_edges->assign(igraph_vector_e_ptr(&incident_edges, 0),
-                              igraph_vector_e_ptr(&incident_edges, degree));
+  _cached_neigh_edges->assign(igraph_vector_e_ptr(incident_edges, 0),
+                              igraph_vector_e_ptr(incident_edges, degree));
   #ifdef DEBUG
     cerr << "Number of edges: " << _cached_neigh_edges->size() << endl;
   #endif
 
-
-  igraph_vector_destroy(&incident_edges);
   #ifdef DEBUG
     cerr << "exit void Graph::cache_neighbour_edges(" << v << ", " << mode << ");" << endl;
   #endif
@@ -512,6 +504,9 @@ void Graph::cache_neighbour_edges(size_t v, igraph_neimode_t mode)
 
 vector<size_t> const& Graph::get_neighbour_edges(size_t v, igraph_neimode_t mode)
 {
+  if (!this->is_directed())
+    mode = IGRAPH_ALL; // igraph ignores mode for undirected graphs
+
   switch (mode)
   {
     case IGRAPH_IN:
@@ -539,13 +534,6 @@ vector<size_t> const& Graph::get_neighbour_edges(size_t v, igraph_neimode_t mode
   throw Exception("Incorrect model for getting neighbour edges.");
 }
 
-pair<size_t, size_t> Graph::get_endpoints(size_t e)
-{
-  igraph_integer_t from, to;
-  igraph_edge(this->_graph, e,&from, &to);
-  return make_pair<size_t, size_t>((size_t)from, (size_t)to);
-}
-
 void Graph::cache_neighbours(size_t v, igraph_neimode_t mode)
 {
   #ifdef DEBUG
@@ -556,9 +544,8 @@ void Graph::cache_neighbours(size_t v, igraph_neimode_t mode)
     cerr << "Degree: " << degree << endl;
   #endif
 
-  igraph_vector_t neighbours;
-  igraph_vector_init(&neighbours, degree);
-  igraph_neighbors(this->_graph, &neighbours, v, mode);
+  igraph_vector_t *neighbours = &this->_temp_igraph_vector;
+  igraph_neighbors(this->_graph, neighbours, v, mode);
 
   vector<size_t>* _cached_neighs = NULL;
   switch (mode)
@@ -576,8 +563,7 @@ void Graph::cache_neighbours(size_t v, igraph_neimode_t mode)
       _cached_neighs = &(this->_cached_neighs_all);
       break;
   }
-  _cached_neighs->assign(igraph_vector_e_ptr(&neighbours, 0),igraph_vector_e_ptr(&neighbours, degree));
-  igraph_vector_destroy(&neighbours);
+  _cached_neighs->assign(igraph_vector_e_ptr(neighbours, 0),igraph_vector_e_ptr(neighbours, degree));
 
   #ifdef DEBUG
     cerr << "Number of edges: " << _cached_neighs->size() << endl;
@@ -590,6 +576,9 @@ void Graph::cache_neighbours(size_t v, igraph_neimode_t mode)
 
 vector< size_t > const& Graph::get_neighbours(size_t v, igraph_neimode_t mode)
 {
+  if (!this->is_directed())
+    mode = IGRAPH_ALL; // igraph ignores mode for undirected graphs
+
   switch (mode)
   {
     case IGRAPH_IN:
@@ -637,7 +626,7 @@ size_t Graph::get_random_neighbour(size_t v, igraph_neimode_t mode, igraph_rng_t
   if (this->degree(v, mode) <= 0)
     throw Exception("Cannot select a random neighbour for an isolated node.");
 
-  if (igraph_is_directed(this->_graph) && mode != IGRAPH_ALL)
+  if (this->is_directed() && mode != IGRAPH_ALL)
   {
     if (mode == IGRAPH_OUT)
     {
@@ -714,60 +703,61 @@ Graph* Graph::collapse_graph(MutableVertexPartition* partition)
   #ifdef DEBUG
     cerr << "Graph* Graph::collapse_graph(vector<size_t> membership)" << endl;
   #endif
-  size_t m = this->ecount();
 
   #ifdef DEBUG
     cerr << "Current graph has " << this->vcount() << " nodes and " << this->ecount() << " edges." << endl;
     cerr << "Collapsing to graph with " << partition->n_communities() << " nodes." << endl;
   #endif
 
-  vector< map<size_t, double> > collapsed_edge_weights(partition->n_communities());
-
-  igraph_integer_t v, u;
-  for (size_t e = 0; e < m; e++)
-  {
-    double w = this->edge_weight(e);
-    igraph_edge(this->_graph, e, &v, &u);
-    size_t v_comm = partition->membership((size_t)v);
-    size_t u_comm = partition->membership((size_t)u);
-    if (collapsed_edge_weights[v_comm].count(u_comm) > 0)
-      collapsed_edge_weights[v_comm][u_comm] += w;
-    else
-      collapsed_edge_weights[v_comm][u_comm] = w;
-  }
-
-  // Now create vector for edges, first determined the number of edges
-  size_t m_collapsed = 0;
   size_t n_collapsed = partition->n_communities();
+  vector<vector<size_t> > community_memberships = partition->get_communities();
 
-  for (vector< map<size_t, double> >::iterator itr = collapsed_edge_weights.begin();
-       itr != collapsed_edge_weights.end(); itr++)
-  {
-      m_collapsed += itr->size();
-  }
+  vector<double> collapsed_weights;
+  double total_collapsed_weight = 0.0;
 
+  vector<double> edge_weight_to_community(n_collapsed, 0.0);
+  vector<bool> neighbour_comm_added(n_collapsed, false);
+
+  // collapsed edges for new graph
   igraph_vector_t edges;
-  vector<double> collapsed_weights(m_collapsed, 0.0);
-  double total_collapsed_weight = 0.0;
+  igraph_vector_init(&edges, 0);
+
+  for (size_t v_comm = 0; v_comm < n_collapsed; v_comm++) {
+    vector<size_t> neighbour_communities;
+    for (size_t v : community_memberships[v_comm]) {
+        for (size_t e : this->get_neighbour_edges(v, IGRAPH_OUT)) {
+            size_t from, to;
+            this->edge(e, from, to);
+
+            if ((size_t) from != v) {
+                // need to skip because IGRAPH_OUT is ignored for undirected graphs
+                continue;
+            }
+
+            size_t u_comm = partition->membership(to);
+
+            double w = this->edge_weight(e);
+            // Self loops appear twice here if the graph is undirected, so divide by 2.0 in that case.
+            if (from == to && !this->is_directed())
+                w /= 2.0;
+
+            if (!neighbour_comm_added[u_comm]) {
+                neighbour_comm_added[u_comm] = true;
+                neighbour_communities.push_back(u_comm);
+            }
+            edge_weight_to_community[u_comm] += w;
+        }
+    }
 
-  igraph_vector_init(&edges, 2*m_collapsed); // Vector or edges with edges (edge[0], edge[1]), (edge[2], edge[3]), etc...
+    for (size_t u_comm : neighbour_communities) {
+        igraph_vector_push_back(&edges, v_comm);
+        igraph_vector_push_back(&edges, u_comm);
+        collapsed_weights.push_back(edge_weight_to_community[u_comm]);
+        total_collapsed_weight += edge_weight_to_community[u_comm];
 
-  size_t e_idx = 0;
-  for (size_t v = 0; v < n_collapsed; v++)
-  {
-    for (map<size_t, double>::iterator itr = collapsed_edge_weights[v].begin();
-         itr != collapsed_edge_weights[v].end(); itr++)
-    {
-      size_t u = itr->first;
-      double w = itr->second;
-      VECTOR(edges)[2*e_idx] = v;
-      VECTOR(edges)[2*e_idx+1] = u;
-      collapsed_weights[e_idx] = w;
-      total_collapsed_weight += w;
-      if (e_idx >= m_collapsed)
-        throw Exception("Maximum number of possible edges exceeded.");
-      // next edge
-      e_idx += 1;
+        // reset edge_weight_to_community to all 0.0 and neighbour_comm_added to all false
+        edge_weight_to_community[u_comm] = 0.0;
+        neighbour_comm_added[u_comm] = false;
     }
   }
 
diff --git a/src/MutableVertexPartition.cpp b/src/MutableVertexPartition.cpp
index bf910ed2..74ed346e 100644
--- a/src/MutableVertexPartition.cpp
+++ b/src/MutableVertexPartition.cpp
@@ -148,6 +148,12 @@ void MutableVertexPartition::init_admin()
   this->_current_node_cache_community_from = n + 1; this->_cached_weight_from_community.resize(this->_n_communities, 0);
   this->_current_node_cache_community_to = n + 1;   this->_cached_weight_to_community.resize(this->_n_communities, 0);
   this->_current_node_cache_community_all = n + 1;  this->_cached_weight_all_community.resize(this->_n_communities, 0);
+  this->_cached_neigh_comms_all.resize(n);
+
+  if (this->get_graph()->is_directed()) {
+    this->_cached_neigh_comms_from.resize(n);
+    this->_cached_neigh_comms_to.resize(n);
+  }
 
   this->_empty_communities.clear();
 
@@ -164,9 +170,8 @@ void MutableVertexPartition::init_admin()
   size_t m = graph->ecount();
   for (size_t e = 0; e < m; e++)
   {
-    pair<size_t, size_t> endpoints = this->graph->get_endpoints(e);
-    size_t v = endpoints.first;
-    size_t u = endpoints.second;
+    size_t v, u;
+    this->graph->edge(e, v, u);
 
     size_t v_comm = this->_membership[v];
     size_t u_comm = this->_membership[u];
@@ -247,10 +252,107 @@ void MutableVertexPartition::renumber_communities()
 {
   vector<MutableVertexPartition*> partitions(1);
   partitions[0] = this;
-  this->set_membership(MutableVertexPartition::renumber_communities(partitions));
+  vector<size_t> new_comm_id = MutableVertexPartition::rank_order_communities(partitions);
+  this->relabel_communities(new_comm_id);
+}
+
+/****************************************************************************
+ Renumber the communities according to the new labels in new_comm_id.
+
+ This adjusts the internal bookkeeping as required, avoiding the more costly
+ setup required in init_admin(). In particular, this avoids recomputation of
+ weights in/from/to each community by simply assigning the previously
+ computed values to the new, relabeled communities.
+
+ For instance, a new_comm_id of <1, 2, 0> will change the labels such that
+ community 0 becomes 1, community 1 becomes 2, and community 2 becomes 0.
+*****************************************************************************/
+void MutableVertexPartition::relabel_communities(vector<size_t> const& new_comm_id) {
+  if (this->_n_communities != new_comm_id.size()) {
+    throw Exception("Problem swapping community labels. Mismatch between n_communities and new_comm_id vector.");
+  }
+
+  size_t n = this->graph->vcount();
+
+  for (size_t i = 0; i < n; i++)
+    this->_membership[i] = new_comm_id[this->_membership[i]];
+
+  this->update_n_communities();
+  size_t nbcomms = this->n_communities();
+
+  vector<double> new_total_weight_in_comm(nbcomms, 0.0);
+  vector<double> new_total_weight_from_comm(nbcomms, 0.0);
+  vector<double> new_total_weight_to_comm(nbcomms, 0.0);
+  vector<size_t> new_csize(nbcomms, 0);
+  vector<size_t> new_cnodes(nbcomms, 0);
+
+  for (size_t c = 0; c < new_comm_id.size(); c++) {
+    size_t new_c = new_comm_id[c];
+    if (this->_csize[c] > 0) {
+      new_total_weight_in_comm[new_c] = this->_total_weight_in_comm[c];
+      new_total_weight_from_comm[new_c] = this->_total_weight_from_comm[c];
+      new_total_weight_to_comm[new_c] = this->_total_weight_to_comm[c];
+      new_csize[new_c] = this->_csize[c];
+      new_cnodes[new_c] = this->_cnodes[c];
+    }
+  }
+
+  this->_total_weight_in_comm = new_total_weight_in_comm;
+  this->_total_weight_from_comm = new_total_weight_from_comm;
+  this->_total_weight_to_comm = new_total_weight_to_comm;
+  this->_csize = new_csize;
+  this->_cnodes = new_cnodes;
+
+  this->_empty_communities.clear();
+  for (size_t c = 0; c < nbcomms; c++) {
+    if (this->_csize[c] == 0) {
+      this->_empty_communities.push_back(c);
+    }
+  }
+
+  // invalidate cached weight vectors
+  this->_current_node_cache_community_from = n + 1; this->_cached_weight_from_community.resize(nbcomms, 0);
+  this->_current_node_cache_community_to = n + 1;   this->_cached_weight_to_community.resize(nbcomms, 0);
+  this->_current_node_cache_community_all = n + 1;  this->_cached_weight_all_community.resize(nbcomms, 0);
+
+  #ifdef DEBUG
+    if (this->_csize.size() < this->_n_communities ||
+        this->_cnodes.size() < this->_n_communities ||
+        this->_total_weight_in_comm.size() < this->_n_communities ||
+        this->_total_weight_to_comm.size() < this->_n_communities ||
+        this->_total_weight_from_comm.size() < this->_n_communities ||
+        this->_cached_weight_from_community.size() < this->_n_communities ||
+        this->_cached_weight_to_community.size() < this->_n_communities ||
+        this->_cached_weight_all_community.size() < this->_n_communities) {
+      cerr << "ERROR: MutableVertexPartition bookkeeping is too small after rearrange_community_labels." << endl;
+    }
+
+    this->init_admin();
+
+    for (size_t c = 0; c < this->_n_communities; c++) {
+      if (fabs(new_total_weight_in_comm[c] - this->_total_weight_in_comm[c]) > 1e-6 ||
+          fabs(new_total_weight_from_comm[c] - this->_total_weight_from_comm[c]) > 1e-6 ||
+          fabs(new_total_weight_to_comm[c] - this->_total_weight_to_comm[c]) > 1e-6 ||
+          new_csize[c] != this->_csize[c] ||
+          new_cnodes[c] != this->_cnodes[c]) {
+        cerr << "ERROR: MutableVertexPartition bookkeeping is incorrect after rearrange_community_labels." << endl;
+        cerr << "Community c has " << endl
+             << "total_weight_in_comm=" << new_total_weight_in_comm[c]
+             << " (should be " << this->_total_weight_in_comm[c] << ")" << endl
+             << "total_weight_from_comm=" << new_total_weight_from_comm[c]
+             << " (should be " << this->_total_weight_from_comm[c] << ")" << endl
+             << "total_weight_to_comm=" << new_total_weight_to_comm[c]
+             << " (should be " << this->_total_weight_to_comm[c] << ")" << endl
+             << "csize=" << new_csize[c]
+             << " (should be " << this->_csize[c] << ")" << endl
+             << "cnodes=" << new_cnodes[c]
+             << " (should be " << this->_cnodes[c] << ")" << endl;
+      }
+    }
+  #endif
 }
 
-vector<size_t> MutableVertexPartition::renumber_communities(vector<MutableVertexPartition*> partitions)
+vector<size_t> MutableVertexPartition::rank_order_communities(vector<MutableVertexPartition*> partitions)
 {
   size_t nb_layers = partitions.size();
   size_t nb_comms = partitions[0]->n_communities();
@@ -296,11 +398,7 @@ vector<size_t> MutableVertexPartition::renumber_communities(vector<MutableVertex
     delete[] csizes[i];
   }
 
-  vector<size_t> membership(n, 0);
-  for (size_t i = 0; i < n; i++)
-    membership[i] = new_comm_id[partitions[0]->_membership[i]];
-
-  return membership;
+  return new_comm_id;
 }
 
 
@@ -308,7 +406,7 @@ vector<size_t> MutableVertexPartition::renumber_communities(vector<MutableVertex
  Renumber the communities using the original fixed membership vector. Notice
  that this doesn't ensure any property of the community numbers.
 *****************************************************************************/
-vector<size_t> MutableVertexPartition::renumber_communities(map<size_t, size_t> const& membership)
+void MutableVertexPartition::renumber_communities(map<size_t, size_t> const& membership)
 {
 
   #ifdef DEBUG
@@ -317,7 +415,7 @@ vector<size_t> MutableVertexPartition::renumber_communities(map<size_t, size_t>
 
   // Skip whole thing if there are no fixed nodes for efficiency
   if (membership.size() == 0)
-    return _membership;
+    return;
 
   // The number of communities does not depend on whether some are fixed
   size_t nb_comms = n_communities();
@@ -358,17 +456,7 @@ vector<size_t> MutableVertexPartition::renumber_communities(map<size_t, size_t>
     }
   }
 
-  // Set the new communities
-  vector<size_t> new_membership(this->graph->vcount());
-  for (size_t i = 0; i < this->graph->vcount(); i++)
-  {
-    #ifdef DEBUG
-      cerr << "Setting membership of node " << i << " from" << _membership[i] << " to " << new_comm_id[_membership[i]] << endl;
-    #endif
-    new_membership[i] = new_comm_id[_membership[i]];
-  }
-
-  return new_membership;
+  this->relabel_communities(new_comm_id);
 }
 
 void MutableVertexPartition::renumber_communities(vector<size_t> const& membership)
@@ -694,48 +782,6 @@ void MutableVertexPartition::from_partition(MutableVertexPartition* partition)
   this->init_admin();
 }
 
-/****************************************************************************
- Calculate what is the total weight going from a node to a community.
-
-    Parameters:
-      v      -- The node which to check.
-      comm   -- The community which to check.
-*****************************************************************************/
-double MutableVertexPartition::weight_to_comm(size_t v, size_t comm)
-{
-  if (this->_current_node_cache_community_to != v)
-  {
-    this->cache_neigh_communities(v, IGRAPH_OUT);
-    this->_current_node_cache_community_to = v;
-  }
-
-  if (comm < this->_cached_weight_to_community.size())
-    return this->_cached_weight_to_community[comm];
-  else
-    return 0.0;
-}
-
-/****************************************************************************
- Calculate what is the total weight going from a community to a node.
-
-    Parameters:
-      v      -- The node which to check.
-      comm   -- The community which to check.
-*****************************************************************************/
-double MutableVertexPartition::weight_from_comm(size_t v, size_t comm)
-{
-  if (this->_current_node_cache_community_from != v)
-  {
-    this->cache_neigh_communities(v, IGRAPH_IN);
-    this->_current_node_cache_community_from = v;
-  }
-
-  if (comm < this->_cached_weight_from_community.size())
-    return this->_cached_weight_from_community[comm];
-  else
-    return 0.0;
-}
-
 void MutableVertexPartition::cache_neigh_communities(size_t v, igraph_neimode_t mode)
 {
   // TODO: We can probably calculate at once the IN, OUT and ALL
@@ -777,7 +823,6 @@ void MutableVertexPartition::cache_neigh_communities(size_t v, igraph_neimode_t
 
   // Reset cached neighbours
   _cached_neighs_comms->clear();
-  _cached_neighs_comms->reserve(degree);
   for (size_t idx = 0; idx < degree; idx++)
   {
     size_t u = neighbours[idx];
@@ -811,6 +856,9 @@ void MutableVertexPartition::cache_neigh_communities(size_t v, igraph_neimode_t
 
 vector<size_t> const& MutableVertexPartition::get_neigh_comms(size_t v, igraph_neimode_t mode)
 {
+  if (!this->get_graph()->is_directed())
+    mode = IGRAPH_ALL; // igraph ignores mode for undirected graphs
+
   switch (mode)
   {
     case IGRAPH_IN:
diff --git a/src/Optimiser.cpp b/src/Optimiser.cpp
index 9754a6eb..e1ae565e 100644
--- a/src/Optimiser.cpp
+++ b/src/Optimiser.cpp
@@ -349,13 +349,13 @@ double Optimiser::optimise_partition(vector<MutableVertexPartition*> partitions,
   // where r is the number of communities. The exception is fixed
   // nodes which should keep the numbers of the original communities
   q = 0.0;
-  vector<size_t> membership = MutableVertexPartition::renumber_communities(partitions);
-  partitions[0]->set_membership(membership);
-  membership = partitions[0]->renumber_communities(original_fixed_memberships);
+  partitions[0]->renumber_communities();
+  partitions[0]->renumber_communities(original_fixed_memberships);
+  vector<size_t> const& membership = partitions[0]->membership();
   // We only renumber the communities for the first graph,
   // since the communities for the other graphs should just be equal
   // to the membership of the first graph.
-  for (size_t layer = 0; layer < nb_layers; layer++)
+  for (size_t layer = 1; layer < nb_layers; layer++)
   {
     partitions[layer]->set_membership(membership);
     q += partitions[layer]->quality()*layer_weights[layer];
@@ -522,12 +522,14 @@ double Optimiser::move_nodes(vector<MutableVertexPartition*> partitions, vector<
   // (2) - The quality function should be exactly the same value after
   //       aggregating/collapsing the graph.
 
+  vector<bool> comm_added(partitions[0]->n_communities(), false);
+  vector<size_t> comms;
+
   // As long as the queue is not empty
   while(!vertex_order.empty())
   {
     size_t v = vertex_order.front(); vertex_order.pop();
-
-    set<size_t> comms;
+    comms.clear();
     Graph* graph = NULL;
     MutableVertexPartition* partition = NULL;
     // What is the current community of the node (this should be the same for all layers)
@@ -539,9 +541,9 @@ double Optimiser::move_nodes(vector<MutableVertexPartition*> partitions, vector<
       {
         for (size_t layer = 0; layer < nb_layers; layer++)
         {
-          if (partitions[layer]->cnodes(comm) > 0)
+          if (partitions[layer]->cnodes(comm) > 0 && !comm_added[comm])
           {
-            comms.insert(comm);
+            comms.push_back(comm);
             break; // Break from for loop in layer
           }
         }
@@ -553,21 +555,26 @@ double Optimiser::move_nodes(vector<MutableVertexPartition*> partitions, vector<
       /****************************ALL NEIGH COMMS*****************************/
       for (size_t layer = 0; layer < nb_layers; layer++)
       {
-        vector<size_t> const& neigh_comm_layer = partitions[layer]->get_neigh_comms(v, IGRAPH_ALL);
-        comms.insert(neigh_comm_layer.begin(), neigh_comm_layer.end());
+        for (size_t u : partitions[layer]->get_graph()->get_neighbours(v, IGRAPH_ALL)) {
+            size_t comm = partitions[layer]->membership(u);
+            if (!comm_added[comm]) {
+                comms.push_back(comm);
+                comm_added[comm] = true;
+            }
+        }
       }
     }
     else if (consider_comms == RAND_COMM)
     {
       /****************************RAND COMM***********************************/
-      comms.insert( partitions[0]->membership(graphs[0]->get_random_node(&rng)) );
+      comms.push_back( partitions[0]->membership(graphs[0]->get_random_node(&rng)) );
     }
     else if (consider_comms == RAND_NEIGH_COMM)
     {
       /****************************RAND NEIGH COMM*****************************/
       size_t rand_layer = get_random_int(0, nb_layers - 1, &rng);
       if (graphs[rand_layer]->degree(v, IGRAPH_ALL) > 0)
-        comms.insert( partitions[0]->membership(graphs[rand_layer]->get_random_neighbour(v, IGRAPH_ALL, &rng)) );
+        comms.push_back( partitions[0]->membership(graphs[rand_layer]->get_random_neighbour(v, IGRAPH_ALL, &rng)) );
     }
 
     #ifdef DEBUG
@@ -576,11 +583,10 @@ double Optimiser::move_nodes(vector<MutableVertexPartition*> partitions, vector<
 
     size_t max_comm = v_comm;
     double max_improv = 0.0;
-    for (set<size_t>::iterator comm_it = comms.begin();
-         comm_it!= comms.end();
-         comm_it++)
+    for (size_t comm : comms)
     {
-      size_t comm = *comm_it;
+      // reset comm_added to all false
+      comm_added[comm] = false;
       double possible_improv = 0.0;
 
       // Consider the improvement of moving to a community for all layers
@@ -617,6 +623,7 @@ double Optimiser::move_nodes(vector<MutableVertexPartition*> partitions, vector<
           // that is has also been added to the other layers
           for (size_t layer = 1; layer < nb_layers; layer++)
               partitions[layer]->add_empty_community();
+          comm_added.push_back(false);
         }
 
         double possible_improv = 0.0;
@@ -701,8 +708,9 @@ double Optimiser::move_nodes(vector<MutableVertexPartition*> partitions, vector<
   }
 
   partitions[0]->renumber_communities();
-  vector<size_t> const& membership = partitions[0]->renumber_communities(original_fixed_memberships);
-  for (size_t layer = 0; layer < nb_layers; layer++)
+  partitions[0]->renumber_communities(original_fixed_memberships);
+  vector<size_t> const& membership = partitions[0]->membership();
+  for (size_t layer = 1; layer < nb_layers; layer++)
   {
     partitions[layer]->set_membership(membership);
     #ifdef DEBUG
@@ -901,8 +909,9 @@ double Optimiser::merge_nodes(vector<MutableVertexPartition*> partitions, vector
   }
 
   partitions[0]->renumber_communities();
-  vector<size_t> const& membership = partitions[0]->renumber_communities(original_fixed_memberships);
-  for (size_t layer = 0; layer < nb_layers; layer++)
+  partitions[0]->renumber_communities(original_fixed_memberships);
+  vector<size_t> const& membership = partitions[0]->membership();
+  for (size_t layer = 1; layer < nb_layers; layer++)
   {
     partitions[layer]->set_membership(membership);
     #ifdef DEBUG
@@ -1177,6 +1186,9 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
 
   vector< vector<size_t> > constrained_comms = constrained_partition->get_communities();
 
+  vector<bool> comm_added(partitions[0]->n_communities(), false);
+  vector<size_t> comms;
+
   // For each node
   for (vector<size_t>::iterator it = vertex_order.begin();
        it != vertex_order.end(); it++)
@@ -1188,7 +1200,7 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
 
     if (partitions[0]->cnodes(v_comm) == 1)
     {
-      set<size_t> comms;
+      comms.clear();
       MutableVertexPartition* partition = NULL;
 
       if (consider_comms == ALL_COMMS)
@@ -1201,7 +1213,10 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
           {
             size_t u = *u_constrained_comm_it;
             size_t u_comm = partitions[0]->membership(u);
-            comms.insert(u_comm);
+            if (!comm_added[u_comm]) {
+              comm_added[u_comm] = true;
+              comms.push_back(u_comm);
+            }
           }
       }
       else if (consider_comms == ALL_NEIGH_COMMS)
@@ -1209,8 +1224,15 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
           /****************************ALL NEIGH COMMS*****************************/
           for (size_t layer = 0; layer < nb_layers; layer++)
           {
-            set<size_t> neigh_comm_layer = partitions[layer]->get_neigh_comms(v, IGRAPH_ALL, constrained_partition->membership());
-            comms.insert(neigh_comm_layer.begin(), neigh_comm_layer.end());
+            for (size_t u : partitions[layer]->get_graph()->get_neighbours(v, IGRAPH_ALL)) {
+              if (constrained_partition->membership(v) == constrained_partition->membership(u)) {
+                size_t comm = partitions[layer]->membership(u);
+                if (!comm_added[comm]) {
+                  comms.push_back(comm);
+                  comm_added[comm] = true;
+                }
+              }
+            }
           }
       }
       else if (consider_comms == RAND_COMM)
@@ -1218,7 +1240,7 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
         /****************************RAND COMM***********************************/
           size_t v_constrained_comm = constrained_partition->membership(v);
           size_t random_idx = get_random_int(0, constrained_comms[v_constrained_comm].size() - 1, &rng);
-          comms.insert(constrained_comms[v_constrained_comm][random_idx]);
+          comms.push_back(constrained_comms[v_constrained_comm][random_idx]);
       }
       else if (consider_comms == RAND_NEIGH_COMM)
       {
@@ -1239,7 +1261,7 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
             if (get_random_int(0, k, &rng) > 0)
             {
               size_t random_idx = get_random_int(0, k - 1, &rng);
-              comms.insert(all_neigh_comms_incl_dupes[random_idx]);
+              comms.push_back(all_neigh_comms_incl_dupes[random_idx]);
             }
           }
       }
@@ -1250,11 +1272,10 @@ double Optimiser::merge_nodes_constrained(vector<MutableVertexPartition*> partit
 
       size_t max_comm = v_comm;
       double max_improv = 0.0;
-      for (set<size_t>::iterator comm_it = comms.begin();
-           comm_it!= comms.end();
-           comm_it++)
+      for (size_t comm : comms)
       {
-        size_t comm = *comm_it;
+        // reset comm_added to all false
+        comm_added[comm] = false;
         double possible_improv = 0.0;
 
         // Consider the improvement of moving to a community for all layers
diff --git a/tests/test_Optimiser.py b/tests/test_Optimiser.py
index ec4d4d2c..24066c3a 100644
--- a/tests/test_Optimiser.py
+++ b/tests/test_Optimiser.py
@@ -74,14 +74,38 @@ def test_optimiser_with_fixed_nodes(self):
               initial_membership=[2, 1, 0])
       # Equivalent to setting initial membership
       #partition.set_membership([2, 1, 2])
-      opt = leidenalg.Optimiser()
       fixed_nodes = [True, False, False]
-      opt.optimise_partition(partition, fixed_nodes=fixed_nodes)
+      original_quality = partition.quality()
+      diff = self.optimiser.optimise_partition(partition, fixed_nodes=fixed_nodes)
+      self.assertAlmostEqual(partition.quality() - original_quality, diff, places=10,
+                             msg="Optimisation with fixed nodes returns inconsistent quality")
       self.assertListEqual(
             partition.membership, [2, 2, 2],
             msg="After optimising partition with fixed nodes failed to recover initial fixed memberships"
             )
 
+  def test_optimiser_fixed_nodes_large_labels(self):
+    G = ig.Graph.Erdos_Renyi(n=100, p=5./100, directed=True, loops=True)
+
+    membership = list(range(G.vcount()))
+    partition = leidenalg.RBConfigurationVertexPartition(G, initial_membership=membership)
+
+    # large enough to force nonconsecutive labels in the final partition
+    fixed_node_idx = 90
+    fixed_nodes = [False] * G.vcount()
+    fixed_nodes[fixed_node_idx] = True
+
+    original_quality = partition.quality()
+    diff = self.optimiser.optimise_partition(partition, fixed_nodes=fixed_nodes)
+
+    self.assertLess(len(set(partition.membership)), len(partition),
+                    msg="Optimisation with fixed nodes yielded too many communities")
+    self.assertAlmostEqual(partition.quality() - original_quality, diff, places=10,
+                           msg="Optimisation with fixed nodes returned inconsistent quality")
+    self.assertEqual(partition.membership[fixed_node_idx], fixed_node_idx,
+                     msg="Optimisation with fixed nodes failed to keep the associated community labels fixed")
+
+
   def test_neg_weight_bipartite(self):
     G = ig.Graph.Full_Bipartite(50, 50);
     G.es['weight'] = -0.1;