better cooordinates for non-tree-child networks

src/phylonetworksPlots.jl

@@ -28,24 +28,47 @@ function getEdgeNodeCoordinates(net::HybridNetwork, useEdgeLength::Bool)
     node_yB = zeros(Float64,net.numNodes) # min (B=begin) and max (E=end)
     node_yE = zeros(Float64,net.numNodes) #   of at children's nodes
     nexty = ymax # first tips at the top, last at bottom
-    for i=length(net.node):-1:1 # post-order traversal in major tree
+    # set node_y of leaves: follow cladewise order
+    for i=length(net.node):-1:1
         ni = net.cladewiseorder_nodeIndex[i]
         if net.node[ni].leaf
             node_y[ni] = nexty
             nexty -= 1.0
-        else
-            node_yB[ni]=ymax; node_yE[ni]=ymin;
-            for e in net.node[ni].edge
-                if net.node[ni] == (e.isChild1 ? e.node[2] : e.node[1]) # if e = child of node
-                    if (!e.isMajor) continue; end
+        end
+    end
+    # set node_y of internal nodes: follow post-order
+    for i=length(net.node):-1:1
+        nn = net.nodes_changed[i]
+        !nn.leaf || continue # previous loop took care of leaves
+        ni = findfirst(net.node, nn)
+        node_yB[ni]=ymax; node_yE[ni]=ymin;
+        minor_yB  = ymax; minor_yE  = ymin;
+        nomajorchild=true
+        for e in nn.edge
+            if nn == PhyloNetworks.getParent(e) # if e = child of node
+                if e.isMajor
+                    nomajorchild = false # we found a child edge that is a major edge
                     yy = node_y[findfirst(net.node, PhyloNetworks.getChild(e))]
                     yy!=0 || error("oops, child has not been visited and its y value is 0.")
                     node_yB[ni] = min(node_yB[ni], yy)
                     node_yE[ni] = max(node_yE[ni], yy)
+                elseif nomajorchild # e is minor edge, and no major found so far
+                    yy = node_y[findfirst(net.node, PhyloNetworks.getChild(e))]
+                    yy!=0 || error("oops, child has not been visited and its y value is 0.")
+                    minor_yB = min(minor_yB, yy)
+                    minor_yE = max(minor_yE, yy)
                 end
-                node_y[ni] = (node_yB[ni]+node_yE[ni])/2
             end
         end
+        if nomajorchild # if all children edges are minor hybrid edges. If so: not level 1, not tree-child
+            if minor_yB == minor_yE # one single child. jitter by 0.1 to make the plot readable
+                minor_yB += (minor_yB < (ymax+ymin)/2 ? 0.1 : -0.1)
+                minor_yE = minor_yB
+            end
+            node_yB[ni] = minor_yB
+            node_yE[ni] = minor_yE
+        end
+        node_y[ni] = (node_yB[ni]+node_yE[ni])/2
     end
 
     # setting branch lengths for plotting

test/test_phylonetworkPlots.jl

@@ -53,4 +53,31 @@
     min=[false,false,true,false,false,false,false,false,false],
     x=collect(6.:14),y=collect(26.:34)))
 
+  # example with level-2 network, non-tree child:
+  # one hybrid node ends up as a leaf in the major tree.
+  # no major child edge to follow to set coordinates
+  net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);")
+  @test_nowarn plot(net, :R, showNodeNumber=true, showGamma=true);
+  @test PhyloPlots.getEdgeNodeCoordinates(net, false) == (
+    [5.,4.,1.,3.,3.,3.,2.,4.,4.,2.,1.],
+    [6.,5.,4.,6.,6.,4.,3.,5.,6.,4.,2.],
+    [2.,2.1,2.275,4.,3.,2.1,3.05,2.,1.,1.5,2.275],
+    [2.,2. ,2.1,  4.,3.,2.1,3.05,2.,1.,1.5,2.275],
+    [6.,5.,4., 6.,6.,3.,  6.,4., 2.,   1.],
+    [2.,2.,2.1,4.,3.,3.05,1.,1.5,2.275,2.275],
+    [0.,2.,2.1,0.,0.,2.1,0.,1.,1.5, 2.275],
+    [0.,2.,2.1,0.,0.,4. ,0.,2.,3.05,2.275],
+    1.,6.,1.,4.)
+  net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2)S1,(#H1,A)S2)S3)S4);")
+  @test_nowarn plot(net, :R, showNodeNumber=true, showGamma=true);
+  @test PhyloPlots.getEdgeNodeCoordinates(net, false) == (
+    [5.,4.,1.,3.,3.,3.,2.,4.,4.,2.,1.],
+    [6.,5.,4.,6.,6.,4.,3.,5.,6.,4.,2.],
+    [2.,2.1,2.1,4.,3.,3.5,3.5,2.,1.,1.5,2.5],
+    [2.,2. ,2.1,4.,3.,2.1,3.5,2.,1.,1.5,2.5],
+    [6.,5.,4., 6.,6.,3., 6.,4., 2., 1.],
+    [2.,2.,2.1,4.,3.,3.5,1.,1.5,2.5,2.3],
+    [0.,2.,2.1,0.,0.,3.,0.,1.,1.5,2.1],
+    [0.,2.,2.1,0.,0.,4.,0.,2.,3.5,2.5],
+    1.,6.,1.,4.0)
 end

test/test_rexport.jl

@@ -32,8 +32,8 @@ class(tree2r) = "phylo"
 @rput net2
 @test convert(Bool, R"isTRUE(all.equal(net2, tree2r))")
 
-# network, h=1, some missing gamma values
-s = "(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);";
+# network, h=1, missing gamma values
+s = "(((A:4.0,(B:1.0)#H1:1.1):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);";
 phy1 = PhyloPlots.rexport(readTopology(s));
 net2 = readTopology(s);
 if useape # needs ape version > 4.1 with read.evonet (not in 4.1)
@@ -73,7 +73,9 @@ storage.mode(tree1r$edge) <- "integer"
 class(tree1r) = "phylo"
 phy2r = list(edge=matrix(c(5,5,6,6,7,8,8,9, 6,4,8,7,3,1,9,2), 8,2),
              Nnode=5L, edge.length=c(3,5,.5,1,.6,4,1.1,NA),
-             reticulation=matrix(c(7L,9L),1,2), tip.label=c("A","B","C","D"))
+             reticulation=matrix(c(7L,9L),1,2),
+             reticulation.gamma=0.1,
+             tip.label=c("A","B","C","D"))
 storage.mode(phy2r$edge) <- "integer"
 class(phy2r) = c("evonet","phylo")
 """;