bug fix: when hybrid edge deletion leads simplified parallel edges at…

… hybrid with 2 children (#207)

src/deleteHybrid.jl

@@ -193,7 +193,7 @@ function deleteHybrid!(node::Node,net::HybridNetwork,minor::Bool, blacklist::Boo
             end
         end
         hybindex = findfirst([e.hybrid for e in other2.edge]);
-        hybindex != nothing || error("didn't find hybrid edge in other2")
+        isnothing(hybindex) && error("didn't find hybrid edge in other2")
         if(hybindex == 1)
             treeedge1 = other2.edge[2];
             treeedge2 = other2.edge[3];
@@ -387,6 +387,7 @@ function deletehybridedge!(net::HybridNetwork, edge::Edge,
         # below: won't delete n1, delete edge instead
     end
 
+    formernumhyb = net.numHybrids
     # next: delete n1 recursively, or delete edge and delete n2 recursively.
     # keep n2 if it has 4+ edges (or if nofuse). 1 edge should never occur.
     #       If root, would have no parent: treat network as unrooted and change the root.
@@ -408,6 +409,9 @@ function deletehybridedge!(net::HybridNetwork, edge::Edge,
                     simplify=simplify, unroot=unroot, multgammas=multgammas,
                     keeporiginalroot=keeporiginalroot)
     end
+    if net.numHybrids != formernumhyb # deleteleaf! does not update containRoot
+        allowrootbelow!(net)
+    end
     return net
 end
 
@@ -429,7 +433,7 @@ function undoPartition!(net::HybridNetwork, hybrid::Node, edgesInCycle::Vector{E
                 p = splice!(net.partition,i)
                 @debug "after splice, p partition has edges $([e.number for e in p.edges]) and cycle $(p.cycle)"
                 ind = findfirst(isequal(hybrid.number), p.cycle)
-                ind != nothing || error("hybrid not found in p.cycle")
+                isnothing(ind) && error("hybrid not found in p.cycle")
                 deleteat!(p.cycle,ind) #get rid of that hybrid number
                 cycles = vcat(cycles,p.cycle)
                 edges = vcat(edges,p.edges)

src/manipulateNet.jl

@@ -901,8 +901,8 @@ tree edges with γ<1, or with reticulations in which the two parent
 `keeporiginalroot`: if true, keep the root even if it is of degree one
 (forcing `unroot` to be false).
 
-Warning: does **not** update attributes related to level-1 networks,
-such as inCycle, partition, gammaz, etc.
+Warning: does **not** update edges' `containRoot` nor attributes
+related to level-1 networks such as inCycle, partition, gammaz, etc.
 Does not require branch lengths, and designed to work on networks
 of all levels.
 """
@@ -993,8 +993,8 @@ function deleteleaf!(net::HybridNetwork, nodeNumber::Integer;
             any(getchild(e) ≡ cn && e !== e1 && e !==e2 for e in cn.edge) &&
                 error("root has 2 hybrid edges, but their common child has an extra parent")
             removeEdge!(cn,e1); removeEdge!(cn,e2)
-            removeHybrid!(net,cn) # removes n1 from net.hybrid, updates net.numHybrids
-            cn.hybrid = false
+            removeHybrid!(net,cn) # removes cn from net.hybrid, updates net.numHybrids
+            cn.hybrid = false # !! allowrootbelow! not called: would require correct isChild1
             empty!(e1.node); empty!(e2.node)
             deleteEdge!(net,e1,part=false); deleteEdge!(net,e2,part=false)
             empty!(nodei.edge)
@@ -1032,11 +1032,13 @@ function deleteleaf!(net::HybridNetwork, nodeNumber::Integer;
             pn  = getparent(e1)
             if pn ≡ getparent(e2)
                 # e1 and e2 have same child and same parent. Remove e1.
-                e2.hybrid=false;
-                e2.isMajor=true;
+                e2.hybrid = false # assumes bicombining at cn: no third hybrid parent
+                e2.isMajor = true
                 e2.gamma = addBL(e1.gamma, e2.gamma)
                 removeEdge!(pn,e1); removeEdge!(cn,e1)
                 deleteEdge!(net,e1,part=false)
+                removeHybrid!(net,cn) # removes cn from net.hybrid, updates net.numHybrids
+                cn.hybrid = false # !! allowrootbelow! not called: would require correct isChild1
                 # call recursion again because pn and/or cn might be of degree 2 (or even 1).
                 deleteleaf!(net, cn.number; nofuse = nofuse, simplify=simplify, unroot=unroot,
                             multgammas=multgammas, keeporiginalroot=keeporiginalroot)
@@ -1234,6 +1236,22 @@ function allowrootbelow!(n::Node, pe::Edge)
     end
     return nothing
 end
+"""
+    allowrootbelow!(net::HybridNetwork)
+
+Set `containRoot` to `true` for each edge below the root node, then
+traverses `net` in preorder to update `containRoot` of all edges (stopping
+at hybrid nodes): see the other methods.
+Assumes correct `isChild1` edge field.
+"""
+function allowrootbelow!(net::HybridNetwork)
+    rn = net.node[net.root]
+    for e in rn.edge
+        if e.containRoot
+            allowrootbelow!(getchild(e), e)
+        end
+    end
+end
 
 """
     unzip_canonical!(net::HybridNetwork)

test/test_compareNetworks.jl

@@ -84,6 +84,14 @@ net = deepcopy(net0)
 @test writeTopology(deletehybridedge!(net, net.edge[5]), round=true) == "((#H1:1.0::0.1,b:1.0):1.0,c:1.0,(a:1.0)#H1:3.0::0.9);"
 @test writeTopology(deletehybridedge!(net0, net0.edge[5],false,true,false,false), round=true) ==
   "((#H2:1.0::0.2,((a:1.0)#H1:1.0::0.9)#H2:3.0::0.8):1.0,(#H1:1.0::0.1,b:1.0):1.0,c:1.0);"
+
+# level-2 degree-2 blob simplifying to parallel edges, + polytomy below
+net0 = readTopology("(africa_east:0.003,((#H3:0::0.003,(non_africa_west:0.2,non_africa_east:0.2)#H1:0.3::1):0.3,(#H1:0::0)#H3:0::0.997)H2:0);");
+PhyloNetworks.deletehybridedge!(net0, net0.edge[2], false, false, false, true, false)
+@test all(!n.hybrid for n in net0.node)
+@test all(e.containRoot for e in net0.edge)
+@test writeTopology(net0) == "(africa_east:0.003,(non_africa_west:0.2,non_africa_east:0.2)H1:0.0);"
+
 end # of testing deletehybridedge!
 
 @testset "testing deleteleaf! and hardwiredClusterDistance" begin