diff --git a/src/main/cpp/ds/tree/IntRootedForest.hpp b/src/main/cpp/ds/tree/IntRootedForest.hpp
index 7b26ef5..6c556f8 100644
--- a/src/main/cpp/ds/tree/IntRootedForest.hpp
+++ b/src/main/cpp/ds/tree/IntRootedForest.hpp
@@ -12,6 +12,13 @@
 //================================================================================
 //    Macros
 //================================================================================
+
+// validation: on
+// #define VALIDATE(s) s
+
+// validation: off
+#define VALIDATE(s)
+
 #define FOR_EACH_CHILD(c, p) for (auto c = nodes_[(p)].first_child; (c) != NOT_AVAILABLE; (c) = nodes_[(c)].right)
 
 namespace ds {
@@ -109,19 +116,19 @@ class IntRootedForest {
   //================================================================================
   // array subscript operator for writing
   Node& operator[](std::size_t index) {
-    if (!is_valid(index)) throw std::invalid_argument("invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("invalid index"));
     return nodes_[index];
   }
 
   // array subscript operator for reading
   Node const& operator[](std::size_t index) const {
-    if (!is_valid(index)) throw std::invalid_argument("invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("invalid index"));
     return nodes_[index];
   }
 
   // get nodes in the one-level lower
   std::vector<int> get_children(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("get_children: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("get_children: invalid index"));
 
     std::vector<int> ret;
     FOR_EACH_CHILD(c, index) ret.push_back(c);
@@ -141,7 +148,7 @@ class IntRootedForest {
    *         leaving:  traverse from node to parent
    */
   std::vector<std::pair<int, int>> dfs_preorder_edges(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("dfs_preorder_edges: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("dfs_preorder_edges: invalid index"));
 
     std::vector<std::pair<int, int>> ret, stack;
     stack.push_back({index, false});
@@ -162,7 +169,7 @@ class IntRootedForest {
   }
 
   std::vector<std::pair<int, int>> dfs_reverse_preorder_edges(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("dfs_preorder_edges: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("dfs_preorder_edges: invalid index"));
 
     std::vector<std::pair<int, int>> ret, stack;
     stack.push_back({index, false});
@@ -183,7 +190,7 @@ class IntRootedForest {
   }
 
   std::vector<int> bfs_nodes(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("bfs_nodes: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("bfs_nodes: invalid index"));
 
     std::vector<int> ret;
     std::queue<int> q;
@@ -199,7 +206,7 @@ class IntRootedForest {
   }
 
   std::vector<int> dfs_preorder_nodes(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("dfs_preorder_nodes: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("dfs_preorder_nodes: invalid index"));
 
     std::vector<int> ret, stack;
     stack.push_back(index);
@@ -216,7 +223,7 @@ class IntRootedForest {
   }
 
   std::vector<int> dfs_reverse_preorder_nodes(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("dfs_reverse_preorder_nodes: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("dfs_reverse_preorder_nodes: invalid index"));
 
     std::vector<int> ret, stack;
     stack.push_back(index);
@@ -238,7 +245,7 @@ class IntRootedForest {
    * @return std::vector<int>
    */
   std::vector<int> get_leaves(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("get_leaves: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("get_leaves: invalid index"));
 
     std::vector<int> ret;
     for (auto x : dfs_reverse_preorder_nodes(index)) {
@@ -248,7 +255,7 @@ class IntRootedForest {
   }
 
   std::vector<int> get_ancestors(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("get_ancestors: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("get_ancestors: invalid index"));
 
     std::vector<int> ret;
     for (auto p = nodes_[index].parent; p != NOT_AVAILABLE; p = nodes_[p].parent) ret.push_back(p);
@@ -256,7 +263,7 @@ class IntRootedForest {
   }
 
   int get_root(int index) const {
-    if (!is_valid(index)) throw std::invalid_argument("get_root: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("get_root: invalid index"));
 
     int ret = NOT_AVAILABLE;
     for (auto p = index; p != NOT_AVAILABLE; p = nodes_[p].parent) ret = p;
@@ -311,9 +318,9 @@ class IntRootedForest {
   }
 
   void remove(int index) {
-    if (!is_valid(index)) throw std::invalid_argument("remove: invalid index");
+    VALIDATE(if (!is_valid(index)) throw std::invalid_argument("remove: invalid index"));
     detach(index);
-    if (!nodes_[index].is_leaf()) throw std::invalid_argument("remove: must be a leaf");
+    VALIDATE(if (!nodes_[index].is_leaf()) throw std::invalid_argument("remove: must be a leaf"));
 
     --num_live_nodes_;
     nodes_[index].alive = false;
@@ -325,12 +332,16 @@ class IntRootedForest {
   //================================================================================
  private:
   void add_child(int parent, int child) {
-    if (!is_valid(parent)) throw std::invalid_argument("add_child: parent invalid index");
-    if (!is_valid(child)) throw std::invalid_argument("add_child: child invalid index");
+    VALIDATE({
+      if (!is_valid(parent)) throw std::invalid_argument("add_child: parent invalid index");
+      if (!is_valid(child)) throw std::invalid_argument("add_child: child invalid index");
+    });
     auto& p = nodes_[parent];
     auto& c = nodes_[child];
 
-    if (!c.is_root()) throw std::invalid_argument("add_child: child must be a root");
+    VALIDATE({
+      if (!c.is_root()) throw std::invalid_argument("add_child: child must be a root");
+    });
 
     if (p.has_child()) {
       nodes_[p.first_child].left = child;
@@ -344,7 +355,9 @@ class IntRootedForest {
 
  public:
   void detach(int index) {
-    if (!is_valid(index)) throw std::invalid_argument("detach: invalid index");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("detach: invalid index");
+    });
 
     auto& node = nodes_[index];
     if (node.parent != NOT_AVAILABLE) nodes_[node.parent].num_children--;
@@ -358,10 +371,11 @@ class IntRootedForest {
   };
 
   void swap(int a, int b) {
-    if (!is_valid(a)) throw std::invalid_argument("swap: a invalid index");
-    if (!is_valid(b)) throw std::invalid_argument("swap: b invalid index");
-    if (get_root(a) == get_root(b)) throw std::invalid_argument("swap: a and b must belong to different trees");
-
+    VALIDATE({
+      if (!is_valid(a)) throw std::invalid_argument("swap: a invalid index");
+      if (!is_valid(b)) throw std::invalid_argument("swap: b invalid index");
+      if (get_root(a) == get_root(b)) throw std::invalid_argument("swap: a and b must belong to different trees")
+    });
     // if (a == b) return; // never happens
 
     auto& na = nodes_[a];
@@ -386,20 +400,25 @@ class IntRootedForest {
    * @param replace_by not to replace
    */
   void replace(int index, int replace_by) {
-    if (!is_valid(index)) throw std::invalid_argument("replace: invalid index");
-    if (!is_valid(replace_by)) throw std::invalid_argument("replace: replace_by invalid index");
-    if (index == replace_by) throw std::invalid_argument("replace: replace_by must differ from index");
-    if (util::contains(get_ancestors(index), replace_by))
-      throw std::invalid_argument("replace: replace_by cannot be an ancestor of index");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("replace: invalid index");
+      if (!is_valid(replace_by)) throw std::invalid_argument("replace: replace_by invalid index");
+      if (index == replace_by) throw std::invalid_argument("replace: replace_by must differ from index");
+      if (util::contains(get_ancestors(index), replace_by)) {
+        throw std::invalid_argument("replace: replace_by cannot be an ancestor of index");
+      }
+    });
 
     detach(replace_by);
     swap(index, replace_by);
   }
 
   void move_to(int index, int new_parent) {
-    if (!is_valid(index)) throw std::invalid_argument("move_to: invalid index");
-    if (!is_valid(new_parent)) throw std::invalid_argument("move_to: new_parent invalid index");
-    if (index == new_parent) throw std::invalid_argument("move_to: index and new_parent cannot be the same");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("move_to: invalid index");
+      if (!is_valid(new_parent)) throw std::invalid_argument("move_to: new_parent invalid index");
+      if (index == new_parent) throw std::invalid_argument("move_to: index and new_parent cannot be the same");
+    });
 
     detach(index);
     add_child(new_parent, index);
@@ -410,12 +429,15 @@ class IntRootedForest {
    * @param node node that is not a root
    */
   void move_to_before(int index, int target) {
-    if (!is_valid(index)) throw std::invalid_argument("move_to_after: invalid index");
-    if (!is_valid(target)) throw std::invalid_argument("move_to_after: target invalid index");
-    if (nodes_[target].is_root()) throw std::invalid_argument("move_to_before: target must not be a root");
-    if (index == target) throw std::invalid_argument("move_to_after: index and target cannot be the same");
-    if (util::contains(get_ancestors(target), index))
-      throw std::invalid_argument("replace: index cannot be an ancestor of target");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("move_to_after: invalid index");
+      if (!is_valid(target)) throw std::invalid_argument("move_to_after: target invalid index");
+      if (nodes_[target].is_root()) throw std::invalid_argument("move_to_before: target must not be a root");
+      if (index == target) throw std::invalid_argument("move_to_after: index and target cannot be the same");
+      if (util::contains(get_ancestors(target), index)) {
+        throw std::invalid_argument("replace: index cannot be an ancestor of target");
+      }
+    });
 
     detach(index);
 
@@ -437,12 +459,15 @@ class IntRootedForest {
    * @param node node that is not a root
    */
   void move_to_after(int index, int target) {
-    if (!is_valid(index)) throw std::invalid_argument("move_to_after: invalid index");
-    if (!is_valid(target)) throw std::invalid_argument("move_to_after: target invalid index");
-    if (nodes_[target].is_root()) throw std::invalid_argument("move_to_after: target must not be a root");
-    if (index == target) throw std::invalid_argument("move_to_after: index and target cannot be the same");
-    if (util::contains(get_ancestors(target), index))
-      throw std::invalid_argument("move_to_after: index cannot be an ancestor of target");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("move_to_after: invalid index");
+      if (!is_valid(target)) throw std::invalid_argument("move_to_after: target invalid index");
+      if (nodes_[target].is_root()) throw std::invalid_argument("move_to_after: target must not be a root");
+      if (index == target) throw std::invalid_argument("move_to_after: index and target cannot be the same");
+      if (util::contains(get_ancestors(target), index)) {
+        throw std::invalid_argument("move_to_after: index cannot be an ancestor of target");
+      }
+    });
 
     detach(index);
 
@@ -462,7 +487,9 @@ class IntRootedForest {
    * @brief Moves this node to the first among all its siblings.
    */
   void make_first_child(int index) {
-    if (!is_valid(index)) throw std::invalid_argument("make_first_child: invalid index");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("make_first_child: invalid index");
+    });
 
     if (nodes_[index].is_root() || nodes_[index].is_first_child()) return;  // do nothing
 
@@ -474,10 +501,13 @@ class IntRootedForest {
    * @param node node whose children are to be added to this node's children
    */
   void add_children_from(int index, int target) {
-    if (!is_valid(index)) throw std::invalid_argument("add_children_from: invalid index");
-    if (!is_valid(target)) throw std::invalid_argument("add_children_from: target invalid index");
-    if (util::contains(get_ancestors(index), target))
-      throw std::invalid_argument("add_children_from: target cannot be an ancestor of index");
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("add_children_from: invalid index");
+      if (!is_valid(target)) throw std::invalid_argument("add_children_from: target invalid index");
+      if (util::contains(get_ancestors(index), target)) {
+        throw std::invalid_argument("add_children_from: target cannot be an ancestor of index");
+      }
+    });
 
     if (index == target) return;  // do nothing
 
@@ -504,12 +534,17 @@ class IntRootedForest {
    * This original node will be detached form its tree but not removed.
    */
   void replace_by_children(int index) {
-    if (!is_valid(index)) throw std::invalid_argument("replace_by_children: invalid index");
-    auto& node = nodes_[index];
-    if (node.is_root()) throw std::invalid_argument("replace_by_children:this must not be a root");
-
-    auto ch = get_children(index);
-    for (auto c : ch) move_to_before(c, index);
+    VALIDATE({
+      if (!is_valid(index)) throw std::invalid_argument("replace_by_children: invalid index");
+      auto& node = nodes_[index];
+      if (node.is_root()) throw std::invalid_argument("replace_by_children:this must not be a root");
+    });
+
+    for (auto c = nodes_[index].first_child; c != NOT_AVAILABLE;) {
+      auto nxt = nodes_[c].right;
+      move_to_before(c, index);
+      c = nxt;
+    }
     detach(index);
   }
 
diff --git a/src/main/cpp/modular-bench.cpp b/src/main/cpp/modular-bench.cpp
index 9321d0b..5110aa4 100644
--- a/src/main/cpp/modular-bench.cpp
+++ b/src/main/cpp/modular-bench.cpp
@@ -10,7 +10,13 @@ int main(int argc, char* argv[]) {
   auto graph = readwrite::read_edge_list(std::cin);
 
   // run algorithm
+#if PROFILE_ON
+  util::Profiler prof;
+  auto result = modular::modular_decomposition_time(graph, true, &prof);
+  prof.print();
+#else
   auto result = modular::modular_decomposition_time(graph, true);
+#endif
 
   // output result
   printf("%d\n", result.first.modular_width());
diff --git a/src/main/cpp/modular/MDTree.cpp b/src/main/cpp/modular/MDTree.cpp
index 77fc29b..7095a82 100644
--- a/src/main/cpp/modular/MDTree.cpp
+++ b/src/main/cpp/modular/MDTree.cpp
@@ -8,9 +8,9 @@ std::ostream &operator<<(std::ostream &os, MDNode const &node) { return os << no
 
 MDTree modular_decomposition(ds::graph::Graph const &graph, bool sorted) { return MDTree(graph, sorted); }
 
-std::pair<MDTree, double> modular_decomposition_time(ds::graph::Graph const &graph, bool sorted) {
+std::pair<MDTree, double> modular_decomposition_time(ds::graph::Graph const &graph, bool sorted, util::Profiler *prof) {
   auto time_start = std::chrono::system_clock::now();
-  auto ret = MDTree(graph, false);
+  auto ret = MDTree(graph, false, prof);
   auto time_finish = std::chrono::system_clock::now();
   auto elapsed = time_finish - time_start;
   auto elapsed_sec = std::chrono::duration_cast<std::chrono::nanoseconds>(elapsed).count() * 1e-9;
diff --git a/src/main/cpp/modular/MDTree.hpp b/src/main/cpp/modular/MDTree.hpp
index 4711ba0..dac1ec0 100644
--- a/src/main/cpp/modular/MDTree.hpp
+++ b/src/main/cpp/modular/MDTree.hpp
@@ -1,5 +1,7 @@
 #pragma once
 
+#include "util/profiler.hpp"
+
 #include "compute/MDSolver.hpp"
 
 namespace modular {
@@ -47,8 +49,8 @@ class MDTree {
 
  public:
   MDTree() : root_(-1){};
-  MDTree(ds::graph::Graph const &graph, bool sorted = false) {
-    auto result = compute::MDSolver::compute(graph);
+  MDTree(ds::graph::Graph const &graph, bool sorted = false, util::Profiler *prof = nullptr) {
+    auto result = compute::MDSolver::compute(graph, prof);
     if (result.second >= 0) {
       *this = MDTree(result.first, result.second);
       if (sorted) this->sort();
@@ -160,5 +162,6 @@ class MDTree {
 
 MDTree modular_decomposition(ds::graph::Graph const &graph, bool sorted = false);
 
-std::pair<MDTree, double> modular_decomposition_time(ds::graph::Graph const &graph, bool sorted = false);
+std::pair<MDTree, double> modular_decomposition_time(ds::graph::Graph const &graph, bool sorted = false,
+                                                     util::Profiler *prof = nullptr);
 }  // namespace modular
diff --git a/src/main/cpp/modular/compute/MDComputeNode.hpp b/src/main/cpp/modular/compute/MDComputeNode.hpp
index cc0bce5..d24c512 100644
--- a/src/main/cpp/modular/compute/MDComputeNode.hpp
+++ b/src/main/cpp/modular/compute/MDComputeNode.hpp
@@ -19,7 +19,9 @@ class MDComputeNode {
   VertexID vertex;  // pivot for problem node
   int comp_number;
   int tree_number;
-  // int num_marks;
+  int num_marks;
+  int num_left_split_children;   // number of the children with split LEFT or MIXED; used only in refinement
+  int num_right_split_children;  // number of the children with split RIGHT or MIXED; used only in refinement
   bool active;
   bool connected;
 
@@ -30,7 +32,9 @@ class MDComputeNode {
         vertex(-1),
         comp_number(-1),
         tree_number(-1),
-        // num_marks(0),
+        num_marks(0),
+        num_left_split_children(0),
+        num_right_split_children(0),
         active(false),
         connected(false) {}
 
@@ -41,7 +45,9 @@ class MDComputeNode {
         vertex(node.vertex),
         comp_number(node.comp_number),
         tree_number(node.tree_number),
-        // num_marks(0),
+        num_marks(node.num_marks),
+        num_left_split_children(0),   // initialize to zero, not copying from the original
+        num_right_split_children(0),  // initialize to zero, not copying from the original
         active(node.active),
         connected(node.connected) {}
 
@@ -67,13 +73,15 @@ class MDComputeNode {
   bool is_operation_node() const { return node_type == NodeType::OPERATION; }
   bool is_problem_node() const { return node_type == NodeType::PROBLEM; }
 
-  // bool is_marked() { return num_marks > 0; }
-  // void add_mark() { ++num_marks; }
-  // void clear_marks() { num_marks = 0; }
+  bool is_marked() const { return num_marks > 0; }
+  void add_mark() { ++num_marks; }
+  int number_of_marks() const { return num_marks; }
+  void clear_marks() { num_marks = 0; }
 
   bool is_split_marked(SplitDirection split_type) {
     return this->split_type == SplitDirection::MIXED || this->split_type == split_type;
   }
+
   void set_split_mark(SplitDirection split_type) {
     if (this->split_type == split_type) {
       // already set
@@ -84,11 +92,34 @@ class MDComputeNode {
     }
   }
 
+  void increment_num_split_children(SplitDirection split_type) {
+    if (split_type == SplitDirection::LEFT) {
+      ++num_left_split_children;
+    } else {
+      ++num_right_split_children;
+    }
+  }
+
+  void decrement_num_split_children(SplitDirection split_type) {
+    if (split_type == SplitDirection::LEFT) {
+      --num_left_split_children;
+    } else {
+      --num_right_split_children;
+    }
+  }
+
+  int get_num_split_children(SplitDirection split_type) const {
+    return split_type == SplitDirection::LEFT ? num_left_split_children : num_right_split_children;
+  }
+
+  void clear_num_split_children() { num_left_split_children = num_right_split_children = 0; }
+
   void clear() {
     comp_number = -1;
     tree_number = -1;
     // num_marks = 0;
     split_type = SplitDirection::NONE;
+    clear_num_split_children();
   }
 
   std::string to_string() const {
diff --git a/src/main/cpp/modular/compute/MDSolver.cpp b/src/main/cpp/modular/compute/MDSolver.cpp
index a0b2ad7..53ac3e4 100644
--- a/src/main/cpp/modular/compute/MDSolver.cpp
+++ b/src/main/cpp/modular/compute/MDSolver.cpp
@@ -11,10 +11,10 @@ int compute(ds::graph::Graph const &graph, CompTree &tree, int main_prob, util::
   int n = graph.number_of_nodes();
   int current_prob = main_prob;
 
-  VI alpha_list[n];
+  VVV alpha_list(n);
+  VVV fp_neighbors(n);  // used only for assembly -> compute_fact_perm_edges()
   bool visited[n];
   for (int i = 0; i < n; ++i) {
-    alpha_list[i] = VI();
     visited[i] = false;
   }
   ds::FastSet vset(n);
@@ -35,14 +35,16 @@ int compute(ds::graph::Graph const &graph, CompTree &tree, int main_prob, util::
 
     auto &fc = tree[cp.first_child];
 
-    if (fc.data.node_type != NodeType::PROBLEM) {
+    if (!fc.data.is_problem_node()) {
       // first, needs to solve subproblems
       visited[cp.first_child] = true;  // set visited
 
       if (cp.has_only_one_child()) {
         // base case
         PROF(util::pcount(prof, "solve(): base case"));
+        PROF(util::pstart(prof, "process_neighbors()"));
         process_neighbors(graph, tree, alpha_list, visited, cp.first_child, current_prob, -1);
+        PROF(util::pstop(prof, "process_neighbors()"));
       } else {
         // pivot at the first child
         PROF(util::pstart(prof, "do_pivot()"));
@@ -65,11 +67,12 @@ int compute(ds::graph::Graph const &graph, CompTree &tree, int main_prob, util::
       PROF(util::pstop(prof, "remove_layers()"));
 
       PROF(util::pstart(prof, "complete_alpha_lists()"));
-      complete_alpha_lists(tree, alpha_list, vset, current_prob);
+      std::vector<int> leaves = tree.get_leaves(current_prob);
+      complete_alpha_lists(tree, alpha_list, vset, current_prob, leaves);
       PROF(util::pstop(prof, "complete_alpha_lists()"));
 
       PROF(util::pstart(prof, "refine()"));
-      refine(tree, alpha_list, current_prob);
+      refine(tree, alpha_list, current_prob, leaves, prof);
       PROF(util::pstop(prof, "refine()"));
 
       PROF(util::pstart(prof, "promote()"));
@@ -77,18 +80,20 @@ int compute(ds::graph::Graph const &graph, CompTree &tree, int main_prob, util::
       PROF(util::pstop(prof, "promote()"));
 
       PROF(util::pstart(prof, "assemble()"));
-      assemble(tree, alpha_list, current_prob);
+      assemble(tree, alpha_list, current_prob, fp_neighbors, vset, prof);
       PROF(util::pstop(prof, "assemble()"));
 
-      PROF(util::pstart(prof, "merge_components()"));
-      merge_components(tree, current_prob, extra_components);
-      PROF(util::pstop(prof, "merge_components()"));
-
       // clear all but visited
+      PROF(util::pstart(prof, "clear all but visited"));
       for (auto c: tree.dfs_reverse_preorder_nodes(tree[current_prob].first_child)) {
         if (tree[c].is_leaf()) alpha_list[c].clear();
         tree[c].data.clear();
       }
+      PROF(util::pstop(prof, "clear all but visited"));
+
+      PROF(util::pstart(prof, "merge_components()"));
+      merge_components(tree, current_prob, extra_components);
+      PROF(util::pstop(prof, "merge_components()"));
     }
 
     result = tree[current_prob].first_child;
diff --git a/src/main/cpp/modular/compute/MDSolver.hpp b/src/main/cpp/modular/compute/MDSolver.hpp
index 97e70ae..e1b9127 100644
--- a/src/main/cpp/modular/compute/MDSolver.hpp
+++ b/src/main/cpp/modular/compute/MDSolver.hpp
@@ -27,7 +27,7 @@ int compute(ds::graph::Graph const &graph, CompTree &tree, int main_prob, util::
 void process_neighbors(             //
     ds::graph::Graph const &graph,  //
     CompTree &tree,                 //
-    VI alpha_list[],                //
+    VVV &alpha_list,                //
     bool const visited[],           //
     VertexID pivot,                 //
     int current_prob,               //
@@ -35,19 +35,19 @@ void process_neighbors(             //
 );
 int do_pivot(ds::graph::Graph const &graph,  //
              CompTree &tree,                 //
-             VI alpha_list[],                //
+             VVV &alpha_list,                //
              bool const visited[],           //
              int prob,                       //
              VertexID pivot                  //
 );
 int remove_extra_components(CompTree &tree, int prob);
 void remove_layers(CompTree &tree, int prob);
-void complete_alpha_lists(CompTree &tree, VI alpha_list[], ds::FastSet &vset, int prob);
+void complete_alpha_lists(CompTree &tree, VVV &alpha_list, ds::FastSet &vset, int prob, std::vector<int> &leaves);
 void merge_components(CompTree &tree, int problem, int new_components);
 
-void refine(CompTree &tree, VI const alpha_list[], int prob);
+void refine(CompTree &tree, VVV const &alpha_list, int prob, std::vector<int> &leaves, util::Profiler *prof);
 void promote(CompTree &tree, int prob);
-void assemble(CompTree &tree, VI const *alpha_list, int prob);
+void assemble(CompTree &tree, VVV const &alpha_list, int prob, VVV &fp_neighbors, ds::FastSet &vset, util::Profiler *prof);
 }  // namespace impl
 
 class MDSolver {
diff --git a/src/main/cpp/modular/compute/assembly.cpp b/src/main/cpp/modular/compute/assembly.cpp
index b61ed6f..68bf2d8 100644
--- a/src/main/cpp/modular/compute/assembly.cpp
+++ b/src/main/cpp/modular/compute/assembly.cpp
@@ -39,7 +39,7 @@ static std::vector<bool> determine_right_comp_fragments(CompTree const &tree, VI
  * @note Read: MDNode::tree_number for each root and leaf
  *             MDLeaf::alpha       for each leaf
  */
-static std::vector<bool> determine_right_layer_neighbor(CompTree const &tree, VI const *alpha_list, VI const &ps, int pivot_index) {
+static std::vector<bool> determine_right_layer_neighbor(CompTree const &tree, VVV const &alpha_list, VI const &ps, int pivot_index) {
   std::vector<bool> ret(ps.size());
   for (int i = pivot_index + 1; i < static_cast<int>(ps.size()); ++i) {
     int current_tree = ps[i];
@@ -64,48 +64,56 @@ static std::vector<bool> determine_right_layer_neighbor(CompTree const &tree, VI
 /**
  * @brief Computes the edges between factorized-permutation elements.
  *
+ * @param fp_neighbors vector of neighboring factorized-permutation indices
+ * 
  * @note Read  : MDLeaf::alpha       for each leaf
  *       Update: MDLeaf::comp_number for each leaf
- * @return vector of neighboring factorized-permutation indices
+ * 
  */
-static VVV compute_fact_perm_edges(CompTree &tree, VI const *alpha_list, VI const &ps) {
+static void compute_fact_perm_edges(CompTree &tree, VVV const &alpha_list, VI const &ps, int pivot_index,
+                                    ds::FastSet &vset, VVV &fp_neighbors) {
   // TRACE("start: %s\n", to_string().c_str());
   int k = static_cast<int>(ps.size());
-  VVV neighbors(k);
-  std::vector<int> element_sizes(k);
+  for (int i = 0; i < pivot_index; ++i) fp_neighbors[i].clear();
+  std::vector<std::vector<int>> leaves(pivot_index);
 
   // initialize
   for (int i = 0; i < k; ++i) {
-    for (auto leaf : tree.get_leaves(ps[i])) {
-      tree[leaf].data.comp_number = i;  // reset the comp number to index
-      ++element_sizes[i];               // increment the element size
+    if (i < pivot_index) {
+      leaves[i] = tree.get_leaves(ps[i]);
+      for (auto leaf : leaves[i]) tree[leaf].data.comp_number = i;  // reset the comp number to index
+    } else {
+      for (auto leaf : tree.get_leaves(ps[i])) tree[leaf].data.comp_number = i;  // reset the comp number to index
     }
   }
 
-  // find joins
-  for (int i = 0; i < k; ++i) {
-    std::vector<int> candidates;
-    std::vector<int> marks(k);
+  // we need the neighbors only up to pivot_index
+  for (int i = 0; i < pivot_index; ++i) {
+    vset.clear();  // used for making candidates unique
 
     // enumerate all edges
-    for (auto leaf : tree.get_leaves(ps[i])) {
+    bool done = false;
+    for (auto leaf : leaves[i]) {
       for (auto a : alpha_list[leaf]) {
         int j = tree[a].data.comp_number;
-        candidates.push_back(j);
-        ++marks[j];
-      }
-    }
 
-    for (auto j : candidates) {
-      if (element_sizes[i] * element_sizes[j] == marks[j]) {
-        // found a join between i and j
-        neighbors[i].push_back(j);
-        marks[j] = 0;  // reset marks so that there will be no duplicates
+        if (!vset.get(j)) {
+          fp_neighbors[i].push_back(j);
+          vset.set(j);
+          if (static_cast<int>(fp_neighbors[i].size()) == k - pivot_index) {  // reached the max possible
+            done = true;
+            break;
+          }
+        }
       }
+      if (done) break;
     }
   }
+  // TODO: Prove:
+  //    If v <- alpha_list[u] with u in component i and v in component j,
+  //    then components i and j must form a full join.
+
   // TRACE("finish: %s\n", to_string().c_str());
-  return neighbors;
 }
 
 //================================================================================
@@ -311,7 +319,7 @@ static void remove_degenerate_duplicates(CompTree &tree, int index) {
 //================================================================================
 //    Main process
 //================================================================================
-void assemble(CompTree &tree, VI const *alpha_list, int prob) {
+void assemble(CompTree &tree, VVV const &alpha_list, int prob, VVV &fp_neighbors, ds::FastSet &vset, util::Profiler *prof) {
   if (tree[prob].is_leaf()) throw std::invalid_argument("roots must not be empty");
 
   // build permutation
@@ -326,14 +334,31 @@ void assemble(CompTree &tree, VI const *alpha_list, int prob) {
   if (pivot_index < 0) throw std::invalid_argument("roots must include a pivot");
 
   // main logic
+  // PROF(util::pstart(prof, "assemble: determine"));
   auto lcocomp = determine_left_cocomp_fragments(tree, ps, pivot_index);
   auto rcomp = determine_right_comp_fragments(tree, ps, pivot_index);
   auto rlayer = determine_right_layer_neighbor(tree, alpha_list, ps, pivot_index);
-  auto neighbors = compute_fact_perm_edges(tree, alpha_list, ps);
-  auto mu = compute_mu(tree, ps, pivot_index, neighbors);
+  // PROF(util::pstop(prof, "assemble: determine"));
+
+  // PROF(util::pstart(prof, "assemble: fact perm"));
+  compute_fact_perm_edges(tree, alpha_list, ps, pivot_index, vset, fp_neighbors);
+  // PROF(util::pstop(prof, "assemble: fact perm"));
+
+  // PROF(util::pstart(prof, "assemble: compute mu"));
+  auto mu = compute_mu(tree, ps, pivot_index, fp_neighbors);
+  // PROF(util::pstop(prof, "assemble: compute mu"));
+
+  // PROF(util::pstart(prof, "assemble: delineate"));
   auto boundaries = delineate(pivot_index, lcocomp, rcomp, rlayer, mu);
+  // PROF(util::pstop(prof, "assemble: delineate"));
+
+  // PROF(util::pstart(prof, "assemble: assemble tree"));
   auto root = assemble_tree(tree, ps, pivot_index, boundaries);
+  // PROF(util::pstop(prof, "assemble: assemble tree"));
+
+  // PROF(util::pstart(prof, "assemble: remove dup"));
   remove_degenerate_duplicates(tree, root);
+  // PROF(util::pstop(prof, "assemble: remove dup"));
 
   // replace the problem node with the result
   tree.replace_children(prob, root);
diff --git a/src/main/cpp/modular/compute/misc.cpp b/src/main/cpp/modular/compute/misc.cpp
index 0a6973e..c990ffd 100644
--- a/src/main/cpp/modular/compute/misc.cpp
+++ b/src/main/cpp/modular/compute/misc.cpp
@@ -30,9 +30,11 @@ int remove_extra_components(CompTree &tree, int prob) {
 void remove_layers(CompTree &tree, int prob) {
   TRACE("start: %s", tree.to_string(prob).c_str());
 
-  for (auto c : tree.get_children(prob)) {
-    tree.replace_by_children(c);
+  for (auto c = tree[prob].first_child; tree.is_valid(c);) {
+    auto nxt = tree[c].right;
+    tree.replace(c, tree[c].first_child);  // should have only one child
     tree.remove(c);
+    c = nxt;
   }
 
   TRACE("finish: %s", tree.to_string(prob).c_str());
@@ -41,28 +43,35 @@ void remove_layers(CompTree &tree, int prob) {
 /**
  * @brief Makes alpha lists in this subproblem symmetric and irredundant.
  */
-void complete_alpha_lists(CompTree &tree, VI alpha_list[], ds::FastSet &vset, int prob) {
+void complete_alpha_lists(CompTree &tree, VVV &alpha_list, ds::FastSet &vset, int prob, std::vector<int> &leaves) {
   TRACE("start: %s", tree.to_string(prob).c_str());
 
   // complete the list
-  for (auto v : tree.get_leaves(prob)) {
+  for (auto v : leaves) {
     assert(v >= 0);
-    for (int a : alpha_list[v]) {
+    for (auto a : alpha_list[v]) {
       assert(a >= 0);
       alpha_list[a].push_back(v);
     }
   }
-  // remove duplicate entries
-  for (auto v : tree.get_leaves(prob)) {
+
+  // remove duplicate entries (in-place)
+  for (auto v : leaves) {
+    auto &vs = alpha_list[v];
     vset.clear();
-    std::vector<int> result;
-    for (auto a : alpha_list[v]) {
-      if (!vset.get(a)) {
-        result.push_back(a);
+
+    std::size_t len = vs.size();
+    for (std::size_t i = 0; i < len;) {
+      auto a = vs[i];
+      if (vset.get(a)) {
+        // found a duplicate; swap with the last element
+        vs[i] = vs[--len];
+        vs.pop_back();
+      } else {
         vset.set(a);
+        ++i;
       }
     }
-    alpha_list[v] = result;
   }
 }
 
diff --git a/src/main/cpp/modular/compute/pivot.cpp b/src/main/cpp/modular/compute/pivot.cpp
index f7f4958..0d1feae 100644
--- a/src/main/cpp/modular/compute/pivot.cpp
+++ b/src/main/cpp/modular/compute/pivot.cpp
@@ -6,10 +6,10 @@ namespace compute {
 namespace impl {
 
 static bool is_pivot_layer(CompTree& tree, int index) {
-  auto &node = tree[index];
+  auto& node = tree[index];
 
   if (!tree.is_valid(node.parent)) return false;
-  auto &p = tree[node.parent];
+  auto& p = tree[node.parent];
 
   return p.data.is_problem_node() && p.data.vertex == node.first_child;
 }
@@ -43,7 +43,7 @@ static void pull_forward(CompTree& tree, VertexID v) {
 void process_neighbors(             //
     ds::graph::Graph const& graph,  //
     CompTree& tree,                 //
-    VI alpha_list[],                //
+    VVV& alpha_list,                //
     bool const visited[],           //
     VertexID pivot,                 //
     int current_prob,               //
@@ -72,7 +72,7 @@ void process_neighbors(             //
  */
 int do_pivot(ds::graph::Graph const& graph,  //
              CompTree& tree,                 //
-             VI alpha_list[],                //
+             VVV& alpha_list,                //
              bool const visited[],           //
              int prob,                       //
              VertexID pivot                  //
diff --git a/src/main/cpp/modular/compute/refinement.cpp b/src/main/cpp/modular/compute/refinement.cpp
index 2da548b..52ab19d 100644
--- a/src/main/cpp/modular/compute/refinement.cpp
+++ b/src/main/cpp/modular/compute/refinement.cpp
@@ -4,6 +4,9 @@ namespace modular {
 namespace compute {
 namespace impl {
 
+// split directions: left or right
+std::vector<SplitDirection> const DIRS = {SplitDirection::LEFT, SplitDirection::RIGHT};
+
 //================================================================================
 //    Set up
 //================================================================================
@@ -44,7 +47,7 @@ static void number_by_tree(CompTree &tree, int prob) {
 //================================================================================
 //    Utilities
 //================================================================================
-static bool is_root_operator(CompTree &tree, int index) {
+static bool is_root_operator(CompTree const &tree, int index) {
   return tree[index].is_root() || !tree[tree[index].parent].data.is_operation_node();
 }
 
@@ -59,12 +62,27 @@ static bool is_root_operator(CompTree &tree, int index) {
  * @param split_type split type
  * @param should_recurse true if it needs to recurse to children
  */
-static void add_split_mark(CompTree &tree, int index, SplitDirection split_type, bool should_recurse) {
-  tree[index].data.set_split_mark(split_type);
+static void add_split_mark(CompTree &tree, int index, SplitDirection split_type, bool should_recurse, util::Profiler *prof) {
+  if (!tree[index].data.is_split_marked(split_type)) {
+    auto p = tree[index].parent;  // must be a valid node
+    // increment the counter if the parent is an operation node
+    if (tree[p].data.is_operation_node()) tree[p].data.increment_num_split_children(split_type);
+    tree[index].data.set_split_mark(split_type);
+  }
+
+  if (!should_recurse || tree[index].data.op_type != Operation::PRIME) return;
 
-  if (should_recurse && tree[index].data.op_type == Operation::PRIME) {
-    for (auto c : tree.get_children(index)) {
-      if (!tree[c].data.is_split_marked(split_type)) tree[c].data.set_split_mark(split_type);
+  // split type is already set to all children
+  if (tree[index].number_of_children() == tree[index].data.get_num_split_children(split_type)) {
+    // PROF(util::pcount(prof, "add_split_mark(): early return"));
+    return;
+  }
+
+  // PROF(util::pcount(prof, "add_split_mark(): proc child"));
+  for (auto c = tree[index].first_child; tree.is_valid(c); c = tree[c].right) {
+    if (!tree[c].data.is_split_marked(split_type)) {
+      tree[index].data.increment_num_split_children(split_type);
+      tree[c].data.set_split_mark(split_type);
     }
   }
 }
@@ -73,16 +91,27 @@ static void add_split_mark(CompTree &tree, int index, SplitDirection split_type,
  * @brief Adds the given mark to all of this node's ancestors.
  * @param split_type mark to be added
  */
-static void mark_ancestors_by_split(CompTree &tree, int index, SplitDirection split_type) {
-  for (auto p : tree.get_ancestors(index)) {
-    if (tree[p].data.is_problem_node()) break;  // passed the operation root
-    add_split_mark(tree, p, split_type, true);
+static void mark_ancestors_by_split(CompTree &tree, int index, SplitDirection split_type, util::Profiler *prof) {
+  for (auto p = tree[index].parent;; p = tree[p].parent) {
+    if (tree[p].data.is_problem_node()) break;
+    if (tree[p].data.is_split_marked(split_type)) {
+      // split type is already set to p but we need to take care of its children
+      add_split_mark(tree, p, split_type, true, prof);
+      break;
+    }
+    add_split_mark(tree, p, split_type, true, prof);
   }
 }
 
 //================================================================================
 //    Get max subtrees
 //================================================================================
+static bool is_parent_fully_charged(CompTree const &tree, int x) {
+  if (is_root_operator(tree, x)) return false;
+  auto p = tree[x].parent;
+  return tree[p].number_of_children() == tree[p].data.number_of_marks();
+}
+
 /**
  * @brief Finds the set of maximal subtrees such that the leaves of
  * each subtree are subsets of the given leaf set.
@@ -91,194 +120,94 @@ static void mark_ancestors_by_split(CompTree &tree, int index, SplitDirection sp
  * @return std::list<NodeP>
  */
 static std::vector<int> get_max_subtrees(CompTree &tree, VI const &leaves) {
-  // charging
-  std::map<int, int> num_charges;
-  std::stack<int> stack;
-  std::unordered_set<int> fully_charged;
-  for (auto v : leaves) {
-    fully_charged.insert(v);
-    stack.push(v);
-  }
+  std::vector<int> full_charged(leaves);
+  std::vector<int> charged;
 
-  while (!stack.empty()) {
-    auto x = stack.top();
-    stack.pop();
-    if (!is_root_operator(tree, x)) {
-      auto p = tree[x].parent;
-      ++num_charges[p];
-
-      if (tree[p].number_of_children() == num_charges[p]) {
-        // fully charged
-        fully_charged.insert(p);
+  // charging
+  for (std::size_t i = 0; i < full_charged.size(); ++i) {
+    auto x = full_charged[i];
+    if (is_root_operator(tree, x)) continue;
 
-        // discharge children
-        for (auto c: tree.get_children(p)) fully_charged.erase(c);
+    auto p = tree[x].parent;
+    if (!tree[p].data.is_marked()) charged.push_back(p);
+    tree[p].data.add_mark();
 
-        stack.push(p);
-      }
+    if (tree[p].data.num_marks == tree[p].number_of_children()) {
+      // fully charged
+      full_charged.push_back(p);
     }
   }
 
-  return std::vector<int>(fully_charged.begin(), fully_charged.end());
+  // discharging
+  std::vector<int> ret;
+  for (auto x : full_charged) {
+    if (!is_parent_fully_charged(tree, x)) ret.push_back(x);
+  }
+  for (auto x : charged) tree[x].data.clear_marks();
+  return ret;
 }
 
-// -----------------------original--------------------------------------------
-// static std::list<int> get_max_subtrees(CompTree &tree, VI const &leaves) {
-//   std::list<int> discharged;
-
-//   // charging process
-//   std::stack<int> stack;
-//   for (auto v : leaves) stack.push(v);
-
-//   while (!stack.empty()) {
-//     auto x = stack.top();
-//     stack.pop();
-//     if (!is_root_operator(tree, x)) {
-//       auto p = tree[x].parent;
-//       tree[p].data.add_mark();
-//       if (tree[p].data.num_marks == tree[p].number_of_children()) stack.push(p);
-//     }
-//     discharged.push_back(x);
-//   }
-
-//   // remove marks on all nodes
-//   for (auto it = discharged.begin(); it != discharged.end(); ++it) {
-//     tree[*it].data.clear_marks();
-//     if (!is_root_operator(tree, *it)) {
-//       auto p = tree[*it].parent;
-//       if (tree[p].data.num_marks == tree[p].number_of_children()) {
-//         it = discharged.erase(it);
-//         --it;
-//       } else {
-//         tree[p].data.clear_marks();
-//       }
-//     }
-//   }
-//   return discharged;
-// }
-// -----------------------original--------------------------------------------
-
 //================================================================================
 //    Group sibling nodes
 //================================================================================
 std::vector<std::pair<int, bool>> group_sibling_nodes(CompTree &tree, std::vector<int> const &nodes) {
-  std::map<int, int> num_marks;
   std::vector<int> parents;
-
-  for (auto node : nodes) {
-    ++num_marks[node];
-
-    if (!is_root_operator(tree, node)) {
-      tree.make_first_child(node);
-      auto p = tree[node].parent;
-
-      if (!util::contains(num_marks, p)) parents.push_back(p);
-      ++num_marks[p];
-    }
-  }
-
   std::vector<std::pair<int, bool>> sibling_groups;
 
-  // (1) the roots of trees
   for (auto node : nodes) {
     if (is_root_operator(tree, node)) {
-      num_marks.erase(node);
+      // (1) the roots of trees
       sibling_groups.push_back({node, false});
+    } else {
+      tree.make_first_child(node);
+      auto p = tree[node].parent;
+
+      if (!tree[p].data.is_marked()) parents.push_back(p);
+      tree[p].data.add_mark();
     }
   }
 
   for (auto p : parents) {
     // there must be at least one mark
     auto c = tree[p].first_child;
+    auto num_marks = tree[p].data.number_of_marks();
 
-    if (num_marks[p] == 1) {
+    if (num_marks == 1) {
       // (2) the non-root nodes without siblings
       sibling_groups.push_back({c, false});
     } else {
       // (3) group sibling nodes as children of a new node inserted in their place.
       auto grouped_children = tree.create_node(tree[p].data);
+      tree[grouped_children].data.clear_marks();
+
+      for (auto st : DIRS) {
+        if (tree[grouped_children].data.is_split_marked(st)) tree[p].data.increment_num_split_children(st);
+      }
 
       auto c = tree[p].first_child;
-      for (int i = 0; i < num_marks[p]; ++i) {
+      for (int i = 0; i < num_marks; ++i) {
         auto nxt = tree[c].right;
         tree.move_to(c, grouped_children);
+
+        for (auto st : DIRS) {
+          if (tree[c].data.is_split_marked(st)) {
+            tree[p].data.decrement_num_split_children(st);
+            tree[grouped_children].data.increment_num_split_children(st);
+          }
+        }
         c = nxt;
       }
       tree.move_to(grouped_children, p);
 
       sibling_groups.push_back({grouped_children, tree[grouped_children].data.op_type == Operation::PRIME});
     }
+    tree[p].data.clear_marks();
   }
 
   // TRACE("return: %s\n", cstr(sibling_groups));
-  // util::pstop(prof, "group_sibling_nodes()");
   return sibling_groups;
 }
 
-// -----------------------original--------------------------------------------
-// std::vector<std::pair<int, bool>> group_sibling_nodes(CompTree &tree, std::list<int> const &nodes) {
-//   // precondition: nodes/parents do not have marks
-//   for (auto node : nodes) {
-//     if (tree[node].data.is_marked()) throw std::invalid_argument("found unclean marks");
-//     if (tree[node].has_parent() && tree[tree[node].parent].data.is_marked()) {
-//       throw std::invalid_argument("found unclean marks");
-//     }
-//   }
-
-//   std::vector<int> parents;
-//   for (auto node : nodes) {
-//     tree[node].data.add_mark();
-//     if (!is_root_operator(tree, node)) {
-//       tree.make_first_child(node);
-//       auto p = tree[node].parent;
-//       if (!tree[p].data.is_marked()) parents.push_back(p);
-//       tree[p].data.add_mark();
-//     }
-//   }
-
-//   std::vector<std::pair<int, bool>> sibling_groups;
-
-//   // (1) the roots of trees
-//   for (auto node : nodes) {
-//     if (is_root_operator(tree, node)) {
-//       tree[node].data.clear_marks();
-//       sibling_groups.push_back({node, false});
-//     }
-//   }
-
-//   for (auto p : parents) {
-//     // there must be at least one mark
-//     auto c = tree[p].first_child;
-
-//     if (tree[p].data.num_marks == 1) {
-//       // (2) the non-root nodes without siblings
-//       tree[p].data.clear_marks();
-//       tree[c].data.clear_marks();
-//       sibling_groups.push_back({c, false});
-//     } else {
-//       // (3) group sibling nodes as children of a new node inserted in their place.
-//       int nm = tree[p].data.num_marks;
-//       tree[p].data.clear_marks();
-//       auto grouped_children = tree.create_node(tree[p].data);
-
-//       auto c = tree[p].first_child;
-//       for (int i = 1; i < nm; ++i) {
-//         auto nxt = tree[c].right;
-//         tree.move_to(c, grouped_children);
-//         c = nxt;
-//       }
-//       tree.move_to(grouped_children, p);
-
-//       sibling_groups.push_back({grouped_children, tree[grouped_children].data.op_type == Operation::PRIME});
-//     }
-//   }
-
-//   // TRACE("return: %s\n", cstr(sibling_groups));
-//   // util::pstop(prof, "group_sibling_nodes()");
-//   return sibling_groups;
-// }
-// -----------------------original--------------------------------------------
-
 //================================================================================
 //    Subroutine
 //================================================================================
@@ -289,7 +218,7 @@ static SplitDirection get_split_type(CompTree &tree, int index, VertexID refiner
   return current < pivot_tn || refiner_tn < current ? SplitDirection::LEFT : SplitDirection::RIGHT;
 }
 
-static void refine_one_node(CompTree &tree, int index, SplitDirection split_type, bool new_prime) {
+static void refine_one_node(CompTree &tree, int index, SplitDirection split_type, bool new_prime, util::Profiler *prof) {
   TRACE("refining tree=%s, index=%d, split_type=%d, new_prime=%d", tree.to_string(index).c_str(), index, split_type, new_prime);
   if (is_root_operator(tree, index)) return;
 
@@ -297,12 +226,18 @@ static void refine_one_node(CompTree &tree, int index, SplitDirection split_type
   int new_sibling = -1;
 
   if (is_root_operator(tree, p)) {
+    // PROF(util::pstart(prof, "refine_one_node: root", 0));
     // parent is a root; must split there
     if (split_type == SplitDirection::LEFT) {
       tree.move_to_before(index, p);
     } else {
       tree.move_to_after(index, p);
     }
+
+    for (auto st : DIRS) {
+      if (tree[index].data.is_split_marked(st)) tree[p].data.decrement_num_split_children(st);
+    }
+
     new_sibling = p;
 
     if (tree[p].has_only_one_child()) {
@@ -310,52 +245,91 @@ static void refine_one_node(CompTree &tree, int index, SplitDirection split_type
       tree.remove(p);
       new_sibling = -1;
     }
+    // PROF(util::pstop(prof, "refine_one_node: root", 0));
   } else if (tree[p].data.op_type != Operation::PRIME) {
+    // PROF(util::pstart(prof, "refine_one_node: non-root", 0));
     // parent is not a root or PRIME
     auto replacement = tree.create_node(tree[p].data);
     tree.replace(p, replacement);
     tree.move_to(index, replacement);
     tree.move_to(p, replacement);
     new_sibling = p;
+
+    for (auto st : DIRS) {
+      if (tree[index].data.is_split_marked(st)) {
+        tree[p].data.decrement_num_split_children(st);
+        tree[replacement].data.increment_num_split_children(st);
+      }
+      if (tree[p].data.is_split_marked(st)) tree[replacement].data.increment_num_split_children(st);
+    }
+    // PROF(util::pstop(prof, "refine_one_node: non-root", 0));
   }
 
-  add_split_mark(tree, index, split_type, new_prime);
-  mark_ancestors_by_split(tree, index, split_type);
+  // PROF(util::pstart(prof, "add_split_mark()", 1));
+  add_split_mark(tree, index, split_type, new_prime, prof);
+  // PROF(util::pstop(prof, "add_split_mark()", 1));
+
+  // PROF(util::pstart(prof, "mark_ancestors_by_split()"));
+  mark_ancestors_by_split(tree, index, split_type, prof);
+  // PROF(util::pstop(prof, "mark_ancestors_by_split()"));
 
   if (new_sibling >= 0) {
     // non-prime or a new root; safe to set should_recurse=true
-    add_split_mark(tree, new_sibling, split_type, true);
+    // PROF(util::pstart(prof, "add_split_mark()", 2));
+    add_split_mark(tree, new_sibling, split_type, true, prof);
+    // PROF(util::pstop(prof, "add_split_mark()", 2));
   }
 }
 
-static void refine_with(CompTree &tree, VI const alpha_list[], VertexID refiner, VertexID pivot) {
+static void refine_with(CompTree &tree, VVV const &alpha_list, VertexID refiner, VertexID pivot, util::Profiler *prof) {
+  // PROF(util::pstart(prof, "get_max_subtrees()"))
   auto subtree_roots = get_max_subtrees(tree, alpha_list[refiner]);
+  // PROF(util::pstop(prof, "get_max_subtrees()"))
+
+  // PROF(util::pstart(prof, "group_sibling_nodes()"))
   auto sibling_groups = group_sibling_nodes(tree, subtree_roots);
+  // PROF(util::pstop(prof, "group_sibling_nodes()"))
 
   TRACE("alpha[%d]: %s", refiner, util::to_string(alpha_list[refiner]).c_str());
   TRACE("subtree_roots: %s", util::to_string(subtree_roots).c_str());
   TRACE("sibling_groups: %s, tree=%s", util::to_string(sibling_groups).c_str(), tree.to_string(tree[pivot].parent).c_str());
 
+  // PROF(util::pstart(prof, "refine_with: main loop"))
   for (auto &x : sibling_groups) {
+    // PROF(util::pstart(prof, "get_split_type()"))
     auto split_type = get_split_type(tree, x.first, refiner, pivot);
-    refine_one_node(tree, x.first, split_type, x.second);
+    // PROF(util::pstop(prof, "get_split_type()"))
+
+    // PROF(util::pstart(prof, "refine_one_node()"))
+    refine_one_node(tree, x.first, split_type, x.second, prof);
+    // PROF(util::pstop(prof, "refine_one_node()"))
   }
+  // PROF(util::pstop(prof, "refine_with: main loop"))
 }
 
 //================================================================================
 //    Refinement
 //================================================================================
-void refine(CompTree &tree, VI const alpha_list[], int prob) {
+void refine(CompTree &tree, VVV const &alpha_list, int prob, std::vector<int> &leaves, util::Profiler *prof) {
   TRACE("start: %s", tree.to_string(prob).c_str());
+  // PROF(util::pstart(prof, "refinement:refine()"))
 
+  // PROF(util::pstart(prof, "refinement:number_by_comp()"))
   number_by_comp(tree, prob);
+  // PROF(util::pstop(prof, "refinement:number_by_comp()"))
+
+  // PROF(util::pstart(prof, "refinement:number_by_tree()"))
   number_by_tree(tree, prob);
+  // PROF(util::pstop(prof, "refinement:number_by_tree()"))
 
-  for (auto v : tree.get_leaves(prob)) {
-    refine_with(tree, alpha_list, v, tree[prob].data.vertex);
+  // PROF(util::pstart(prof, "refinement:refine_with()"));
+  for (auto v : leaves) {
+    refine_with(tree, alpha_list, v, tree[prob].data.vertex, prof);
     TRACE("refined at %d: tree=%s", v, tree.to_string(prob).c_str())
   }
+  // PROF(util::pstop(prof, "refinement:refine_with()"));
 
+  // PROF(util::pstop(prof, "refinement:refine()"))
   TRACE("finish: %s", tree.to_string(prob).c_str());
 }
 
diff --git a/src/main/cpp/util/profiler.hpp b/src/main/cpp/util/profiler.hpp
index 7a76058..2ffae35 100644
--- a/src/main/cpp/util/profiler.hpp
+++ b/src/main/cpp/util/profiler.hpp
@@ -110,11 +110,13 @@ class Profiler {
         auto count = kv.second.frequency;
         double ms = std::chrono::duration_cast<std::chrono::nanoseconds>(kv.second.accumulated_time).count();
         double bs = std::chrono::duration_cast<std::chrono::nanoseconds>(kv.second.best_time).count();
-        auto percall = std::string(ms / count < 1e4 ? format("%10.0f ns /call", ms / count)
-                                                    : format("%10.6f sec/call", ms / count * 1e-9));
-        auto mincall = std::string(format("%8.6f s", bs * 1e-9));
-        fprintf(stderr, "%-30s (%5d): %10.3f sec %9d calls [%s; (min)%s]\n", kv.first.first.c_str(), kv.first.second,
-                ms * 1e-9, count, percall.c_str(), mincall.c_str());
+
+        fprintf(stderr, "%-30s (%5d): %10.3f sec %9d calls ", kv.first.first.c_str(), kv.first.second, ms * 1e-9, count);
+        if (ms / count < 1e4) {
+          fprintf(stderr, "[%10.0f ns /call; (min)%10.6f s]\n", ms / count, bs * 1e-9);
+        } else {
+          fprintf(stderr, "[%10.6f sec/call; (min)%10.6f s]\n", ms / count * 1e-9, bs * 1e-9);
+        }
       }
     }
 
diff --git a/src/main/python/modular/MDTree.py b/src/main/python/modular/MDTree.py
index 5783998..4b230e0 100644
--- a/src/main/python/modular/MDTree.py
+++ b/src/main/python/modular/MDTree.py
@@ -16,13 +16,13 @@ class MDTree:
     Modular decomposition tree for undirected graphs.
     """
 
-    def __init__(self, G: nx.Graph, solver: str = 'linear') -> None:
+    def __init__(self, G: nx.Graph, solver: str = 'linear', verify: bool = False) -> None:
         assert len(G) > 0, 'graph cannot be empty'
 
         if solver == 'linear':
-            tree, root, vertices = MDSolver.compute(G)
+            tree, root, vertices = MDSolver.compute(G, verify=verify)
         elif solver == 'naive':
-            tree, root, vertices = MDNaiveSolver.compute(G)
+            tree, root, vertices = MDNaiveSolver.compute(G, verify=verify)
         else:
             raise ValueError(f'Unknown solver: {solver}')
 
@@ -40,7 +40,7 @@ def sort(self) -> None:
         for node in reversed(level_order):
             if node.is_leaf():
                 min_label[node] = node.data.vertex
-            
+
             if node.parent is not None:
                 min_label[node.parent] = min(min_label[node.parent], min_label[node]) if node.parent in min_label else min_label[node]
 
@@ -73,12 +73,12 @@ def __repr__(self) -> str:
         return repr(self.root)
 
 
-def modular_decomposition(G: nx.Graph, sorted: bool = False, solver: str = 'linear') -> Optional[MDTree]:
+def modular_decomposition(G: nx.Graph, sorted: bool = False, solver: str = 'linear', verify: bool = False) -> Optional[MDTree]:
     """Alias to the constructor."""
 
     if len(G) == 0:
         return None
-    ret = MDTree(G, solver=solver)
+    ret = MDTree(G, solver=solver, verify=verify)
     if sorted:
         ret.sort()
     return ret
diff --git a/src/main/python/modular/compute/MDComputeNode.py b/src/main/python/modular/compute/MDComputeNode.py
index a3cb5aa..8f7e5dd 100644
--- a/src/main/python/modular/compute/MDComputeNode.py
+++ b/src/main/python/modular/compute/MDComputeNode.py
@@ -34,6 +34,9 @@ def __init__(
         self.vertex = vertex
         self.comp_number = comp_number
         self.tree_number = tree_number
+        self.num_marks = 0  # used only in refinement
+        self.num_left_split_children = 0  # used only in refinement
+        self.num_right_split_children = 0  # used only in refinement
         self.active = active
         self.connected = connected
 
@@ -62,6 +65,12 @@ def is_vertex_node(self) -> bool: return self.node_type == NodeType.VERTEX
     def is_problem_node(self) -> bool: return self.node_type == NodeType.PROBLEM
     def is_operation_node(self) -> bool: return self.node_type == NodeType.OPERATION
 
+    def add_mark(self) -> None:
+        self.num_marks += 1
+
+    def clear_marks(self) -> None:
+        self.num_marks = 0
+
     def is_split_marked(self, split_type: SplitDirection):
         return self.split_type in [split_type, SplitDirection.MIXED]
 
@@ -75,10 +84,36 @@ def set_split_mark(self, split_type: SplitDirection):
         else:
             self.split_type = SplitDirection.MIXED
 
+    def increment_num_split_children(self, split_type: SplitDirection) -> None:
+        assert split_type in [SplitDirection.LEFT, SplitDirection.RIGHT]
+
+        if split_type == SplitDirection.LEFT:
+            self.num_left_split_children += 1
+        else:
+            self.num_right_split_children += 1
+
+    def decrement_num_split_children(self, split_type: SplitDirection) -> None:
+        assert split_type in [SplitDirection.LEFT, SplitDirection.RIGHT]
+
+        if split_type == SplitDirection.LEFT:
+            self.num_left_split_children -= 1
+        else:
+            self.num_right_split_children -= 1
+
+    def get_num_split_children(self, split_type: SplitDirection) -> int:
+        assert split_type in [SplitDirection.LEFT, SplitDirection.RIGHT]
+
+        return self.num_left_split_children if split_type == SplitDirection.LEFT else self.num_right_split_children
+
+    def clear_num_split_children(self) -> None:
+        self.num_left_split_children = 0
+        self.num_right_split_children = 0
+
     def clear(self) -> None:
         self.comp_number = -1
         self.tree_number = -1
         self.split_type = SplitDirection.NONE
+        self.clear_num_split_children()
 
     def __str__(self) -> str:
         if self.node_type == NodeType.VERTEX:
diff --git a/src/main/python/modular/compute/MDSolver.py b/src/main/python/modular/compute/MDSolver.py
index 747203c..1fcde54 100644
--- a/src/main/python/modular/compute/MDSolver.py
+++ b/src/main/python/modular/compute/MDSolver.py
@@ -69,8 +69,6 @@ def compute(
             # trace(f'after promote: current={current_prob}')
             assemble(tree, vertex_nodes, alpha_list, current_prob)
             # trace(f'after assemble: current={current_prob}')
-            merge_components(tree, current_prob, extra_components)
-            # trace(f'after merge_components: current={current_prob}')
 
             # clear all but visited
             assert current_prob.first_child is not None
@@ -79,6 +77,9 @@ def compute(
                     del alpha_list[c.data.vertex]
                 c.data.clear()
 
+            merge_components(tree, current_prob, extra_components)
+            # trace(f'after merge_components: current={current_prob}')
+
         # trace(f'solve finish: current ({t}): {current_prob}, tree={tree}, visited={visited}, alpha={dict(alpha_list)}')
         result = current_prob.first_child
         current_prob = current_prob.parent if current_prob.is_last_child() else current_prob.right
@@ -97,7 +98,7 @@ def compute(
 
 class MDSolver:
     @staticmethod
-    def compute(G: nx.Graph) -> tuple[RootedForest[MDNode], Node[MDNode], list[VertexId]]:
+    def compute(G: nx.Graph, verify: bool = False) -> tuple[RootedForest[MDNode], Node[MDNode], list[VertexId]]:
         n = len(G)
         assert n > 0, 'empty graph'
 
@@ -107,7 +108,7 @@ def compute(G: nx.Graph) -> tuple[RootedForest[MDNode], Node[MDNode], list[Verte
         labels = nx.get_node_attributes(G, 'label')
 
         # build computation tree
-        tree: RootedForest[MDComputeNode] = RootedForest()
+        tree: RootedForest[MDComputeNode] = RootedForest(verify=verify)
 
         # create the main problem
         main_prob = tree.create_node(MDComputeNode.new_problem_node(connected=False))
diff --git a/src/main/python/modular/compute/assembly.py b/src/main/python/modular/compute/assembly.py
index 45f926e..db9672c 100644
--- a/src/main/python/modular/compute/assembly.py
+++ b/src/main/python/modular/compute/assembly.py
@@ -30,7 +30,7 @@ def assemble(
     lcocomp = determine_left_cocomp_fragments(ps, pivot_index)
     rcomp = determine_right_comp_fragments(ps, pivot_index)
     rlayer = determine_right_layer_neighbor(vertex_nodes, alpha_list, ps, pivot_index)
-    neighbors = compute_fact_perm_edges(vertex_nodes, alpha_list, ps)
+    neighbors = compute_fact_perm_edges(vertex_nodes, alpha_list, ps, pivot_index)
     mu = compute_mu(ps, pivot_index, neighbors)
     boundaries = delineate(pivot_index, lcocomp, rcomp, rlayer, mu)
     # print(f'[TRACE] boundaries={boundaries}')
@@ -97,7 +97,8 @@ def f(t: Node[MDComputeNode]) -> bool:
 def compute_fact_perm_edges(
     vertex_nodes: dict[VertexId, Node[MDComputeNode]],
     alpha_list: dict[VertexId, set[VertexId]],
-    ps: list[Node[MDComputeNode]]
+    ps: list[Node[MDComputeNode]],
+    pivot_index: int
 ) -> list[list[int]]:
     k = len(ps)
     neighbors: list[list[int]] = [[] for _ in range(k)]
@@ -110,14 +111,15 @@ def compute_fact_perm_edges(
             elem_size[i] += 1  # increment the element size
 
     # find joins
-    for i, p in enumerate(ps):
-        candidates = []
-        marks = [0] * k
+    marks = [0] * k
+    for i in range(pivot_index):  # we need the neighbors only up to pivot_index
+        p = ps[i]
+        candidates = set()
 
         for leaf in p.get_leaves():
             for a in alpha_list[leaf.data.vertex]:
                 j = vertex_nodes[a].data.comp_number
-                candidates += [j]
+                candidates.add(j)
                 marks[j] += 1
 
         for j in candidates:
@@ -125,7 +127,8 @@ def compute_fact_perm_edges(
                 # found a join between i and j
                 neighbors[i] += [j]
                 marks[j] = 0  # reset marks so that there will be no duplicates
-
+            else:
+                assert False, 'never happens?'
     return neighbors
 
 
diff --git a/src/main/python/modular/compute/refinement.py b/src/main/python/modular/compute/refinement.py
index e4926e9..81025e0 100644
--- a/src/main/python/modular/compute/refinement.py
+++ b/src/main/python/modular/compute/refinement.py
@@ -75,54 +75,77 @@ def is_root_operator(node: Node[MDComputeNode]) -> bool:
 # ===============================================================================
 def mark_ancestors_by_split(node: Node[MDComputeNode], split_type: SplitDirection) -> None:
     """Adds the given mark to all of the node's ancestors."""
-    for p in node.get_ancestors():
+    p = node.parent
+    while p is not None:
         if p.data.is_problem_node():
             break  # passed the operation root
+        if p.data.is_split_marked(split_type):
+            add_split_mark(p, split_type, should_recurse=True)
+            break
         add_split_mark(p, split_type, should_recurse=True)
+        p = p.parent
 
 
 def add_split_mark(node: Node[MDComputeNode], split_type: SplitDirection, should_recurse: bool) -> None:
     """Add the given mark to the node and possibly its children."""
-    node.data.set_split_mark(split_type)
+    if not node.data.is_split_marked(split_type):
+        if node.parent is not None and node.parent.data.is_operation_node():
+            node.parent.data.increment_num_split_children(split_type)
+        node.data.set_split_mark(split_type)
+
+    if not should_recurse or node.data.op_type != OperationType.PRIME:
+        return
+
+    if node.number_of_children() == node.data.get_num_split_children(split_type):
+        return  # split type is already set to all children
 
-    if should_recurse and node.data.op_type == OperationType.PRIME:
-        for c in node.get_children():
-            if not c.data.is_split_marked(split_type):
-                c.data.set_split_mark(split_type)
+    for c in node.get_children():
+        if not c.data.is_split_marked(split_type):
+            node.data.increment_num_split_children(split_type)
+            c.data.set_split_mark(split_type)
 
 
 # ===============================================================================
 #    Get max subtrees
 # ===============================================================================
+def is_parent_fully_charged(x: Node[MDComputeNode]) -> bool:
+    if is_root_operator(x):
+        return False
+    p = x.parent
+    assert p is not None
+    return p.number_of_children() == p.data.num_marks
+
+
 def get_max_subtrees(leaves: list[Node[MDComputeNode]]) -> list[Node[MDComputeNode]]:
     """
     Finds the set of maximal subtrees such that the leaves of each subtree are
     subsets of the given leaf set.
     """
-
-    num_charges: dict[Node[MDComputeNode], int] = defaultdict(int)
+    full_charged = list(leaves)
+    charged: list[Node[MDComputeNode]] = []
 
     # charging
-    fully_charged: set[Node[MDComputeNode]] = set(leaves)
-    st = [x for x in leaves]
-
-    while st:
-        x = st.pop()
+    idx = 0
+    while idx < len(full_charged):
+        x = full_charged[idx]
         if not is_root_operator(x):
-            p = x.parent
-            assert p is not None
-            num_charges[p] += 1
-            if p.number_of_children() == num_charges[p]:
-                # fully charged
-                fully_charged.add(p)
+            p = x.get_parent()
+            if p.data.num_marks == 0:
+                charged += [p]
+            p.data.add_mark()
 
-                # discharge children
-                for c in p.get_children():
-                    fully_charged.remove(c)
+            if p.number_of_children() == p.data.num_marks:
+                # fully charged
+                full_charged += [p]
+        idx += 1
 
-                st += [p]
+    # discharging
+    ret: list[Node[MDComputeNode]] = [x for x in full_charged if not is_parent_fully_charged(x)]
 
-    return list(fully_charged)
+    # cleaning
+    for x in charged:
+        x.data.clear_marks()
+    return ret
 
 
 # ===============================================================================
@@ -132,48 +155,56 @@ def group_sibling_nodes(
     tree: RootedForest[MDComputeNode],
     nodes: list[Node[MDComputeNode]]
 ) -> list[tuple[Node[MDComputeNode], bool]]:
-    num_marks: dict[Node[MDComputeNode], int] = defaultdict(int)
+    """
+    Precondition: `nodes` must be maximal subtrees (no node is an ancestor of another).
+    """
     parents: list[Node[MDComputeNode]] = []
+    sibling_groups: list[tuple[Node[MDComputeNode], bool]] = []
 
     for node in nodes:
-        num_marks[node] += 1
-
-        if not is_root_operator(node):
+        if is_root_operator(node):
+            # (1) roots of trees
+            sibling_groups += [(node, False)]
+        else:
             tree.make_first_child(node)
             p = node.parent
             assert p is not None
 
-            if p not in num_marks:
+            if p.data.num_marks == 0:
                 parents += [p]
-            num_marks[p] += 1
-
-    sibling_groups: list[tuple[Node[MDComputeNode], bool]] = []
-
-    # (1) roots of trees
-    for node in nodes:
-        if is_root_operator(node):
-            del num_marks[node]
-            sibling_groups += [(node, False)]
+            p.data.add_mark()
 
     for p in parents:
         # there must be at least one mark
         c = p.first_child
         assert c is not None
 
-        if num_marks[p] == 1:
+        if p.data.num_marks == 1:
             # (2) non-root nodes without siblings
             sibling_groups += [(c, False)]
         else:
             # (3) group sibling nodes as the children of a newly inserted node
             grouped_children = tree.create_node(p.data.copy())
 
+            for split_type in [SplitDirection.LEFT, SplitDirection.RIGHT]:
+                if grouped_children.data.is_split_marked(split_type):
+                    p.data.increment_num_split_children(split_type)
+
             # FIX FROM ORIGINAL CODE (RecSubProblem.java line 971); bug when a child of `p` is also in `parents`
-            for c in p.get_children()[:num_marks[p]]:
+            for c in p.get_children()[:p.data.num_marks]:
+                for split_type in [SplitDirection.LEFT, SplitDirection.RIGHT]:
+                    if c.data.is_split_marked(split_type):
+                        p.data.decrement_num_split_children(split_type)
+                        grouped_children.data.increment_num_split_children(split_type)
+
                 tree.move_to(c, grouped_children)
             tree.move_to(grouped_children, p)
 
             sibling_groups += [(grouped_children, grouped_children.data.op_type == OperationType.PRIME)]
 
+        # clean marks
+        p.data.clear_marks()
+
     return sibling_groups
 
 
@@ -212,6 +243,10 @@ def refine_one_node(
             tree.move_to_before(node, p)
         else:
             tree.move_to_after(node, p)
+        for st in [SplitDirection.LEFT, SplitDirection.RIGHT]:
+            if node.data.is_split_marked(st):
+                p.data.decrement_num_split_children(st)
+
         new_sibling = p
 
         if p.has_only_one_child():
@@ -232,6 +267,13 @@ def refine_one_node(
         tree.move_to(p, replacement)
         new_sibling = p
 
+        for st in [SplitDirection.LEFT, SplitDirection.RIGHT]:
+            if node.data.is_split_marked(st):
+                p.data.decrement_num_split_children(st)
+                replacement.data.increment_num_split_children(st)
+            if p.data.is_split_marked(st):
+                replacement.data.increment_num_split_children(st)
+
     add_split_mark(node, split_type, should_recurse=new_prime)
     mark_ancestors_by_split(node, split_type)
 
@@ -257,3 +299,17 @@ def refine_with(
         # trace(f'before: refiner={refiner}, node={current}, tree={tree}')
         refine_one_node(tree, current, split_type, new_prime)
         # trace(f'after : refiner={refiner}, node={current}, tree={tree}')
+
+        # sanity check
+        # for r in tree.roots:
+        #     for node in r.dfs_reverse_preorder_nodes():
+        #         if not node.data.is_operation_node():
+        #             assert node.data.num_left_split_children == 0
+        #             assert node.data.num_right_split_children == 0
+        #         else:
+        #             lact = node.data.num_left_split_children
+        #             ract = node.data.num_right_split_children
+        #             lexp = sum(1 if c.data.is_split_marked(SplitDirection.LEFT) else 0 for c in node.get_children())
+        #             rexp = sum(1 if c.data.is_split_marked(SplitDirection.RIGHT) else 0 for c in node.get_children())
+        #             assert lact == lexp, f'node={node}, expect={lexp}, actual={lact}'
+        #             assert ract == rexp, f'node={node}, expect={rexp}, actual={lact}'
diff --git a/src/main/python/modular/compute_naive/MDNaiveSolver.py b/src/main/python/modular/compute_naive/MDNaiveSolver.py
index d34f6f8..afd6abf 100644
--- a/src/main/python/modular/compute_naive/MDNaiveSolver.py
+++ b/src/main/python/modular/compute_naive/MDNaiveSolver.py
@@ -59,7 +59,7 @@ def compute_implication_matrix(G: nx.Graph) -> np.ndarray:
         return M + M.T
 
     @staticmethod
-    def compute(G: nx.Graph) -> tuple[RootedForest[MDNode], Node[MDNode], list[VertexId]]:
+    def compute(G: nx.Graph, verify: bool = False) -> tuple[RootedForest[MDNode], Node[MDNode], list[VertexId]]:
         # Implements the algorithm described in
         # "A Fast Algorithm for the Decomposition of Graphs and Posets"
         # Hermann Buer and Rolf H. Möhring (1983)
@@ -73,7 +73,7 @@ def compute(G: nx.Graph) -> tuple[RootedForest[MDNode], Node[MDNode], list[Verte
         n = len(G)
         vertices = list(G.nodes())
 
-        tree: RootedForest[MDNode] = RootedForest()
+        tree: RootedForest[MDNode] = RootedForest(verify=verify)
         root = tree.create_node(MDNode(None, None, 0, n))
 
         q: deque[Node[MDNode]] = deque()
diff --git a/src/main/python/modular/tree/RootedForest.py b/src/main/python/modular/tree/RootedForest.py
index 9ede8d9..162f764 100644
--- a/src/main/python/modular/tree/RootedForest.py
+++ b/src/main/python/modular/tree/RootedForest.py
@@ -45,6 +45,11 @@ def has_child(self) -> bool: return self.first_child is not None
     def has_only_one_child(self) -> bool: return self.num_children == 1
     def number_of_children(self) -> int: return self.num_children
 
+    def get_parent(self) -> Node[T]:
+        ret = self.parent
+        assert ret is not None
+        return ret
+
     def get_children(self) -> list[Node[T]]:
         """Returns a list of the (direct) children of this node."""
         x = self.first_child
@@ -155,9 +160,14 @@ class RootedForest(Generic[T]):
     Generic rooted forest (disjoint set of rooted trees).
     """
 
-    def __init__(self) -> None:
+    def __init__(self, verify: bool = False) -> None:
+        """
+        Args:
+            verify: True if verification is enabled
+        """
         self.roots: set[Node[T]] = set()
         self.size = 0
+        self.verify = verify
 
     def __len__(self) -> int: return self.size
     def __bool__(self) -> bool: return bool(self.roots)
@@ -213,7 +223,8 @@ def move_to(self, node: Node[T], new_parent: Node[T]) -> None:
         self.roots.remove(node)
 
     def swap(self, a: Node[T], b: Node[T]) -> None:
-        assert a.get_root() != b.get_root(), 'swap: a and b must belong to different trees'
+        if self.verify:
+            assert a.get_root() != b.get_root(), 'swap: a and b must belong to different trees'
 
         for x, y in [(a, b), (b, a)]:
             if x.is_first_child():
@@ -234,7 +245,9 @@ def swap(self, a: Node[T], b: Node[T]) -> None:
     def replace(self, node: Node[T], replace_by: Node[T]) -> None:
         """Replaces the node and its subtree with the given node."""
         assert node != replace_by, 'replace: replace_by must differ from node'
-        assert replace_by not in node.get_ancestors(), 'replace: replace_by cannot be an ancestor of node'
+
+        if self.verify:
+            assert replace_by not in node.get_ancestors(), 'replace: replace_by cannot be an ancestor of node'
 
         self.detach(replace_by)
         self.swap(node, replace_by)
@@ -243,7 +256,9 @@ def move_to_before(self, node: Node[T], target: Node[T]) -> None:
         """Moves the node to the left sibling of the given target."""
         assert target.parent is not None, 'move_to_before: target must not be a root'
         assert node != target, 'move_to_before: node and target must differ'
-        assert node not in target.get_ancestors(), 'move_to_before: target cannot be an ancestor of node'
+
+        if self.verify:
+            assert node not in target.get_ancestors(), 'move_to_before: target cannot be an ancestor of node'
 
         self.detach(node)
         self.roots.remove(node)
@@ -263,7 +278,9 @@ def move_to_after(self, node: Node[T], target: Node[T]) -> None:
         """Moves the node to the right sibling of the given target."""
         assert target.parent is not None, 'move_to_after: target must not be a root'
         assert node != target, 'move_to_after: node and target must differ'
-        assert node not in target.get_ancestors(), 'move_to_after: target cannot be an ancestor of node'
+
+        if self.verify:
+            assert node not in target.get_ancestors(), 'move_to_after: target cannot be an ancestor of node'
 
         self.detach(node)
         self.roots.remove(node)
@@ -285,7 +302,9 @@ def make_first_child(self, node: Node[T]) -> None:
 
     def add_children_from(self, node: Node[T], target: Node[T]) -> None:
         """Moves all the children of the target to node."""
-        assert target not in node.get_ancestors(), 'add_children_from: target cannot be an ancestor of node'
+
+        if self.verify:
+            assert target not in node.get_ancestors(), 'add_children_from: target cannot be an ancestor of node'
 
         if node == target:
             return  # do nothing
diff --git a/src/test/python/__init__.py b/src/test/python/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/test/python/modular/__init__.py b/src/test/python/modular/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/test/python/modular/compute/test_refinement.py b/src/test/python/modular/compute/test_refinement.py
index b5edaa0..dd1fce3 100644
--- a/src/test/python/modular/compute/test_refinement.py
+++ b/src/test/python/modular/compute/test_refinement.py
@@ -145,3 +145,29 @@ def setup():
         sg = group_sibling_nodes(tree, [op0])[0]
         refine_one_node(tree, sg[0], SplitDirection.LEFT, False)
         self.assertEqual(repr(prob), '(C-(P-(P-(J-(P-(P-(0-)(1-)(2-))(3-))(4-))(5-))(6-)(7-)(8-)))')
+
+        def setup2():
+            n = 16
+            tree: RootedForest[MDComputeNode] = RootedForest()
+            vs = list(reversed([tree.create_node(MDComputeNode.new_vertex_node(n - 1 - i)) for i in range(n)]))
+            prob = tree.create_node(MDComputeNode.new_problem_node(False))
+            ops = [tree.create_node(MDComputeNode.new_operation_node(OperationType.PRIME)) for _ in range(n - 1)]
+
+            tree.move_to(ops[0], prob)
+            for i in range(7):
+                tree.move_to(ops[i * 2 + 2], ops[i])
+                tree.move_to(ops[i * 2 + 1], ops[i])
+            for i in range(8):
+                tree.move_to(vs[i * 2 + 1], ops[7 + i])
+                tree.move_to(vs[i * 2], ops[7 + i])
+            return tree, prob, vs, ops
+
+        tree, prob, vs, ops = setup2()
+        refine_one_node(tree, vs[0], SplitDirection.RIGHT, False)
+        self.assertEqual(repr(prob), '(C-(P>(P>(P>(P>(0>)(1>))(P>(2-)(3-)))(P>(P-(4-)(5-))(P-(6-)(7-))))(P>(P-(P-(8-)(9-))(P-(10-)(11-)))(P-(P-(12-)(13-))(P-(14-)(15-))))))')
+
+        refine_one_node(tree, vs[15], SplitDirection.RIGHT, False)
+        self.assertEqual(repr(prob), '(C-(P>(P>(P>(P>(0>)(1>))(P>(2-)(3-)))(P>(P-(4-)(5-))(P-(6-)(7-))))(P>(P>(P-(8-)(9-))(P-(10-)(11-)))(P>(P>(12-)(13-))(P>(14>)(15>))))))')
+
+        refine_one_node(tree, vs[5], SplitDirection.RIGHT, False)
+        self.assertEqual(repr(prob), '(C-(P>(P>(P>(P>(0>)(1>))(P>(2-)(3-)))(P>(P>(4>)(5>))(P>(6-)(7-))))(P>(P>(P-(8-)(9-))(P-(10-)(11-)))(P>(P>(12-)(13-))(P>(14>)(15>))))))')
diff --git a/src/test/python/modular/test_MDTree.py b/src/test/python/modular/test_MDTree.py
index b8cc86a..fab2a64 100644
--- a/src/test/python/modular/test_MDTree.py
+++ b/src/test/python/modular/test_MDTree.py
@@ -11,15 +11,15 @@ class TestMDTree(unittest.TestCase):
     def check_property(self, G: nx.Graph):
         n = len(G)
 
-        t_naive = modular_decomposition(G, sorted=True, solver='naive')
-        t_linear = modular_decomposition(G, sorted=True, solver='linear')
+        t_naive = modular_decomposition(G, sorted=True, solver='naive', verify=True)
+        t_linear = modular_decomposition(G, sorted=True, solver='linear', verify=True)
 
         self.assertEqual(str(t_naive), str(t_linear), msg=f'n={len(G)}, edges={G.edges()}')
         self.assertEqual(t_naive.modular_width(), t_linear.modular_width(), msg=f'n={len(G)}, edges={G.edges()}')
 
         # complement graph
-        c_naive = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
-        c_linear = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        c_naive = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
+        c_linear = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
 
         self.assertEqual(str(c_naive), str(c_linear), msg=f'n={len(G)}, edges={G.edges()}')
         self.assertEqual(str(c_naive), str(t_naive).replace('J', 'X').replace('U', 'J').replace('X', 'U'), msg=f'n={len(G)}, edges={G.edges()}')
@@ -30,8 +30,8 @@ def check_property(self, G: nx.Graph):
         # vertex permutation
         mapping = {i: n - 1 - i for i in range(n)}
         H = nx.relabel_nodes(G, mapping, copy=True)
-        p_naive = modular_decomposition(H, sorted=True, solver='naive')
-        p_linear = modular_decomposition(H, sorted=True, solver='linear')
+        p_naive = modular_decomposition(H, sorted=True, solver='naive', verify=True)
+        p_linear = modular_decomposition(H, sorted=True, solver='linear', verify=True)
 
         self.assertEqual(t_naive.modular_width(), p_naive.modular_width(), msg=f'n={len(G)}, edges={G.edges()}')
         self.assertEqual(t_naive.modular_width(), p_linear.modular_width(), msg=f'n={len(G)}, edges={G.edges()}')
@@ -168,20 +168,20 @@ def test_modular_decomposition(self):
         G = nx.Graph()
         G.add_nodes_from(['a', 'b', 'c', 'd'])
         G.add_edges_from([('a', 'b'), ('b', 'c'), ('c', 'a'), ('b', 'd')])
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 0)
         self.assertEqual(str(t), "(J(U(J('a')('c'))('d'))('b'))")
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 0)
         self.assertEqual(str(t), "(J(U(J('a')('c'))('d'))('b'))")
 
         G = nx.hypercube_graph(3)
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 8)
         self.assertEqual(str(t), '(P((0, 0, 0))((0, 0, 1))((0, 1, 0))((0, 1, 1))((1, 0, 0))((1, 0, 1))((1, 1, 0))((1, 1, 1)))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 8)
         self.assertEqual(str(t), '(P((0, 0, 0))((0, 0, 1))((0, 1, 0))((0, 1, 1))((1, 0, 0))((1, 0, 1))((1, 1, 0))((1, 1, 1)))')
 
@@ -189,19 +189,19 @@ def test_modular_decomposition_2(self):
         G = nx.empty_graph(6)
         G.add_edges_from([(0, 2), (2, 4), (4, 3)])
 
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(U(P(0)(2)(3)(4))(1)(5))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(J(P(0)(2)(3)(4))(1)(5))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(U(P(0)(2)(3)(4))(1)(5))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(J(P(0)(2)(3)(4))(1)(5))')
 
@@ -209,19 +209,19 @@ def test_modular_decomposition_3(self):
         G = nx.empty_graph(7)
         G.add_edges_from([(5, 2), (5, 0), (5, 6), (1, 3)])
 
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 0)
         self.assertEqual(str(t), '(U(J(U(0)(2)(6))(5))(J(1)(3))(4))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 0)
         self.assertEqual(str(t), '(J(U(J(0)(2)(6))(5))(U(1)(3))(4))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 0)
         self.assertEqual(str(t), '(U(J(U(0)(2)(6))(5))(J(1)(3))(4))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 0)
         self.assertEqual(str(t), '(J(U(J(0)(2)(6))(5))(U(1)(3))(4))')
 
@@ -229,19 +229,19 @@ def test_modular_decomposition_4(self):
         G = nx.empty_graph(13)
         G.add_edges_from([(1, 9), (2, 8), (2, 11), (3, 5), (3, 7), (4, 11), (5, 9), (6, 12), (7, 10)])
 
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(U(0)(P(1)(3)(5)(7)(9)(10))(P(2)(4)(8)(11))(J(6)(12)))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(J(0)(P(1)(3)(5)(7)(9)(10))(P(2)(4)(8)(11))(U(6)(12)))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(U(0)(P(1)(3)(5)(7)(9)(10))(P(2)(4)(8)(11))(J(6)(12)))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(J(0)(P(1)(3)(5)(7)(9)(10))(P(2)(4)(8)(11))(U(6)(12)))')
 
@@ -249,23 +249,23 @@ def test_modular_decomposition_5(self):
         G = nx.empty_graph(8)
         G.add_edges_from([(0, 1), (1, 2), (2, 3), (4, 5), (5, 6), (6, 7)])
 
-        t1 = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
-        t2 = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t1 = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
+        t2 = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t1), str(t2))
 
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(U(P(0)(1)(2)(3))(P(4)(5)(6)(7)))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(J(P(0)(1)(2)(3))(P(4)(5)(6)(7)))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(U(P(0)(1)(2)(3))(P(4)(5)(6)(7)))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 4)
         self.assertEqual(str(t), '(J(P(0)(1)(2)(3))(P(4)(5)(6)(7)))')
 
@@ -273,19 +273,19 @@ def test_modular_decomposition_6(self):
         G = nx.empty_graph(6)
         G.add_edges_from([(0, 1), (1, 2), (2, 3), (2, 4), (4, 5)])
 
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
 
@@ -295,19 +295,19 @@ def test_modular_decomposition_7(self):
 
         # self.check_property(G)
 
-        t = modular_decomposition(G, sorted=True, solver='naive')
+        t = modular_decomposition(G, sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='naive', verify=True)
         self.assertEqual(t.modular_width(), 6)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
         self.assertEqual(t.modular_width(), 6)
 
-        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear')
+        t = modular_decomposition(nx.complement(G), sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5))')
         self.assertEqual(t.modular_width(), 6)
 
@@ -318,7 +318,7 @@ def test_modular_decomposition_8(self):
             (6, 10), (6, 13), (7, 8), (7, 11), (9, 13), (11, 13)
         ])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13))')
         self.assertEqual(t.modular_width(), 14)
 
@@ -328,7 +328,7 @@ def test_modular_decomposition_9(self):
         G = nx.empty_graph(9)
         G.add_edges_from([(0, 7), (1, 4), (2, 4), (2, 7), (2, 8), (4, 5)])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(U(P(0)(U(1)(5))(2)(4)(7)(8))(3)(6))')
         self.assertEqual(t.modular_width(), 6)
 
@@ -338,7 +338,7 @@ def test_modular_decomposition_10(self):
         G = nx.empty_graph(12)
         G.add_edges_from([(0, 8), (1, 6), (1, 7), (4, 8), (5, 7), (6, 8), (6, 9), (8, 9), (9, 11)])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(U(P(U(0)(4))(1)(5)(6)(7)(8)(9)(11))(2)(3)(10))')
         self.assertEqual(t.modular_width(), 8)
 
@@ -348,7 +348,7 @@ def test_modular_decomposition_11(self):
         G = nx.empty_graph(11)
         G.add_edges_from([(0, 5), (1, 3), (1, 8), (3, 8), (4, 9), (7, 8), (8, 9)])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(U(J(0)(5))(P(U(J(1)(3))(7))(4)(8)(9))(2)(6)(10))')
         self.assertEqual(t.modular_width(), 4)
 
@@ -361,7 +361,7 @@ def test_modular_decomposition_12(self):
             (4, 7), (5, 9), (5, 10), (9, 10), (11, 13)
         ])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(U(P(0)(1)(U(2)(11))(3)(5)(9)(10)(13))(J(4)(7))(6)(8)(12))')
         self.assertEqual(t.modular_width(), 8)
 
@@ -371,7 +371,7 @@ def test_modular_decomposition_13(self):
         G = nx.empty_graph(8)
         G.add_edges_from([(0, 3), (0, 7), (1, 3), (1, 6), (2, 3), (2, 4), (2, 5), (3, 4), (3, 6), (3, 7), (4, 5), (4, 6)])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(J(0)(7))(1)(2)(3)(4)(5)(6))')
         self.assertEqual(t.modular_width(), 7)
 
@@ -384,7 +384,7 @@ def test_modular_decomposition_14(self):
             (3, 8), (3, 10), (3, 11), (3, 12), (4, 11), (5, 10), (6, 11), (6, 12),
             (7, 12), (8, 9), (8, 12), (9, 10), (9, 12), (10, 11)
         ])
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(U(0)(7))(1)(2)(3)(4)(5)(6)(8)(9)(10)(11)(12))')
         self.assertEqual(t.modular_width(), 12)
 
@@ -398,7 +398,7 @@ def test_modular_decomposition_15(self):
             (5, 12), (6, 13), (7, 8), (8, 12), (9, 11), (9, 13)
         ])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13))')
         self.assertEqual(t.modular_width(), 14)
 
@@ -409,12 +409,12 @@ def test_modular_decomposition_15(self):
     def test_modular_decomposition_16(self):
         G = nx.empty_graph(25)
         G.add_edges_from([
-            (0, 1),(0, 2),(0, 3),(1, 4),(1, 5),(1, 6),(2, 7),(2, 8),(2, 9),
-            (3, 10),(3, 11),(3, 12),(4, 13),(4, 14),(4, 15),(5, 16),(5, 17),(5, 18),
-            (6, 19),(6, 20),(6, 21),(7, 22),(7, 23),(7, 24),
+            (0, 1), (0, 2), (0, 3), (1, 4), (1, 5), (1, 6), (2, 7), (2, 8), (2, 9),
+            (3, 10), (3, 11), (3, 12), (4, 13), (4, 14), (4, 15), (5, 16), (5, 17), (5, 18),
+            (6, 19), (6, 20), (6, 21), (7, 22), (7, 23), (7, 24),
         ])
 
-        t = modular_decomposition(G, sorted=True, solver='linear')
+        t = modular_decomposition(G, sorted=True, solver='linear', verify=True)
         self.assertEqual(str(t), '(P(0)(1)(2)(3)(4)(5)(6)(7)(U(8)(9))(U(10)(11)(12))(U(13)(14)(15))(U(16)(17)(18))(U(19)(20)(21))(U(22)(23)(24)))')
         self.assertEqual(t.modular_width(), 14)
         self.check_property(G)
@@ -423,7 +423,7 @@ def test_modular_decomposition_random(self):
         rand = Random(12345)
         min_n = 5
         max_n = 20
-        num_iterations = 5
+        num_iterations = 2
         # num_iterations = 1000
         def ps(n): return [1 / n, 0.1, 0.2, 0.3]
 
@@ -437,15 +437,15 @@ def test_modular_decomposition_random_mw_bounded(self):
         rand = Random(12345)
         min_n = 30
         max_n = 50
-        step_n = 5
-        num_iterations = 10
+        step_n = 10
+        num_iterations = 2
         mw = 5
         p = 0.5
 
         for n in range(min_n, max_n + 1, step_n):
             for _ in range(num_iterations):
                 G = self.generate_mw_bounded_graph(n, mw, p, rand)
-                t = modular_decomposition(G, solver='linear')
+                t = modular_decomposition(G, solver='linear', verify=True)
                 self.assertLessEqual(t.modular_width(), mw)
                 self.check_property(G)
 
@@ -454,11 +454,11 @@ def test_modular_decomposition_random_mw_bounded_large(self):
         rec_lim = sys.getrecursionlimit()
         rand = Random(12345)
 
-        for _ in range(10):
-            G = self.generate_mw_bounded_graph(500, 4, 0.5, rand)
+        for _ in range(2):
+            G = self.generate_mw_bounded_graph(100, 4, 0.5, rand)
 
             sys.setrecursionlimit(70)
-            t = modular_decomposition(G, solver='linear')
+            t = modular_decomposition(G, solver='linear', verify=True)
             sys.setrecursionlimit(rec_lim)
 
             self.assertLessEqual(t.modular_width(), 4)
diff --git a/src/test/python/modular/tree/test_RootedForest.py b/src/test/python/modular/tree/test_RootedForest.py
index 7e7b7a4..c2013b6 100644
--- a/src/test/python/modular/tree/test_RootedForest.py
+++ b/src/test/python/modular/tree/test_RootedForest.py
@@ -7,7 +7,7 @@ def parse_bool(s: str): return [c != '0' for c in s]
 
 
 def create_forest():
-    tree = RootedForest()  # type: RootedForest[int]
+    tree = RootedForest(verify=True)  # type: RootedForest[int]
     nodes = [tree.create_node(i) for i in range(20)]
 
     #  0    3       10   13