Handle self-loops in SPQR tree output

.gitignore

@@ -10,3 +10,4 @@ diff.py
 dist.py
 out
 .DS_Store
+*.spqr

src/main.cpp

@@ -18306,6 +18306,9 @@ namespace solver
                 ccNodes[component[v]].push_back(v);
             }
 
+            // Track which edges were output already to avoid duplicating self-loops
+            std::unordered_set<spqr_compat::edge> self_loops_output;
+
             const int kBlockProgressStep = 256;
             const int kLogEachBlockIfLeq = 50;
             const int kLargeBlockNodes = 200000; 
@@ -18340,12 +18343,23 @@ namespace solver
 
                 for (spqr_compat::node vOrig : ccNodes[ccIdx])
                 {
+                    self_loops_output.clear();
                     C.G.forEachAdj(vOrig, [&](node /*neighbor*/, edge e) {
                         if (C.G.source(e) != vOrig) // Only process from source side
                             return;
 
                         spqr_compat::node src = C.G.source(e);
                         spqr_compat::node tgt = C.G.target(e);
+
+                        // Self-loops cannot be deduplicated by the above condition about the source side.
+                        // Hence, we deduplicate them via a set.
+                        if (src == tgt) {
+                            if (self_loops_output.count(e)) {
+                                return;
+                            } else {
+                                self_loops_output.insert(e);
+                            }
+                        }
 
                         spqr_compat::node srcCc = origToCc[src];
                         spqr_compat::node tgtCc = origToCc[tgt];
@@ -18392,6 +18406,10 @@ namespace solver
                 std::vector<spqr_compat::node> tmpNodes;
                 tmpNodes.reserve(1024);
 
+                // We use this to track which self-loops have been output.
+                // We want to output each self-loop exactly once, assigned to an arbitrary block or SPQR node.
+                self_loops_output.clear();
+
                 for (spqr_compat::node bNode : bc.bcTree().nodes)
                 {
                     if (bc.typeOfBNode(bNode) != spqr_compat::BCTree::BNodeType::BComp)
@@ -18530,6 +18548,30 @@ namespace solver
                             std::string eName = makeId("E", edgeCtr++);
                             writeE(eName, blockName, C.node2name[uOrig], C.node2name[vOrig]);
                         }
+
+                        // The BC tree contains no self-loops, hence we need to output self-loops manually.
+                        for (spqr_compat::node vCc : tmpNodes) {
+                            spqr_compat::node vOrig = ccToOrig[vCc];
+
+                            C.G.forEachAdj(vOrig, [&](node /*neighbor*/, edge e) {
+                                spqr_compat::node src = C.G.source(e);
+                                spqr_compat::node tgt = C.G.target(e);
+                                if (src != tgt) // Only process self loops
+                                    return;
+
+                                // Self-loops are yielded twice by the iterator and may occur in multiple blocks and SPQR nodes,
+                                // hence we need to deduplicate them manually.
+                                if (self_loops_output.count(e)) {
+                                    return;
+                                } else {
+                                    self_loops_output.insert(e);
+                                }
+
+                                std::string eName = makeId("E", edgeCtr++);
+                                writeE(eName, blockName, C.node2name[vOrig], C.node2name[vOrig]);
+                            });
+                        }
+
                         continue;
                     }
 
@@ -18731,6 +18773,31 @@ namespace solver
                                     C.node2name[uOrig],
                                     C.node2name[vOrig]);
                             }
+
+                            // The SPQR tree contains no self-loops, hence we need to output self-loops manually.
+                            for (spqr_compat::node skelNode : skel.getGraph().nodes) {
+                                spqr_compat::node vB = skel.original(skelNode);
+                                spqr_compat::node vCc = blockToCC[vB];
+                                spqr_compat::node vOrig = ccToOrig[vCc];
+
+                                C.G.forEachAdj(vOrig, [&](node /*neighbor*/, edge e) {
+                                    spqr_compat::node src = C.G.source(e);
+                                    spqr_compat::node tgt = C.G.target(e);
+                                    if (src != tgt) // Only process self loops
+                                        return;
+
+                                    // Self-loops are yielded twice by the iterator and may occur in multiple blocks and SPQR nodes,
+                                    // hence we need to deduplicate them manually.
+                                    if (self_loops_output.count(e)) {
+                                        return;
+                                    } else {
+                                        self_loops_output.insert(e);
+                                    }
+
+                                    std::string eName = makeId("E", edgeCtr++);
+                                    writeE(eName, spqrNodeNames[treeNode], C.node2name[vOrig], C.node2name[vOrig]);
+                                });
+                            }
                         }
 
                         if (logThisBlock)

test_graphs/self_loops.gfa

@@ -0,0 +1,16 @@
+H	VN:Z:1.0
+S	a	*
+S	b	*
+S	c	*
+S	d	*
+S	e	*
+S	f	*
+S	g	*
+L	a	+	a	-	0M
+L	b	+	b	+	0M
+L	c	-	d	-	0M
+L	c	+	c	+	0M
+L	e	+	f	+	0M
+L	f	+	g	+	0M
+L	g	+	e	+	0M
+L	e	+	e	+	0M