Flip Hybrid Edge Move (JuliaPhylo#139)

add new structure optimization move to flip a hybrid edge's direction with careful attention to new root selection if hybrid edge to flip is below the root
Co-authored-by: Cecile Ane <[email protected]>

docs/src/lib/public.md

@@ -50,6 +50,7 @@ blobDecomposition
 treeedgecomponents
 getlabels
 nparams
+mapindividuals
 nni!
 ```
 
@@ -140,5 +141,4 @@ maxParsimonyNet
 nstates
 stationary
 empiricalDNAfrequencies
-mapindividuals
 ```

src/PhyloNetworks.jl

@@ -57,6 +57,7 @@ module PhyloNetworks
         readTableCF,
         readTableCF!,
         writeTableCF,
+        mapAllelesCFtable,
         readInputTrees,
         readNexusTrees,
         summarizeDataCF,
@@ -77,7 +78,6 @@ module PhyloNetworks
         rotate!,
         setLength!,
         setGamma!,
-        mapAllelesCFtable,
         deleteHybridThreshold!,
         displayedTrees,
         majorTree,
@@ -90,6 +90,8 @@ module PhyloNetworks
         biconnectedComponents,
         blobDecomposition!,
         blobDecomposition,
+        mapindividuals,
+        nni!,
         checkroot!,
         treeedgecomponents,
         ## Network Bootstrap
@@ -153,8 +155,6 @@ module PhyloNetworks
         readfastatodna,
         stationary,
         empiricalDNAfrequencies,
-        nni!,
-        mapindividuals,
         # phyLiNC,
         # neighbor joining
         nj

src/auxiliary.jl

@@ -99,7 +99,7 @@ function setNode!(edge::Edge, node::Node)
         else
           if edge.hybrid
             if node.hybrid
-                @debug (edge.node[1].hybrid ? "hybrid edge $(edge.number) has two hybrid nodes" : "")
+                # @debug (edge.node[1].hybrid ? "hybrid edge $(edge.number) has two hybrid nodes" : "")
                 edge.isChild1 = false;
 	        else
 	            edge.node[1].hybrid || error("hybrid edge $(edge.number) has no hybrid nodes");

src/manipulateNet.jl

@@ -444,6 +444,9 @@ end
 
 Removes `i`th node in net.node, if it is of degree 2.
 The parent and child edges of this node are fused.
+If either of the edges is hybrid, the hybrid edge is retained. Otherwise, the
+edge with the lower edge number is retained.
+
 Reverts the action of breakedge!.
 
 returns the fused edge.
@@ -612,8 +615,10 @@ compatible with the direction of existing hybrid edges.
 Relies on hybrid nodes having exactly 1 major hybrid parent edge,
 but checks for that if checkMajor=true.
 
-Warning: Assumes that isChild1 is correct on hybrid edges
+Warnings:
+1. Assumes that isChild1 is correct on hybrid edges
 (to avoid changing the identity of which nodes are hybrids and which are not).
+2. Does not check for cycles (to maintain a network's DAG status)
 
 Returns the network. Throws a 'RootMismatch' Exception if the root was found to
 conflict with the direction of any hybrid edge.

src/moves_semidirected.jl

@@ -326,7 +326,6 @@ function nni!(net::HybridNetwork, e::Edge, nohybridladder::Bool=true, no3cycle::
             return nothing
         end
     end
-    hybparent = getParent(e).hybrid
     nnirange = 0x01:nnimax(e) # 0x01 = 1 but UInt8 instead of Int
     nnis = Random.shuffle(nnirange)
     for nummove in nnis # iterate through all possible NNIs, but in random order
@@ -896,3 +895,188 @@ function moveroot!(net::HybridNetwork, constraints=TopologyConstraint[]::Vector{
     # if we get here: none of the root positions worked
     return nothing
 end
+
+"""
+    fliphybrid!(net::HybridNetwork, minor=true::Bool, nohybridladder=false::Bool,
+                constraints=TopologyConstraint[]::Vector{TopologyConstraint})
+
+Cycle through hybrid nodes in random order until an admissible flip is found.
+At this hybrid node, flip the indicated hybrid parent edge (minor or major).
+
+If an admissible flip is found, return the tuple: newhybridnode, flippededge, oldchildedge.
+Otherwise, return nothing.
+
+The flip can be undone with
+`fliphybrid!(net, newhybridnode, minor, constraints)`, or
+`fliphybrid!(net, newhybridnode, !flippededge.isMajor, constraints)`
+more generally, such as if the flipped edge had its γ modified after the
+original flip.
+
+*Warnings*
+- if the root needed to be reset and if the original root was of degree 2,
+  then a node of degree 2 remains in the modified network.
+- undoing the flip may not recover the original root in case
+  the root position was modified during the original flip.
+"""
+function fliphybrid!(net::HybridNetwork, minor=true::Bool,
+                     nohybridladder=false::Bool,
+                     constraints=TopologyConstraint[]::Vector{TopologyConstraint})
+    hybridindex = Random.shuffle(1:length(net.hybrid)) # indices in net.hybrid
+    while !isempty(hybridindex) # all minor edges
+        i = pop!(hybridindex)
+        undoinfo = fliphybrid!(net, net.hybrid[i], minor, nohybridladder, constraints)
+        if !isnothing(undoinfo) # if it failed, the network was not changed
+            return undoinfo
+        end # else, continue until explored all i in hybridindex
+    end
+    return nothing
+end
+
+"""
+    fliphybrid!(net::HybridNetwork, hybridnode::Node, minor=true::Bool,
+                nohybridladder=false::Bool,
+                constraints=TopologyConstraint[]::Vector{TopologyConstraint})
+
+Flip the direction of a single hybrid edge:
+the minor parent edge of `hybridnode` by default,
+or the major parent edge if `minor` is false.
+The parent node of the hybrid edge becomes the new hybrid node.
+The former hybrid edge partner is converted to a tree edge (with γ=1),
+and `hybridnode` becomes a tree node.
+
+For the flip to be admissible, the new network must be a semi-directed
+phylogenetic network: with a root such that the rooted version is a DAG.
+If `nohybridladder` is false (default), the flip may create a hybrid ladder
+If `nohybridladder` is true and if the flip would create a hybrid ladder,
+then the flip is not admissible.
+A hybrid ladder is when a hybrid child of another hybrid.
+
+The new hybrid partner is an edge adjacent to the new hybrid node,
+such that the flip is admissible (so it must be a tree edge).
+The flipped edge retains its original γ.
+The new hybrid edge is assigned inheritance 1-γ.
+
+Output:
+`(newhybridnode, flippededge, oldchildedge)` if the flip is admissible,
+`nothing` otherwise.
+
+The network is unchanged if the flip is not admissible.
+If the flip is admissible, the root position may be modified, and
+the direction of tree edges (via `isChild1`) is modified accordingly. If the
+root needs to be modified, then the new root is set to the old hybrid node.
+
+The index of the new hybrid node in `net.hybrid` is equal to that of the
+old `hybridnode`.
+
+Warning: Undoing this move may not recover the original root if
+the root position was modified.
+"""
+function fliphybrid!(net::HybridNetwork, hybridnode::Node, minor=true::Bool,
+                     nohybridladder=false::Bool,
+                     constraints=TopologyConstraint[]::Vector{TopologyConstraint})
+    #= for species constraints, there is nothing to check, because hybrids cannot point into or come out of the group
+    for con in constraints
+        if con.type == # 0x02/0x03 for types 2 and 3, need to consider more cases
+            return nothing
+        end
+    end
+    =#
+    runDirectEdges = false
+    edgetoflip, edgetokeep = minor ?
+        (getMinorParentEdge(hybridnode), getMajorParentEdge(hybridnode)) :
+        (getMajorParentEdge(hybridnode), getMinorParentEdge(hybridnode))
+    oldchildedge = getChildEdge(hybridnode)
+    newhybridnode = getParent(edgetoflip)
+    if newhybridnode.hybrid # already has 2 parents: cannot had a third.
+        return nothing
+    end
+    ## choose newhybridedge ##
+    p2 = getParent(edgetokeep) # parent node
+    isdesc = Bool[]
+    for e in newhybridnode.edge # is p2 undirected descendant of nhn via this edge?
+        isp2desc = false
+        #= isp2desc = true means:
+        there is a semi-directed path  newhybridnode -e-> neighbor --...-> p2
+        so: flipping hybridedge without making e the new partner would:
+        - make e a child of newhybridnode, and
+        - create a directed cycle: p2 -hybridedge-> newhybridnode -e-> neibr --...-> p2
+        so we HAVE to make e the new partner hybrid edge if its "isp2desc" is true
+        =#
+        if e !== edgetoflip # e cannot be hybrid --e-> nhn because earlier check
+            neibr = getOtherNode(e, newhybridnode) # neighbor of new hybrid node via e
+            if !neibr.leaf && (p2 === neibr ||
+                isdescendant_undirected(p2, neibr, e))
+                isp2desc = true
+            end
+        end
+        push!(isdesc, isp2desc)
+    end
+    sum_isdesc = sum(isdesc)
+    if sum_isdesc > 1 # if 2+ edges have a semi-directed path to p2,
+        # flipping either would create a semi-directed cycle via the other(s)
+        return nothing
+    end
+    if sum_isdesc == 0 # no risk to create directed cycle.
+        # choose the current parent edge of nhn to: keep its current direction
+        # and keep reaching all nodes from a single root
+        newhybridedge = getMajorParentEdge(newhybridnode)
+    else
+        newhybridedge = newhybridnode.edge[findfirst(isdesc)]
+        if newhybridedge.hybrid # cannot flip another hybrid edge, but we needed
+            return nothing      # to calculate its isp2desc for total # of true's
+        end
+    end
+    # @debug "edgetoflip is $edgetoflip, edgetokeep is $edgetokeep, newhybridedge is $newhybridedge"
+    if nohybridladder
+        # case when edgetoflip would be bottom rung of ladder
+        if getOtherNode(newhybridedge, newhybridnode).hybrid
+            return nothing
+        end
+        # case when edgetoflip would be the top rung: happens when
+        # newhybridnode = center point of a W structure initially
+        for e in newhybridnode.edge
+            if e !== edgetoflip && e.hybrid # newhybridedge is NOT hybrid, so far
+                return nothing
+            end
+        end
+    end
+    # change hybrid status and major status for nodes and edges
+    hybridnode.hybrid = false
+    edgetokeep.hybrid = false
+    newhybridnode.hybrid = true
+    newhybridedge.hybrid = true
+    newhybridedge.isMajor = minor
+    edgetokeep.isMajor = true
+    # update node order to keep isChild1 attribute of hybrid edges true
+    edgetoflip.isChild1 = !edgetoflip.isChild1 # just switch
+    if getChild(newhybridedge) !== newhybridnode # includes the case when the newhybridnode was the root
+        # then flip newhybridedge too: make it point towards newhybridnode
+        newhybridedge.isChild1 = !newhybridedge.isChild1
+        net.root = findfirst(n -> n === hybridnode, net.node)
+        runDirectEdges = true # many edges are likely to need to have directions flipped here
+    end
+    # give new hybridnode a name
+    newhybridnode.name = hybridnode.name
+    hybridnode.name = "" # remove name from former hybrid node
+    # update gammas
+    newhybridedge.gamma = edgetokeep.gamma
+    edgetokeep.gamma = 1.0
+    # update hybrids in network (before directEdges!)
+    hybridindex = findfirst(n -> n === hybridnode, net.hybrid)
+    net.hybrid[hybridindex] = newhybridnode
+    if runDirectEdges # direct edges only if root is moved
+        directEdges!(net)
+    else # if not, only update containRoot attributes
+        # norootbelow in child edges of newhybridedge
+        for ce in newhybridnode.edge
+            ce !== newhybridedge || continue # skip e
+            getParent(ce) === newhybridnode || continue # skip edges that aren't children of cn
+            norootbelow!(ce)
+        end
+        # allow root for edges below old hybrid node
+        if edgetokeep.containRoot #else, already forbidden below due to hybrid above
+            allowrootbelow!(hybridnode, edgetokeep) # old hybrid node
+        end
+    end
+    return newhybridnode, edgetoflip, oldchildedge
+end