v0.15.2: redirect edges if re-rooting failed. new: shrinkedge! (#184)

Project.toml

@@ -1,7 +1,7 @@
 name = "PhyloNetworks"
 uuid = "33ad39ac-ed31-50eb-9b15-43d0656eaa72"
 license = "MIT"
-version = "0.15.1"
+version = "0.15.2"
 
 [deps]
 BioSequences = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59"

docs/src/man/bootstrap.md

@@ -35,7 +35,7 @@ bootnet = bootsnaq(net0, bootTrees, hmax=1, nrep=10, runs=3,
 
 The bootstrap networks are saved in the `boostrap.out` file, so they
 can be read in a new session with
-`bootnet = readMultiTopology("bootsnap.out")`. To save the bootstrap networks to
+`bootnet = readMultiTopology("bootsnaq.out")`. To save the bootstrap networks to
 a different file (perhaps after having re-rooted them with an
 outgroup), we could do this: `writeMultiTopology(bootnet, "bootstrapNets.tre")`.
 
@@ -225,7 +225,7 @@ We can plot the bootstrap values of the 2 hybrid edges in the best network:
 ```@example bootstrap
 R"svg(name('boot_net_net.svg'), width=4, height=4)" # hide
 R"par"(mar=[0,0,0,0]) # hide
-plot(net1, edgelabel=BSe[!,[:edge,:BS_hybrid_edge]]);
+plot(net1, edgelabel=BSe[:,[:edge,:BS_hybrid_edge]]);
 R"dev.off()" # hide
 nothing # hide
 ```
@@ -255,7 +255,7 @@ select!(tmp, [:edge,:BS_hybrid_edge])            # select 2 columns only
 rename!(tmp, :BS_hybrid_edge => :proportion)     # rename those columns, to match names in BSe_tree
 rename!(tmp, :edge => :edgeNumber)
 tmp = vcat(BSe_tree, tmp)
-plot(net1, edgelabel=tmp, nodelabel=BSn[!, [:hybridnode,:BS_hybrid_samesisters]])
+plot(net1, edgelabel=tmp, nodelabel=BSn[:, [:hybridnode,:BS_hybrid_samesisters]])
 ```
 
 ### Who are the hybrids in bootstrap networks?

docs/src/man/netmanipulation.md

@@ -28,7 +28,8 @@ To modify a network, for example:
 - [`deleteleaf!`](@ref),
   [`deleteaboveLSA!`](@ref) (the "least stable ancestor" may be different from the root)
 - [`deleteHybridThreshold!`](@ref) to simplify a network by deleting edges with small γ's
-- [`removedegree2nodes!`](@ref), [`shrink3cycles!`](@ref), [`shrink2cycles!`](@ref)
+- [`removedegree2nodes!`](@ref), [`shrink3cycles!`](@ref), [`shrink2cycles!`](@ref),
+  [`PhyloNetworks.shrinkedge!`](@ref)]
 - [`PhyloNetworks.addleaf!`](@ref),
   [`PhyloNetworks.deletehybridedge!`](@ref), [`PhyloNetworks.addhybridedge!`](@ref)
 - [`nni!`](@ref) (nearest neighbor interchange),

src/auxiliary.jl

@@ -260,8 +260,9 @@ end
 
 Return child edge below a hybrid node.
 
-Warning: Does not check that the node is a hybrid.
-If not a hybrid, returns the first child edge.
+Warning: Does not check that the node is a hybrid, or that the node has a single
+child edge. If not a hybrid or if there's a polytomy (hybrid with 2 children),
+returns the first child edge.
 """
 @inline function getChildEdge(node::Node)
     for e in node.edge
@@ -1459,6 +1460,49 @@ function hashybridladder(net::HybridNetwork)
     return false
 end
 
+"""
+    shrinkedge!(net::HybridNetwork, edge::Edge)
+
+Delete `edge` from net, provided that it is a non-external tree edge.
+Specifically: delete its child node (as determined by `isChild1`) and connect
+all edges formerly incident to this child node to the parent node of `edge`,
+thus creating a new polytomy, unless the child was of degree 2.
+
+Warning: it's best for `isChild1` to be in sync with the root for this. If not,
+the shrinking may fail (if `edge` is a tree edge but its "child" is a hybrid)
+or the root may change arbitrarily (if the child of `edge` is the root).
+
+Output: true if the remaining node (parent of `edge`) becomes a hybrid node with
+more than 1 child after the shrinking; false otherwise (e.g. no polytomy was
+created, or the new polytomy is below a tree node)
+"""
+function shrinkedge!(net::HybridNetwork, edge2shrink::Edge)
+    edge2shrink.hybrid && error("cannot shrink hybrid edge number $(edge2shrink.number)")
+    cn = getChild(edge2shrink)
+    cn.hybrid && error("cannot shrink tree edge number $(edge2shrink.number): its child node is a hybrid. run directEdges! ?")
+    pn = getParent(edge2shrink)
+    (cn.leaf || pn.leaf) &&
+      error("won't shrink edge number $(edge2shrink.number): it is incident to a leaf")
+    removeEdge!(pn,edge2shrink)
+    empty!(edge2shrink.node) # should help gc
+    for ee in cn.edge
+        ee !== edge2shrink || continue
+        cn_index = findfirst(x -> x === cn, ee.node)
+        ee.node[cn_index] = pn # ee.isChild1 remains synchronized
+        push!(pn.edge, ee)
+    end
+    pn.hasHybEdge = any(e -> e.hybrid, pn.edge)
+    empty!(cn.edge) # should help to garbage-collect cn
+    deleteEdge!(net, edge2shrink; part=false)
+    deleteNode!(net, cn)
+    badpolytomy = false
+    if pn.hybrid # count the number of pn's children, without relying on isChild1 of tree edges
+        nc = sum((!e.hybrid || getChild(e) !== pn) for e in pn.edge)
+        badpolytomy = (nc > 1)
+    end
+    return badpolytomy
+end
+
 @doc raw"""
     shrink2cycles!(net::HybridNetwork, unroot=false::Bool)
 
@@ -1691,15 +1735,14 @@ end
     adjacentedges(centeredge::Edge)
 
 Vector of all edges that share a node with `centeredge`.
-Warning: assumes that there aren't "duplicated" edges, that is, no 2-cycles.
 """
 function adjacentedges(centeredge::Edge)
     n = centeredge.node
     length(n) == 2 || error("center edge is connected to $(length(n)) nodes")
     @inbounds edges = copy(n[1].edge) # shallow copy, to avoid modifying the first node
     @inbounds for ei in n[2].edge
         ei === centeredge && continue # don't add the center edge again
-        # a second edge between nodes n[1] and n[2] would appear twice
+        getOtherNode(ei, n[2]) === n[1] && continue # any parallel edge would have been in `edges` already
         push!(edges, ei)
     end
     return edges

src/compareNetworks.jl

@@ -55,7 +55,7 @@ function tree2Matrix(T::HybridNetwork, S::Union{Vector{String},Vector{Int}}; roo
 end
 
 """
-`hardwiredClusters(net::HybridNetwork, S::Union{Vector{String},Vector{Int}})`
+`hardwiredClusters(net::HybridNetwork, S::Union{AbstractVector{String},AbstractVector{Int}})`
 
 
 Returns a matrix describing all the hardwired clusters in a network.
@@ -70,7 +70,7 @@ Both parent hybrid edges to a given hybrid node only contribute a single row (th
 - next columns: 0/1. 1=descendant of edge, 0=not a descendant, or missing taxon.
 - last column:  10/11 values. 10=tree edge, 11=hybrid edge
 """
-function hardwiredClusters(net::HybridNetwork, S::Union{Vector{String},Vector{Int}})
+function hardwiredClusters(net::HybridNetwork, S::Union{AbstractVector{String},AbstractVector{Int}})
     ne = length(net.edge)-net.numTaxa # number of internal branch lengths
     ne -= length(net.hybrid)          # to remove duplicate rows for the 2 parent edges of each hybrid
     if (net.node[net.root].leaf)      # root is leaf: the 1 edge stemming from the root is an external edge
@@ -84,8 +84,8 @@ function hardwiredClusters(net::HybridNetwork, S::Union{Vector{String},Vector{In
     return M
 end
 
-function hardwiredClusters!(node::Node, edge::Edge, ie::Vector{Int}, M::Matrix{Int},
-                            S::Union{Vector{String},Vector{Int}})
+function hardwiredClusters!(node::Node, edge::Edge, ie::AbstractVector{Int}, M::Matrix{Int},
+                            S::Union{AbstractVector{String},AbstractVector{Int}})
     child = getOtherNode(edge,node)
 
     !child.leaf || return 0 # do nothing if child is a leaf.
@@ -133,15 +133,14 @@ end
 
 
 """
-    hardwiredCluster(edge::Edge,taxa::Union{Vector{String},Vector{Int}})
-    hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}})
-    hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}},
-                      visited::Vector{Int})
+    hardwiredCluster(edge::Edge, taxa::Union{AbstractVector{String},AbstractVector{Int}})
+    hardwiredCluster!(v::Vector{Bool}, edge::Edge, taxa)
+    hardwiredCluster!(v::Vector{Bool}, edge::Edge, taxa, visited::Vector{Int})
 
-Calculate the hardwired cluster of `node`, coded a vector of booleans:
-true for taxa that are descendent of nodes, false for other taxa (including missing taxa).
+Calculate the hardwired cluster of `edge`, coded as a vector of booleans:
+true for taxa that are descendent of the edge, false for other taxa (including missing taxa).
 
-The node should belong in a rooted network for which `isChild1` is up-to-date.
+The edge should belong in a rooted network for which `isChild1` is up-to-date.
 Run `directEdges!` beforehand. This is very important, otherwise one might enter an infinite loop,
 and the function does not test for this.
 
@@ -170,16 +169,16 @@ julia> hardwiredCluster(net5.edge[12], taxa) # descendants of 12th edge = CEF
  0
 ```
 """
-function hardwiredCluster(edge::Edge,taxa::Union{Vector{String},Vector{Int}})
+function hardwiredCluster(edge::Edge,taxa::Union{AbstractVector{String},AbstractVector{Int}})
     v = zeros(Bool,length(taxa))
     hardwiredCluster!(v,edge,taxa)
     return v
 end
 
-hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}}) =
+hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{AbstractVector{String},AbstractVector{Int}}) =
     hardwiredCluster!(v,edge,taxa,Int[])
 
-function hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}},
+function hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{AbstractVector{String},AbstractVector{Int}},
                            visited::Vector{Int})
     n = getChild(edge)
     if n.leaf
@@ -193,7 +192,7 @@ function hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String}
     end
     push!(visited, n.number)
     for ce in n.edge
-        if n == getParent(ce)
+        if n === getParent(ce)
             hardwiredCluster!(v,ce,taxa,visited)
         end
     end
@@ -760,7 +759,7 @@ function hardwiredClusterDistance_unrooted!(net1::HybridNetwork, net2::HybridNet
        sister to a hybrid edge, when a the leaf edge is the donor. =#
     for i in length(net1roots):-1:1 # reverse order, to delete some of them
         try
-            rootatnode!(net1, net1roots[i]; verbose=false)
+            rootatnode!(net1, net1roots[i])
             # tricky: rootatnode adds a degree-2 node if i is a leaf,
             #         and delete former root node if it's of degree 2.
         catch e
@@ -771,7 +770,7 @@ function hardwiredClusterDistance_unrooted!(net1::HybridNetwork, net2::HybridNet
     net2roots = [n.number for n in net2.node if !n.leaf]
     for i in length(net2roots):-1:1
         try
-            rootatnode!(net2, net2roots[i]; verbose=false)
+            rootatnode!(net2, net2roots[i])
         catch e
             isa(e, RootMismatch) || rethrow(e)
             deleteat!(net2roots, i)
@@ -780,9 +779,9 @@ function hardwiredClusterDistance_unrooted!(net1::HybridNetwork, net2::HybridNet
     bestdissimilarity = typemax(Int)
     bestns = missing
     for n1 in net1roots
-        rootatnode!(net1, n1; verbose=false)
+        rootatnode!(net1, n1)
         for n2 in net2roots
-            rootatnode!(net2, n2; verbose=false)
+            rootatnode!(net2, n2)
             diss = hardwiredClusterDistance(net1, net2, true) # rooted = true now
             if diss < bestdissimilarity
                 bestns = (n1, n2)

src/manipulateNet.jl

@@ -224,9 +224,9 @@ end
 
 
 """
-    rootatnode!(HybridNetwork, nodeNumber::Integer; index=false::Bool, verbose=true::Bool)
-    rootatnode!(HybridNetwork, Node; verbose=true)
-    rootatnode!(HybridNetwork, nodeName::AbstractString; verbose=true)
+    rootatnode!(HybridNetwork, nodeNumber::Integer; index=false::Bool)
+    rootatnode!(HybridNetwork, Node)
+    rootatnode!(HybridNetwork, nodeName::AbstractString)
 
 Root the network/tree object at the node with name 'nodeName' or
 number 'nodeNumber' (by default) or with index 'nodeNumber' if index=true.
@@ -243,9 +243,8 @@ Warnings:
   the edge adjacent to the leaf. This might add a new node.
 - If the desired root placement is incompatible with one or more hybrids, then
 
-  * a RootMismatch error is thrown; use `verbose=false` to silence
-    the root mismatch info printed before the error is thrown.
-  * the input network will still have some attributes modified.
+  * the original network is restored with its old root and edges' direction.
+  * a RootMismatch error is thrown.
 
 See also: [`rootonedge!`](@ref).
 """
@@ -263,7 +262,7 @@ function rootatnode!(net::HybridNetwork, nodeName::AbstractString; kwargs...)
     rootatnode!(net, tmp[1]; kwargs..., index=true)
 end
 
-function rootatnode!(net::HybridNetwork, nodeNumber::Integer; index=false::Bool, verbose=true::Bool)
+function rootatnode!(net::HybridNetwork, nodeNumber::Integer; index=false::Bool)
     ind = nodeNumber # good if index=true
     if !index
       try
@@ -279,9 +278,9 @@ function rootatnode!(net::HybridNetwork, nodeNumber::Integer; index=false::Bool,
         length(net.node[ind].edge)==1 || error("leaf has $(length(net.node[ind].edge)) edges!")
         pn = getOtherNode(net.node[ind].edge[1], net.node[ind])
         if length(pn.edge) <= 2 # if parent of leaf has degree 2, use it as new root
-            rootatnode!(net, pn.number; verbose=verbose)
+            rootatnode!(net, pn.number)
         else # otherwise, create a new node between leaf and its parent
-            rootonedge!(net,net.node[ind].edge[1]; verbose=verbose)
+            rootonedge!(net,net.node[ind].edge[1])
         end
     else
         rootsaved = net.root
@@ -290,10 +289,11 @@ function rootatnode!(net::HybridNetwork, nodeNumber::Integer; index=false::Bool,
           directEdges!(net)
         catch e
           if isa(e, RootMismatch) # new root incompatible with hybrid directions: revert back
-            verbose && println("RootMismatch: reverting to old root position.")
             net.root = rootsaved
+            directEdges!(net)
           end
-          rethrow(e)
+          throw(RootMismatch("""the desired root is below a reticulation,
+                                reverting to old root position."""))
         end
         if (net.root != rootsaved && length(net.node[rootsaved].edge)==2)
             fuseedgesat!(rootsaved,net) # remove old root node if degree 2
@@ -304,8 +304,8 @@ end
 
 
 """
-    rootonedge!(HybridNetwork, edgeNumber::Integer; index=false::Bool, verbose=true::Bool)
-    rootonedge!(HybridNetwork, Edge; verbose=true::Bool)
+    rootonedge!(HybridNetwork, edgeNumber::Integer; index=false::Bool)
+    rootonedge!(HybridNetwork, Edge)
 
 Root the network/tree along an edge with number `edgeNumber` (by default)
 or with index `edgeNumber` if `index=true`.
@@ -322,7 +322,7 @@ function rootonedge!(net::HybridNetwork, edge::Edge; kwargs...)
     rootonedge!(net, edge.number, index=false; kwargs...)
 end
 
-function rootonedge!(net::HybridNetwork, edgeNumber::Integer; index=false::Bool, verbose=true::Bool)
+function rootonedge!(net::HybridNetwork, edgeNumber::Integer; index=false::Bool)
     ind = edgeNumber # good if index=true
     if !index
       try
@@ -340,12 +340,12 @@ function rootonedge!(net::HybridNetwork, edgeNumber::Integer; index=false::Bool,
       directEdges!(net)
     catch e
       if isa(e, RootMismatch) # new root incompatible with hybrid directions: revert back
-        verbose && println("RootMismatch: reverting to old root position.")
         fuseedgesat!(net.root,net) # reverts breakedge!
         net.root = rootsaved
         directEdges!(net)
       end
-      rethrow(e)
+      throw(RootMismatch("""the desired root is below a reticulation,
+                                reverting to old root position."""))
     end
     if (net.root != rootsaved && length(net.node[rootsaved].edge)==2)
         fuseedgesat!(rootsaved,net) # remove old root node if degree 2

src/parsimony.jl

@@ -1055,7 +1055,7 @@ function maxParsimonyNetRun1!(currT::HybridNetwork, tolAbs::Float64, Nfail::Inte
             flag = proposedTop!(move,newT,true, count,10, movescount,movesfail,false) #N=10 because with 1 it never finds an edge for nni
             newTr = deepcopy(newT) ##rooted version only to compute parsimony
             try
-                rootatnode!(newTr,outgroup; verbose=false)
+                rootatnode!(newTr,outgroup)
             catch
                 flag = false
             end

test/test_auxillary.jl

@@ -38,7 +38,15 @@ PhyloNetworks.addhybridedge!(tree, tree.edge[2], tree.edge[1], true)
 @test PhyloNetworks.hashybridladder(tree)
 end # of testing hashybridladder
 
-@testset "shrink 2/3 cycles" begin
+@testset "shrink edges and cycles" begin
+# shrink 1 edge
+net = readTopology("((A:2.0,(((B1,B2):1.0)0.01)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5);")
+@test !PhyloNetworks.shrinkedge!(net, net.edge[10])
+@test_throws Exception PhyloNetworks.shrinkedge!(net, net.edge[6]) # hybrid edge
+@test_throws Exception PhyloNetworks.shrinkedge!(net, net.edge[3]) # external edge
+@test !PhyloNetworks.shrinkedge!(net, net.edge[5])
+@test  PhyloNetworks.shrinkedge!(net, net.edge[4])
+# shrink cycles
 net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
 @test !shrink2cycles!(net)
 @test !shrink3cycles!(net)

test/test_manipulateNet.jl

@@ -107,7 +107,7 @@ net.root = 5
 @test [net.edge[i].containRoot for i in [9,5,18,2]] == [true for i in 1:4]
 # or: "directEdges! didn't correct containRoot of hyb edges."
 @test_logs rootatnode!(net, -10); # or: rootatnode! complained, node -10
-@test_throws PhyloNetworks.RootMismatch rootatnode!(net, "M"; verbose=false);
+@test_throws PhyloNetworks.RootMismatch rootatnode!(net, "M");
 # println("the rootmismatch about node 5 is good and expected.")
 @test_logs rootonedge!(net, 9); # or: rootonedge! complained, edge 9
 @test_logs PhyloNetworks.fuseedgesat!(27, net);
@@ -140,9 +140,9 @@ net.root=15; # node number -5 (clau: previously -4)
 @test_throws PhyloNetworks.RootMismatch directEdges!(net);
 # occursin(r"non-leaf node 9 had 0 children",e.msg))
 @test_logs rootatnode!(net, -13); # or: rootatnode complained...
-@test_throws PhyloNetworks.RootMismatch rootatnode!(net, -5); # verbose = true this time
+@test_throws PhyloNetworks.RootMismatch rootatnode!(net, -5);
 # occursin(r"non-leaf node 9 had 0 children", e.msg))
-@test_throws PhyloNetworks.RootMismatch rootatnode!(net,"H2"; verbose=false); #try rethrow();
+@test_throws PhyloNetworks.RootMismatch rootatnode!(net,"H2"); #try rethrow();
 # occursin(r"hybrid edge 17 conflicts", e.msg))
 # earlier: """node 12 is a leaf. Will create a new node if needed, to set taxon "10" as outgroup."""
 @test_logs rootatnode!(net,"10");